Larvae and plant material
The Glanville fritillary butterfly, Melitaea cinxia, L. 1758, is a widely distributed species across Eurasia and North Africa. Over the last three decades, many aspects of the ecology, life history, demography and eco-evolutionary dynamics of the Finnish population inhabiting the Åland islands, the Baltic Sea, have been intensively studied [60-62]. Our laboratory stock population was originally assembled by collecting three individual larvae from 38 (F0) families across the Åland metapopulation in September 2015, as previously described [30]. Although larvae were not genotyped, our sampling strategy insured that we worked on representative individuals from the genetic diversity of the metapopulation [63, 64]. Larvae were reared in family groups under optimum conditions (during periods of growth: Day:Night (D:N), 27°C:10°C, 12h:12h; during diapause: D:N, 5°C:5°C, 12h:12h) over two generations at the Lammi Biological station, University of Helsinki, Finland. Butterflies were mated with non-siblings partners in the laboratory, and mated females were individually isolated in small cages to lay eggs on Plantago lanceolata host plants. The host plants were checked daily for egg clutches, which were carefully transferred into individually labeled petri dishes.
On the day of emergence from the eggs, L1 (F2) larvae from each of 27 selected families (NTotal=2,160) were divided between four treatments in groups of 20 individuals each in a full factorial design (Family x Treatment), and reared in the laboratory until 3rd instar (L3), as described below. In parallel, we also reared post-diapause (F1) larvae from five families under optimum laboratory conditions (above), and collected fresh frass every day once the larvae reached 7th instar (L7). Sample sizes for each treatment and experiment are described in Table 1.
Plantago lanceolata (N = 120) was used as the larval food throughout the experiment (see below). The plants were collected as seeds across the Åland islands in 2015, and grown in optimum laboratory greenhouse conditions at the Lammi Biological Research station (D:N, 27°C:10°C, 12h:12h). Plants were watered every 3rd day. Plant leaves were only harvested for the experiment, thus preserving all natural defensive metabolites and original microbial load of the plants for the experiment. We also harvested and froze in liquid nitrogen some extra leaves to provide controls to the experiments described below.
Treatments
The 2,160 larvae from 27 families were equally divided between four treatments. Each group of 20 larvae was given daily a freshly harvested 1.7cm2 piece of randomly collected P. lanceolata leaf [40], which was supplemented differentially according to treatment (Fig 1):
- (Control): 200µl of sterile water was left to dry on the leaves before being provided to the larvae, from day1 (L1) until the larvae molt into L3.
- (Antibiotic): 200µl of the antibiotic solution was left to dry on the leaves before being provided to the larvae, from day 1 (L1) the larvae molt into L3.
- (Reinfection): 200µl of antibiotic solution was left to dry on the leaves before being provided to the larvae, from day1 to day3 (L1). On day4, 200µl of the antibiotic solution supplemented with 5% of L7’s frass was left to dry on the leaves before being provided to the larvae. From day5 and until the larvae molted into L3, 200µl of sterile water was left to dry on the leaves before being provided to the larvae.
- (Antibiotic during Reinfection): 200µl of the antibiotic solution was left to dry on the leaves before being provided to the larvae, from day1 to 3 (L1). On day4, 200µl of the antibiotic solution was supplemented with 5% of L7’s frass and left to dry on the leaves before being provided to the larvae. From day5, 200µl of sterile water was left to dry on the leaves before being provided to the larvae until the larvae molt into L3.
The antibiotic solution was prepared by mixing three anti-bacterial agents (2x10-4g.ml-1 of neomycin sulfate, with 1x10-3g.ml-1 of aureomycin, 6x10-5g.ml-1 of streptomycin) and two antifungal agents (8x10-4g.ml-1 of methyl paraben, and 6x10-4g.ml-1 of sorbic acid) as described by Chung et al. [65].
Larval performance: Development and Survival
For each larval group, transition to the 2nd larval instar was checked every day over a 13 day long period, while survival until 3rd larval instar within each group was estimated every third day. On the day the surviving larvae reached the 3rd instar, they were frozen in liquid nitrogen and stored at -80°C until further manipulated. As the larvae were not starved before being killed, the gut content of most larvae may still include material from the diet. Due to the large number of larval families and treatment groups, not all larval groups could be reared during the exact same days, instead, the emergence dates of the larvae from the eggs spread over eight successive days. The larvae from the same family all emerged on the same day. Larvae were reared at 23°C with lights on between 8:00-10:00 am, and between 3:00pm and 5:00pm, at 28°C with lights on between 10:00 am and 3:00 pm, and at 18°C in the dark between 5:00 pm and 8:00 am.
Metabarcoding of the gut microbiota
We surface sterilized three L3 larvae from each of the four treatments for 13 families before individually dissecting their gut out under a microscope in a sterile laminar flow hood. All larval carcasses were preserved to perform the metabolomics analyses described below (see Metabolomics section). We individually extracted the DNA from the gut of the 156 larvae under sterile conditions. The DNA was extracted using a Qiagen DNeasy Blood and Tissue kit (Qiagen, Germany) following an optimized protocol as described by Minard, Tran [66]. Three additional extractions were carried out on sterile water to control for environmental contamination during the procedure. We amplified the hypervariable V5-V6 bacterial region of the rrs gene using the primers 784F (5’-AGGATTAGATACCCTGGTA) and 1061R (5′-CRRCACGAGCTGACGAC) [67], and the LSU region of the ITS gene of Ascomycota fungi using the primers LSU200A-F (5’-AACKGCGAGTGAAGCRG) and LSU476A-R (5’-CSATCACTSTACTTGTKC) [68]. Each sample was amplified in duplicate and using 3µl of the DNA extract for each PCR reaction [66], and the duplicates for each sample were pooled in sterile condition after amplification. Sequencing was performed by the Institute for Molecular Medecine Finland (FIMM, Finland) on a Miseq v.3. Sequencing platform (Illumina, USA) using both reverse and forward primers.
We analyzed the libraries using Mothur v.1.37.6 (http://www.mothur.org/wiki/MiSeq_SOP). [69]. We selected all 250-350bp-long sequences, with less than eight homopolymers, no ambiguous position, and which aligned to the rrs Silva v.123 database. Chimeric sequences were removed using UCHIME implemented in Mothur [70]. Sequences were clustered within Operational Taxonomic Units (OTUs) according to average neighbor method with 3% distance maximum within each OTU. All OTUs showing at least a 10x higher proportion in any given sample than in the negative controls were considered as contaminant and removed from our dataset using an in-house R script [40].
Metabolomic analysis of the larval carcasses
We used the carcasses of the larvae used in the microbiota assays described above for the metabolomics analyses. This allowed us to provide information of the metabolites from the larvae without contamination from the plant diet. After dissection, the three larval carcasses from each family (Ntotal=13 families) were pooled, and crushed in liquid nitrogen using a sterile pestle. Similarly, we also crushed two samples of 30mg of P. lanceolata leaves each and two samples of 30mg of larval frass in liquid nitrogen using sterile pestle, to use as controls for diet metabolite that might still contaminate our larval samples. All 55 samples (52 pooled larval carcasses, one host-plant and two frass controls) were then freeze-dried for 48 hours in a freeze dryer (MechaTech Systems Ltd, UK). Dry samples were weighted on an analytical balance (d=0.1mg, Fisher Scientific, UK). Metabolites were then extracted using the protocol described by Kim et al. [71]. In brief, for each sample, we placed 10mg of dry material in 350µl of CD3OD (VWR Chemicals, Belgium) and 350µl of KH2PO4 (Sigma-Aldrich, Germany) buffer mixed in D2O (pH6) containing 0.05% (wt/wt) of TSP (sodium trimethylsilylpropionic acid) (Sigma-Aldrich, USA). Samples were then sonicated for 20min, and centrifuged 10min at 17,000g. We transferred 600µl of the clear supernatants into individual 5mm diameter NMR tube (Wilmad, USA), and the metabolite content of each sample was analyzed using a Bruker 850MHz Advance III HD NMR spectrometer equipped with a TCI Cryoprobe (Bruker, USA) at the Finnish Biological NMR Center, the University of Helsinki, Finland.
Proton Nuclear Magnetic Resonance (1H-NMR) spectra were acquired at 298K and recorded using 1D pre-saturation pulse sequence (zgpr). For each 1H spectrum, 256 transients were collected into 32K time domain points using a 60° flip angle, spectral width of 10.2kHz, relaxation delay of 5.0s, an acquisition time of 1.6s, and a mixing time of 5ms. Fourier transformation of the free-induction decay was applied with zero filling to give 65K frequency domain data points. The preliminary treatments of the 1H-NMR spectra were performed using the software MNOVA v.10.0.2 (Mestrelab research S.L., Spain). Standard solutions containing 1mg of one of the five antibiotics used for the treatments were measured individually in order to enable their identification, quantification, and trimming from the antibiotic-treated larval samples.
Immune gene expression
We individually sampled up to 12 larvae from each of the 12 families, and flash-froze them in liquid nitrogen once they had reached the third larval instar. For five of the 12 families, we included three L3 larvae for each treatment, while for the remaining seven families, we only had three larvae in at least two of the treatment groups, but only one or two larvae for the other treatments due to mortality during development (Ntotal= 133; Table 1).
The RNA from each larva was individually extracted following a protocol described by Woestmann et al. [58] using TRIzol reagent (Life Technologies Corporation, Carlbad, USA), acid-phenol:chloroform:isoamyl alcohol (25:24:1, pH=5) and Chloroform. The RNA was then precipitated using isopropanol, washed in 75% ethanol, air-dried in a flow hood, and re-suspended in 50µL MQ water. Potential genomic DNA contaminants were removed using DNase I (Thermo Fisher Scientific Inc., UK). The RNA was reverse-transcribed using an iScriptTM cDNA Synthesis Kit (Bio-Rad Laboratories, Hercules, USA) following the manufacturer’s protocol.
The qPCR assays were performed with three technical replicates for each sample, and one negative control and plate control (same sample across all plates) for each 384-well plate used, in a 10µL volume, on a C1000TM Thermal Cycler (Bio-Rad Laboratories, USA). We amplified each of the seven immune genes (lysozyme C, prophenoloxidase, Attacin, peptidoglycan recognition protein LC, ß-1,3-glucan recognition protein, serpin 3a, and pelle) and three housekeeping genes (histone variant H2A.Z, and mitochondrial ribosomal protein L37 and S24) using primers and appropriated protocols described by Woestmann et al. [58]. For each qPCR reaction, we mixed 1μL of the 1/5 diluted cDNA, with 5μL of SYBR® Green containing master mix (iQ™ SYBR® Green Supermix, Bio-Rad Laboratories, USA), 3μL of nuclease-free water, and 0.5μL of the forward and reverse primers (10μm). Non-reverse transcribed samples were used as controls for the lack of genomic DNA contamination.
Statistical analysis
All the statistical analysis were performed with the software R v3.3.1 [72].
Larval performance: We first tested correlations among the variables using linear models (lm), from the package lmer4 [73]. As larval group size at L2 and treatment are highly correlated, only the treatment variable is used in the following models. The development time to L2 was log-transformed prior to analysis. The development time to L2 was compared among larvae from the 27 families using a linear mixed model including the ‘treatment’ as an explanatory variable and the ‘family’ as a random variable. The survival rate at day 13 of the L3 larvae from the 27 families were compared using a general linear mixed model assuming a Gamma distribution of the data, with ‘treatment’ as an explanatory variable and the ‘family’ as a random factor. We used the packages lme4 [73] and MASS [74] for the mixed model analyses. Interclass correlation coefficients (ICC) were calculated based on variance of the random factor and residual to estimate how much of the variance was explained by the random factor ‘family’ in each model.
Microbiota: We used VEGAN [75] in R to compute a Bray-Curtis dissimilarity matrix, and analyze bacterial composition variations among samples using non metric multidimensional scaling (NMDS) or distance based redundancy analysis (dbRDA) [76]. The α-diversity (diversity of the microbiota within each samples) of the microbiota was estimated through the Shannon index while the β-diversity (dissimilarity among samples) of the microbiota was estimated through the Bray-Curtis index. For the α-diversity comparisons, a linear model was used after a logarithmic transformation of the index. The impact of the treatment, the larval family and their interaction were considered as explanatory variables. Similarly, for the β-diversity comparisons, we used a Permutational analysis of variance (adonis-ANOVA) [77] with the treatment, the larval family and their interaction as explanatory variables.
Metabolomics: We extracted entire spectra values from each sample using the program MestReNova 12 (Mestrelab Research, Spain) [78]. For multivariate analysis, the signals were binned to 0.04 ppm, the TSP, H2O and CD3OD signals were removed, and the integral values were transformed following the formula given below: (See Formula 1 in the Supplemental Files)
with ‘m’ as the exact dry mass of the sample (±0.1mg), ‘δ’ as the 1H chemical shift and ‘9’ as to the number of equivalent 1H atoms contained within the TSP reference molecule. Characteristic signals corresponding to α-glucose (δ 4.59, d, J=7.9Hz), β-glucose (δ 5.19, d, J=3.7Hz), Alanine (δ 1.49, d, J=7.2Hz), Formic acid (δ 8.47, s), Acetic acid (δ 1.91, s), Fumaric acid (δ 6.53, s) and Ethanol (δ 1.19, t, J=7Hz) were annotated based on previously published datasets applying the same protocol [71, 79, 80]. Plantago lanceolata also contains variable quantities of iridoid and phenylpropanoid glycosides: namely aucubin, catalpol and verbascoside [43-45]. Two characteristic peaks were identified for aucubin (δ 6.31, dd, J=1.9Hz) and catalpol (δ 6.40, dd, J=1.9Hz) based on the 1H-NMR profiles of standard compounds.
We calculated a first principal component analysis (PCA) including all signal bin values from 49 of the 52 larval samples. Three larval samples (family 1: treatment A and AR, and family 19: treatment R) were removed prior PCA as they clearly appeared as outliers driving most of the variation from the dataset. We considered only seven principal components (PCs) as they showed Eigen-values>3, and together represented over 65% of the observed variation in the dataset. We analyzed the seven PCs using a linear mixed model (lmer) [73] after log transformation. Interclass correlation coefficients (ICC) from each model were calculated based on the proportion of the variance explained by the random family factor. We included ‘larval treatment’ as an explanatory variable and ‘family’ as a random factor in each model. Tukey tests, after correction for the family effect using the glht function, were used as posthoc tests to explore paired comparison between treatments. The resulting P-values were corrected for multiple testing using a Bonferroni correction (a=0.025).
We tested the relationship between each of the first seven PCs and the development time to L2 and survival to L3 between the treatment groups using a linear mixed model (lmer), with ‘treatment’ and the ‘PC’ of interest as fixed factors, and ‘family’ as a random factor to each model. Tukey tests after correction for the family effect with the glht function, were used as posthoc tests to explore paired comparison between treatments. The resulting P-values were corrected for multiple testing using a Bonferroni correction (a=0.025).
Immune gene expression: We calculated the mean immune gene expression from the three technical replicates (with exception of few outliers) considering the geometric mean of the three reference genes, for each family but one. For unknown reason, the control samples for the S24 house-keeping gene was not expressed for the family16, thus, the immune gene expression for family16 was calculated based on the geometric mean of the two remaining house-keeping genes only. The immune genes expression (Log2) were compared among larvae from 12 families using generalized linear models [73, 77] including the ‘larval treatment’ and ‘gene’ as fixed factors (including interaction term), and ‘family’ as a random factor. We performed a post-hoc analysis using the lsmeans function with Tukey’s HSD adjustment for pairwise comparisons [81], to explore paired comparison between treatments and genes, and corrected resulting P-values for multiple testing using a Bonferroni correction (a=0.025).
Finally, we tested whether development time to L2 (corrected for family effect) and survival to L3 (corrected for family effect) were differently affected by variation in the expression levels (Log2) of the Attacin immune gene (an antimicrobial peptide active against Gram-negative bacteria [82]) from the different larval treatment groups, including the ‘immune genes expression levels’ and ‘treatment’ as explanatory variables. Tukey tests, after correction for the family effect using the glht function, were used as posthoc tests to explore paired comparison between treatments. The resulting P-values were corrected for multiple testing using a Bonferroni correction (a=0.025).