Terrestrial herbivory drives adaptive evolution in an aquatic community via indirect effects
Data files
Oct 01, 2024 version files 639.48 MB
-
README.md
3.01 KB
-
ScriptsAndData.zip
639.48 MB
Abstract
Indirect ecological interactions are thought to be ubiquitous in nature and may be important in shaping evolutionary processes in multitrophic communities. However, direct evidence of indirect interactions driving evolutionary processes is still lacking. Here, we characterized how the indirect effects of terrestrial insect (aphids) herbivory on macrophytes (duckweed) affected the evolution of an aquatic community using real-time evolutionary experiments in a two-year outdoor pond experiment. Aphid herbivory reduced macrophyte growth and increased the abundance of phytoplankton, which in turn increased the abundance of zooplankton consumers. Pool sequencing and phenotypic assays showed that the aphid herbivory altered the genetic compositions of Daphnia magna populations, a key member of the zooplankton community. In the second year, transplant experiments further indicated that evolutionary changes in D. magna driven by the aphid herbivory on the macrophytes were adaptive and that changes in the aquatic community altered aphid-macrophyte interactions. These results suggest that indirect ecological interactions can shape eco-evolutionary dynamics in multitrophic communities.
README: Terrestrial herbivory drives adaptive evolution in an aquatic community via indirect effects
Here, we tested whether changes in the intensity of terrestrial herbivory on macrophytes can drive the evolution of species in an aquatic community via indirect effects, and whether changes in the aquatic community can alter interactions between macrophytes and terrestrial insect herbivores. To this end, we established a replicated semi-natural aquatic community using outdoor experimental ponds. We grew the macrophyte Spirodela polyrhiza (giant duckweed) in the experimental ponds and manipulated the intensity of herbivory by waterlily aphids (Rhopalosiphum nymphaeae) on the duckweed. We then measured the ecosystem changes in the aquatic community and assessed the evolutionary responses of Daphnia magna, a key consumer of phytoplankton. By performing transplant experiments, we quantified the consequences of changes in the aquatic community on the fitness of aphids and duckweed and assessed the eco-evolutionary dynamics mediated by indirect interactions.
Description of the data and file structure
The data folder contains the results from different experiments.
ADD-Monitoring_v2.xlsx: Results from ecosystem monitoring on the ponds. Each sheet contains different parameters. The units are provided in each sheet. NA refers to data not available.
DaphinaTransPlant.xlsx: Results from Daphina transplant experiments. Each sheet contains the results from one round of experiments. The number refers to the number of individuals.
FeedbackExperiment2022.xlsx: Results from duckweed and aphid transplant experiments. Each sheet contains the results from one round of experiments. The number refers to the total frond weight (in mg).
filt.cmh2021.txt: Results from CMH test from the 2021 pool-seq data. Columns are: chromosome name, position, P-value, SNP-ID
filt.cmh2022.txt: Results from CMH test from the 2022 pool-seq data. Columns are: chromosome name, position, P-value, SNP-ID
final.dp20_400.AD.txt: Final SYNC file obtained from the BAM files. This file is used to calculate allele frequencies and selection.
final.dp20_400.recode.vcf.gz: The VCF file contains the SNP information from 2021 and 2022 pool-seq data.
glmmTMB.rds: Results from the glmmTMB tests.
PondInf.txt: Information on each pond and its treatment.
resistotypes_binary.csv: Results from parasite-resistance assays. R: resistance (no parasite attachment); S: susceptible (with parasite attachment); NA: data not available.
treatments.csv: Treatment information for parasite-resistance assays.
SequenceID_info.txt: Sample information from the sequencing data.
snplist.txt: A list of SNPs is included for generating glmmTMB results.
SRRInfo: Pond and treatment information for the short reads accession IDs.
Code/Software
The R markdown file (EAWAG_2024.Rmd) contains the scripts for generating figures.
The bash scripts for genomic analysis are deposited in the GitHub (see Related Works section)
Methods
Establishing outdoor experimental ponds
In 2021, we established 16 outdoor experimental ponds at the Swiss Federal Institute of Aquatic Science and Technology (EAWAG) in Dübendorf, Switzerland, using a similar setup as reported by Narwani et al. Each pond has a size of approximately 4m ´ 4m, with a maximum depth of 1.5m and a water volume of ~15m3. At the beginning of the experiment, each pond was cleaned with a high-pressure cleaner. Then, we added 220 L of leaf litter collected from a forest nearby (next to Maur/Volketswil, 47°24'33.5"N, 8°40'05.5"E) to each pond to provide nutrients and shelter for different aquatic organisms. After filling with tap water (calendar week 22), we inoculated the ponds with an initial plankton and microbe community collected from Lake Greifensee (47°20'59.592"N, 8°40'40.736"E), located near the experimental site (calendar week 24). In each pond, we added 27 L of lake water that was collected from 5-6 m depth after removing large organisms using a net with a pore size of 95µm. Additionally, we added 2.8 L concentrated plankton community (30-95 µm fraction), each collected from filtering approximately 1,750 L of lake water. The 30-95 µm plankton fraction was collected by sampling water columns from 0-10 m depth with a 30 µm mesh and pouring the concentrate subsequently through a 95 µm mesh to exclude larger organisms. Since our focal plant species, Spirodela polyrhiza, is known to produce turions (resting stages) in low phosphate environments, we additionally added 40.8g of KH2PO4 per pond (ca. 20 µM; calendar week 24).
Establishing the aphid-duckweed-daphnia populations in the ponds
In calendar week 25 of 2021, we added into each pond approximately 1,300 Daphnia magna individuals from a mix of 122 genotypes with additional algae (Tetradesmus spp. and Nanochlorpsis spp.) that were used to feed D. magna in the lab, approximately 5,000 Spirodela polyrhiza fronds from a single genotype (genotype SP102, originally collected in Switzerland and pre-cultivated in the lab before the field experiment) covering approximately 1% of the water surface and 18 great pond snails (Lymnaea stagnalis). All D. magna genotypes had originally been collected in Lake Aegelsee near Frauenfeld (47°33.48'N, 8°51.66' E), Switzerland. These clones were produced from wild-caught females in the years 2014 to 2019 and kept in the laboratory by clonal propagation. Each clone has a fully sequenced genome (Illumina short reads with an average coverage of about 20x) following the sequencing protocol.
All D. magna clones were propagated in 400-mL jars (6 jars per clone) with abundant food (green algae Tetradesmus sp.), at a 16:8 hours light:dark cycle and a temperature of 20 °C. When all clones were at their approximate carrying capacity, we mixed all clones into four 60-L tanks. The tanks were topped up to 45 L volume. This procedure was replicated 4 times. Subsamples were counted, indicating an estimated 7,300 Daphnia per 60-L tank. The following day, the four tanks were transported to the field site. The populations in the tanks were well mixed, and 20 subsamples of equal size (2 L volume) were collected from each of the four tanks, resulting in 20 samples of 8 L volume with an estimated 1,300 D. magna individuals. Sixteen sub-samples were used to inoculate the 16 experimental ponds. The initial inoculum for each mesocosm was large enough to minimize possible stochastic effects, influencing SNP and resistance phenotype frequencies. The remaining four samples were frozen in liquid nitrogen for later genomics analyses (see below), which were considered as samples from timepoint 0. From the timepoint-0 samples, we also collected 150 individuals of D. magna to estimate parasite attachment phenotype (see below).
In each pond, we introduced 200 waterlily aphids (Rhopalosiphum nymphaeae) from a single genotype into eight ponds while we kept the other eight ponds free from aphids as control. The aphid populations were developed from a single individual that was collected at the University of Münster campus (51°57'40.7"N, 7°36'55.8"E) in September 2020. The ponds were arranged in four blocks of four ponds, each containing two ponds with aphids and two control ponds. The populations of aphid, duckweed, and plankton communities were randomized.
To minimize aphid migration between ponds, we installed pond covers made from mosquito nets, and sticky flags and sticky tape were placed on the ground between ponds. Control ponds were regularly monitored for aphids. Aphids that accidentally invaded control ponds were manually removed as long as feasible.
The experiment was run for two years. The pond covers were removed over the winter season to avoid damage due to snowfall and reinstalled in March 2022. We re-introduced the aphids in calendar week 25 in 2022 to the same ponds that already received aphids the year before since the aphids do not overwinter within the ponds.
Due to an unexpected accidental extinction of pond snails in pond 6A in 2021, we removed this replicate from all analyses. We also excluded the 2022 data from ponds 1D and 3D from the aphid herbivory treatment because the aphid populations didn’t successfully establish in these two ponds in 2022.
Monitoring population dynamics in the experimental ponds
Duckweed growth was assessed based on an estimation of the surface coverage of the ponds with duckweed. Values above 100% represent multi-layered growth (via close inspections). Additionally, at each timepoint overview pictures of the ponds were taken with a camera (GoPro Hero 7 and Hero 9) that was attached to a long stick. Aphid abundance was calculated based on close-up pictures that were taken with a mobile phone from 8 random positions of each pond. Within a subsection of four pictures per pond, all aphids and fronds were counted to calculate the number of aphids per frond.
We monitored the population dynamics of plankton organisms using a previously established protocol with minor modifications. Briefly, we sampled three vertical profiles of the water column using a Leibold sampler, which consists of a 180 cm long PVC tube (5 cm diameter) that can be closed via a wire with a plug at the bottom. From each pond, we collected 10 L of water per sampling. From this, we took a 1 L subsample for nutrient analyses and 1 L for chlorophyll-a and phytoplankton analysis. During this subsampling, zooplankton was prevented from leaving the beaker with a 150 µm mesh. The remaining water (with zooplankton corresponding to 10 L pond water) was filtered through a 150 µm mesh to collect the retained zooplankton.
We estimated the total phytoplankton abundance in the pond by measuring chlorophyll-a content in the water. To this end, each sample was filtered using a 25 mm GF/F filter on a 50 mL syringe. Samples were filtered until a noticeable backpressure occurred or the filter turned green. The volume of filtered water was noted for later normalization. Afterward, excess water was removed from the filters by pushing another syringe filled with air through the filter. Filters were then removed from the support and shortly dried on a paper towel before they were placed in aluminum-covered tubes on ice and subsequently stored at -20°C until further processing. We extracted chlorophyll-a from the filters with 5 mL 90% EtOH after vortexing, sonication, and overnight incubation at 4°C in the dark. Subsequently, the extracts were passed through 0.2 μm cellulose acetate syringe filters before they were analyzed on a Hitachi 2000UV/VIS Spectrophotometer LLG-uniSPEC4 based on the absorption at 750 nm and 665 nm and an external chlorophyll-a standard. The chlorophyll-a content was calculated using the following formula:
((E665 - E750 - E665cuvette) × v) ÷ (CHL factor × V × d)
where E665 refers to extinction at 665 nm, E750 refers to extinction at 750 nm, E665cuvette refers to extinction at 665 nm of a cuvette filed with 90% EtOH, v refers to the volume of 90% EtOH, CHL factor refers to a correction factor (0.000082), V refers to the filtrated volume and d refers to the cuvette length (5cm).
To determine the taxonomic composition of phytoplankton, the water samples were preserved by adding each 5 mL of Lugol’s solution into 50 mL of pond water. Samples were stored at 4 °C until further processing. For analysis, samples were placed in a Utermöhl counting chamber for 24 hours. Afterward, the tube was removed, a cover slide was placed on top of the counting chamber, and the sedimented plankton was analyzed under a Zeiss Axiovert 135 inverted microscope. To this end, phytoplanktons were identified (mostly to genus level) and counted in 40 randomly chosen fields of view (out of 1,686) with 320X magnification and 40 fields of view (out of 6,900) with 640X magnification. If one taxon was highly abundant (>100 cells in 10 fields of view), we counted this taxon in only 10 fields of view.
The zooplankton collected from the 10 L of pond water via a 150 μm mesh was eluted from the mesh with 100 % EtOH each into a 50 mL Kautex bottle. Samples were stored at room temperature until further analyses. Samples were analyzed under a Leica M205 C stereomicroscope at 10-50X magnification, and counts were converted to individuals per liter of pond water. Adult Cladocera were identified at the species level, juvenile Cladocera and Copepoda at the genus level, and other zooplankton at higher taxonomic levels.
Monitoring nutrients and abiotic environment in the experimental ponds
Orthophosphate-phosphorus (PO43--P) was measured photometrically using an Agilent Cary 60 spectrophotometer after reaction to the phosphorus molybdenum blue complex. Total phosphorous (P) was measured after decomposition by autoclaving at 121°C with subsequent determination of PO43--P as above. Ammonium-nitrogen (NH4+-N) was measured photometrically using an Agilent Cary 60 spectrophotometer after reaction to a blue complex by Berthelot-reaction. Total organic carbon (TOC) was measured on a Shimadzu TOC-L CSH after catalytic combustion at 720°C with subsequent measurement of CO2 via an infrared detector.
The pH and oxygen concentrations of the pond water were regularly monitored using handheld sondes (Multi 3630 IDS with an FDO 925 oxygen sensor and a SenTix 940-3 pH sensor, all from WTW). The dipping probes were placed 75-80 cm below the water surface. For a continuous analysis of the water temperature and light penetration into the pond water, we deployed a logger in the middle of each pond 30-35 cm below the water surface. Six ponds (3 per treatment) were provided with a temperature logger (HOBO pendant temperature 64K, UA-001-64), and 10 ponds (5 per treatment) were provided with a temperature logger with a light sensor (HOBO Pendant Temperature+Light 64K, UA-002-64). The data were recorded every 15 minutes. The temperature was expressed as the weekly average temperature. Light data (in LUX) were summed up for each day, and weekly averages of the daily sum of light were calculated.
Assessing parasite attachment and genomic changes in the D. magna populations
At defined time points during the experiment, samples were collected to assess parasite attachment and genomic changes using a pool-seq approach. Samples at timepoint 0 were collected from the four samples from the original inoculum, as described above. Additional samples were collected from the experimental ponds. For this, we used hand-held nets with a mesh width of about 0.2 mm. Care was taken to ensure that every depth and every area of the ponds were swept with the net. The samples were brought to a laboratory and were inspected using stereo microscopes. We sorted the samples to include only D. magna. From each sample, about 250 adult D. magna were collected and frozen with little water in liquid nitrogen. Sampling dates for the pool-seq samples were September 22nd, 2021, and July 28th, 2022. Samples for the assessment of parasite attachment phenotypes were transported in 2-L bottles to the University of Basel laboratory, where 100 animals from each pond were placed individually in a 100-mL jar. Animals were tested for parasite resistance as described below. Sampling dates roughly coincided with the sampling dates for the pool-seq samples but were more spread out, because the workload associated with testing 16 times 100 animals was too high to be performed in a short time interval. Therefore, samples were taken across a time period of about 4 weeks (from September 22nd to October 16th, 2021, and from July 3rd to July 28th, 2022). Within each sampling campaign, samples were randomized (treatment and blocks were exactly balanced), and the assessment of attachment phenotypes was blinded with regard to the treatment a sample came from.
For each clone we assessed its susceptibility to attachment of five isolates of the common bacterial pathogen, Pasteuria ramosa, using the attachment assay. In short, fluorescently labeled parasite endospores were exposed to D. magna individuals. After 20 to 60 minutes, animals were assayed with a fluorescence microscope searching for endospores attached to the foregut or the hindgut. The ability of the parasite to attach to the gut epithelium is a prerequisite for infection and correlates very strongly with successful infection (6). The loci determining attachment are Mendelian but interact epistatically with each other. The diversity of our panel of D. magna clones was enriched for rare attachment phenotypes.
For the pool-seq samples, we extracted DNA from pooled individuals from each pond using the CTAB method. The total DNA was sequenced with 150bp paired-end reads at ~480X coverage for each sample using an Illumina NovaSeq 6000 instrument by Novogene in Beijing, China. Raw data was quality-checked and trimmed using TrimGalore v0.6.1, and reads were mapped towards the D. magna reference genome using BWA and SAMtools. Signatures of selection in D. magna populations from both herbivory and control ponds were detected using CLEAR based on time-series data. We then used SAMtools, PoPoolation2, and poolfstat to estimate genome-wide FST between and within treatments and sample time points, with a minimum read coverage of 20 and a maximum coverage of 400 as the cutoff. The principle component analysis was performed using the R-package “pcadapt”. The heatmap showing the similarity among each pool was created using a pairwise FST matrix estimated by poolfstat. To identify allele frequency changes for each SNP within each time point, we used the Cochran‐Mantel‐Haenszel (CMH) test. We used Bonferroni’s corrected P-value of 0.05 as the cutoff to identify significantly differentiated genomic regions. To identify genomic regions that were under selection during the two years, we used the beta-binomial mixed-effects model (from R-package glmmTMB v1.1.9) with sampling time and pond block as random factors. The genomic position of the resistance locus was identified using the sequence from Fredericksen et al..
We used the pool-seq data to test for the presence of P. ramosa in our samples. Using SAMtools, we selected the sequences not aligning with the reference genome of D. magna. These data were then used as input for Kraken2 to characterize the Daphnia-associated microbial community using a custom reference dataset that included the bacterial and fungal RefSeqs from NCBI and, additionally, the reference genome of P. ramosa. This analysis did not identify any sequence of the parasite in any of our samples.
Transplant experiments
To assess the effects of the aquatic community on the D. magna growth, we used a 140 cm long vertical PVC tube (5 cm diameter). The tube was closed at the bottom and contained 4 equally distributed openings (12 ´ 6 cm) that were covered with a 150 µm mesh to allow for water and nutrient exchange. At the beginning of the transplant experiments, we added 50 similar-sized D. magna females (without eggs) originating from the ponds into each tube. The D. magna were sorted using a stereomicroscope. The D. magna were either kept in the same pond, moved to another pond from the same treatment, or moved to another pond with a different treatment. At the end of the experiment, the D. magna inside each tube were collected on a 1 mm mesh, from which they were kept in 100% EtOH. Samples were stored in a 50 mL Kautex bottle at room temperature until further analyses. We counted D. magna individuals using a Leica M205 C stereomicroscope. We performed the experiments twice. The first experiment lasted for 13 days (between calendar weeks 23 and 25), and the second experiment lasted for 12 days (between calendar weeks 28 and 30). The logistic growth rate was calculated as log(Nend) - log(Nstart), where Nend and Nstart, refer to the number of D. magna at the end and beginning of experiments, respectively.
To assess the effects of the aquatic environment on the fitness of the duckweed and aphid, we performed growth assays using the floating boxes, which allowed for the exchange of water and plankton while preventing the escape of the duckweed and aphids via mesh. The floating boxes were constructed using plastic boxes with lids (Eurobox, Auer Packaging, inside measures: w ´ d ´ h 17 cm ´ 12 cm ´ 11.5 cm). We cut holes in the bottom (15.3 cm ´ 10.3 cm) and lid (16.8 cm ´ 11.8 cm) and covered them with steel meshes (bottom: mesh size 0.63 mm, wire diameter 0.22 mm, material 1.4301 steel; lid: mesh size 0.14 mm, wire diameter 0.112 mm, material 1.4301 steel). Each box was fitted with a styrodur frame on the outside to make it float on the pond surface.
We performed the growth assays using four duckweed genotypes (SP102, SP56, SP58, and SP65) that were originally collected in different areas in Europe. At the beginning of the experiment, we used approximately 20 duckweed individuals (fronds) from one of those genotypes for each floating box. For each pond and each genotype, we used two boxes: in one, we added five aphids, while the other one was left without aphids (control). Approximately two weeks later, we collected both duckweeds and aphids and froze them in liquid nitrogen. The dry weight was measured after freeze-drying. To estimate the number of individuals, we collected and counted up to 100 duckweed individuals for each box and measured their dry weight. Then, dry weight per frond was calculated and used to estimate the total number of duckweed individuals for each box based on their total dry weight. To estimate the number of aphids in each box, we also determined the number of aphids per frond from those up to 100 fronds and multiplied it by the total number of fronds per box. The logistic growth rate was calculated as log(Nend) - log(Nstart), where Nend and Nstart refer to the number of duckweed or aphids at the end and beginning of experiments, respectively. The resistance of duckweed was estimated by calculating the differences in growth rate between control and aphid-treated boxes for each genotype inside each pond.
Anthocyanin analysis in the duckweed populations
Samples from the duckweed populations in each pond were collected, cleaned from aphids and other large contaminants, washed with fresh tap water, dried using paper towels, and shock frozen in liquid nitrogen. Samples were kept until further processing at -20°C. Anthocyanins were analyzed by LC-MS and HPLC-PDA similar to those described by Malacrinò et al (24). In brief, approximately 10 mg of fine-ground, freeze-dried plant tissue was extracted with acidified methanol. Fresh weight was estimated based on the fresh weight to dry weight ratio of each sample. Cyanidin-3-O-glucoside was analyzed by LC-MS. Before analysis, the sample extracts were diluted 1:100 in an aqueous mix of isotope-labeled amino acids (algal amino acid mixture-13C-15N; Sigma-Aldrich). Chromatographic separation was done on a Shimadzu Nexera X3 equipped with an Agilent 1290 infinity II inline filter (0.3 µm) and a ZORBAX RRHD Eclipse XDB-C18 column (3´50 mm, 1.8 µm; Agilent Technologies). Separation was done in gradient mode with 0.05 % formic acid (Fisher Chemical), 0.1% acetonitrile (Fisher Chemical) in water as Solvent A and methanol (Fisher Chemical) as Solvent B, with the column oven set to 42 °C and a flow rate of 500 µL/min using following solvent settings (Time [min]/B[%]): 0.0/2, 1.5/2, 3.5/100, 4.5/100, 5.0/2, 6.0/2. The LC-system was coupled to a Shimadzu LCMS-8060 mass spectrometer equipped with an ESI source, which was operated in positive ionization mode and the following settings: Nebulizing Gas Flow: 3 L/min; Heating Gas Flow: 10 L/min; Drying Gas Flow: 10 L/min; Interface Temperature: 300 °C; DL Temperature 250 °C; Heat Block Temperature: 400 °C; CID Gas: 270 kPa; Q1 Resolution: Unit; Q3 Resolution: Unit. Analysis was done in the multi-reaction-monitoring (MRM) modus using following precursor to product ion transitions: Cyanidin-3-O-glucoside (quantifier), 449,10 -> 137,10; Cyanidin-3-O-glucoside (qualifier), 449,10 -> 213,15; 13C9, 15N1-Phenylalanine, 176,11 -> 129,25. The Dwell time [ms], collision energy [V], Q1Pre Bias [V] and Q3 Pre Bias [V] were: 32/32/247, -54/-54/-15, -20/-20/-20 and -14/-22/-20, respectively for Cyanidin-3-O-glucoside (quantifier)/Cyanidin-3-O-glucoside (qualifier)/13C9, 15N1-Phenylalanine. The retention time of Cyanidin-3-O-glucoside was 3.38 min and for 13C9, 15N1-Phenylalanine 2.52 min. Relative quantification of Cyanidin-3-O-glucoside was done based on the 13C9, 15N1-Phenylalanine internal standard from the isotope-labeled amino acid mixture and after normalization to the fresh weight of extracted plant material.
For absolute quantification of Cyanidin-3-O-glucoside and additional analysis of Cyanidin-3-O-(6-O-malonyl-beta-glucoside) the pure sample extracts were further analyzed via HPLC-PDA, which was done on a Shimadzu Nexera XR equipped with a photodiode array detector, a Nucleodur Sphinx RP column (250´4.6 mm, 5µm, Macherey-Nagel) and an EC 4/3 Nucleodur Sphinx RP pre-column (5 µm, Macherey-Nagel). Separation was done in gradient mode with 0.2% formic acid (Fisher Chemical), 0.1% acetonitrile (Fisher Chemical) in water as Solvent A and acetonitrile (Fisher Chemical) as Solvent B, with the column oven set to 20 °C and a flow rate of 1300 µL/min using following solvent settings (Time [min]/B[%]): 0.0/10, 8.0/21, 18.0/49, 18.1/100, 19.0/100, 19.1/10, 24.0/10. Measurement was performed with a PDA detector. Cyanidin-3-O-glucoside and Cyanidin-3-O-(6-O-malonyl-beta-glucoside) were analyzed at an absorption wavelength of 517 nm and their retention times of 7.055 min and 9.410 min, respectively. Quantification was done based on the comparison with the molar quantity of an external Cyanidin-3-O-glucoside standard curve and after normalization to the fresh weight of extracted plant material.