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).