Intergenerational effects, also known as parental effects in which the offspring phenotype is influenced by the parental phenotype, can occur in response to factors that occur not only in early but also in late parental life. However, little is known about how these parental life stage-specific environments interact with each other and with the offspring environment to influence offspring phenotypes, particularly in organisms that realize distinct niches across ontogeny. We examined the effects of parental larval starvation and adult reproductive environment on offspring traits under matching or mismatching offspring larval starvation conditions using the holometabolous, haplodiploid insect Athalia rosae (turnip sawfly). We show that parental larval starvation had trait-dependent intergenerational effects on both life-history and consumption traits of offspring larvae, partly in interaction with offspring conditions, while there was no significant effect of parental adult reproductive environment. In addition, while offspring larval starvation led to numerous gene- and pathway-level expression differences, parental larval starvation impacted far fewer genes and only the ribosomal pathway. Our findings reveal that parental starvation evokes complex intergenerational effects on offspring life-history traits, consumption patterns as well as gene expression, although the effects are less pronounced than those of offspring starvation.

Set-Up of Insect Rearing and Plant Cultivation

Adults of A. rosae (F0) were collected in May 2019 at two locations (population A: 52°02'48.0"N 8°29'17.7"E, population B: 52°03'54.9"N 8°32'22.2"E). These individuals were reared for further two generations for details of breeding design see S1). White mustard (Sinapis alba) was provided for oviposition and Chinese cabbage (Brassica rapa var. pekinensis) as food plant. Adults of the F2 generation were kept individually in Petri dishes and provided with a honey:water mixture (1:50). Athalia rosae is haplodiploid, i.e. virgin females produce male offspring (Naito and Suzuki 1991). Thus, mated as well as virgin females were placed individually into boxes (25 x 15 x 10 cm) to increase the likelihood of gaining similar numbers of females and males. Females were supplied with middle-aged leaves of non-flowering cabbage plants for oviposition and a honey:water mixture, which were replenished daily. Females were removed from the boxes after one week and their offspring used to set up the experimental generations, which were kept in a climate chamber (20 °C:16 °C, 16 h: 8 h light:dark, 70% r.h.).

Plants of S. alba and B. rapa were grown from seeds (Kiepenkerl, Bruno Nebelung GmbH, Konken, Germany) in a greenhouse (20 °C, 16 h: 8 h light:dark, 70% r.h.) and a climate chamber (20 °C, 16 h: 8 h light:dark, 70% r.h.). Plants of Ajuga reptans used for clerodanoid supply were grown from seeds (RHS Enterprise Ltd, London, UK) in the greenhouse and transferred outside in late spring when about 2 months old. Middle-aged leaves of 8 months old plants were offered to adults.

Experimental Overview and Measurements of Life-History Traits

We conducted a fully factorial experiment where we manipulated the parental larval environment (parental larval starvation), the parental reproductive environment (exposure to clerodanoids) and the offspring environment (offspring larval starvation) (Fig. 1). We thus tested the effects of parental and offspring larval starvation, and parental adult reproductive environment on the offspring phenotype. On the day of hatching, larvae of the parental generation were individually placed into Petri dishes (5.5 cm diameter) with moistened filter paper and B. rapa leaf discs cut from 7-10 week old plants, which were replaced daily. Per maternal line, larvae were split equally between one of two larval starvation treatments, no starvation (N) or starvation (S). For starvation, individuals were starved twice for 24 h, the day after moulting into 2^nd instar and the day of moulting into 4^th instar, thereby mimicking the food deprivation larvae may experience when occurring in high densities (Riggert 1939). Larval instars were tracked by checking daily for the presence of exuviae.

Eonymphs were placed in soil for pupation. Adults were kept individually in Petri dishes and provided with honey:water mixture. From the parental generation, pairs of non-sib females and males reared under the same larval starvation treatments were assigned to one of two reproductive environment treatments, where both parents either had clerodanoids (C+) or did not (C-). C+ individuals were exposed to a leaf section (1 cm²) of A. reptans for 48 h prior to mating, giving individuals time to take up clerodanoids (preprint Paul et al. 2021). Mated females (2-9 days old) from each of these four treatments (NC-, NC+, SC-, SC+) and C- virgin females (NC-, SC-) were then placed in individual breeding boxes. Their offspring were distributed to offspring starvation treatments (N or S) that either matched the parental starvation treatment or differed from it (mismatch) (Fig. 1). In both generations, larval, pupal, and total development time (from larva to adult) as well as the adult body mass at the day of emergence (Sartorius AZ64, M-POWER Series Analytical Balance, Göttingen, Germany) were recorded for each individual. In total 358 larvae of the parental generation (177 females, 181 males) and 484 larvae of the offspring generation (282 females, 202 males) reached adulthood, out of the 607 and 688 larvae, respectively, that were initially reared in each generation.

Consumption Assays

To test effects of parental and offspring larval starvation experience on offspring consumption, assays were performed with larvae of the offspring generation at the start of the 3^rd instar (directly after the first starvation event), measuring the relative growth rate, relative consumption rate, and efficiency of conversion of ingested food. Each larva was weighed at the beginning of the consumption assays (= initial body mass) (ME36S, accuracy 0.001 mg; Sartorius, Göttingen, Germany) and provided with four discs of middle-aged cabbage leaves (surface area of 230.87 mm² per disc) on moistened filter paper. After 24 hours, larvae were weighed again (= final body mass) and the leaf disc remains scanned (Samsung SAMS M3375FD, resolution 640 x 480). In that period, leaf discs showed no signs of wilting. The total area of leaf consumed (mm²) was then calculated as the difference between the average leaf area and the remaining leaf area.

Statistical Analyses

All data were analyzed using R 4.0.2 (2020-06-22). We set α = 0.05 for all tests and checked model residuals for normality and variance homogeneity: All All linear mixed effects models (LMMs) and generalized linear mixed models (GLMM) were run in lme4 using maximum likelihood. Stepwise backwards deletion using X² ratio tests (package:MASS; version 7.3-53.1) for the life history traits and, due to the much smaller sample sizes (Luke 2017), conditional F-tests with df correction using Satterthwaite method (package lmerTest; version 3.1-3; Kuznetsova et al. 2017) for the consumption analyses were employed to reach the minimum adequate models (Crawley 2012). Posthoc analyses were carried out using the package ‘multcomp’ (version 1.4-13; Hothorn et al. 2008). Post data entry, raw data were visually inspected thrice, all variables plotted and outliers and possible anomalies in the data (e.g. strings of similar values) interrogated (package:pointblank; version 0.6.0). Intragenerational effects of starvation on life-history traits of individuals of the parental generation were tested as described in S1.

In A. rosae, usually 6 instars for female and 5 instars for males are found (Sawa et al. 1989). Due to observations made during the experiment (no a priori hypothesis), we tested in the offspring generation whether the likelihood of an additional larval instar (7 for females and 6 for males) differed based on offspring larval starvation treatment (independent of parental treatments) using a binomial GLMM (package: lme4), where the predictor was offspring starvation treatment and parental pair was included as a random effect. The effects of the parental larval starvation treatment, parental clerodanoid exposure, offspring starvation treatment and their interaction on larval, pupal and total developmental time as well as adult mass of offspring individuals were assessed in separate LMMs, with parental pair included as a random effect (controlling for non-independence of sibling larvae). Female and male data were analyzed separately to enable model convergence.

Relative growth rate, relative consumption rate, and food conversion efficiency of larvae of the offspring generation were analyzed using LMMs. We excluded parental clerodanoid exposure as a predictor variable from the consumption assay analysis due to the low number of individuals in certain treatments (Fig. 1) and analyzed male and female data separately. To assess variation in relative growth rate, the change in larval body mass [final mass — initial mass] was used as the response variable and initial larval body mass, parental starvation treatment, offspring starvation treatment, and their interactions were used as predictors. To assess relative consumption rate, we used the total area of consumed leaf material as the response variable and initial larval body mass, parental starvation treatment, offspring starvation treatment, and their interactions as the predictors. Finally, for food conversion efficiency the change in larval body mass was taken as the response variable and the predictor variables were total area of consumed leaf material, parental starvation treatment, offspring starvation treatment, and their interactions. Parental pair was included in all three models as a random effect.

Sample Collection and Sequencing

A total of 24 male larvae (4^th instar, 9 d old), comprising six biological replicates per treatment level (parent/offspring treatment: N/N, N/S, S/N, S/S, all C-, Fig. 2), were collected to investigate the effects of parental and offspring starvation treatment on gene expression in offspring larvae. Individuals were chosen in a way that maximized the equal spread of siblings across treatments and frozen at -80 °C prior to extraction. RNA was extracted with an Invitrogen PureLink™ RNA Mini Kit (ThermoFisher Scientific, Germany), including a DNase step (innuPREP DNase I Kit, analyticJena, Jena, Germany). RNA quality was assessed on a bioanalyzer 2100 (Agilent, CA, United States) and Xpose (PLT SCIENTIFIC SDN. BHD, Malaysia). Library preparation (Ribo-Zero for rRNA removal) and sequencing (NovaSeq6000 and S4 Flowcells, Illumina, CA, United States) were provided by Novogene (Cambridge, UK). Sequence quality before and after trimming was assessed using FastQC (v. 0.11.9; Andrews 2010). After short (< 75 bp), low quality (Q < 4, 25 bp sliding window), and adapter sequences were removed using Trimmomatic (Bolger et al. 2014), more than 98% of reads remained.

Differential Expression Analysis

Cleaned reads were mapped to the annotated genome of A. rosae, version AROS v.2.0 (GCA_000344095.2) with RSEM v1.3.1 (Li and Dewey 2011), implementing mapping with STAR v2.7.1a (Dobin et al. 2013). Analysis of differential gene expression was conducted with DESeq2 (version 1.28.1; Love et al. 2014). The results of mapping with RSEM were passed to DESeq2 for gene-level analyses using Tximport (version 1.16.1; Soneson et al. 2015). Prior to analysis, genes with zero counts in all samples and those with low counts (< 10) in less than a quarter of samples were excluded. Model fitting was assessed by plotting dispersion estimates of individual gene models and outlier samples were inspected using principle component analysis of all expressed genes and pairwise-distance matrices between samples. Expression was modelled based on the entire dataset, with the four levels representing the combination of parental.offspring starvation treatments (N.N, N.S, S.N, and S.S). Differential expression was assessed in four pairwise comparisons between the treatments: 1) N.N vs N.S and 2) S.N vs S.S were used to assess the effects of offspring starvation treatment for individuals whose parents experienced the same starvation treatment, whereas 3) N.N vs S.N and 4) N.S vs S.S were used to assess the effects of differing parental starvation treatment on individuals that experienced the same offspring starvation treatment. Significance was based on a Wald Test and shrunken log fold-change values with apeglm (version 1.10.1; Zhu et al. 2019). P-values were adjusted for multiple testing using a Benjamini-Hochberg procedure (Benjamini and Hochberg 1995); we retained significantly differentially expressed genes (relatively up- or downregulated) if this adjusted P-value was < 0.05. A Venn diagram (Venn.Diagram, version 1.6.20) was used to depict the overlap of significantly differentially expressed genes for each pairwise comparison. The expression of the significantly differentially expressed genes of each comparison was visualized in a heatmap using normalized counts scaled per gene (scale=“row”) (pheatmap, version 1.0.12).

