The tobacco hawkmoth Manduca sexta was reared and staged as described previously (Schwartz and Truman 1983). The ISMs are a very discrete cells and can be cleanly isolated without contamination from other tissues. The lateral intersegmental muscles (ISMs) were dissected under ice-cold saline, flash frozen on dry ice, and stored in liquid nitrogen until used for RNA isolation.

Some animals were injected on day 17 of pupal-adult development with 25 µg of 20-hydroxyecdysone (20E) (Sigma) in 10% isopropanol to delay ISM death (Schwartz and Truman 1984) and then the ISMs removed prior to the normal time of eclosion on day 18.

The ISMs of three to four animals per developmental time point (eight time points in total: days 13, 14, 15, 16, 17, 18, and 1-h post-eclosion (PE); plus 20-hydroxyecdysone injection on day 17: 20E) were homogenized and total RNA was isolated using a mirVana RNA Isolation kit (Life Technologies).

For small RNA-seq library construction, 50 µg of total RNA was fractionated by 15% urea polyacrylamide gel electrophoresis and the 18–30 nt fraction extracted for library construction. (While some longer piRNAs may have been missed with this size cut off, it prevented contamination of the library by a ∼32 nt 2S ribosomal RNA that is present in some insect species (Tautz et al., 1988)). 3′ and 5′ adaptors were ligated to the small RNA and the cDNA reverse transcribed and PCR amplified. The libraries were purified by polyacrylamide gel electrophoresis, and subjected to 50 nt single-end sequencing on an Illumina HiSeq™ 2000 (San Diego, CA, United States) by Beijing Genomics Institute (Hong Kong).

For RNA-seq, directional and random primed cDNA libraries were constructed with poly(A)+ RNA, analyzed with a Bioanalyzer (Agilent Technologies; Santa Clara, CA, United States) and 50 nt single-end sequencing was performed as above. For each of the eight time points, we prepared three biological replicates of RNA-seq libraries.

All the sequencing libraries are accessible from Gene Expression Omnibus (GEO) (accession number GSE80830).

piRNAs are 2′-O-methylated at their 3′ termini, which renders them resistant to oxidation (Vagin et al., 2006). In contrast, while oxidation does not destroy miRNAs or mRNA degradation products, it does render them unclonable (Zhang et al., 2014). Therefore, we oxidized small RNAs as outlined in (Gunawardane et al., 2007), then cloned the resulting piRNAs as described above. We sequenced the oxidized small RNA library at the Deep Sequencing Core Facility at the UMass Medical School.

