2.Material and Methods
2.1. Host plants
This study is based on two switch experiments that differed in the host plants involved. The host plants for the first experiments (Switch A) comprised Urtica dioica (stinging nettle), Salix caprea(sallow) and Ulmus glabra (elm). For the second switch experiment (Switch B) Urtica dioica , Salix caprea and Ribes uva-crispa (gooseberry) were used. Leaf material of the host plants was collected in the close surrounding area of Stockholm University, Stockholm, Sweden (N 59° 21’ 47”, E 18° 03’ 38”).
2.2. Butterfly rearing and host switch
Mated females of Polygonia c-album were collected in Åkersberga (N 59° 28’ 56”, E 18° 21’ 18”) and Hansta (N 59° 26’ 11”, E 17° 54’ 14”) near Stockholm. For oviposition, females were individually placed into separate cages (50x50x40cm) with a transparent roof and a net in the front. The cages contained a cut Urtica branch and sugar water for feeding. Moistened paper towels were distributed at the floor of each oviposition cage. All cages were placed below a light source (LD 8:16). Eggs were collected every morning.
Immediately after hatching, larvae were placed into separate cups (volume: 1.25 oz) containing one leaf of either Urtica, Salix orUlmus (Switch A), or Urtica, Salix or Ribes (Switch B) respectively. To prevent the leaves from drying out quickly, wet cotton was wrapped around their petioles. The quality of the leaves was checked regularly. Leaf material that was too dry, or of too small an amount was replaced. Larvae were kept at 18C (L:D 18:6) during their development.
One day after molting into the fourth instar, larvae were switched to another plant (Figure 1). For this, larvae were first allowed to feed for one hour on the original host (Host 1) – to ensure larval feeding this was preceded by a 6-hour fasting phase. After feeding on Host 1, larvae had to fasten again for one more hour before being transferred to the new plant (Host 2; Figure 1b). Host 2 could either be the same plant species (Control) or one of the alternative host plants (Switch). This resulted in a total of 9 different host combinations per switch experiment.
To measure the transcriptional response to these host switches, the midguts of the larvae were dissected under phosphate saline solution (pH 7.4, 10 mM) 2 or 17 hours after they initiated feeding on the second host (Host2). The gut was emptied and Malpighian tubules were removed. Dissected gut tissue was snap frozen using liquid nitrogen and stored at -80C. Sampling gut tissue at two different time points allowed us to assess how quickly caterpillars can adjust to their new host environment. Moreover, it also enabled a further distinction between transcriptional changes due to the switch itself (i.e. possible stress responses) and “true” host plant effects. In total, midguts were taken from 72 caterpillars (1 larva x 4 families x 9 different host combinations x 2 different time points per treatment). In the same way, the second switch experiment was conducted using Urtica ,Salix and Ribes . Here, midgut tissue from 72 caterpillars (2 larvae x 2 families x 9 different host combinations x 2 different time points per treatment) was collected for the differential gene expression analysis. In both experiments, one sample had to be discarded due to low quality.
The effect of host switches on the larval performance was investigated, but only for Switch B. For this, 1-2 caterpillars from 5 different families were reared at 18C (LD: 12:12) on either Urtica ,Salix or Ribes (n=72). Larvae were moved to a new plant (same or alternative host species) upon reaching the fifth instar. Their weight was measured daily until pupation. The daily growth rate was calculated as:
Growth rate =\(\frac{\log(\text{mass})}{\text{developmental}\ \text{time}}\)
2.3.RNA extraction
The frozen gut tissue was first transferred into 300µl Tri-reagent (ZymoResearch) and homogenized with a clean pistil. The tissue continued to be homogenized after adding 700µl of additional Tri-reagent and samples were incubated for 5 minutes at room temperature. After the incubation, 100µl of bromo-chloro propane (BCP) (Sigma-Aldrich) per milliliter Tri-reagent were added. The samples were mixed thoroughly and incubated for another 10 minutes. Subsequently, the samples were centrifuged at 12,000g for 13 minutes at 4C. The aqueous supernatant was then transferred into a fresh tube and 500µl isopropanol (Sigma-Aldrich) per milliliter Tri-reagent were added. Samples were mixed and incubated for 8 minutes at room temperature. After centrifuging for 8 minutes at room temperature and 12,000g, the supernatant was discarded. 1 milliliter of 75% Ethanol (KiiltoClean) was added per milliliter Tri-reagent and samples were centrifuged at 7500g for 5 minutes at room temperature. The ethanol was fully removed and the remaining pellets were air dried. Dry pellets were subsequently dissolved in 90µl pre-heated (55C) RNA-storage solution (Invitrogen). DNA was removed using the DNA-free ™ DNA removal Kit (Invitrogen) following the manufacturer’s protocol. The quality of the extracts was assessed using Bioanalyzer (Agilent RNA 6000 Nano). RNA concentration was quantified using a Qubit 2.0 fluorometer (Invitrogen). Library preparation (with Poly-A-selection) and sequencing was done at the NGI SciLifelab in Stockholm. RNA libraries were sequenced with 2x101bp paired-end sequencing on an Illumina TrueSeq platform.
2.4. Analysis
2.4.1. Read alignment and quantification
The quality of the raw sequence data was checked using FastQC (version 0.11.9; Andrews 2010) and MultiQC (version 1.10; Ewels et a l. 2016). Adapter and quality trimming was performed using Trimmomatic (version 0.36; -phred 33, SLIDINGWINDOW:4:5, LEADING:5, TRAILING:5, MINLEN25; Bolger et al. 2014). The trimmed sequences were again checked with FastQC and MultiQC. Trimmed reads were mapped against an inhouse created reference genome of Polygonia c-album(Celorio-Mancera et al. 2021) using a STAR 2-pass approach (version 2.7.2b; Dobin et al. 2013). Mapping success was checked with samtools (version 1.12; Danecek et al. 2021) and MultiQC. Abundance estimations were done using featureCounts from Subread (version 2.0.0; Liao et al. 2014).
2.4.2. Differential gene expression
Differential gene expression analysis was performed using the R-package edgeR (version 3.40.0; Robinson et al. 2010). For this, the count matrix including 21,761 genes was first split into the two time points (2 vs.17 hours) and separate gene expression lists were created. Filtering was done using the filterByExpression-function of edgeR (default setting) with information about the host combination as a group-factor, which retained 13,126 (Switch A, 2hrs) 12,811 (Switch A, 17hrs), 12,622 (Switch B, 2hrs) and 12,910 (Switch B, 17hrs) genes. Filtered data were then normalized using CalcNormFactors. Differential gene expression was analysed using QLFit (robust=TRUE) and glmQLFTest (FDR <0.05). Here, an ANOVA-like framework was used to investigate the main effect of Host 1 and Host2, as well as their interaction. We further used pairwise comparisons to assess the differences between treatments. To account for potential differences due to relatedness between samples, family was included in the models. Significantly differentially expressed genes were identified as those with an adjusted P-value (FDR, Benjamini-Hochberg correction) less than 0.05.
We generated a functional annotation using the web-based version of eggNOG-emapper (Huerta-Cepas et al. 2019; Cantalapiedra et al. 2021), which relies on HIMMER (Eddy 2011), DIAMOND (Buchfinket al. 2021) and MMSEQS2 (Steinegger and Söding 2017) when searching for sequence alignments. Briefly, genes were converted to amino acid sequences using gffread (Pertea and Pertea 2020) and submitted to eggNOG emapper using default parameters, including automatic identification of taxonomic scope (i.e., functional annotations were not limited to those confirmed in arthropods). We obtained functional annotations for 12,469 of 17,606 genes in theP. c-album annotation, however only 8034 of these genes were annotated with gene ontology (GO) terms. Enrichment of GO terms among genes that differed significantly in their expression between the treatments was tested using topGo (Alexa and Rahnenfuhrer 2016). Significant enrichment was tested using Fisher’s exact tests corrected using topGO’s parent-child algorithm. This algorithm evaluates enrichment of a term in relation to its parent terms and so reduces false positives due to the “inheritance” problem – where more specific terms are likely to be falsely enriched if their parent terms are enriched (Grossmann et al. 2007). Here we focused on enrichment of biological process terms and only considered GO-terms that were annotated to at least five genes.
2.5. Congruence with earlier results
In order to investigate the repeatability of our observations, and whether host-specific gene expression is primarily a direct response to the host or a down-stream effect, we compared the set of differentially expressed genes per contrast from this study to the sets of co-expressed genes reported under similar host plant treatments but without a switch (Celorio-Mancera et al. 2023). In order to do this, it was necessary to join the gene sets from two different annotation databases. For this purpose, we used DIAMOND (blastp flag, default parameters, 16930 unique queries aligned) to annotate the predicted protein set from the P. c-album genome in the present study, which relied on an ortholog Heliconius melpomene database generated and described by Celorio-Mancera et al. (2023). The output table of annotated sets of genes from P. c-album was used in turn to name the lists of differentially expressed genes obtained from the different contrasts.
2.6 Performance
We tested the effect of host plant switches on larval performance using linear mixed models (lme4, version 1.1.31, Bates et al. 2015), including ID of the mothers and individual-ID as random effects. Data wrangling, analysis and visualization were performed in R with the support of tidyverse packages (ggplot2, version 3.4.0, Wickham 2016).