Junko Tsuji1†
Travis Thomson2,3
Christine Brown4
Subhanita Ghosh3‡
William E. Theurkauf5
Zhiping Weng1*
Lawrence M. Schwartz4*- 1Program in Bioinformatics and Integrative Biology, University of Massachusetts Medical School, Worcester, MA, United States
- 2Program in Molecular Medicine, University of Massachusetts Medical School, Worcester, MA, United States
- 3Department of Neurobiology, University of Massachusetts Medical School, Worcester, MA, United States
- 4Department of Biology, University of Massachusetts, Amherst, MA, United States
- 5Program in Molecular Medicine, University of Massachusetts Medical School, Worcester, MA, United States
PIWI-interacting RNAs (piRNAs) are small single-stranded RNAs that can repress transposon expression via epigenetic silencing and transcript degradation. They have been identified predominantly in the ovary and testis, where they serve essential roles in transposon silencing in order to protect the integrity of the genome in the germline. The potential expression of piRNAs in somatic cells has been controversial. In the present study we demonstrate the expression of piRNAs derived from both genic and transposon RNAs in the intersegmental muscles (ISMs) from the tobacco hawkmoth Manduca sexta. These piRNAs are abundantly expressed, ∼27 nt long, map antisense to transposons, are oxidation resistant, exhibit a 5’ uridine bias, and amplify via the canonical ping-pong pathway. An RNA-seq analysis demonstrated that 19 piRNA pathway genes are expressed in the ISMs and are developmentally regulated. The abundance of piRNAs does not change when the muscles initiate developmentally-regulated atrophy, but are repressed coincident with the commitment of the muscles undergo programmed cell death at the end of metamorphosis. This change in piRNA expression is correlated with the repression of several retrotransposons and the induction of specific DNA transposons. The developmentally-regulated changes in the expression of piRNAs, piRNA pathway genes, and transposons are all regulated by 20-hydroxyecdysone, the steroid hormone that controls the timing of ISM death. Taken together, these data provide compelling evidence for the existence of piRNA in somatic tissues and suggest that they may play roles in developmental processes such as programmed cell death.
Introduction
Small silencing RNAs are powerful regulators of gene expression. They can lead to epigenetic silencing of transcription, transcript degradation, and inhibition of mRNA translation. The best characterized class of small silencing RNAs are the microRNAs (miRNAs)(Grewal and Elgin 2007; Ghildiyal and Zamore 2009; Castel and Martienssen 2013). miRNAs have been shown to be ∼22 nucleotides (nt) long, ubiquitously expressed and can repress their target transcripts through seed-based partial complementarity (Guo et al., 2010; Ha and Kim 2014; Bartel 2018).
The most recently discovered group of small silencing RNAs are PIWI-interacting RNAs (piRNAs). piRNAs are 23–35 nt in length and predominantly expressed in the germline of animals, including humans (Iwasaki et al., 2015; Williams et al., 2015; Czech and Hannon 2016; Czech et al., 2018). piRNAs guide the PIWI clade of Argonaute proteins to silence transposons and other selfish elements and function to protect the integrity of the germline genome (reviewed in Ozata et al., 2019).
The biogenesis and function of piRNAs has been well studied in the fruit fly Drosophila melanogaster (Huang, Fejes Toth et al., 2017). piRNAs are processed from long transcripts that can be up to hundreds of kilobases long that are transcribed from discrete genomic loci termed “piRNA clusters” (Brennecke et al., 2007; Malone et al., 2009; Thomson and Lin 2009). piRNAs are amplified via reciprocal target cleavages by PIWI proteins, a mechanism known as the ping-pong cycle (Brennecke et al., 2007; Gunawardane et al., 2007). Because PIWI proteins cleave the phosphodiester bond between the nucleotides in the target RNA that pair with the 10th and the 11th nucleotides of the guide piRNA, and 3′ cleavage products is subsequently made into another piRNA, there is an enrichment of piRNAs that perfectly reverse complement each other in their first 10 nucleotides, the hallmark of the ping-pong cycle. This process typically creates piRNAs with a uridine residue as the first nucleotide of the primary piRNA, and complementarity over the first 10 nt of post-transcriptionally amplified piRNAs (Brennecke et al., 2007; Cora et al., 2014; Wang et al., 2014; Stein et al., 2019).
piRNAs predominantly target transposons, retroviruses and other “selfish” genetic elements. In the absence of piRNAs, transposons can mobilize, resulting in double-stranded DNA breaks in the germline genome leading to infertility (Klattenhoff and Theurkauf 2008; Weick and Miska 2014). Consequently, the expression and function of piRNAs has been most extensively studied in germ cells, gonadal somatic cells, and certain progenitor cells (Brennecke et al., 2007; Rojas-Rios and Simonelig 2018; Sato and Siomi 2020).
While there have been several reports demonstrating the presence of both piRNAs and the associated protein machinery in non-gonadal cells, their role in cellular processes has been a controversial subject (Pandey et al., 2017; Huang et al., 2019; Perera et al., 2019; Onishi et al., 2020; Teefy et al., 2020). Somatic piRNAs have been observed broadly in arthropods, where they may provide genome defense against transposable elements and viruses (Lewis et al., 2018). In agreement with this hypothesis, piRNAs from fat body and midgut cells elicits antiviral response against nucleopolyhedrovirus in the silkmoth Bombyx mori (Feng et al., 2020). In addition, piRNAs derived from endogenous viral elements system in the mosquito Aedes aegypti help maintain long-lasting adaptive immunity (Whitfield et al., 2017; Tassetto et al., 2019).
In the current study we provide substantial data demonstrating the presence of piRNAs and their synthetic machinery in the intersegmental muscles (ISMs) from the tobacco hawkmoth Manduca sexta. We further demonstrate that transposon expression becomes deregulated coincident with the commitment of the muscles to initiate programmed cell death at the end of metamorphosis.
The ISMs are a classical model system for skeletal muscle atrophy and death (Schwartz and Truman 1982; Schwartz 2018; Sheel et al., 2020). These muscles are composed of giant syncytial cells, each of which is about 5 mm long and up to 1 mm in diameter. The ISMs are used by the larvae to crawl and by the developing adult moth to eclose (emerge) from the overlying pupal cuticle at the end of metamorphosis. Three days before adult eclosion (day 15 of the normal 18-day period of pupal-adult development), the ISMs initiate a program of atrophy that results in a ∼40% loss of muscle mass prior to eclosion. This atrophy is non-pathological and the muscles retain normal physiological properties such as resting potential and force per cross-sectional area (Schwartz and Ruff 2002). Late on day 18, coincident with adult eclosion, the ISMs initiate programmed cell death (PCD) and die during the subsequent 30 h. In fact, the term PCD was coined by Lockshin and Williams in 1965 to describe the death of these specific cells (Lockshin and Williams 1965). ISM PCD is a fundamentally different program than atrophy and results in the complete destruction of the contractile apparatus, loss of the resting potential, and enhanced autophagy (Schwartz and Ruff 2002; Tsuji et al., 2020). The dying cells are not phagocytosed and instead activate the molecular machinery required for both cellular destruction and nutrient recycling (Jones and Schwartz 2001; Tsuji et al., 2020). The developmental timing of both atrophy and death are controlled by circadian declines in the circulating titer of the insect molting hormone 20-hydroxyecdysone (20E) (Schwartz and Truman 1983). Judiciously timed administration of exogenous 20E can prevent atrophy or death, and the differential gene expression associated with these programs, but once either of these programs has been initiated, it cannot be altered or delayed by 20E treatment (Tsuji et al., 2020).
Several studies have demonstrated that ISM PCD requires de novo gene expression, and numerous death-associated genes have been identified (Schwartz et al., 1990; Schwartz et al., 1990; Tsuji et al., 2020). During the transition from atrophy to death, some ISM transcripts display significant changes in both their stability and translatability that is mediated by their 3′ untranslated regions (UTRs) (Cascone and Schwartz 2001). In subsequent studies, it was shown that specific microRNAs can regulate the translation of death-associated transcripts, establishing a role for this class of small RNAs (Tsuji et al., 2020).
To extend this analysis, we performed RNA-seq with the small RNAs isolated from the ISMs each day of development, ranging from before the initiation of atrophy (day 13) until when the muscles were committed to die (day 18) (Tsuji et al., 2020). We also analyzed ISMs from day 18 animals that had been injected the previous day with 20E to delay cell death. In addition to miRNAs, we also found that the ISMs also contain high levels of piRNAs, an unexpected result. Detailed analysis demonstrated that these piRNAs are: ∼27 nt long, map antisense to transposons, are oxidation resistant, exhibit a uridine bias at their first nucleotide, and amplify via the canonical ping-pong pathway. In addition, the ISMs express the genes required for piRNA synthesis and activity. When the ISMs become committed to die, there is a concurrent loss of piRNAs and an apparent deregulation of transposable elements, with repression of some retrotransposons and the induction of DNA transposons. The expression of piRNAs, piRNA regulatory machinery and transposable elements expression are all under the control of 20E. Thus, piRNAs are expressed in somatic tissues where they may regulate developmental processes such as PCD.
Materials and Methods
Animals
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.
RNA Isolation, Library Construction and Sequencing
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).
Oxidized Small RNA-Seq Library
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.
Genomic Sequence and Annotation Data
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.
Sequence Extraction and Annotation of piRNAs
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).
piRNA Ping-Pong Signature
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.
Relative Mapping Position of piRNAs on Genes
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.
Gene Expression and Differentially Expressed Genes
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.
Results
piRNAs in the Intersegmental Muscles Prior to Atrophy
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).
FIGURE 1. Length and mapping locations of Manduca piRNAs. Length distribution of the small RNAs in day 13 unoxidized (top) and oxidized (bottom) libraries. Genomic annotations of the locations where the reads mapped are summarized in pie charts (right).
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.
piRNA Expression Changes During the ISM Development
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).
FIGURE 2. Characteristics of piRNAs expressed in the ISM. (A) Changes in piRNA abundance during the ISM development (left); the frequency of 5′-5′ 10 nt overlaps between piRNA pairs (middle); and the corresponding ping-pong Z-scores during development (right). The shaded bars on the left panels indicate the fraction of piRNAs that are amplified via the ping-pong amplification loop. Almost all of the Z-scores were statistically significant (red dots). (B) Size distribution and strand bias of piRNAs mapped to transposons (left), and the change in piRNAs expression during the ISM development (right). Blue reflects sense strand bias while red is indicative of antisense. As examples, piRNAs mapped to LTR/Gypsy are in the upper panel while those mapping to DNA/PIF-Harbinger are on the bottom panel.
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).
piRNAs Predominantly Map Antisense to Transposons
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.
piRNAs Expressed in the ISM Amplify via the Ping-Pong Amplification Loop
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.
Genic piRNAs Preferentially Map to 5′UTRs
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−2 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).
FIGURE 3. Relative mapping positions of genic piRNAs. The relative mapping position of piRNAs on highly expressed genes. piRNAs mapping to the sense and antisense strands of genes are highlighted in blue and red respectively.
(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).
piRNA Biogenesis Pathway Factors Are Differentially Expressed During ISM Atrophy and the Commitment to Die
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).
FIGURE 4. Expression of piRNA biogenesis pathway factors. (A) Gene expression of ago3 and aub/piwi based on RNA-seq analyses (top) (B) Downregulated (zuc, qin, spn-E, hen1, brother of Yb, tejes and papi), mixed (shu), and upregulated (mael) piRNA pathway factors in the ISMs prior to the initiation of cell death. The difference between each developmental time point (days 13–17) were statistically different from the day 18 value, although we only highlighted the comparison of day 13–day 18 (red and blue circles). Red and blue circles on day 18 indicate up- and down-regulation of the factors compared to expression on day 13. Those circles in the 20E lane indicate statistically significant changes in ISM gene expression animals injected on day 17 with 20E to delay cell death on day 18. Means are presented plus/minus standard deviations calculated from three replicates of the mRNA-seq libraries.
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−9, fc = −2.6), qin (q = 2.9 × 10−7, fc = −2.1), zuc (q = 3.9 × 10−3, fc = −2.1), and spn-E (q = 1.7 × 10−9, 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−19, 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.
Transposon Expression Becomes Deregulated when the ISM Become Committed to Die
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).
FIGURE 5. Transposon expression in atrophy and ping-pong piRNAs. Changes in transposon expression during ISM development relative to day 13 (from day 14 to day 18, and 20E). Up-regulated and down-regulated transposons are highlighted in red and blue respectively.
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−7, 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−13, fc = 4.0), hAT-Pegasus (q = 4.1 × 10−5, fc = 2.6), TcMar-Tc1 (q = 6.4 × 10−12, fc = 2.2), and general DNA mobile elements (q = 4.6 × 10−11, 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−4, fc = −2.0), LOA (LINE; q = 1.1 × 10−8, fc = −2.0), Penelope (LINE; q = 3.1 × 10−6, fc = −2.6), and CMC-Transib (DNA transposon; q = 4.7 × 10−4, fc = −2.3). Interestingly, pre-treatment with the steroid 20E, which delays ISM death, muted differential transposon expression of transposons.
Discussion
In this report we present substantial evidence supporting the presence of developmentally-regulated piRNAs in striated skeletal muscle, a highly differentiated somatic cell. Like piRNAs characterized from ovary and testis, Manduca ISM piRNAs: are ∼27 nucleotides long, have a strong 5’ uridine bias, are oxidation resistant, amplify via a ping-pong mechanism, map to transposable elements, and are predominantly antisense. (We were unable to determine if these piRNAs physically associated with Aub as none of the anti-Drosophila Aub antibodies that we tested recognized Manduca Aub/PIWI (data not shown)). While piRNAs from Manduca ovary have not been examined, ovarian piRNAs in the related Lepidopteran Bombyx mori are also ∼27 nt long, map antisense to transposons, and amplify via ping-pong mechanism (Kawaoka et al., 2011).
ISM piRNAs represent a substantial percentage of the small RNA pool within the tissue and are neither low abundance sequences nor potential contaminants from mRNA degradation, as oxidized mRNA degradation products cannot be cloned (Zhang et al., 2014). The ISMs are a very large and discrete tissue that can be isolated cleanly from the animal at each stage of development without contamination from other piRNA-rich tissues such as gonads (Schwartz 2018). (Indeed, gonad-specific transcripts, such as the vitellogenin receptor and caudal, which are restricted to the ovary, are expressed at 0.01–0.04 RPKM during the 6 days of development examined in our study, which are orders of magnitude lower than the majority of the piRNA pathway genes that we analyzed). As well, unlike mammalian muscle, the ISMs are composed of only a single fiber type and do not contain regenerative stem cells like satellite cells, endothelial cells or pericytes that might complicate the analysis (Beaulaton and Lockshin 1977).
The expression of piRNAs and the machinery required for their synthesis and action in somatic cells is controversial and the subject of debate (Pandey et al., 2017; Lewis et al., 2018; Onishi et al., 2020). Some studies have suggested that piRNAs may act in the nervous system and influence transposon activity as well as aspects of learning and memory. For example, in Drosophila, the piRNA pathway proteins Aub and Ago3 regulate transposon expression and mobilization in the mushroom bodies, brain structures that regulate memory formation and cognitive function (Perrat et al., 2013). These proteins are also found in specialized structures within glial cells in the adult Drosophila brain where they appear to repress transposon activity (Tindell et al., 2020). In the sea slug Aplysia, a specific piRNA has been shown to modulate synaptic plasticity and memory (Lee et al., 2011; Rajasethupathy et al., 2012).
piRNAs have been also reported in cancer cells (Meister, 2013) and regenerative stem-like cells in invertebrates (Kim et al., 2020). For example, PIWI-like proteins and piRNAs have been identified in somatic cells of the Cnidarian Hydra and the flatworm Macrostomum, although it is thought that their expression is restricted to stem-like progenitor cells that have the capacity to regenerate all cell types of the animal, including gonad (Juliano et al., 2014; Zhou et al., 2015). Loss of somatic piRNA pathway in Drosophila fat body disrupts metabolic homeostasis by depleting lipids and stored metabolites leading shortened lifespan (Jones et al., 2016). In the silkmoth Bombyx mori, the primary sex-determining factor is regulated by a single piRNA originated from the sex-determining genomic locus of W chromosome (Kiuchi et al., 2014). In the nematode C. elegans, the piRNA pathway components, including PRDE-1 and PRG-1, are expressed in neurons and their inhibition facilitates sensory neuron regeneration following injury in a cell autonomous manner (Kim et al., 2018).
The best characterized system for the analysis of somatic piRNAs is follicle cells in the Drosophila ovary (Olivieri et al., 2010; Ishizu et al., 2015; Onishi et al., 2020). Follicle cells are epithelial cells that are needed for the survival and maturation of the underlying germline cells, and disruption of the piRNA pathway in these cells leads to enhanced transposon activity and sterility (Olivieri et al., 2010). While follicle cells are derived somatically, they are physically connected to the germline via large ring canals, so it is not clear how much interplay may exist between them. Nevertheless, piRNA production can take place in the follicle-derived OSS (ovary somatic sheet) cell line, supporting the hypothesis that these specialized somatic cells can produce piRNAs (Pelisson et al., 2007; Li et al., 2009; Saito et al., 2009; Sato and Siomi 2020). However, piRNA production in these cells differs from that seen in ovarian cells in several key ways (Li et al., 2009; Malone et al., 2009). First, follicle cells lack both Aub and Ago3, and consequently do not display ping-pong amplification (Malone et al., 2009). Secondly, these piRNAs are expressed predominantly from uni-strand genomic piRNA clusters, while ovarian cells produce piRNAs from both single and double-stranded transposable element clusters (Malone et al., 2009; Yashiro et al., 2018). In contrast to follicle cells, the ISMs express both Aub and Ago3 and display efficient ping-pong amplification with high Z-scores. These data suggest that while the ISMs are clearly somatic cells, their piRNA pathway more closely represents the hallmarks of germline piRNA production.
Sequencing and mapping analyses have demonstrated that Manduca piRNAs are generated from both transposable elements and protein coding genes. In the few instances where genic piRNAs have been analyzed, they have been shown to be derived predominantly from the 3′ UTR sequences within mRNAs (Robine et al., 2009; Saito et al., 2009; Jones et al., 2016; Sokolova et al., 2019). In contrast, genic piRNAs analyzed from Manduca appear to be derived primarily from the start of translation in the 5′UTR. The possible regulatory role that these sequences might serve is unclear.
The ISMs are a classic model system for the study of skeletal muscle atrophy and death (Finlayson 1956; Lockshin and Williams 1965; Schwartz et al., 2016). Data presented here demonstrate that the expression of ISM piRNAs is developmentally-regulated, with low levels during the window when the ISMs undergo developmentally-regulated atrophy and peak expression on day 17, prior to the time when the ISMs become committed to die late on day 18. Consequently, we examined the expression of transposable elements since they are the primary target of piRNAs. Transposon expression appears to be tightly controlled in the ISMs since there was almost no variation their abundance until day 17, at which point they became deregulated. Within a matter of hours, the expression of several DNA transposons increased significantly, while some elements, most notably retrotransposons, were concurrently repressed. The reason for this disparity is not known, but DNA elements integrate directly into the genome via a cut-and-paste mechanism, while retrotransposons may integrate via a copy-and-paste mechanism (Sultana et al., 2017). The ability of the ISMs to initiate PCD occurs when the circulating levels of the steroid hormone 20E decline below a specific threshold on day 17 (Schwartz and Truman 1983). Hormone replacement with exogenous 20E on day 17 not only delays ISM death, it also prevented the developmental changes in the expression of piRNAs, piRNA protein machinery and transposon. This data supports the hypothesis that the piRNA pathway is also under hormonal control. It should be noted that the changes in transposon expression occur well in advance of the initiation of cell death and therefore do not appear to be a secondary consequence of cellular suicide. In the germline, transposon expression is repressed in order to protect the genome from insertional mutagenesis and subsequent catastrophe. However, this same process might confer an advantage for the organism by ensuring that cell death is indeed an irreversible process. For example, during apoptosis, genome destruction is insured by the activation of endogenous nucleases that cleave chromosomal DNA into nucleosome sized fragments (Wyllie et al., 1980). However, the ISMs die by an autophagic process that does not include DNA fragmentation (Beaulaton and Lockshin 1977; Schwartz et al., 1993; Schwartz et al., 2016). One intriguing hypothesis is that transposon expression and reintegration might serve a similar role for cells undergoing non-apoptotic forms of cell death where it insures that condemned cells are truly “dead” by fragmenting the genome and depleting the cell of beta-nicotinamide adenine dinucleotide and ATP (Berger 1985). Indeed, some support for this hypothesis comes from the observation that transposon expression is elevated in certain neurodegenerative disorders (Li et al., 2012; Tan et al., 2012; Jonsson et al., 2020), a phenomenon that has been called a “transposon storm” (Reilly et al., 2013; Dubnau 2018). To date, none of these studies have examined the possible expression of piRNAs in the diseased tissues. We attempted to directly test the hypothesis that reintegration of transposons in Manduca would result in double stranded breaks in genomic DNA by sequencing genomic libraries generated from ISMs taken before and after adult eclosion. Unfortunately, laboratory reared Manduca are not isogenic, and the individual-to-individual variability in transposon copy number precluded a direct test of the hypothesis (data not shown). As well, we sought to use antibodies against phosphorylated histone H2Av since it is a marker of double stranded DNA breaks (Madigan et al., 2002), but the antibodies directed against Drosophila phosphorylated H2Av did not cross react with the Manduca protein, and a Lepidopteran-specific antibody has not been identified (data not shown). In the absence of data demonstrating that the elevated expression of transposons is associated with enhanced genomic integration, we can only speculate about their possible role(s) in the cell death process.
Taken together, we have presented compelling evidence that piRNAs are expressed in highly differentiated somatic cells. Their expression is developmentally regulated and controlled by the steroid hormone 20E. piRNA expression is repressed coincident with the commitment of the ISMs to die, which is correlated with the deregulation of transposon expression. These data verifies the expression of somatic piRNAs and suggest the intriguing possibility that they impact the process of cell death during development, and perhaps pathogenesis.
Author Contributions
JT performed the bioinformatics data analyses, wrote the materials and methods and results sections, and contributed substantially to the discussion. TT helped with the library constructions, data analysis and writing. CB purified mRNA, performed quality control analyses of the libraries, and general lab analyses. SG contributed to the data analysis and discussion. WT and ZW helped design the study and analyses, and provided supervision. LS performed all of the dissections, library constructions, wrote the background and discussion and oversaw the study.
Funding
This work was supported by a Life Sciences Moment Fund grant from the University of Massachusetts Center for Clinical and Translational Science to LS, ZW, and WT, the Eugene M. and Ronnie Isenberg Professorship Endowment (LS), NIH grant P01HD078253 to WT and ZW, and National Institute of Child Health and Human Development Grant P01HD078253 to WT and ZW.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
We thank Wendy Smith for the provision of animals and Leonid Pobezinsky for helpful discussions.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.775369/full#supplementary-material
References
Beaulaton, J., and Lockshin, R. A. (1977). Ultrastructural Study of the Normal Degeneration of the Intersegmental Muscles ofAntheraea polyphemus andManduca Sexta (Insecta, Lepidoptera) with Particular Reference to Cellular Autophagy. J. Morphol. 154 (1), 39–57. doi:10.1002/jmor.1051540104
Berger, N. A. (1985). Poly(ADP-ribose) in the Cellular Response to DNA Damage. Radiat. Res. 101 (1), 4–15. doi:10.2307/3576299
Brennecke, J., Aravin, A. A., Stark, A., Dus, M., Kellis, M., Sachidanandam, R., et al. (2007). Discrete Small RNA-Generating Loci as Master Regulators of Transposon Activity in Drosophila. Cell 128 (6), 1089–1103. doi:10.1016/j.cell.2007.01.043
Cao, X., He, Y., Hu, Y., Wang, Y., Chen, Y.-R., Bryant, B., et al. (2015). The Immune Signaling Pathways of Manduca Sexta. Insect Biochem. Mol. Biol. 62, 64–74. doi:10.1016/j.ibmb.2015.03.006
Cascone, P. J., and Schwartz, L. M. (2001). Post-transcriptional Regulation of Gene Expression During the Programmed Death of Insect Skeletal Muscle. Dev. Genes Evol. 211 (8-9), 397–405. doi:10.1007/s004270100172
Castel, S. E., and Martienssen, R. A. (2013). RNA Interference in the Nucleus: Roles for Small RNAs in Transcription, Epigenetics and Beyond. Nat. Rev. Genet. 14 (2), 100–112. doi:10.1038/nrg3355
Cora, E., Pandey, R. R., Xiol, J., Taylor, J., Sachidanandam, R., McCarthy, A. A., et al. (2014). The MID-PIWI Module of Piwi Proteins Specifies Nucleotide- and Strand-Biases of piRNAs. RNA 20 (6), 773–781. doi:10.1261/rna.044701.114
Czech, B., and Hannon, G. J. (2016). One Loop to Rule Them All: The Ping-Pong Cycle and piRNA-Guided Silencing. Trends Biochem. Sci. 41 (4), 324–337. doi:10.1016/j.tibs.2015.12.008
Czech, B., Munafò, M., Ciabrelli, F., Eastwood, E. L., Fabry, M. H., Kneuss, E., et al. (2018). piRNA-Guided Genome Defense: From Biogenesis to Silencing. Annu. Rev. Genet. 52, 131–157. doi:10.1146/annurev-genet-120417-031441
Dubnau, J. (2018). The Retrotransposon Storm and the Dangers of a Collyer's Genome. Curr. Opin. Genet. Develop. 49, 95–105. doi:10.1016/j.gde.2018.04.004
Feng, M., Kolliopoulou, A., Zhou, Y. H., Fei, S. G., Xia, J. M., Swevers, L., et al. (2020). The piRNA Response to BmNPV Infection in the Silkworm Fat Body and Midgut. Insect Sci., 662–679. doi:10.1111/1744-7917.12796
Finlayson, L. H. (1956). Normal and Induced Degeneration of Abdominal Muscles During Metamorphosis in the Lepidoptera. J. Cel Sci. s3-97 (38), 215–233. doi:10.1242/jcs.s3-97.38.215
Ghildiyal, M., and Zamore, P. D. (2009). Small Silencing RNAs: An Expanding Universe. Nat. Rev. Genet. 10 (2), 94–108. doi:10.1038/nrg2504
Grewal, S. I. S., and Elgin, S. C. R. (2007). Transcription and RNA Interference in the Formation of Heterochromatin. Nature 447 (7143), 399–406. doi:10.1038/nature05914
Grimson, A., Srivastava, M., Fahey, B., Woodcroft, B. J., Chiang, H. R., King, N., et al. (2008). Early Origins and Evolution of microRNAs and Piwi-Interacting RNAs in Animals. Nature 455 (7217), 1193–1197. doi:10.1038/nature07415
Gunawardane, L. S., Saito, K., Nishida, K. M., Miyoshi, K., Kawamura, Y., Nagami, T., et al. (2007). A Slicer-Mediated Mechanism for Repeat-Associated siRNA 5' End Formation in Drosophila. Science 315 (5818), 1587–1590. doi:10.1126/science.1140494
Guo, H., Ingolia, N. T., Weissman, J. S., and Bartel, D. P. (2010). Mammalian microRNAs Predominantly Act to Decrease Target mRNA Levels. Nature 466 (7308), 835–840. doi:10.1038/nature09267
Ha, M., and Kim, V. N. (2014). Regulation of microRNA Biogenesis. Nat. Rev. Mol. Cel Biol 15 (8), 509–524. doi:10.1038/nrm3838
Huang, S., Ichikawa, Y., Igarashi, Y., Yoshitake, K., Kinoshita, S., Omori, F., et al. (2019). Piwi-interacting RNA (piRNA) Expression Patterns in Pearl Oyster (Pinctada Fucata) Somatic Tissues. Sci. Rep. 9 (1), 247. doi:10.1038/s41598-018-36726-0
Huang, X., Fejes Tóth, K., and Aravin, A. A. (2017). piRNA Biogenesis in Drosophila melanogaster. Trends Genet. 33 (11), 882–894. doi:10.1016/j.tig.2017.09.002
Ishizu, H., Iwasaki, Y. W., Hirakata, S., Ozaki, H., Iwasaki, W., Siomi, H., et al. (2015). Somatic Primary piRNA Biogenesis Driven by Cis-Acting RNA Elements and Trans-acting Yb. Cel Rep. 12 (3), 429–440. doi:10.1016/j.celrep.2015.06.035
Iwasaki, Y. W., Siomi, M. C., and Siomi, H. (2015). PIWI-interacting RNA: Its Biogenesis and Functions. Annu. Rev. Biochem. 84, 405–433. doi:10.1146/annurev-biochem-060614-034258
Jones, B. C., Wood, J. G., Chang, C., Tam, A. D., Franklin, M. J., Siegel, E. R., et al. (2016). A Somatic piRNA Pathway in the Drosophila Fat Body Ensures Metabolic Homeostasis and Normal Lifespan. Nat. Commun. 7, 13856. doi:10.1038/ncomms13856
Jones, M., and Schwartz, L. M. (2001). Not All Muscles Meet the Same Fate When They Die. Cel Biol. Int. 25 (6), 539–545. doi:10.1006/cbir.2000.0669
Jönsson, M. E., Garza, R., Johansson, P. A., and Jakobsson, J. (2020). Transposable Elements: A Common Feature of Neurodevelopmental and Neurodegenerative Disorders. Trends Genet. 36 (8), 610–623. doi:10.1016/j.tig.2020.05.004
Juliano, C. E., Reich, A., Liu, N., Götzfried, J., Zhong, M., Uman, S., et al. (2014). PIWI Proteins and PIWI-Interacting RNAs Function in Hydra Somatic Stem Cells. Proc. Natl. Acad. Sci. USA 111 (1), 337–342. doi:10.1073/pnas.1320965111
Kawaoka, S., Arai, Y., Kadota, K., Suzuki, Y., Hara, K., Sugano, S., et al. (2011). Zygotic Amplification of Secondary piRNAs During Silkworm Embryogenesis. RNA 17 (12), 1401–1407. doi:10.1261/rna.2709411
Kim, D., Pertea, G., Trapnell, C., Pimentel, H., Kelley, R., and Salzberg, S. L. (2013). TopHat2: Accurate Alignment of Transcriptomes in the Presence of Insertions, Deletions and Gene Fusions. Genome Biol. 14 (4), R36. doi:10.1186/gb-2013-14-4-r36
Kim, I. V., Riedelbauch, S., and Kuhn, C.-D. (2020). The piRNA Pathway in Planarian Flatworms: New Model, New Insights. Biol. Chem. 401 (10), 1123–1141. doi:10.1515/hsz-2019-0445
Kim, K. W., Tang, N. H., Andrusiak, M. G., Wu, Z., Chisholm, A. D., and Jin, Y. (2018). A Neuronal piRNA Pathway Inhibits Axon Regeneration in C. elegans. Neuron 97 (3), 511–519. doi:10.1016/j.neuron.2018.01.014
Kiuchi, T., Koga, H., Kawamoto, M., Shoji, K., Sakai, H., Arai, Y., et al. (2014). A Single Female-specific piRNA Is the Primary Determiner of Sex in the Silkworm. Nature 509 (7502), 633–636. doi:10.1038/nature13315
Klattenhoff, C., and Theurkauf, W. (2008). Biogenesis and Germline Functions of piRNAs. Development 135 (1), 3–9. doi:10.1242/dev.006486
Langmead, B., Trapnell, C., Pop, M., and Salzberg, S. L. (2009). Ultrafast and Memory-Efficient Alignment of Short DNA Sequences to the Human Genome. Genome Biol. 10 (3), R25. doi:10.1186/gb-2009-10-3-r25
Lee, E. J., Banerjee, S., Zhou, H., Jammalamadaka, A., Arcila, M., Manjunath, B. S., et al. (2011). Identification of piRNAs in the Central Nervous System. RNA 17 (6), 1090–1099. doi:10.1261/rna.2565011
Lewis, S. H., Quarles, K. A., Yang, Y., Tanguy, M., Frézal, L., Smith, S. A., et al. (2018). Pan-arthropod Analysis Reveals Somatic piRNAs as an Ancestral Defence Against Transposable Elements. Nat. Ecol. Evol. 2 (1), 174–181. doi:10.1038/s41559-017-0403-4
Li, C., Vagin, V. V., Lee, S., Xu, J., Ma, S., Xi, H., et al. (2009). Collapse of Germline piRNAs in the Absence of Argonaute3 Reveals Somatic piRNAs in Flies. Cell 137 (3), 509–521. doi:10.1016/j.cell.2009.04.027
Li, W., Jin, Y., Prazak, L., Hammell, M., and Dubnau, J. (2012). Transposable Elements in TDP-43-Mediated Neurodegenerative Disorders. PLoS One 7 (9), e44099. doi:10.1371/journal.pone.0044099
Lockshin, R. A., and Williams, C. M. (1965). Programmed Cell Death-I. Cytology of Degeneration in the Intersegmental Muscles of the Pernyi Silkmoth. J. Insect Physiol. 11, 123–133. doi:10.1016/0022-1910(65)90099-5
Madigan, J. P., Chotkowski, H. L., and Glaser, R. L. (2002). DNA Double-Strand Break-Induced Phosphorylation of Drosophila Histone Variant H2Av Helps Prevent Radiation-Induced Apoptosis. Nucleic Acids Res. 30 (17), 3698–3705. doi:10.1093/nar/gkf496
Malone, C. D., Brennecke, J., Dus, M., Stark, A., McCombie, W. R., Sachidanandam, R., et al. (2009). Specialized piRNA Pathways Act in Germline and Somatic Tissues of the Drosophila Ovary. Cell 137 (3), 522–535. doi:10.1016/j.cell.2009.03.040
Meister, G. (2013). Argonaute Proteins: Functional Insights and Emerging Roles. Nat. Rev. Genet. 14 (7), 447–459. doi:10.1038/nrg3462
Olivieri, D., Sykora, M. M., Sachidanandam, R., Mechtler, K., and Brennecke, J. (2010). An In Vivo RNAi Assay Identifies Major Genetic and Cellular Requirements for Primary piRNA Biogenesis in Drosophila. EMBO J. 29 (19), 3301–3317. doi:10.1038/emboj.2010.212
Onishi, R., Sato, K., Murano, K., Negishi, L., Siomi, H., and Siomi, M. C. (2020). Piwi Suppresses Transcription of Brahma-dependent Transposons via Maelstrom in Ovarian Somatic Cells. Sci. Adv. 6 (50), 7420. doi:10.1126/sciadv.aaz7420
Ozata, D. M., Gainetdinov, I., Zoch, A., O’Carroll, D., and Zamore, P. D. (2019). PIWI-interacting RNAs: Small RNAs with Big Functions. Nat. Rev. Genet. 20 (2), 89–108. doi:10.1038/s41576-018-0073-3
Pandey, R. R., Homolka, D., Chen, K.-M., Sachidanandam, R., Fauvarque, M.-O., and Pillai, R. S. (2017). Recruitment of Armitage and Yb to a Transcript Triggers its Phased Processing into Primary piRNAs in Drosophila Ovaries. Plos Genet. 13 (8), e1006956. doi:10.1371/journal.pgen.1006956
Pélisson, A., Sarot, E., Payen-Groschêne, G., and Bucheton, A. (2007). A Novel Repeat-Associated Small Interfering RNA-Mediated Silencing Pathway Downregulates Complementary Sense Gypsy Transcripts in Somatic Cells of the Drosophila Ovary. J. Virol. 81 (4), 1951–1960. doi:10.1128/JVI.01980-06
Perera, B. P. U., Tsai, Z. T.-Y., Colwell, M. L., Jones, T. R., Goodrich, J. M., Wang, K., et al. (2019). Somatic Expression of piRNA and Associated Machinery in the Mouse Identifies Short, Tissue-specific piRNA. Epigenetics 14 (5), 504–521. doi:10.1080/15592294.2019.1600389
Perrat, P. N., DasGupta, S., Wang, J., Theurkauf, W., Weng, Z., Rosbash, M., et al. (2013). Transposition-driven Genomic Heterogeneity in the Drosophila Brain. Science 340 (6128), 91–95. doi:10.1126/science.1231965
Rajasethupathy, P., Antonov, I., Sheridan, R., Frey, S., Sander, C., Tuschl, T., et al. (2012). A Role for Neuronal piRNAs in the Epigenetic Control of Memory-Related Synaptic Plasticity. Cell 149 (3), 693–707. doi:10.1016/j.cell.2012.02.057
Reilly, M. T., Faulkner, G. J., Dubnau, J., Ponomarev, I., and Gage, F. H. (2013). The Role of Transposable Elements in Health and Diseases of the Central Nervous System. J. Neurosci. 33 (45), 17577–17586. doi:10.1523/jneurosci.3369-13.2013
Robine, N., Lau, N. C., Balla, S., Jin, Z., Okamura, K., Kuramochi-Miyagawa, S., et al. (2009). A Broadly Conserved Pathway Generates 3′UTR-Directed Primary piRNAs. Curr. Biol. 19 (24), 2066–2076. doi:10.1016/j.cub.2009.11.064
Rojas-Ríos, P., and Simonelig, M. (2018). piRNAs and PIWI Proteins: Regulators of Gene Expression in Development and Stem Cells. Development 145 (17), 786. doi:10.1242/dev.161786
Saito, K., Inagaki, S., Mituyama, T., Kawamura, Y., Ono, Y., Sakota, E., et al. (2009). A Regulatory Circuit for Piwi by the Large Maf Gene Traffic Jam in Drosophila. Nature 461 (7268), 1296–1299. doi:10.1038/nature08501
Sato, K., and Siomi, M. C. (2020). The piRNA Pathway in Drosophila Ovarian Germ and Somatic Cells. Proc. Jpn. Acad. Ser. B: Phys. Biol. Sci. 96 (1), 32–42. doi:10.2183/pjab.96.003
Schwartz, L. M. (2018). Skeletal Muscles Do Not Undergo Apoptosis During Either Atrophy or Programmed Cell Death-Revisiting the Myonuclear Domain Hypothesis. Front. Physiol. 9, 1887. doi:10.3389/fphys.2018.01887
Schwartz, L. M., Brown, C., McLaughlin, K., Smith, W., and Bigelow, C. (2016). The Myonuclear Domain Is Not Maintained in Skeletal Muscle During Either Atrophy or Programmed Cell Death. Am. J. Physiology-Cell Physiol. 311 (4), C607–C615. doi:10.1152/ajpcell.00176.2016
Schwartz, L. M., Kosz, L., and Kay, B. K. (1990). Gene Activation Is Required for Developmentally Programmed Cell Death. Proc. Natl. Acad. Sci. 87 (17), 6594–6598. doi:10.1073/pnas.87.17.6594
Schwartz, L. M., Myer, A., Kosz, L., Engelstein, M., and Maier, C. (1990). Activation of Polyubiquitin Gene Expression During Developmentally Programmed Cell Death. Neuron 5 (4), 411–419. doi:10.1016/0896-6273(90)90080-y
Schwartz, L. M., and Ruff, R. L. (2002). Changes in Contractile Properties of Skeletal Muscle During Developmentally Programmed Atrophy and Death. Am. J. Physiology-Cell Physiol. 282 (6), C1270–C1277. doi:10.1152/ajpcell.01275.2000
Schwartz, L. M., Smith, S. W., Jones, M. E., and Osborne, B. A. (1993). Do all Programmed Cell Deaths Occur via Apoptosis. Proc. Natl. Acad. Sci. 90 (3), 980–984. doi:10.1073/pnas.90.3.980
Schwartz, L. M., and Truman, J. W. (1984). Hormonal Control of Muscle Atrophy and Degeneration in the Moth Antheraea polyphemus. J. Exp. Biol. 111, 13–30. doi:10.1242/jeb.111.1.13
Schwartz, L. M., and Truman, J. W. (1983). Hormonal Control of Rates of Metamorphic Development in the Tobacco Hornworm Manduca Sexta. Develop. Biol. 99 (1), 103–114. doi:10.1016/0012-1606(83)90257-9
Schwartz, L. M., and Truman, J. W. (1982). Peptide and Steroid Regulation of Muscle Degeneration in an Insect. Science 215 (4538), 1420–1421. doi:10.1126/science.6278594
Sheel, A., Shao, R., Brown, C., Johnson, J., Hamilton, A., Sun, D., et al. (2020). Acheron/Larp6 Is a Survival Protein that Protects Skeletal Muscle from Programmed Cell Death During Development. Front. Cel Dev. Biol. 8, 622. doi:10.3389/fcell.2020.00622
Sokolova, O. A., Ilyin, A. A., Poltavets, A. S., Nenasheva, V. V., Mikhaleva, E. A., Shevelyov, Y. Y., et al. (2019). Yb Body Assembly on the Flamenco piRNA Precursor Transcripts Reduces Genic piRNA Production. MBoC 30 (12), 1544–1554. doi:10.1091/mbc.e17-10-0591
Stein, C. B., Genzor, P., Mitra, S., Elchert, A. R., Ipsaro, J. J., Benner, L., et al. (2019). Decoding the 5′ Nucleotide Bias of PIWI-Interacting RNAs. Nat. Commun. 10 (1), 828. doi:10.1038/s41467-019-08803-z
Sultana, T., Zamborlini, A., Cristofari, G., and Lesage, P. (2017). Integration Site Selection by Retroviruses and Transposable Elements in Eukaryotes. Nat. Rev. Genet. 18 (5), 292–308. doi:10.1038/nrg.2017.7
Tan, H., Qurashi, A., Poidevin, M., Nelson, D. L., Li, H., and Jin, P. (2012). Retrotransposon Activation Contributes to Fragile X Premutation rCGG-Mediated Neurodegeneration. Hum. Mol. Genet. 21 (1), 57–65. doi:10.1093/hmg/ddr437
Tassetto, M., Kunitomi, M., Whitfield, Z. J., Dolan, P. T., Sánchez-Vargas, I., Garcia-Knight, M., et al. (2019). Control of RNA Viruses in Mosquito Cells Through the Acquisition of vDNA and Endogenous Viral Elements. Elife 8, 44. doi:10.7554/eLife.41244
Tautz, D., Hancock, J. M., Webb, D. A., Tautz, C., and Dover, G. A. (1988). Complete Sequences of the rRNA Genes of Drosophila melanogaster. Mol. Biol. Evol. 5 (4), 366–376. doi:10.1093/oxfordjournals.molbev.a040500
Teefy, B. B., Siebert, S., Cazet, J. F., Lin, H., and Juliano, C. E. (2020). PIWI-piRNA Pathway-Mediated Transposable Element Repression in Hydra Somatic Stem Cells. RNA 26 (5), 550–563. doi:10.1261/rna.072835.119
Thomson, T., and Lin, H. (2009). The Biogenesis and Function of PIWI Proteins and piRNAs: Progress and Prospect. Annu. Rev. Cel Dev. Biol. 25, 355–376. doi:10.1146/annurev.cellbio.24.110707.175327
Tindell, S. J., Rouchka, E. C., and Arkov, A. L. (2020). Glial Granules Contain Germline Proteins in the Drosophila Brain, Which Regulate Brain Transcriptome. Commun. Biol. 3 (1), 699. doi:10.1038/s42003-020-01432-z
Tsuji, J., Thomson, T., Chan, E., Brown, C. K., Oppenheimer, J., Bigelow, C., et al. (2020). High-resolution Analysis of Differential Gene Expression During Skeletal Muscle Atrophy and Programmed Cell Death. Physiol. Genomics 52 (10), 492–511. doi:10.1152/physiolgenomics.00047.2020
Vagin, V. V., Sigova, A., Li, C., Seitz, H., Gvozdev, V., and Zamore, P. D. (2006). A Distinct Small RNA Pathway Silences Selfish Genetic Elements in the Germline. Science 313 (5785), 320–324. doi:10.1126/science.1129333
Wang, W., Yoshikawa, M., Han, B. W., Izumi, N., Tomari, Y., Weng, Z., et al. (2014). The Initial Uridine of Primary piRNAs Does Not Create the Tenth Adenine that Is the Hallmark of Secondary piRNAs. Mol. Cel 56 (5), 708–716. doi:10.1016/j.molcel.2014.10.016
Weick, E.-M., and Miska, E. A. (2014). piRNAs: From Biogenesis to Function. Development 141 (18), 3458–3471. doi:10.1242/dev.094037
Whitfield, Z. J., Dolan, P. T., Kunitomi, M., Tassetto, M., Seetin, M. G., Oh, S., et al. (2017). The Diversity, Structure, and Function of Heritable Adaptive Immunity Sequences in the Aedes aegypti Genome. Curr. Biol. 27 (22), 3511–3519. doi:10.1016/j.cub.2017.09.067
Williams, Z., Morozov, P., Mihailovic, A., Lin, C., Puvvula, P. K., Juranek, S., et al. (2015). Discovery and Characterization of piRNAs in the Human Fetal Ovary. Cel Rep. 13 (4), 854–863. doi:10.1016/j.celrep.2015.09.030
Wyllie, A. H., Kerr, J. F. R., and Currie, A. R. (1980). Cell Death: The Significance of Apoptosis. Int. Rev. Cytol. 68, 251–306. doi:10.1016/s0074-7696(08)62312-8
Yashiro, R., Murota, Y., Nishida, K. M., Yamashiro, H., Fujii, K., Ogai, A., et al. (2018). Piwi Nuclear Localization and its Regulatory Mechanism in Drosophila Ovarian Somatic Cells. Cel Rep. 23 (12), 3647–3657. doi:10.1016/j.celrep.2018.05.051
Zhang, Z., Wang, J., Schultz, N., Zhang, F., Parhad, S. S., Tu, S., et al. (2014). The HP1 Homolog Rhino Anchors a Nuclear Complex that Suppresses piRNA Precursor Splicing. Cell 157 (6), 1353–1363. doi:10.1016/j.cell.2014.04.030
Zhang, Z., Xu, J., Koppetsch, B. S., Wang, J., Tipping, C., Ma, S., et al. (2011). Heterotypic piRNA Ping-Pong Requires Qin, a Protein with Both E3 Ligase and Tudor Domains. Mol. Cel 44 (4), 572–584. doi:10.1016/j.molcel.2011.10.011
Keywords: Manduca sexta, development, ping-pong amplification, transposon, RNA interference, small RNAs
Citation: Tsuji J, Thomson T, Brown C, Ghosh S, Theurkauf WE, Weng Z and Schwartz LM (2021) Somatic piRNAs and Transposons are Differentially Expressed Coincident with Skeletal Muscle Atrophy and Programmed Cell Death. Front. Genet. 12:775369. doi: 10.3389/fgene.2021.775369
Received: 13 September 2021; Accepted: 30 November 2021;
Published: 22 December 2021.
Edited by:
Chiara Gamberi, Coastal Carolina University, United StatesReviewed by:
Thomas Francis Duchaine, McGill University, CanadaTamás I. Orbán, Hungarian Academy of Sciences (MTA), Hungary
Copyright © 2021 Tsuji, Thomson, Brown, Ghosh, Theurkauf, Weng and Schwartz. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Lawrence M. Schwartz, LMA@umass.edu; Zhiping Weng, zhiping.weng@umassmed.edu
†Current address: Broad Institute, Cambridge, MA
‡Current address: MRC London Institute of Medical Sciences, Hammersmith Hospital Campus, London, United Kingdom