Discussion
Diapause is a complex, plastic phenotype that requires the sequential
coordination of many physiological processes, making miRNAs interesting
candidates due to their regulation of mRNA expression patterns by
targeting many mRNAs at once. Here, we identified miRNAs across the
genome of P. napi, and found that there are a diverse range of
miRNA expression patterns that occur throughout diapause. Ultimately, we
identified several miRNAs that may play an important role in diapause
timing, notably miR-14-5p, which could have a key role in diapause
termination by targeting genes in the ecdysone synthesis pathway.
We identified a number of miRNAs in P. napi, including 129 novel
miRNAs that we describe for the first time here. These genes were
located throughout the genome, with 32% being located within intron of
genes, which is lower than the average of 52% across 80 species
(Guerra-Assunção & Enright, 2012), but higher than the 12% reported in
zebrafish (Thatcher, Bond, Paydar, & Patton, 2008). miRNAs that are
located in introns can act as functional miRNAs from the spliced-out
intron, although these function differently than canonical miRNAs and
are hard to detect using normal miRNA detection tools (Rorbach, Unold,
& Konopka, 2018). Nevertheless, finding these intronic miRNAs suggests
that there is a substantial potential for the expressed gene harboring
the intronic miRNA to be functionally related to the miRNA targets
(Olena & Patton, 2010).
There is also a high concentration of miRNA on the Z chromosome compared
to the autosomes, which is notable since this is the sex chromosome in
Lepidoptera. The high concentration on the Z chromosome appears to be
driven by several known clusters of miRNAs. These clusters include the
miR-2 cluster, of which we identified 4 members of the miR-2 family
(miR-2a, miR-13a, miR-13b, and miR-2c), along with miR-71, which flanks
the miR-2 family in many other insects, and is heavily involved in
neural development and maintenance (Marco, Hooks, & Griffiths-Jones,
2012). Finally, we note that our miRNA sequencing of two body parts
composed of diverse tissues, across several developmental time points,
provides a rather deep sampling of the miRNA, at least for these body
parts. Nevertheless, the numbers of miRNAs we have identified are
roughly comparable with similar efforts in other Lepidopteran species
(Kozomara et al., 2019; Quah, Hui, & Holland, 2015; Reynolds et al.,
2017), suggesting that the number of identified miRNAs is probably quite
representative of detectable miRNAs in this butterfly.
We next quantified miRNA expression patterns across tissues and
developmental trajectories, finding that their general expression
patterns closely match expression patterns found in previous metabolomic
and transcriptomic (mRNA and alternative splicing expression)
assessments (Lehmann et al., 2018; Pruisscher et al., 2022; Steward et
al., in review ). Using multidimensional scaling (MDS) for a
global assessment of the head tissues, there is a clear alternative
trajectory mostly in dimension 1 away from the direct developing
individuals until the latest stages of diapause. The main axis of
variation in total expression accounts for diapause, and has a larger
impact on miRNA expression than development, where 40% of variance is
accounted for, likely due to the high activity of neuroendocrine
regulation occurring throughout diapause, which is highly pleiotropic
(Williams, 1952). Despite the overall lower expression of miRNAs
patterns in abdomen tissue (there is a 7-fold higher overall expression
in the head tissue), these patterns were nearly identical.
In an initial analysis containing all targets of a miRNA expression
cluster, the largest GO terms were often the most enriched sets, but by
reducing targets to include only genes targeted by multiple miRNAs, we
were able to find much more specific and informative GO terms with
potentially interesting biological roles in diapause. Given our interest
in the termination of diapause, we were primarily interested in clusters
with peaks between 6 and 114 days in diapause, which includes H4, H6, A3
and A4. Cluster H2 was enriched for GO terms related to lipid
metabolism, since the expression of this cluster decreases in the first
few days of pupation, this could indicate that there is a suppression of
lipid metabolism upon entry in diapause, and a cessation of miRNA
regulation as pupae shift to heavily relying on lipid metabolism in
diapause. However, previous lipidomic studies of these timepoints inP. napi shows little decrease in storage lipids through diapause,
and only displays a significant decrease after diapause has terminated
(Lehmann et al., 2016). Even though there was no detectable lipid store
depletion found, there is a decrease in the respiratory quotient to
below 0.7 ten days after pupation, indicating they are primarily
metabolizing lipid, matching the miRNA expression patterns of cluster H2
(Lehmann et al., 2018). Alternatively, it could be related to the
extensive remodeling of both cell membrane lipid and lipid pool
composition that can occur in response to cold in insects (Enriquez &
Teets, 2023; V. r. Koštál, Berková, & Šimek, 2003).Cluster H4 is
enriched for both regulation of Wnt signaling and imaginal disc pattern
formation. Wnt signaling has been identified as an important regulator
of diapause termination (V. Koštál et al., 2017; Ragland, Egan, Feder,
Berlocher, & Hahn, 2011), and expression in cluster H4 peaks after
diapause is terminated. Since Wnt signaling has been found to have
higher expression just before diapause is terminated, perhaps the
expression of cluster H4 after diapause is terminated is a signal
related to cessation of the need for an upregulation of Wnt signaling as
a signal to terminate diapause (V. Koštál et al., 2017). As for imaginal
discs patterning, it is a critical component of pupal development, and
in some species the lack of imaginal disc development is used as an
indicator of diapause (V. Koštál, Šimůnková, Kobelková, & Shimada,
2009). The upregulation of cluster H4 after termination therefor may
reflect a difference between diapause and post-diapause quiescence, in
that during quiescence the pupae is able to develop if conditions allow,
but diapausing pupae cannot (V. Koštál, 2006; Lehmann et al., 2017;
Süess et al., 2022). Cluster H6 is enriched for the GO term
“glucocorticoid metabolic processes”, which notably is a class of
steroid hormones that are involved in mammalian stress response, and
absent in insects. Insects do respond to cortisol, which is a
glucocorticoid, but cannot produce it, instead insects use ecdysteroids
in a similar manner as mammals use glucocorticoids (Gawienowski,
Kessler, Tan, & Yin, 1987; Ivanovic, 2018). Regardless, the enrichment
of a process that is absent in the sequenced taxa warrants further
study, particularly due to the importance of hormonal signaling in the
regulation of diapause timing.
In order to gain physiological informative insights, we further parsed
these results to only focus upon miRNAs that were differentially
expressed across termination time points and within the aforementioned
clusters, in order to gain biological insight via first principals
complexity reduction. This narrowing of focus identified interesting
candidates playing a role in diapause regulation through their
interaction with the PTTH/20E axis, a main candidate being miR-14-5p.
The PTTH/20E axis is a critical regulator of determining the timing of
transition from diapause to post-diapause quiescence in P. napi ,
where it is hypothesized that the presence of PTTH starts a cascade of
ecdysone release, which leads to termination of diapause (Süess et al.,
2022). miR-14, both 3p and 5p have been found to delay larval
development by decreasing ecdysone signaling in Bombyx mori , by
targeting both ecdysone receptors (Liu et al., 2018) and targeting genes
in the ecdysone-signaling pathway (He et al., 2019). In the beetleGaleruca daurica there is further evidence that both miR-14-3p
and miR-14-5p impact the production of ecdysone, with neither miRNA
expressed in response to ecdysone, as when larvae were exposed to
ecdysone these miR-14s were not among the miRNAs that were
differentially expressed (Jin et al., 2020). Our target identification
analysis identified a miR-14-5p binding site in the 3’-UTR of the genephantom, which is in the ecdysone-synthesis pathway. The
ecdysone-synthesis pathway converts cholesterol to 20E, meaning that
downregulating a gene in the pathway would likely delay diapause
termination. In B. mori miR-14-5p targets two ecdysone signaling
genes, neither of which are phantom , in addition to several
ecdysone response genes (He et al., 2019). While our results do not
fully align with those in B. mori , it is encouraging to find the
same miRNA targeting the same pathway. The expression of miR-14-5p is
highest in the deepest phase of diapause, and is down-regulated as
diapause termination occurs, matching the timing of our expectation of
ecdysone presence. This provides a putative mechanism for a miRNA to
have a large role in the timing of a critical transition in diapause.
In P. napi, the complementary sequence from miR-14-5p, miR-14-3p
was found to target a receptor of ecdysone, but is only expressed early
in diapause, which is a period when pupae are not sensitive to ecdysone
to terminate diapause (Süess et al., 2022). In B. mori miR-14-3p
targets the mRNA of two 20E receptors, ecdysone receptor andecdysone-induced protein 75 (E75) (Liu et al., 2018), while we
found that miR-14-3p in P. napi targets E75 . However,
miR-14-3p is not differentially expressed during diapause termination,
instead it is downregulated in the first three days of pupation. ThatP. napi is insensitive to ecdysone during the first few weeks of
diapause (Süess et al., 2022), could therefore partially be due to a
lack of available ecdysone receptors (a dimer consisting of the two
proteins ecdysone receptor (EcR) and ultraspiracle (USP) (Rinehart,
Cikra-Ireland, Flannagan, & Denlinger, 2001)). In sum, it is notable
that we are finding connections between the miRNAs, their phenotypically
relevant expression, their identified targets in our focal species, and
previous literature making similar connections.
The other interesting miRNA that is differentially expressed in diapause
termination is miR-2a-3p, a miRNA that has been found to be responsive
to ecdysone. miR-2a-3p is involved in metamorphosis by targeting the
juvenile hormone response gene krüppel-homolog 1 (Lozano, Montañez, &
Belles, 2015). In B. mori , miR-2a-3p is more highly expressed in
embryogenic cells that are treated with 20E than in unexposed
embryogenic cells (Jin et al., 2020). Being expressed in response to 20E
in P. napi would both match our expectation based on the timing
of ecdysone sensitivity in diapause, and support that miR-14-5p is
involved in the regulation 20E synthesis, due to the expression patterns
of miR-2a-3p being inverse of miR-14-5p. However, in a hemimetabolous
insect, Nilaparvata lugens , miR-2a-3p is transcriptionally
repressed by an early 20E response gene, broad-complex (BR-C)(Chen, Liang, Liang, Pang, & Zhang, 2013), which if that were the case
in P. napi , would lead to an inverse expression pattern instead
of the one we observed. We cannot disentangle the role of miR-2a-3p in
diapause ecdysone signaling in P. napi from this data, but it
provides an interesting candidate for a miRNA involved in diapause
regulation and a target for future research.
When we compare differentially expressed miRNAs identified in previous
diapause studies, we find an overlap with several miRNA, but many more
are missing. Of the 16 miRNA identified in the literature the only
overlap with DE miRNAs we identified in diapause termination is
miR-14-5p. miR-14-5p is found to be downregulated in diapausing pharate
larvae of the mosquito Aedes albopictus (Batz et al., 2017),
which would only match our findings as long as the diapause timepoint is
taken after day114, when it is lower than day 6 of direct developing
individuals. We also see patterns in miRNAs that match existing
literature in expression direction, but only if taken at the right
timepoints. miR-92-3p was found to be downregulated in diapausing pupae
of the flesh fly Sarcophaga bullata , which match our findings,
but when looking at the time course of expression there is a large
upregulation at day 144. Additionally, the expression of miR-275-3p
matches expression direction in flesh flies, but this depends on the
time of sampling the direct developing individual (Figure 7). Overall,
the extent that the results from this study overlap with existing work
depends on the timepoint used for comparison. These discrepancies of
expression patterns at specific times during diapause highlight the
dynamic nature of diapause and the importance of sampling a time course
to identify patterns that would be otherwise missed.