Multiple stressor effects of insecticide exposure and increased fine sediment deposition on the gene expression profiles of two freshwater invertebrate species

Freshwater ecosystem degradation and biodiversity decline are strongly associated with intensive agricultural practices. Simultaneously occurring agricultural stressors can interact in complex ways, preventing an accurate prediction of their combined effects on aquatic biota. Here, we address the limited mechanistic understanding of multiple stressor effects of two globally important stressors, an insecticide (chlorantraniliprole), and increased fine sediment load and assessed their impact on the transcriptomic profile of two stream macroinvertebrates: the amphipod Gammarus pulex and the caddisfly Lepidostoma basale. We identified mainly antagonistic stressor interactions at the transcriptional level, presumably because the insecticide adsorbed to fine sediment particles. L. basale, which is phylogenetically more closely related to the insecticide’s target taxon Lepidoptera, exhibited strong transcriptional changes when the insecticide stressor was applied, whereas no clear response patterns were observed in the amphipod G. pulex. These differences in species vulnerability can presumably be attributed to molecular mechanisms determining the cellular affinity toward a stressor as well as differential exposure patterns resulting from varying ecological requirements between L. basale and G. pulex. Interestingly, the transcriptional response induced by insecticide exposure in L. basale was not associated with a disruption of the calcium homeostasis, which is the described mode of action for chlorantraniliprole. Instead, immune responses and alterations of the developmental program appear to play a more significant role. Our study shows how transcriptomic data can be used to identify multiple stressor effects and to explore the molecular mechanisms underlying stressor-induced physiological responses. As such, stressor effects assessed at the molecular level can inform about modes of action of chemicals and their interplay with non-chemical stressors. We demonstrated that stressor effects vary between different organismic groups and that insecticide effects are not necessarily covered by their described mode of action, which has important implications for environmental risk assessment of insecticides in non-target organisms.