We downloaded the genomic assembly of Manduca sexta (Msex1.0) and the transcript and the protein sequences (revised-OGS-June2012) from the Manduca Base (http://agripestbase.org/manduca/) (Tsuji et al., 2020). The genomic assembly contains 20,868 scaffolds, with a median length of 994 bp. Unfortunately, the annotated Manduca genome available lacks the resolution to identify piRNA clusters. We annotated protein-coding genes, transposable elements, low complexity regions, miRNAs, and other non-coding RNAs such as rRNA, tRNA, snoRNA, snRNA etc. The detailed annotation protocol and statistics are described in Supplementary Material (below). All the sequencing libraries are accessible from GSE80830 in Gene Expression Omnibus.

After computationally removing the adaptor sequences, we mapped the extracted sequences to the reference Manduca genome using the Bowtie algorithm (Langmead et al., 2009). We only retained reads that matched the genome perfectly for downstream analysis. To identify potential piRNAs, we selected sequences that were longer than 22 nt and not annotated as miRNAs or other non-coding RNAs (see Supplementary Material for annotation of miRNAs and other non-coding RNAs). The reads of piRNAs that mapped to multiple loci in the genome were apportioned over these loci, and piRNA reads were normalized by the total number of genome mapping reads excluded other non-coding RNAs (rRNA, tRNA etc.). piRNA abundance is quantified in parts per million (ppm).

To determine if Manduca piRNAs amplify via the ping-pong cycle, we computed the Z-score as described in (Li et al., 2009). Briefly, we identified the piRNA pairs that mapped to overlapping genomic positions but on opposite genomic strands. We counted the numbers of such pairs with 5′-5′ overlapping distances from 1 to 20 nts, and calculated Z-score for the 10-nt overlap (the expected overlapping distance due to ping-pong) using the counts of 1–9 nt and 11–19 nt overlaps as the background.

In order to characterize the relative positions in mRNAs that Manduca piRNAs map to, we scaled the 5′ UTRs, coding regions, and 3′ UTRs of mRNAs to 350, 1000, and 800 nts respectively, which we calculated are the median lengths of these regions in annotated M. sexta genes. Using scaled non-overlapping windows which are equivalent to each of the 2,150 nts of the scaled genes, we counted piRNA abundance in RPKM. For measuring the enrichment, we performed Wilcoxon ranksum test on calculated RPKM values across the 5′ UTRs, coding regions, and 3′ UTRs of mRNAs.

We mapped reads in each RNA-seq library to the reference Manduca genome using the TopHat2 algorithm (Kim et al., 2013) allowing two mismatches (“-v 2”). To detect reads mapping to transposons, we allowed reads to map to multiple locations of the genome with the “-g” option in TopHat2. The most abundant transposon in Manduca is SINE (44,487 copies when all subfamilies are combined) (Supplementary Table S2), we ran TopHat2 with “-g 45,000.” Reads were apportioned by the number of times they mapped to the genome.

After mapping, we counted the number of RNA-seq reads for each gene and transposon, expressed in the unit of RPKM (Reads Per Kilobase of transcript per Million mapped reads). To identify differentially expressed genes and transposons between a pair of time points during ISM development, we ran the DESeq2 algorithm (version 1.5.5) implemented in R (Anders and Huber, 2010) using mapped read counts as the input, taking advantage of the three biological replicates per stage. The statistical analysis was performed by DESeq2 (version 1.5.5). Genes with q-value <0.01 and absolute fold-change (fc) >2 (both computed with DESeq2) were considered to be differentially expressed between the two time points.

On day 13 of the pupal-adult development the ISMs are fully functional and have yet to initiate either the atrophy or PCD (Schwartz and Truman 1983). We isolated the 18–30 nt small RNAs at this time point and split the sample into two aliquots. One aliquot was directly cloned and sequenced. The second aliquot was oxidized to render the miRNAs and mRNA degradation products unclonable due to the absence of 2′-O-methylated 3′ termini (Vagin et al., 2006; Grimson et al., 2008). This manipulation preferentially enriches the piRNAs. The small RNA-seq library mapping results are presented in Supplementary Table S1.

We obtained 17.8 million and 11.4 million reads in the control and oxidized day 13 small RNA libraries respectively, of which 71.7 and 75.8% mapped to the Manduca genome. When the mapped reads in the unoxidized small RNA library were partitioned by size, we observed a bimodal distribution, with peaks at 22 and 27 nt (Figure 1). Among the reads shorter than or equal to 23 nt, we annotated 77.7% (78.1% for total reads) as miRNAs (Tsuji et al., 2020). In sharp contrast, among the reads longer than 23 nt, only 0.4% were miRNAs and instead, these RNAs displayed a strong 5′ uridine bias and a weak adenine signal at the 10th position (Supplementary Figure S1).

Oxidation resulted in an almost complete loss of the 22 nt peak and 98.6% of the reads in the oxidized library were >23 nt (Figure 1). The majority of these reads started with a uridine base at the 5′ position (Supplementary Figure S1). Only 0.3% of the reads in this library were identified as miRNAs. In contrast, 53.6% of the reads mapped to genes, 21.4% mapped to transposons, and 23.1% mapped to unannotated regions of the genome (Figure 1). These percentages were much higher than the ≤23 nt reads in the unoxidized library (0.8, 0.3, and 0.6% respectively), and comparable to the reads >23 nt in the unoxidized library (9.7, 3.2, and 5.6% respectively). Based on size and resistance to oxidation the 27 nt peak appears to represent piRNAs.

In addition to the libraries on day 13, we generated and sequenced unoxidized small RNA libraries from seven more time points of ISM development: days 14, 15, 16, 17, 18, 1-h post-eclosion (PE) and 20-hydroxyecdysone (20E) treated animals. Similar to what we observed with the day 13 unoxidized sample, the small RNAs from the other time points also displayed a bimodal distribution of small RNAs, with the majority of the ∼22 nt sequences mapping to miRNAs (Supplementary Figure S1). In all of these samples, there was a strong bias for uridine in the first nucleotide, especially for the piRNAs (reads >23 nt).

piRNAs abundance from each stage was normalized by the sequencing depth and calculated as “parts per million” (ppm). The expression of both genic and transposable element piRNAs gradually declined from day 13 to day 16, increased sharply on day 17, and then declined dramatically on day 18 (Figure 2A). piRNA abundance was elevated in the 1-h post-eclosion (PE) sample. While treatment with 20E on day 17 delays ISM death (Schwartz and Truman 1983), it did not significantly alter piRNA abundance in the muscles relative to the corresponding PE timepoint. Interestingly the number of genic piRNAs was greater than those that map to transposons, the significance of which is unknown (Figure 2A).

Of the 64 known transposon families in Manduca, 34 generated piRNAs that mapped to the genome. Among them, piRNAs for 27 transposon families mirrored the pattern of expression observed for all piRNAs and transposon-mapping piRNAs (e.g. LTR:Copia in Figure 2A) during ISM development. In contrast, the expressions of piRNAs mapped to five transposon families (two DNA transposon families: P and TcMar; and two LINE families: CR1 and I; and one SINE family: tRNA) were reduced gradually during the ISM development (Supplementary Figure S2).

The majority of piRNAs in flies that map to transposons are antisense to the mRNAs in order to guide the PIWI-family proteins Aub and PIWI to repress transposon expression (Czech and Hannon 2016). Consequently, mutations in piRNA pathway proteins that alter this antisense bias, such as Qin, lead to transposon de-repression (Zhang et al., 2011). We examined the strand bias of transposon-mapping piRNAs at each stage of ISM development. Of the 34 transposon families with mapped piRNAs, we found that the majority of them (28 families) displayed strong antisense bias, while the remaining six families displayed a sense bias. The antisense piRNAs to transposons tend to be 27 nt in length (e.g. LTR:Gypsy in Figure 2B), and the sense piRNAs to transposons tend to be either 26 or 27 nt (e.g. DNA:PIF-Harbinger in Figure 2B). Although the relative abundance of piRNAs was regulated during the course of ISM development their length distributions were constant.

Following transcription from piRNA clusters, the primary transcripts are processed and cyclically amplified by a mechanism known as the “ping-pong amplification loop” (Meister 2013; Czech and Hannon 2016). The secondary piRNAs generated via ping-pong tend to have a 10 nt 5′ end overlap with other piRNAs on the opposite strand. This ping-pong activity is quantified by calculating the frequency of 5′-5′ overlaps between piRNA pairs, normalized as a Z-score (Li et al., 2009).

We detected high Z-scores at all stages examined. For example, piRNAs from day 17 had an overall Z-score of 234.6 (Figure 2A). piRNAs from both transposons and genic sequences were amplified by ping-pong as demonstrated by Z-scores of 98.7 and 105.5 respectively. Z-scores were high on day 13, fell precipitously during the next 3 days, and then rose on day 17 and remained high for the rest of development, which agrees well with the abundance of piRNAs in the tissue. It should be noted that these developmental changes in Z-score were not artifacts arising from variations in sequencing depth since these same patterns were observed when we down-sampled to 8 million reads per stage prior to our analysis.

As part of our analysis, we computed the Z-scores for piRNAs that mapped to each transposon class. Out of 34 transposon families with detectable piRNAs, 20 families displayed a high ping-pong signature throughout ISM development (e.g. LTR:Copia in Figure 2A), and 14 families displayed statistically significant Z-scores at least transiently. Intriguingly, in 12 out of 14 transposon families, the ping-pong signature peaked on day 17 in advance of the commitment of the ISMs to die. These data support the hypothesis that ISM piRNAs amplify via the traditional ping-pong amplification loop.

Combining the RepeatMasker result with our own gene annotation, 8,633 (45.9%) of the 18,806 genes in M. sexta contain transposons within their introns. From the mapping results with our oxidized small RNA-seq library, we observed that piRNAs (24.38 ppm per intron in median) fell into the transposon-derived introns from 6,902 genes (80.0% of 8,633 genes containing transposons). piRNAs also mapped to introns that did not contain transposons, although the abundance was far lower (1,731 out of 18,806 genes (9.2%)). Compared to the abundance of piRNA reads on transposon-derived introns, there were very few non-transposon introns (2.29 × 10⁻² ppm per intron in median).

It has been demonstrated that genic piRNAs tend to map to the sense orientation of 3′UTRs in germline and somatic cells (Robine et al., 2009; Jones et al., 2016; Sokolova et al., 2019). However, in the ISMs, we observed that the genic piRNAs tended to map to the sense strand preferentially within the 5′UTRs rather than the 3′ UTRs or coding sequences (CDSs). The strand bias patterns of the genic piRNAs uniquely mapped to 5′UTRs, CDSs, 3′UTRs, and introns were the same as those with all genic piRNA reads. We then focused specifically on the top 25% of the genes that are highly enriched with piRNAs. In the oxidized small RNA-seq library, 5′UTRs tended to be more enriched with piRNAs mapped to the sense strands of genes as compared to other gene domains. In unoxidized small RNA-seq libraries, the overall trends were stably observed through all time points across development.

Considering the result that genic piRNAs tend to fall into the sense strands of 5′UTRs, we investigated the relative mapping positions of piRNAs in genes by examining the oxidized small RNA-seq library from day 13. Interestingly, when the piRNAs that mapped to highly expressed genes were plotted onto the gene map, there was a striking peak in the sense orientation in the 5′UTRs (Figure 3).

(Only a modest number of piRNAs mapped to low abundance transcripts). In unoxidized small RNA-seq libraries, genic piRNAs tended to map to 5′UTRs and to form the peaks of mapped piRNAs (Supplementary Figure S3).

We also sought to determine if genic piRNA expression was correlated with the abundance of their source mRNAs. Toward that end, we examined the top 25% piRNA-generating genes that were either up- or down-regulated between days 17 and 18. (It should be noted that none of the piRNA pathway related genes were on the list, suggesting that genic piRNAs are unlikely to alter the abundance of the piRNA pathway machinery). As shown in Supplementary Figure S4, piRNA generation levels were unrelated to the RKPM levels. The RPKM distribution of the top 25% piRNA generating mRNAs ranged from 0.089 to 13,773.77 RPKM (median = 15.45 RPKM).

We next sought to determine which piRNA pathway protein components are expressed in the ISMs. Based on piRNA pathway genes characterized in Drosophila, we identified all 19 of the genes we sought in the ISMs (ago3, armi, aub, piwi, BoYb, egg, krimp, mael, papi, qin, rhi, shu, spn-E, tej, tud, vas, vret, hen-1, and zuc), plus two small RNA biogenesis pathway factors (ago1 and ago2). In agreement with published reports, we only identified a single Aub/PIWI protein sequence, which is also the case for other Lepidopterans like Tricoplusia ni and Bombyx mori (Cao et al., 2015). Consequently, we refer to this protein as Aub/PIWI.

Next, we examined the mRNA-seq reads from six developmental stages (days 13, 14, 15, 16, 17, 18) plus 20E-treated to determine if these factors are differentially expressed. While the expression levels of ago3 and aub appear to be low, their expression was nevertheless ranked within the top 64.4 and 31.1% respectively of all expressed genes (Figure 4A). The expression levels of other genes in this pathway were also within the top 18–65% of all expressed genes (Supplementary Figure S4).

Computational analysis suggests that the expression of many of the small RNA biogenesis pathway factors were developmentally regulated, with a general trend for stable or increasing expression prior to the initiation of atrophy (days 13–15), a sharp decline on day 17, and a near loss of expression when the ISMs became commitment to die on day 18. Four genes were significantly repressed on day 18: papi (q = 5.8 × 10⁻⁹, fc = −2.6), qin (q = 2.9 × 10⁻⁷, fc = −2.1), zuc (q = 3.9 × 10⁻³, fc = −2.1), and spn-E (q = 1.7 × 10⁻⁹, fc = −2.0) (Figure 4B) and five other demonstrated this general trend (ago3, armi, egg, krimp, and tud), although their changes did not reach statistical significance (Supplementary Figure S5). In contrast, mael mRNA levels were increased with the commitment to die (q = 6.4 × 10⁻¹⁹, fc = 2.5; Figure 4B) as were some members of the shu family (Figure 4B). In all cases, expression of piRNA pathway components were regulated by 20E and displayed their highest levels of expression when exposed to the exogenous steroid, suggesting that these genes are regulated by the same developmental signals that control developmental changes in the ISM, including atrophy and death.

We observed that the abundance of both piRNAs, and the majority of the factors that mediate their biogenesis, declined precipitously coincident with the commitment of the ISMs to die on day 18. Consequently, we sought to test the hypothesis that this loss might be correlated with changes in transposon expression. Using the RNA-seq reads mapped to transposon loci as the input, we computed the expression of each transposable element and the fold change during ISM development relative to day 13 (from day 14 to day 18, and 20E) (Figure 5).

Prior to day 18, the only significant change was the Tourist DNA mobile element, which was up-regulated on day 17 (q = 7.0 × 10⁻⁷, fc = 2.1) compared to day 13. Coincident with the muscle’s commitment to die on day 18, (but prior to the initiation of death later in the day), the patterns of the transposon expression changed dramatically. Four DNA transposable elements were up-regulated compared to those on day 13: Tourist (q = 8.0 × 10⁻¹³, fc = 4.0), hAT-Pegasus (q = 4.1 × 10⁻⁵, fc = 2.6), TcMar-Tc1 (q = 6.4 × 10⁻¹², fc = 2.2), and general DNA mobile elements (q = 4.6 × 10⁻¹¹, fc = 2.4). Concurrently, the expression levels of four transposons were significantly down-regulated on day 18 as compared to day 13: 5S-Deu (SINE; q = 4.7 × 10⁻⁴, fc = −2.0), LOA (LINE; q = 1.1 × 10⁻⁸, fc = −2.0), Penelope (LINE; q = 3.1 × 10⁻⁶, fc = −2.6), and CMC-Transib (DNA transposon; q = 4.7 × 10⁻⁴, fc = −2.3). Interestingly, pre-treatment with the steroid 20E, which delays ISM death, muted differential transposon expression of transposons.