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.