Background
Anthropogenic activities lead to stream degradation worldwide.One of the most pressing factors is intensive land use, which poses a variety of stressors to stream ecosystems and their biodiversity [1].Agricultural streams are exposed to stressors such as increased fine sediment deposition [2], reduced discharge [3], or contamination with organic substances like pesticides [4], which are applied at adjacent fields for crop protection.Predicting the impact of simultaneously present stressors on stream ecosystems is challenging, as their combined effects can strongly deviate from expectations which were derived based on individual effects assuming stressor additivity.In fact, stressors often do not operate independently but rather mitigate (antagonistic interaction) or amplify (synergistic interaction) each other [5].These interactions can arise from a wide variety of mechanisms such as physico-chemical properties of stressors and complex spatial and temporal exposure patterns and might be further modulated indirectly via species interaction and trophic cascades [6,7].
For instance, precipitation and flood events lead to elevated fine sediment levels in streams due to sediment erosion and surface runoff.The subsequent water turbidity negatively affects visual predators [8] and primary production [9], resulting in top-down or bottomup effects in the food web, respectively.Sedimentation further decreases the structural complexity of microhabitats because interstitial spaces in the streambed and hyporheic zone are clogged [10,11].In agricultural areas, increased fine sediment load is accompanied by chemical stressors like pesticides, which also enter the stream via surface runoff, i.e., simultaneously with sediment input [12,13].Since contaminants such as heavy metals or pesticides attach to small particles, sedimentation can either represent a major source of chemical pollution [14,15] or can reduce its bioavailability and mitigate the toxic effects of a pollutant [16].
Although it is crucial to understand the environmental impact of these globally important stressors, few studies have addressed multiple stressor effects of pesticides and fine sediment on freshwater invertebrates [e.g., [17][18][19].Only one focused on a substance relevant for invertebrates, the insecticide chlorpyrifos [19].Based on specimen abundance information, i.e., population level data, all studies identified increased fine sediment load as the major stressor for macroinvertebrate communities, while no or minor effects were attributed to pesticides [17][18][19].However, relying solely on abundance information can be misleading, as organisms can remain at affected sites due to physiological coping mechanisms at the individual level and even a strong physiological stress might not be immediately reflected in population abundance data [20].For instance, although the chlorpyrifos exposure in [19] did not induce mortality of macroinvertebrates, the authors found evidence for indirect effects such as feeding inhibition, reflected in reduced organic matter decomposition [19], emphasizing the need for an integration of physiological endpoints.Because such physiological responses are regulated by gene expression networks, studying stressor effects at the transcriptomic scale can provide significant insights into the mechanistic stress responses of an organism long before abundances change [20].
To the best of our knowledge, no study to date has looked at the impact of the individual and joint application of an insecticide and increased fine sediment load at the transcriptional level within freshwater invertebrates.In the current study, we addressed this research gap and assessed single and combined effects of increased fine sediment deposition and the insecticide chlorantraniliprole on the gene expression profiles of the caddisfly Lepidostoma basale (Kolenati, 1848), and the amphipod crustacean Gammarus pulex (Linnaeus, 1758).
Chlorantraniliprole is an anthranilic diamide and designed for crop protection against lepidopteran pest species [21].It selectively binds to the insect ryanodine receptor, a non-voltage-gated calcium channel controlling the intracellular calcium homeostasis, resulting in uncontrolled calcium release in neurons and muscle cells [22].Although chlorantraniliprole was initially reported to be highly selective for target pest species [22,23], it was shown to be acutely toxic for aquatic macroinvertebrates at concentrations, which are observed for insecticides in agricultural streams [4,13,24].
For the species under investigation, no toxicity information is available, but ecotoxicological studies indicated that caddisflies are more sensitive to chlorantraniliprole than amphipods (for the caddisfly Chimarra aterrima and the amphipod Gammarus pseudolimnaeus LC 50 values of 11.7 µg/L and 35.1 µg/L were determined, respectively [25]).Despite its fast photolysis [26], chlorantraniliprole is persistent in the environment and has been found to accumulate in sediments [24].Therefore, both the dissolved and the adsorbed phase pose a relevant risk to aquatic biota, but the exposure phasedependent physiological responses are largely unknown.Exploring their molecular basis could provide insights into mechanisms underlying the different sensitivity toward chlorantraniliprole between caddisflies and amphipods, which both possess key positions in aquatic food webs [27,28].
We used an innovative mesocosm setup, the ExStream system [29], that allows to integrate ecological complexity when studying the interplay between chlorantraniliprole and increased fine sediment load.L. basale and G. pulex have slightly different habitat preferences which may determine their susceptibility to the applied stressors.For instance, the clinger organism L. basale is strongly associated with wood and typically attaches to woody structures in streams [30].As such, L. basale does not directly rely on a structural heterogeneity of the channel substratum and might thus be less impacted by fine sediment deposition.In contrast, G. pulex tends to use the interstitium as shelter for predation or during low flow periods in natural streams [31].Moreover, sheltering of G. pulex may reduce its exposure to the insecticide stressor.If this avoidance strategy is, however, prevented by fine sediment clogging the interstitium, stressor interactions could arise.
Fine sediment deposition has been shown to induce a downregulation of genes associated with the energy metabolism in the amphipod Gammarus fossarum [20].Since metabolic depression is a common physiological response to environmental stressors in different groups of invertebrates [32,33], our first hypothesis is that (i) both species react with energy allocation and metabolic suppression in response to the fine sediment stressor.
Stressor-induced molecular response pathways, which determine the cellular affinity toward a stressor, are expected to be more conserved between more closely related taxa.Since chlorantraniliprole is designed as pest control agent against lepidopteran species, which represents the sister order of caddisflies (Trichoptera), our second hypothesis is that (ii) the insecticide has a more pronounced effect on the gene expression profile of L. basale than on G. pulex.
Due to the described mode of action of chlorantraniliprole, we further hypothesize that (iii) chlorantraniliprole induces a differential expression of genes involved in maintaining the calcium homeostasis and therefore in neuromuscular processes, at least in L. basale.

Outdoor mesocosm experiment
We aimed to disentangle single and joint exposure effects of the insecticide chlorantraniliprole and increased fine sediment on the transcriptional profile of L. basale and G. pulex under semi-natural conditions.The ExStream system (Fig. 1A; ExStream Systems Ltd., Dunedin New Zealand) was set up next to the the Bieber, a fine substrate dominated siliceous low mountain stream (Hessian Ministry of the Environment) in Hesse, Germany (50°09′38.9"N,9°17′58.6"E;213 m a.s.l).The river belongs to the Rhine-Main-Observatory (https:// deims.org/ 9f9ba 137-342d-4813-ae58-a6091 1c3ab c1), a long-term ecological research (LTER) site [34,35].The experiment was conducted from August 09 to September 19, 2020 and comprised a 21-day colonization period followed by a 21-day stressor period (Fig. 1B).Stream water was pumped continuously into four header tanks, each supplying 16 circular mesocosms (diameter = 25 cm, vol.= 3.5 L, area = 450 cm 2 ) permanently with fresh stream water.The stream water flow within the mesocosms was calibrated daily to 2 L/min.Each mesocosm was filled with substrate collected from the Bieber stream bed: 600 g gravel < 1 cm, 600 g gravel 1-3 cm, 300 g stones > 3 cm, and 3 large flat stones.This composition reflected stream bed areas where flow velocities resembled the ones realized in the mesocosms (~ 10 cm/s).The water temperature was assessed throughout the whole experiment in 5 min intervals using HOBO pendant loggers (Onset, Bourne, United States).Other physico-chemical parameters (oxygen concentration, pH, conductivity) were measured once a week with a MultiLine P4 probe (WTW, Washington, D.C., United States).

Colonization period
Natural drift colonization of the experimental system was possible during the whole experiment for aquatic organisms < 4 mm via the water pumps.This passive drift colonization was supplemented by active sampling of the macroinvertebrate community one week prior to the start of the manipulative period (day -7).Eight kick-net samples, each comprising two pool and two riffle microhabitats, were taken along a 50 m river section 250 m upstream of the experimental setup (50°09′41.0"N,9°18′12.0"E),and specimens obtained from one kick-net sample were randomly distributed across eight mesocosms.After specimens were supplied to each mesocosm once, the procedure was repeated but kick-net samples were taken from a different 50 m river section (50°09′41.0"N,9°18′14.0"E).During the subsequent acclimatization phase (days -7 to 0), all drifting organisms were captured with a sieve (mesh size: 2 mm) placed in the central circular opening (Fig. 1C) and reintroduced in the mesocosms.

Stressor period
During the stress period, the macroinvertebrate community was exposed to three different concentrations of the insecticide chlorantraniliprole and increased fine sediment in a full-factorial design with 8 replicates (i.e., mesocosms) per treatment.Control mesocosms received no stressor treatment.The stressors were applied in a randomized block design, i.e., each header tank supplied stream water to two experimental replicates of each unique treatment combination.Fine sediment (< 2 mm particle size) was collected from the adjacent field and 450 mL were manually added to each mesocosm that was exposed to increased fine sediment levels.This led to approximately 90% fine sediment coverage, which has been shown to be realistic for streams in areas with intensive agriculture [2] and to induce stressor responses of the macroinvertebrate community [36,37].The fine sediment level was constant throughout the whole manipulative period.The height of the sediment layer was measured on day 4 and day 21 at three distinct points in each mesocosm (Additional file 2: Table S1).
We used the commercially available product Coragen (batch no.: MAY19CL13A, DuPont, Wilmington, Delaware, United States) to apply chlorantraniliprole in three different concentration levels, i.e., low, medium and high concentration treatments.Notably, other ingredients of the formulation Coragen likely evoke gene regulatory responses as well, which cannot be separated from the effect of chlorantraniliprole in this study.A dosing pump delivered an insecticide stock solution to three dripper lines.Each dripper line then supplied the insecticide directly to the respective mesocosms.The different insecticide levels were targeted by different pumping pulses of the dosing pump, thereby delivering different volumes of the insecticide stock solution to the dripper lines and finally to the mesocosms.During the first 4 days of the stressor period, we aimed for nominal concentrations of 0.2 µg/L, 2 µg/L, and 20 µg/L, but expected to detect lower concentrations in our experiment due to potential   [24,25].Sublethal effects were reported at concentrations as low as 0.2 µg/L in the caddisfly Sericostoma personatum [38].Since environmental concentrations of chlorantraniliprole of up to 28 µg/L were reported in the past [39], our nominal concentration gradient represented a field realistic exposure scenario with stressor effects expected to range from weak (0.2 µg/L) to strong (20 µg/L).
When fine sediment and insecticides enter the natural stream during rainfall, aquatic organisms are shortly confronted with an intense stressor exposure ('stressor pulse'), followed by a reduced but persistent stressor occurrence ('base exposure').To resemble these natural dynamics, the insecticide concentrations were lowered by the factor of 10 after four days for the remaining 17 days.The insecticide concentration in different treatment combinations was measured from water samples taken at day 3 (Additional file 2: Table S2) and day 21.Water samples (500 mL) were collected from the outflow of eight randomly selected mesocosms, each representing one unique treatment combination.

Sampling
In this study, we focused on two species and sampled the specimens at day 4 to identify the stressor-induced transcriptional response after the stressor pulse.Because the experiment was continued for further 17 days, we avoided any sampling-induced disturbance of the mesocosms.Therefore, we were restricted to sample the specimens directly from the water column within the mesocosm or from the large flat stones.During this sampling procedure, we were not able to sample gammarids in one increased fine sediment mesocosm.From each of the remaining mesocosms, three G. pulex and five L. basale specimens were sampled, irrespective of sex or developmental stage.In total, 189 gammarids (3 × 63 mesocosms) and 320 caddisflies (5 × 64 mesocosms) were sampled, directly preserved on dry ice, and stored at -80 °C.

DNA barcoding
The extracted nucleic acids were used as template for the amplification of the cytochrome c oxidase subunit one barcoding fragment using HCO2198-JJ/LCO1490-JJ (Eurofins Genomics, Konstanz, Germany) primers [41].Per sample, the reaction contained 5 μL DreamTaq (Thermo Fisher Scientific, Waltham, Massachusetts, United States) master mix, 0.03 μL of each HCO2198-JJ and LCO1490-JJ (100 μM), 2.9 μl molecular grade water (Thermo Fisher Scientific, Waltham, Massachusetts, United States), and 1 μL nucleic acid extract.The initial denaturation of DNA was conducted for 3 min at 95 °C, followed by 45 cycles of denaturation at 95 °C for 30 s, annealing at 46 °C for 30 s, elongation at 72 °C for 1 min, and a final elongation at 72 °C for 5 min.Amplicons were purified with 1 unit ExoI and 10 units FastAP (both Thermo Fisher Scientific, Waltham, Massachusetts, United States).Unidirectional sequencing was performed at Eurofins Genomics (Konstanz, Germany) and barcode sequences were queried against the Barcode of Life Data System database [42].Through this, all trichopterans were verified to be L. basale and all amphipods G. pulex.

RNA purification
Nucleic acid extracts were digested with 2 units DNase (Ambion, Waltham, Massachusetts, United States) in order to remove genomic DNA.Then, the samples were mixed with RNA binding buffer [40] and magnetic beads in TE-min buffer (4:1, v/v, alcoholic RNA binding buffer, aqueous solutions, i.e., digested extract and beads in TEmin).The RNA was washed twice with 80% ethanol and finally resolubilized in TE-min buffer.The concentration and quality of the cleaned RNA samples were assessed on a Fragment Analyzer system using the 15 nt RNA kit (both Agilent, Santa Clara, California, United States).
Based on sample quality, RNA extracts of each mesocosm were pooled per species.Accordingly, a sequencing library represented one mesocosm, comprising the pooled RNA extracts of either G. pulex or L. basale specimens.RNA extracts obtained from G. pulex specimens that showed a low quality or quantity were excluded from downstream processing.Therefore, 35 G. pulex libraries contained pooled RNA extracts from all three specimens, 26 libraries included RNA from two specimens, and 2 libraries contained RNA of only one specimen.All L. basale libraries consisted of 5 specimens.Protocols to prepare the GITC lysis buffer and the RNA binding buffer can be found under https:// bomb.bio/ proto cols/ (protocol 8.2).

Library preparation and sequencing
All 127 RNA libraries were sent to the West German Genome Center (WGGC) for further library preparation and sequencing.mRNA enrichment was conducted with 800 ng total RNA as input.The poly-A selection and cDNA library construction were performed using the NEBNext Ultra II Directional RNA kit (New England Biolabs, Ipswich, Massachusetts, United States).Libraries were paired-end sequenced (150 bp) on a NovaSeq 6000 (Illumina, San Diego, California, United States).

Chlorantraniliprole analysis
To qualify chlorantraniliprole concentrations, water samples were pre-concentrated via solid-phase extraction using Oasis HLB 6 cc 500 mg extraction cartridges (Waters Corp., Milford, Massachusetts, United States).After conditioning the cartridges with 5 mL of methanol (Carl Roth, Karlsruhe, Germany) and equilibrating with 10 mL ultrapure water, an aliquot of ~ 500 mL (exact volume was recorded) was loaded to the cartridge at a flow rate of ~ 7 mL/min.Subsequently, the cartridges were dried with nitrogen (Carl Roth, Karlsruhe, Germany) for 2 h.The samples were eluded with 6 mL methanol:ethyl acetate (1:1, v/v, methanol LC grade, ethyl acetate; Carl Roth, Karlsruhe, Germany) and 2 mL methanol.The eluates were then dried to 50 µL at room temperature under a gentle nitrogen flow and reconstituted with 450 µL ultrapure water.The resulting extracts were centrifuged (4000 rpm, 30 min) and the supernatants were used for chemical analysis.To correct for potential evaporation and extraction losses and matrix effects, selected standards in ultrapure water as well as in stream water from the Bieber were treated identically.The samples were quantified using an Exactive (LCHRMS) Orbitrap system (Thermo Fisher Scientific, Waltham, Massachusetts, United States).The chromatographic separation was achieved with an Atlantis T3 5 μm 3.0 × 150 mm column (Waters Corp., Milford, Massachusetts, United States).A calibration row of chlorantraniliprole (PESTANAL analytical standards, Merck, Darmstadt, Hesse, Germany) was linear from 3 µg/L (analytical limit of quantification (LOQ) in water calculated as 0.007 µg/L) to 500 µg/L.Details on gradient settings were given in [43].

Data analysis
Raw sequencing libraries were homopolymer trimmed using a custom C + + [44] program followed by quality and adapter trimming with the Cutadapt v3.2 [45] wrapper script TrimGalore! v0.6.6 (https:// github.com/ Felix Krueg er/ TrimG alore) in paired-end mode, retaining only bases with Phred > 20 and reads with a length ≥ 25.The quality trimmed read data were used for de novo transcriptome assembly with Trinity v2.13 [46] and rnaSPADES v3.15.0 [47] for L. basale and G. pulex, respectively.Assemblers were run in default modes, except for parameters which control memory usage and multi-threading.Details about the assembly evaluation are provided in Additional file 1. Quality trimmed reads were mapped against the transcriptomes using bow-tie2 [48] and transcript abundances were estimated with RSEM v1.3.3 [49].The number of reads obtained after sequencing and quality trimming as well as results of the transcriptome assembly evaluation are presented in Additional file 2: Table S3.
The functional annotation of the generated transcriptomes was performed on protein level.Protein coding sequences were identified with TransDecoder v5.5.0 (https:// github.com/ Trans Decod er/).HMMER v3.3 (http:// hmmer.org/) was used to search for protein domains in the putative protein sequences based on the Pfam database [50].Blastp v2.9.0 [51] searches for the putative protein sequences were conducted (e-value < 1e-5) against the Swissprot/Uniprot database [52] for protein prediction.Gene ontology (GO) terms [53,54] were retrieved from the eggNOG v5.0.2 database [55] with the eggNOG-mapper v2.1.9[56].DIA-MOND v2.0.15.153 [57] searches were performed to query the final protein sequences against the database, using ' Arthropoda' as target taxon.Hits were required to have an e-value ≤ 1e−5, sequence identity ≥ 50% and a bitscore ≥ 50 to be included in the final annotation.For G. pulex, the annotation data were scarce, i.e., less than 20% of all genes tested for differential expression were annotated with GO terms.Scarce annotation data are a general problem for functional profiling of nonmodel organisms and crustaceans have been found to be poorly represented in databases [58,59].Therefore, we iteratively supplemented the annotation of the G. pulex transcriptome with annotations retrieved from the protein sets of (i) the amphipod Hyalella azteca (GenBank accession no.GCA_000764305.4), (ii) the cladoceran D. pulex (https:// wflea base.org), and (iii) the invertebrate Uniprot/Swissprot TrEMBL database.Blastp hits were considered valid if they showed an e-value ≤ 1e−5, sequence identity ≥ 50% and a bitscore ≥ 50 and were only included if no annotations were obtained from the previous annotation round.

Differential expression analyses and clustering of insecticide-responsive L. basale genes
Transcript abundances were summarized to gene level estimates using tximport v1.22.0 [60].Isoform-to-gene relationships were inferred based on assembler assumptions; therefore, 'genes' are approximated, and the term is used loosely here.Only genes with ≥ 20 normalized counts in at least eight samples were included.Gene counts were modeled in DESeq2 v1.34.0 [61] specifying the model design ~ Sediment*Insecticide. Variance stabilized counts were extracted, and a surrogate variable analysis (SVA) was performed with the sva R package v3.42.0 [62] (full model: ~ Sediment*Insecticide; reduced model: ~ 1).Including surrogates in the model allows to account for data variation that is related to unmeasured sources (e.g., expression differences between sexes or developmental stages).All significant surrogates were incorporated as covariates in an updated DESeq2 model.The Wald test was used to identify differentially expressed genes with a | log2FC | > 0 (FDR adjusted p-values < 0.05).Effect sizes were shrunk with the adaptive shrinkage estimator from the R package ashr v2.2.54 [63].
The differential expression analyses revealed, by far, the strongest treatment effect in the expression data set from L. basale exposed to high insecticide concentrations.This is partly related to the low insecticide concentrations achieved in our experiment (Additional file 2: Table S2) but might be also an artifact of the low signal-to-noise ratio inherent to RNA-sequencing data obtained from wild organisms in a semi-natural setting.To explore the expression data in an unsupervised manner, we clustered insecticide-responsive genes in L. basale according to their expression profiles along an increasing insecticide concentration gradient.This co-expression analysis can inform about subtle changes in expression (i.e., small effect sizes) beyond statistical testing.Clustering was performed separately for genes which were differentially expressed in (i) the single stressor high insecticide concentration treatment and in (ii) the treatment combining high insecticide concentration and increased fine sediment using DEGreport v1.30.3 [64].

Stressor interactions
To identify genes showing a stressor interaction, DESeq2 was used to test the interaction terms for being non-zero (FDR adjusted p-values < 0.05, obtained from the Wald test).The reported effect sizes of these contrasts represent an additional log2FC which is specifically attributed to the interaction between stressors.These genes were then classified in genes showing a synergistic and antagonistic stressor interaction when the combined stressor effect was larger and smaller than expected based on the individual stressor effects, respectively.The expectation was derived from the null model of additive stressor effects, i.e., the sum of the individual log2FC.In the gene regulatory context, a positive synergistic interaction (S+) is defined as a stronger upregulation than expected based on the assumption of additivity, whereas a negative synergistic (S−) interaction indicates a stronger downregulation than expected.Vice versa, positive antagonistic (A+) interactions denote 'less upregulated than expected, ' whereas a negative antagonistic (A−) interaction refers to genes, which are 'less downregulated than expected.' For instance, a significant interaction was identified for a given gene, indicating that this gene changed its expression due to the interaction of two stressors.If this gene shows a log2FC = 1 in both individual stressor treatments, the combination of stressors is expected to evoke a log2FC = 2 under the scenario of additivity.If, however, the combined stressor treatment induced a log2FC > 2, a 'stronger upregulation than expected' was observed, and a positive synergistic interaction is reported.If a log2FC < 2 was observed, the gene was 'less upregulated than expected, ' and a positive antagonistic interaction is reported.

Functional enrichment
The functional enrichment analysis aims to identify gene ontology (GO) terms, which are more frequently observed in a specific set of genes than expected by chance.We used the elim algorithm implemented in the R package topGO v2.50.0 [65] and Fisher's exact test (p-value < 0.05) to identify overrepresented biological process terms.
The functional enrichment was performed for each set of differentially expressed genes from the different stressor treatments.All genes that were initially tested for differential expression were used as the gene universe.Upregulated and downregulated genes were tested separately.
In order to identify functional submodules within the insecticide-responsive genes of L. basale, we further performed a functional enrichment for the identified clusters from the high insecticide concentration treatments.The subset of genes belonging to each cluster was tested for overrepresentation of biological process terms against the set of genes from which the cluster was derived, i.e., either genes differentially expressed due to (i) the high insecticide concentration treatment alone or (ii) the combination of high insecticide concentration and increased fine sediment.

Visualization
Variance stabilized gene counts were corrected with a frozen SVA [62] and a principle component analysis (PCA) was performed with PCAtools v2.6.0 [66].Heatmaps were created with ComplexHeatmap v2.14.0 [67] and treatments were clustered based on Euclidean (G.pulex) and Canberra (L.basale) distances.Clusters of co-expressed genes and functional enrichment results were visualized with ggplot2 v3.4.1 [68].To determine how cluster identity changes due to the presence of increased fine sediment, the intersection of genes from (i) and (ii) was visualized in a chord diagram with the circlize library v0.4.15 [69].All statistical analyses and data visualizations were performed in R v4.1.2[70].

Chlorantraniliprole concentration
Measured chlorantraniliprole concentrations in the water samples differed strongly from the nominal concentrations: low insecticide concentrations ranged from < 0.007 μg/L to 0.14 μg/L (nominal concentration: 0.02 μg/L), medium insecticide concentrations from 0.009 μg/L to 0.65 μg/L (nominal concentration: 2 μg/L), and high insecticide concentrations from 0.13 μg/L to 2.73 μg/L (nominal concentration: 20 μg/L) (Additional file 2: Table S2).In a complementary study [71], we performed a similar experiment (i.e., exposing L. basale and G. pulex specimens from the same populations to different concentrations of chlorantraniliprole using the formulation Coragen) in an indoor setting under highly controlled conditions and observed gene expression patterns consistent to the results reported in this study.Given the match between nominal and measured concentrations under highly controlled conditions, we attribute the here observed discrepancies between nominal and measured insecticide concentrations mainly to our application approach using a dosing pump that applied micro-pulses of the insecticide stock solution to the distributing tubes.Compared to the constant concentrations under controlled conditions, this micro-pulse approach led to highly variable short-term concentrations, a limitation inherent to the semi-natural experimental setting.
Accordingly, the insecticide stressor levels cannot be treated quantitatively and constitute only relative exposure levels.However, the signal in the expression data suggests that we achieved different exposure levels that were relatively consistent in their effects within the treatments.Therefore, we retain a qualitative ranking of the insecticide levels, i.e., low, medium, and high exposure concentrations.Additional explanation is given in Additional file 1.

Treatment-induced gene expression patterns in L. basale and G. pulex
We observed a strong physiological response in the caddisfly L. basale obtained from mesocosms in which high chlorantraniliprole concentrations were applied.The PCA revealed a separation of these from the remaining treatments along the first PCA axis, which accounted for 13% variance (Fig. 2A).This differentiation was most pronounced when the insecticide was applied as single stressor and slightly declined when increased fine sediment was added as second stressor.Accordingly, we detected the strongest differential expression response in L. basale exposed to high insecticide concentration alone (5,237 genes), followed by the combined exposure of high insecticide concentration and increased fine sediment (1,080 genes) (Fig. 3A).Consistent to these results, only antagonistic stressor interactions (A + /A−) between increased fine sediment and insecticide exposure were identified in L. basale (Table 1, Additional file 2: Table S5).
The transcriptional response of G. pulex to the applied stressor combinations was less pronounced compared to the caddisfly: the PCA revealed no clear pattern for G. pulex, and the first axis explained only 4% data variance (Fig. 2B).Relatively few genes were differentially expressed, and the detected expression patterns were rather inconsistent.The highest number of differentially expressed genes were found when increased fine sediment was applied as single stressor, inducing a differential expression of 125 genes, of which the majority (73 genes) were downregulated (Fig. 3B).Among these downregulated genes, we identified the mitochondrial ND4 and COX6A2 gene.In G. pulex, 13 genes were stronger downregulated than anticipated based on the single stressor effects (Table 1, Additional file 2: Table S4).Among these genes showing a synergistic interaction, one encoded the NADH dehydrogenase subunit 4 protein.

Functional enrichment of L. basale genes responding to high insecticide concentration treatments
We could only infer robust functional enrichment results in the set of differentially expressed genes obtained from L. basale exposed to high insecticide concentrations.Due to the limited effect of the remaining treatment combinations, only a few overrepresented terms were detected in the corresponding sets of differentially expressed genes, and these terms were often represented by a single gene (Additional file 2: Table S7).Such enrichment results are unlikely to allow a biologically meaningful interpretation.In the G. pulex enrichment results (Additional file 2: Table S6), this problem is aggravated by scarce annotation data.Therefore, we focused on the functional profile of genes that were regulated in L. basale due to high insecticide concentration exposure, applied as single stressor treatment and in combination with increased fine sediment, as well as the functional profiles of the identified clusters.
Genes which were suppressed in the high insecticide concentration treatment applied as single stressor were significantly associated with the cellular machinery of protein biosynthesis (Additional file 2: Table S7): overrepresented biological process terms described the accession, replication and repair of DNA (e.g., 'chromosome organization, ' 'DNA biosynthetic process, ' 'cellular response to DNA damage stimulus'), transcription (e.g., 'mRNA polyadenylation, ' 'alternative mRNA splicing, via spliceosome'), and RNA metabolism (e.g., 'RNA processing, ' 'RNA phosphodiester bond hydrolysis') as well as translation (e.g., 'ribosome assembly, ' 'translational elongation') and post-translational protein modification (e.g., 'protein deneddylation, ' 'N-terminal protein amino acid acetylation').In addition, mitochondrial gene expression (e.g., 'mitochondrial translation') and processes (e.g., 'mitochondrial electron transport, ubiquinol to cytochrome c, ' 'mitochondrial respiratory chain complex assembly') appeared to play a dominant role in the set of downregulated genes.Further, we found evidence for the suppression of the cell cycle (e.g., 'chromosome segregation, ' 'cell cycle G2/M phase transition').Clustering of the genes that were differentially expressed in this treatment resulted in six significant groupings (Fig. 4A-F), of which three clusters represented downregulated genes (clusters 1, 4, and 5).The most prominently enriched terms, i.e., terms linked with protein biosynthesis and cell cycle, were consistently detected in all of these three clusters (Additional file 2: Table S8, Additional file 1: Fig. S2), indicating that the co-expression networks have similar functional profiles.In addition, we detected some terms referring to neuromuscular processes in cluster 4 (Fig. 4D) and cluster 5 (Fig. 4E), but these were annotated to only a few genes and referred in several cases to neuromuscular development (e.g., 'neuromuscular junction development, skeletal muscle fiber'; Additional file 2: Table S8, Additional file 1: Fig. S2).The addition of fine sediment did not significantly alter the functional profile of downregulated genes (Additional file 2: Tables S7, S9, Additional file 1: Fig. S3), and most of the genes remained co-expressed (Fig. 4K).However, the presence of fine sediment slightly modified the average expression trend of the largest cluster (single stressor treatment cluster 1, Fig. 4A) from a linear expression decrease toward an expression pattern characterized by higher gene expression under medium concentration than under low insecticide concentration (combined stressor treatment cluster 1, Fig. 4G).
Genes which were stimulated in response to high insecticide concentration could be broadly categorized into two functional groups.These functional categories were observed in both treatment combinations and all derived clusters (Additional file 2: Tables S7-S9, Additional file 1: Figs.S2, S3): the first group comprised genes involved in molecular response pathways to adverse, exogenous stimuli: specifically, we detected a strong overrepresentation of not only terms describing the insect immune system (e.g., 'Toll signaling pathway, ' 'regulation of innate immune response, ' 'humoral immune response') but also many general terms associated with physiological responses to xenobiotics (e.g., 'response to antibiotic, ' 'response to DDT') or stress conditions (e.g., 'response to starvation, ' 'response to oxidative stress') (Additional file 2: Table S7).The second dominant functional enrichment category was related to the insect larval development (e.g., 'instar larval or pupal development, ' 'determination of adult lifespan, ' 'larval salivary gland morphogenesis'), including neuromuscular development (e.g., 'larval somatic muscle development, ' 'nephrocyte differentiation'), and endocrine metabolism (e.g., 'cellular hormone metabolic process, ' 'juvenile hormone biosynthetic process').In the combined stressor treatment, we observed a similar overrepresentation of terms related to the perception of (adverse) abiotic and biotic stimuli with special focus on the insect immune system, as well as terms referring to the developmental cycle of insects (Additional file 2: Table S7).This suggests again that the addition of fine sediment induced no pronounced shift in the functional perspective of upregulated genes.Further, fine sediment addition induced no change in the expression profiles of most upregulated genes: the largest cluster derived from the single stressor treatment (cluster 2, Fig. 4B) comprised genes which similarly aggregated together in cluster 4 of the combined stressor treatment (Fig. 4K).These genes exhibited a linear increase in expression with increasing insecticide concentration (Fig. 4B, J).

Stressor effects of fine sediment deposition and insecticide exposure
Although increased fine sediment applied as a single stressor induced the largest number of differentially expressed genes in G. pulex, the overall effect in terms of the number of differentially expressed genes due to this treatment was limited in both species.This might be surprising, considering the literature evidence that fine sediment load is often among the most pressing stressors for freshwater organisms [17][18][19]36], but could be explained by the different mechanistic targets of chemical and nonchemical stressors.While chemical stressors interact with receptors, thereby inducing physiological response cascades, increased fine sediment levels do not directly interfere with cells at the molecular level but rather alter physico-chemical conditions (e.g., light, oxygen) of microhabitats [10,11].Therefore, immediate effects on the transcriptome may be more difficult to detect.Similar observations were made in a previous ExStream experiment in which fine sediment levels, flow velocity, and salinity were manipulated for 22 days and macroinvertebrate community responses as well as the transcriptional stress response of the amphipod G. fossarum were assessed [20,36].During this experiment, fine sediment deposition was identified as the most pervasive stressor for the macroinvertebrate community [36], but the strongest transcriptional effects were observed in response to the chemical stressor treatment, i.e., increased ion concentrations [20].These results further suggest that the limited effect of fine sediment observed in this study is not a result of insufficient stressor exposure time.Since the quantity of added fine sediment is comparable to previous ExStream experiments [36,37], we further argue that the amount of added fine sediment was sufficient to act as stressor for stream macroinvertebrates.However, it is possible that the fine sediment stressor did not alter habitat characteristics in the mesocosms that are specifically relevant for the two study species L. basale and G. pulex.Because an ecological dimension shaping stressor effects is inherently part of a semi-natural setup such as the ExStream system, stressor effects must be discussed in the ecological context of both species: fine sediment load directly increases water turbidity, thereby potentially affecting macroinvertebrate taxa that directly rely on primary production or visual predators [8,9].Further, the subsequent deposition leads to a homogenization of the channel substratum [10,11].
For the leaf-shredding organisms L. basale and G. pulex, the food availability is unlikely affected since leaf material can be found in sufficient quantities in the leaf packs within the mesocosms, as has also been reported by [72].Another consequence of particle load in the water is impeded respiration when gills are covered [73].While this may affect G. pulex, in which we detected a slightly stronger fine sediment effect, caddisflies such as L. basale largely rely on cutaneous respiration [74] and might therefore be less impacted by water turbidity.A similar observation was made in a recent multiple stressor study, in which fine sediment addition induced functional shifts of the macroinvertebrate community by favoring organismic groups which rely on integumentary rather than branchial respiration [72].
In line with our expectations, we detected strong transcriptional changes in the caddisfly L. basale when the insecticide stressor was applied, as opposed to G. pulex in which we could not infer a clear physiological response associated with chlorantraniliprole exposure.We can only speculate to which degree the observed differential species vulnerability is the result of differences in toxicant uptake between the species, driven by biological and ecological traits, and varying intrinsic sensitivity.However, we argue that the overall high toxicological sensitivity of L. basale is inherently driven by predispositions present at the molecular and the cellular level, which is in line with ecotoxicological studies reporting higher toxicity of chlorantraniliprole for aquatic insects than for amphipods [25] and our indoor experiment [71]: first, the stressor receptors at the molecular level likely show a higher affinity toward the insecticide in caddisflies than in amphipods.This is, for instance, illustrated by the divergence of the ryanodine receptor amino acid sequences between insects and crustaceans (Additional file 1: Fig. S1).Second, water-borne insecticide exposure at the cellular level might be stronger in caddisfly larvae than in amphipods due to the varying degree of cuticle sclerotization.These differential toxicodynamics might be further influenced by differences in ecological requirements, which may have enhanced insecticide exposure in the mesocosm for L. basale compared to G. pulex: the strong association of L. basale to wood [30] narrows its ecological niche, thereby potentially reducing its ability to find shelter from the insecticide stressor in less affected microhabitats within the mesocosm.As a more flexible and relatively mobile species, G. pulex could have effectively avoided the insecticide stressor by withdrawal to microhabitats where contaminant exposure is minimized.Apparently, fine sediment deposition had no significantly adverse impact on this avoidance behavior.
We conclude that fine sediment addition did not induce strong transcriptional stress responses in G. pulex and L. basale in our experiment.Therefore, we reject our first hypothesis that fine sediment addition induces a metabolic depression in the two invertebrate species.In contrast, applying the insecticide stressor in its highest concentration resulted in a highly pronounced effect on the transcriptome of L. basale, whereas no considerable transcriptional changes were detected in G. pulex.In favor of our second hypothesis, we attribute Fig. 4 Clustering of insecticide-responsive genes according to their expression profile along an increasing insecticide concentration gradient.A-F Differentially expressed genes in L. basale under high insecticide concentration.G-J Differentially expressed genes in L. basale under high insecticide concentration combined with increased fine sediment.Z-scores represent the scaled and centered gene expression value.K Cluster identities of differentially expressed genes in the single (purple arc, right) and the two-stressor treatment (red arc, left).Black lines indicate groups of genes for which the average expression pattern is similar between the two stressor scenarios.The purple and red arcs only contain the genes that are differentially expressed in both treatments the species-specific insecticide sensitivity to the degree of molecular pathway conservation between butterflies, caddisflies, and amphipods and conclude that the magnitude of stressor effects depends on molecular stress receptor mechanisms.These comprise different mechanistic targets of stressors at the molecular level as well as the cellular affinity toward stressors, which differs between evolutionary lineages.In addition, we acknowledge that ecological determinants such as varying stressor exposure or the ability to avoid stressors are likely to shape stressor effects under natural conditions [75], which we aimed to simulate in our mesocosm experiment.

Consistent stressor effects, antagonism and synergism
Identifying stressor interactions is crucial to improve our understanding of the ecological impact of co-occurring stressors in natural systems.At the same time, an identification of consistent stressor effects is similarly important [76], as they provide the complementary perspective on multiple stressor dynamics.Despite the strong discrepancies between the stressor effects on L. basale and G. pulex, we can report the following consistent effects at the inter-and intraspecific level: (i) a limited effect of fine sediment on the transcriptomic response, (ii) a downregulation of mitochondrial genes in response to environmental stressors such as fine sediment deposition (G.pulex) and insecticide pollution (L.basale), and (iii) a predominant antagonistic nature of stressor interactions between increased fine sediment and insecticide exposure.Within L. basale (iv), the suppression of mitochondrial processes as well as the upregulation of developmental and immunity-related genes evoked by insecticide exposure persisted, independent of the expression profile of genes or the mitigating effect by increased fine sediment addition.The associated genes were consistently co-expressed in both treatment scenarios, indicating that these genes are part of key gene regulatory networks, which contribute to the protective pathways that allow L. basale to cope with sublethal insecticide effects.
We detected only antagonistic gene expression changes in L. basale, and fine sediment addition substantially reduced the number of insecticide-responsive genes.These observations suggest that the bioavailability of chlorantraniliprole is reduced by fine sediment addition.While other studies reported food-related uptake as the most relevant source of pesticide exposure for freshwater organisms [77], our findings imply that one of the main exposure sources is dissolved insecticide in the aqueous phase and adsorption to fine sediment particles decreases the insecticide concentration in the water.While this pattern might be specific for a shredder organism like L. basale, other types of stressor interactions could arise for organisms that show a strong ecological association to sediment or soil.

Insecticide-induced molecular response mechanisms in L. basale -a comparison between indoor and outdoor experiments
The insecticide mode of action is clearly described: binding of chlorantraniliprole to the insect ryanodine receptor releases calcium from its intracellular storages, resulting in uncontrolled muscle contraction, paralysis, and eventually the death of the target pest organism [22].In L. basale, several biological process terms referring to calcium homeostasis/signaling and neuromuscular processes were detected as functionally enriched in the high insecticide concentration treatments.However, these terms often refer to (embryonic or larval) neuromuscular development (single stressor treatment: clusters 3-6; two-stressor treatment: clusters 2, 3) (Additional file 1: Figs.S2, S3, Additional file 2: Tables S8-S9) than to processes directly associated to the calcium homeostasis or involved in muscle contraction and were in many instances annotated to only a few genes.Therefore, we conclude that these genes play a subordinate role in the molecular response mechanisms of L. basale exposed to the insecticide stressor and reject our third hypothesis.
Instead, our data suggest that the caddisfly reacts to the highest insecticide concentration treatments with (i) a suppression of the protein biosynthesis, particularly of genes involved in mitochondrial processes, (ii) an increase in expression of genes involved in immunity, and (iii) a stimulation of the genetic developmental program.Since these results are highly consistent with the expression results obtained from [71], in which we exposed L. basale to an insecticide concentration gradient of 0.3-19.5 µg/L (Additional file 1), we argue that the observed effects are not a result of the relatively low insecticide concentrations achieved during this experiment.However, in the indoor experiment, the stimulation of immunity genes was less pronounced than in this study, and we detected additionally an overrepresentation of muscular terms under higher insecticide concentrations [71].Taken together, we propose the following molecular response mechanism: lower insecticide concentrations evoke rather general physiological stress responses, which comprise a metabolic depression (reflected in suppression of the protein biosynthesis, the cell cycle, and reduced mitochondrial activity) or the stimulation of cellular pathways associated with the perception of external stimuli and general stress responses, including the immune system.The expression of innate insect immune genes such as antimicrobial peptides can be induced by non-immune stressors [78] and a similar stimulation of immune gene expression in response to chlorantraniliprole exposure was observed in other non-target insects such as the honeybee Apis mellifera [79] and the fruit fly Drosophila melanogaster [80].Other arthropods such as the amphipod H. azteca showed a similar differential expression of immune genes following insecticide exposure [81], indicating a signaling overlap between immunity and molecular stress perception which is conserved across evolutionary lineages.For instance, reactive oxygen species (ROS) which are excessively generated during mitochondrial stress [82] are linked to the activation of innate immune pathways in vertebrates [83] and invertebrates [78,84].Cellular oxidative stress could further contribute to the activation of unspecific stress responses such as cellular detoxification pathways, which were induced following chlorantraniliprole exposure in freshwater macroinvertebrates such as the non-biting midge Chironomus riparius [85] or the caddisfly S. vittatum [38].Since cellular detoxification mechanisms and (undirected) immune responses are expected to come at the expense of increased energetic maintenance costs, even low insecticide exposure concentrations could affect growth rates or emergence patterns, as has been found for, e.g., C. riparius exposed to chlorantraniliprole [85].
Under higher insecticide concentrations, as the ones achieved in the indoor system, we observed in addition to the suppression of mitochondrial processes an inhibition of muscular processes [71].Considering that muscle tissue is rich in mitochondria, these two mechanisms are not exclusive and could be explained by a disruption of mitochondrial processes specifically in muscle cells.Further, alterations of the developmental program became more prominent under stronger insecticide stress, whereas the induction of immune genes played a subordinate role.Changes in the developmental program could be the result of endocrine disruption, an off-target effect known for many insecticides [86], or represent a physiological coping strategy to mitigate toxic effects [87,88].Due to the importance of calcium as a signaling molecule, an alternative explanation might be an insecticide-induced alteration of intracellular calcium levels (for further discussion, see [71]).However, the disruption of the calcium homeostasis is clearly not reflected in our data and further research is required.

Conclusions
Our results show that multiple stressor effects are shaped by various factors such as molecular mechanisms of stress perception.Ecological requirements, which lead to different stressor exposure patterns, could further contribute to variation of stressor effects in different organismic groups.Although the molecular response pathways clearly diverge between amphipods and caddisflies, we observed a predominantly antagonistic stressor interaction in both species.Furthermore, we found evidence that the downregulation of mitochondrial gene expression represents a conserved mechanism to cope with environmental stress.Suppression of the respiratory chain will decrease the energy budget available to an organism, resulting in negative long-term effects at the population level.Consequently, biotic communities in freshwater systems will be affected differently by stressor exposure, depending on, e.g., community composition or molecular targets of chemical stressors.The latter are not always known since insecticides are complex chemical mixtures and effects are not only induced by the active compound, adding further layers of complexity to multiple stressor dynamics in natural systems.Our findings underpin the need of gene expression analyses to supplement multiple stressor research with endpoints informing about stressor-induced physiological mechanisms and have important implications for results obtained from ecotoxicological assays, which typically test pure substances and are limited in time and biotic complexity.

Fig. 1
Fig.1Overview of the field experiment.A Stream water was pumped into header tanks and supplied to 64 mesocosms.The macroinvertebrate communities living in the mesocosms were exposed to the stressors in a 4 × 2 factorial design (4 pesticide levels, 2 fine sediment levels) with 8 replicates per treatment.B Time course of the experiment.A 21-day colonization period was followed by a 21-day stressor period.For this study, specimens were sampled after the insecticide pulse (day 4, red square).C Stream water entered the mesocosms via the inflow, flowed in a clockwise orientation, and left the system through the central circular opening.The sieve in the outflow allowed sampling of drifting organisms to estimate stressor effects for the macroinvertebrate community

Fig. 2
Fig. 2 Biplot of the first two principal component axes for A L. basale and B G. pulex

Fig. 3
Fig.3Heatmap of differentially expressed genes in L. basale (A) and G. pulex (B).Gray bars represent genes with no significant differences in expression between the treatment and control conditions (adjusted p-value ≥ 0.05).Red and blue numbers represent up-and downregulated genes, respectively.Fold changes are reported at the binary logarithmic scale (Log2FC)

Table 1
Number of genes that showed antagonistic (A + /A−) or synergistic (S + /S−) stressor interactionsPositive synergistic (S+) and positive antagonistic (A+) interactions were reported when the upregulation of a gene was stronger and weaker, respectively, than assumed based on additive stressor effects.Negative synergistic (S−) and negative antagonistic (A−) interactions were reported when the downregulation of a gene was stronger and weaker, respectively, than expected based on additivity