To examine differential expression of genes with known roles in stress and/or starvation response, we searched specifically for target significantly differentially expressed genes, based on their identity (from the genome annotations), using the keywords “heat shock protein” (and “hsp”), “cytochrome P450”, “octopamine” and “tyramine” in the putative gene names.

Pathway-Level Analysis of Differential Expression

We used the KEGG (Kyoto Encyclopedia of Genes and Genomes) database to assign the predicted genes in the A. rosae genome to gene pathways with the KEGG Automatic Annotation Server (KAAS; Moriya et al. 2007) (S3). Of the annotated genes for which we had read counts, 65% were assigned to at least one KEGG pathway. A gene set enrichment analysis was then performed on the normalized counts using GAGE (Luo et al. 2009) and pathview (Luo and Brouwer 2013) in R, applying a false discovery rate-adjusted P-value cut-off of < 0.05 to identify differentially expressed pathways. We derived the non-redundant significant gene set lists, meaning those that did not overlap heavily in their core genes, using esset.grp and a P-value cut-off for the overlap between gene sets of 10e^-10.

Unmapped reads were inspected to check for the differential expression of genes not present in the reference genome. These were extracted from the RSEM-produced BAM file using samtools view (Li et al. 2009) and converted to fastq (using bamtools bamtofastq; Barnett et al. 2011), and de novo assembled using Trinity (default settings) (Grabherr et al. 2011). This assembly was done jointly with reads from additional samples (preprint Paul et al. 2021) to optimize coverage. The unmapped reads were mapped back to this reference and expected read counts extracted using eXpress (Roberts et al. 2011). Transcripts with >10,000 mapped reads were identified using BLASTN (default settings, limited to one match per gene and 1e^-1) against the NCBI nt nucleotide database. Differential expression of unmapped reads in the larval samples was carried out as for mapped reads (see above).