Respiratory syncytial virus infection changes the piwi-interacting RNA content of airway epithelial cells

Piwi-interacting RNAs (piRNAs) are small non-coding RNAs (sncRNAs) of about 26–32 nucleotides in length and represent the largest class of sncRNA molecules expressed in animal cells. piRNAs have been shown to play a crucial role to safeguard the genome, maintaining genome complexity and integrity, as they suppress the insertional mutations caused by transposable elements. However, there is growing evidence for the role of piRNAs in controlling gene expression in somatic cells as well. Little is known about changes in piRNA expression and possible function occurring in response to viral infections. In this study, we investigated the piRNA expression profile, using a human piRNA microarray, in human small airway epithelial (SAE) cells infected with respiratory syncytial virus (RSV), a leading cause of acute respiratory tract infections in children. We found a time-dependent increase in piRNAs differentially expressed in RSV-infected SAE cells. We validated the top piRNAs upregulated and downregulated at 24 h post-infection by RT-qPCR and identified potential targets. We then used Gene Ontology (GO) tool to predict the biological processes of the predicted targets of the most represented piRNAs in infected cells over the time course of RSV infection. We found that the most significant groups of targets of regulated piRNAs are related to cytoskeletal or Golgi organization and nucleic acid/nucleotide binding at 15 and 24 h p.i. To identify common patterns of time-dependent responses to infection, we clustered the significantly regulated expression profiles. Each of the clusters of temporal profiles have a distinct set of potential targets of the piRNAs in the cluster Understanding changes in piRNA expression in RSV-infected airway epithelial cells will increase our knowledge of the piRNA role in viral infection and might identify novel therapeutic targets for viral lung-mediated diseases.

Introduction

Respiratory syncytial virus (RSV) is a single-stranded RNA virus that belongs to the family of Pneumoviridae (Rima et al., 2017) and is the leading cause of acute lower respiratory tract infections in children worldwide, and it also results in severe respiratory infections in the elderly and immunocompromised patients (Shi et al., 2017; Shi et al., 2020). An effective treatment or vaccine is not currently available for RSV infection, and many molecular mechanisms regarding RSV-induced lung disease are still unknown (Mazur et al., 2015; Domachowske et al., 2021). Piwi-interacting RNAs (piRNAs) are a novel class of small non-coding RNAs (sncRNAs) which are between 26 and 32 nucleotides in length and have a 2′-O-methylated residue on their 3′ terminal base (Kirino and Mourelatos, 2007 (Fu and Wang, 2014). piRNAs can be grouped into four classes according to their origin: 1) transposon—derived piRNAs, coming from RNA transcripts of active transposable element copies; 2) piRNAs originated from the three untranslated regions (UTRs) of messenger RNAs (mRNAs); 3) non-coding RNAs derived piRNAs and 4) Caenorhabditis piRNAs found in worms with unique properties (Han and Zamore, 2014; Yang et al., 2016; Ernst et al., 2017), although their biogenesis is still not fully understood. piRNAs interact with the P-element Induced Wimpy Testis (PIWI) proteins, forming RNA-protein complexes known as piwi-interacting RNA complex (piRC). The piRC are involved in multiple cellular functions as gene regulation, chromatin modifications and transposon silencing in different organisms (Iwasaki et al., 2016; Monga and Banerjee, 2019). piRNAs were originally discovered in germ cells, maintaining the genome integrity and repressing mutations via transposable elements (Weick and Miska, 2014). Recently, they have been characterized in multiple organs and somatic cells (Zamore, 2010; Théron et al., 2014), and their unique expression profile is associated with the development of several diseases, including cancer, diabetes, cardiovascular and neurodegenerative disorders (Moyano and Stefani, 2015; Wakisaka and Imai, 2019; Wu et al., 2020; Rayford et al., 2021; Zeng et al., 2021). However, little is known about changes in piRNA expression and possible function occurring in response to viral infections. In our previous published work, we identified piRNAs among the sncRNAs present in extracellular vesicles altered by infection of airway epithelial cells with RSV (Chahar et al., 2018). In this study, we examined the changes in piRNA profile in primary small airway epithelial (SAE) cells infected with RSV, using a human piRNA microarray. We found that RSV infection of SAE cells was associated with a time-dependent increase in differentially expressed piRNAs, when compared to uninfected cells. Among the detected piRNAs, we validated six of the top ten upregulated at 24 h p. i. By RT-qPCR method. We used Gene Ontology (GO) analysis to predict the biological processes, molecular function and cell components of the most represented predicted targets of piRNAs over the three time points and we performed cluster analysis of the expression profiles. Identification of human piRNAs and their targets in airway epithelial cells will help to uncover new roles of this class of sncRNAs following viral infections, and to identify possible biomarkers and/or novel therapeutic targets for viral-mediated lung diseases.

Materials and methods

Small airway epithelial (SAE) culture and Respiratory syncytial virus infection

SAE cells (Lonza Inc., San Diego, CA), derived from terminal bronchioli of two different cadaveric donors (25- and 38-years old donors), were grown in growth medium, containing 7.5 mg/ml bovine pituitary extract (BPE), 0.5 mg/ml hydrocortisone, 0.5 μg/ml hEGF, 0.5 mg/ml epinephrine, 10 mg/ml transferrin, 5 mg/ml insulin, 0.1 μg/ml retinoic acid, 0.5 μg/ml triiodothyronine, 50 mg/ml gentamicin and 50 mg/ml bovine serum albumin. For preparation of RSV stocks, RSV Long strain was grown in Hep-2 cells and purified by centrifugation on discontinuous sucrose gradient as described previously (Ueba, 1978). SAE cells were switched to basal media (not supplements added) prior to RSV infection. At 90%–95% confluence, cell monolayers were infected with RSV at multiplicity of infection (MOI) of 3. An equivalent amount of 30% sucrose solution was added to uninfected SAE cells as a control (mock cells). Cells were collected at 6, 15 and 24 h RSV post-infection (p.i.) for the next analysis. We used two different donors and three replicates per time point.

Total RNA extraction and Piwi-interacting RNAs microarray

Total RNA was extracted using ToTALLY RNA kit from Ambion (Life Technologies). piRNA microarray analysis was performed by Go Beyond RNA Company (Arraystar Inc., Rockville, MD). RNA samples were quantified using a Nanodrop Spectrophotometer (Nanodrop Technologies) and integrity was assessed using standard denaturing agarose gel electrophoresis. The Arraystar HG19 piRNA array was designed for profiling of piRNAs in human cells. Human piRNAs were downloaded from the NCBI database and mapped to the HG19 genome sequence using UCSC Blat. DQ (accession ID) was used to search exact piRNA in the GenBank sequences, NCBI database (Wang et al., 2019). piRNAs with good match were selected, and on this array, probes of about 23,000 piRNAs were designed successfully using a duplex method. piRNA sample labeling was performed using an RNA ligase method. Briefly, 1 μg of each sample was 3′-end-labeled with Cy3 fluorescent label, using T4 RNA ligase using the following procedure: RNA in 2.0 μl of water was combined with 1.0 μl of CIP buffer and CIP (Exiqon). The mixture was incubated for 30 min (min) at 37°C, and was terminated by incubation for 5 min at 95°C. Then 3.0 μl of labeling buffer, 1.5 μl of fluorescent label (Cy3), 2.0 μl of DMSO, 2.0 μl of labeling enzyme were added into the mixture. The labeling reaction was incubated for 1 h at 16°C, and terminated by incubation for 15 min at 65°C. After stopping the labeling procedure, the Cy3-labeled samples were hybridized on the Arraystar Human piRNA Array. Hybridization was performed at 65°C for 18 h in Agilent’s SureHyb Hybridization Chambers. After being washed in an ozone-free environment, the slides were fixed and scanned using the Agilent DNA Microarray Scanner.

Data analysis of the Piwi-interacting RNAs microarray

Agilent Feature Extraction software (version 11.0.1.1) was used to analyze the acquired array images. Quantile normalization and subsequent data processing were performed using the GeneSpring GX v12.1 software package (Agilent Technologies). After quantile normalization of the raw data, piRNAs with Present or Marginal calls (microarray quality control flags) for at least four out of the 24 samples (six replicates for four time points) were chosen for further data analysis. Consistency between the two groups of three replicates at each time point was confirmed by visual inspection of the corresponding scatterplot. Differentially expressed piRNAs were identified through Fold Change filtering, combined with paired t-test across the six replicates. The t-test analysis was performed on log-transformed expression values. The use of logarithmic transformation is justified based on the skewness test (mean skewness of log transformed data k₃ = −0.3, which signifies only moderate deviation from normal distribution). Differentially expressed piRNAs between two groups were defined through a cut-off fold ratio, and significance, equivalent to Volcano Plot filtering. The resulting data were listed as fold ratios and p-values.

Predicted targets of Piwi-interacting RNAs and functional analysis

Prediction of potential targets of the piRNAs was performed by combining piRNA alignments with lists of genomic features, using the bedtools package. To annotate the potential targets of the identified piRNAs, we used the genomic coordinates of the piRNAs, as available from the piRNAQuest resource, http://bicresources.jcbose.ac.in/zhumur/pirnaquest/, (Ghosh et al., 2022). We used the UCSC liftOver tool (Hinrichs et al., 2006) and BEDTools (Quinlan and Hall, 2010) to identify overlaps of piRNA sequences with genomic features, as in the piRNAdb database. We verified the implementation by comparing selected results against piRNAdb. We used Panther version 17.0 to analyze the Gene Ontology (GO) annotations and enrichments within the significant predicted target genes for the differentially active piRNAs between each two timepoints (p value <0.01 for differential expression). The cutoff for calling enrichments in Panther was set as FDR <0.05 for Fisher exact test. We examined the three primary GO categories: biological process, molecular function, and cell component. Relative frequencies were calculated against the entire human genome as background.

Reverse transcription quantitative PCR (RT-qPCR)

Briefly, 600 ng total RNA containing the piRNA was reversed transcribed into cDNA in 1X RT buffer, 50 nM RT primers (Table 1), 250 nM dNTPs, 0.6 U RNase Inhibitor, and 3 U M-MuLV Reverse Transcriptase in a 20 μl reaction volume. The reaction was incubated at 16°C for 30°min, 42°C for 40°min, and 85°C for 5 min 20°μl of the synthesized cDNA was diluted in 100 μL water. The 2 μl cDNA dilution was mixed with 5 μl Arraystar SYBR^® Green Real-time qPCR Master Mix and 0.5 μl of 10 μM Forward and Reverse PCR primers (Table 2) in a 10 μL reaction volume. The qPCR was carried out in QuantStudio5 Real-time PCR System (Applied Biosystems) by one cycle of (95°C for 10 min) and 40 cycles of (95°C for 10°s, 60°C for 1°min, optical reading). U6 RNA was used as the housekeeping gene reference. ΔΔCt method was used to calculate the differential expression. Reverse transcription was performed using the First Strand cDNA Synthesis Kit (Arraystar). Quantitative PCR was conducted using the SYBR Green qPCR Master Mix (Arraystar), according to the manufacturer’s instructions. Primer sequences are available upon request. For the viral replication quantification, RSV N gene was amplified using specific primers, as previously described (Li et al., 2015).

TABLE 1

TABLE 1. Top 10 highly up- (≥2 fold) and down-regulated (≤2 fold) piRNAs by RSV 6 h p.i (p ≤ 0.01). DQ is the accession number of piRNA in the NCBI database.

TABLE 2

TABLE 2. Top 10 highly up- (≥4 fold) and down-regulated (≤4 fold) piRNAs by RSV 15 h p.i. (p ≤ 0.01). DQ is the accession number of piRNA in the NCBI database.

Statistical analysis

The statistical analysis of the piRNA microarray was run using the PDL:Stats and Statistics:PointEstimation modules from CPAN. Differentially expressed piRNAs were identified for every pair of timepoints (6 h versus mock, 15 h versus mock, 24 h versus mock). Fold change of the differentially expressed piRNA by microarray were calculated as linear average or Log2 average ratio. Fold change of RT-qPCR experiments was calculated by 2^-ΔΔCt method and represent mean ± SEM using GraphPad Prism v4 (GraphPad 160 Software). p value <0.05 or 0.01 was considered statistically significant.

Results

RSV infection induces changes in piRNAs profile expression

To investigate the impact of infection on the piRNA expression profile in primary airway epithelial cells, we analyzed piRNAs changes in human SAE cells infected with RSV for 6, 15 and 24 h p. i., using the Arraystar HG19 piRNA array, which probes for 23,677 different piRNAs. Total RNA from each sample was extracted and assessed for quality control. We measured viral repliation by detecting RSV N gene as shown in Supplementary Figure 1. Samples were labeled and hybridized to perform the piRNA microarray according to the company’s protocol. After quantile normalization and subsequent data processing, we detected 548 (6 h p.i.), 897 (15 h p.i.) and 1,644 (24 h p.i.) differentially expressed piRNAs in RSV-infected cells compared to mock (uninfected) cells. A total of 120, 321 and 772 piRNAs were significantly upregulated at 6, 15 and 24 h p. i., respectively. Interestingly, there was a higher number of downregulated piRNAs in response to viral infection, with 428, 576 and 872 differentially expressed piRNAs at 6, 15 and 24 h p. i., respectively (Figure 1). Among the differentially expressed piRNAs, 244 were shared between the time points of 6 and 15 h p.i.; 221 between 6 and 24 h p.i. And 515 between 15 and 24 h p.i., with 157 piRNAs commonly detected over all time points, as represented in the intersecting Venn diagram (Figure 2). The top ten piRNAs upregulated and downregulated (≤2 Fold) for each time point p. i. Are shown in Table 1 (6 h), Table 2 (15 h) and Table 3 (24 h). The complete list of piRNAs whose expression in SAE cells was changed in a time-dependent manner after RSV infection is shown in Supplementary files 1, two and 3.

FIGURE 1

FIGURE 1. Numbers of upregulated (black bars) and downregulated (grey bars) piRNAs in RSV infected SAE cells at 6, 15 and 24 h p. i. Compared to the mock (uninfected) SAE cells. Data is representative of six replicates for each time point with p < 0.01.

FIGURE 2

FIGURE 2. Intersecting Venn diagram showing the numbers of up- and down-regulated piRNAs differentially expressed at the three different time points (6, 15 and 24 h) RSV p. i. Compared to the mock (uninfected) SAE cells. The numbers in overlapping area(s) represent the differential regulated number of piRNAs shared among all the time points. Data is representative of six replicates for each time point with p < 0.01.

TABLE 3

TABLE 3. Top 10 highly up- (≥4 fold) and down-regulated (≤4 fold) piRNAs by RSV 24 h p.i. (p ≤ 0.01). DQ is the accession number of piRNA in the NCBI database.

We then validated the expression of six or eight of the top ten upregulated or downregulated piRNAs at 24 h p.i. Listed in Table 3. Total RNA was isolated from SAE mock or RSV infected cells at 24 h p.i. And subjected to RT-qPCR. We found that piR-54651 (DQ587539), piR-32372 (DQ582260) and piR-61160 (DQ595048) levels increased more than threefold, and the other three piRNAs, piR-61505 (DQ595393), piR-37886 (DQ599820) and piR-50722 (DQ583610) increased one and half fold in viral-infected cells compared to mock (uninfected) cells, as shown in Figure 3A. We observed more than two fold decreased levels of piR-46025 (DQ577913) and piR-38945 (DQ600879), and the other two piRNAs, piR-41627 (DQ573515) and piR-30049 (DQ5699370) showed more than one and half fold of decreased levels in viral-infected cells compared to mock (uninfected) cells Figure 3B. In general, fold change values of piRNAs analyzed by microarray were at least two times higher than the fold change of the same piRNAs validated by RT-qPCR.

FIGURE 3

FIGURE 3. Validation of up (A) or down (B) regulated piRNAs expression in SAE cells infected with RSV at 24 h p. i. RNA extracted from SAE mock (uninfected) and RSV infected cells (24 h) was subjected to piRNAs analysis by RT-qPCR. Fold changes in piRNA expression were determined by 2^−ΔΔCT method and represent mean ± SEM normalized to a small nuclear RNA, U6. * and ** indicates a statistical difference comparing RSV infected versus mock cells (*p value <0.05; **p value <0.01). Data is representative of three independent experiments.

Potential targets and Gene Ontology prediction of the viral-induced up- and down-regulated piRNAs

To characterize potential targets of piRNAs identified in RSV-infected SAE cells for each time point, we combined the piRNA alignments with lists of genomic features, using the bedtools package, and verified the approach by comparing selected results against the piRNAdb database (Piuco and Galante, 2021). At the stricter threshold for piRNA regulation at p < 0.01, we identified a total of 629, 1,044 and 2,750 targets, including 519, 874, and 2,118 targets of down-regulated piRNAs, and 153, 245 and 886 targets of up-regulated piRNAs, all respectively at 6, 15 and 24 h p. i. With RSV. Using the inclusive approach, we found a total of 2,874, 2,341, and 6,967 potential targets of regulated piRNAs at 6, 15 and 24 h p. i. Compared to mock (uninfected cells) with a p value <0.05, respectively. Of these, 1961, 1978, and 5,778 are targets of down-regulated piRNAs, and 1,313, 568, and 2,684 are targets of up-regulated piRNAs at 6, 15 and 24 h p. i. We found that most of the piRNAs had multiple target genes and that several piRNAs shared common gene targets. Potential targets of selected upregulated piRNAs following RSV infection are listed in Table 4. Most of the predicted targets belong to the class of long noncoding RNAs (lncRNAs), including long intergenic non-coding RNAs (lincRNAs) such as LINC01297, LINC01087 and LINC01087. LincRNAs are autonomously transcribed RNAs of more than 200 nucleotides in length that do not overlap protein-coding genes (Ransohoff et al., 2018). Other identified targets of the upregulated piRNAs were the endosome protein Pleckstrin Homology Domain Containing B2 (PLEKHB2), and the cytoplasm proteins the Adaptor Related Protein Complex five Subunit Zeta 1 (AP5Z1), the Probable ATP-dependent RNA helicase (DDX4), and the Signal transducer and activator of transcription 2 (STAT2). PLEKHB2 is involved in retrograde transport of recycling endosomes. AP5Z1 belongs to the family of Clathrin adaptor proteins, known also as adaptin proteins (APs), and participates in the processes of homologous recombination DNA double-strand break repair (HR-DSBR). DDX4, highly expressed in the ovary and testis, has been involved in the production of piRNAs in fetal male germ cells through a ping-pong amplification cycle, and in the secondary piRNAs metabolic process (Li et al., 2010). STAT2 protein is critical to the biological response of type I interferons (IFNs) with a DNA-binding transcription factor activity (Park et al., 2000). The complete lists of the potential target genes for each time point of RSV infection are presented in Supplementary files 1 (6 h p.i.), 2 (15 h p.i.) and 3 (24 h p.i.).

TABLE 4

TABLE 4. Predicted target genes of selected upregulated piRNAs during the course of RSV infection.

To characterize the prevalent functions of the predicted targets of the differentially expressed piRNAs we use the Gene Ontology (GO) functional classification (biological process, cellular component, and molecular function). We found that the most significant groups of targets of regulated piRNAs are related to cytoskeletal or Golgi organization and nucleic acid/nucleotide binding at 15 and 24 h p.i. (p < 0.01). A summary of the GO biological process of the predict targets of piRNAs upregulated at 24 h RSV p.i (p < 0.05) is listed in Table 5. The most enriched biological processes at 24 h p.i. included endosomal, intracellular and nitrogen compound transports with regulation of vesicles-mediated transport. The GO analysis also showed Golgi organelle and cytoskeleton organizations and interestingly immune response mediated by cells involved in somatic diversification process (immunological memory). Complete GO analysis of biological process, cellular component, and molecular function (BP, CC and MF) for all three time points is provided as supplementary files 4, 5 and 6 (6 h p.i.), files 7, 8 and 9 (15 h p. i.), and files 10, 11 and 12 (24 h p.i.).

TABLE 5

TABLE 5. Gene Ontology (GO), biological processes (BP) analysis for the top gene targets of piRNAs (time point, RSV 24 h p.i) with p values <0.01 or <0.05

Clustering of piRNA expression profiles

To identify common patterns of time-dependent responses to infection, we clustered the significantly regulated expression profiles. Specifically, we used the piRNAs with a significant change at p < 0.01 in at least one of the time points, p.i., compared to control. The expression profiles were converted to the log2 scale, averaged over the six replicates, and normalized to levels in uninfected cells and root mean square deviation (rmsd) was equal to 1, rmsd = 1. The preprocessed data were clustered using an in-house implementation of the Hartigan-Wong k-Means algorithm. After conducting an elbow test as shown in Figure 4, we assessed that the most informative clustering was either into three or five groups of temporal profiles. The results of the clustering into five groups are shown in Figure 5. The five basic types of behavior are early (6 h) upregulation, followed by remaining in the upregulated, or uninfected-base level (cluster 0), early downregulation, remaining downregulated (cluster 1), steady downregulation (cluster 2), steady upregulation (cluster 3), and early upregulation, followed by strong downregulation (cluster 4). The complete list of piRNAs for both type of clustering, three or five groups, is included in the supplementary files 13 and 14 named “clustered_list_3. xlsx, and clustered_list_5. xlsx”.

FIGURE 4

FIGURE 4. The elbow test indicates that the most informative representations of the categories of temporal profiles of the differentially regulated piRNAs are classifications into three clusters or into five clusters. X-axis: number of clusters, Y-axis: detrended (residual) Total Within Sum of Squares, a measure of structure unexplained by the clustering.

FIGURE 5

FIGURE 5. The k-means clustering of the temporal profiles of X-axis: time post-infection with RSV (cluster zero represents uninfected cells). Y-axis: normalized expression profile, averaged over the replicates, and shown in the log₂ scale.

Each of the clusters of temporal profiles have a distinct set of potential targets of the piRNAs in the cluster. Moreover, the targets are associated with enrichment of different functional categories that can be expressed in terms of Gene Ontology biological process enrichments. For example, the early upregulated cluster 0 is enriched with targets with a role in protein transport and localization. The targets of early downregulated piRNAs (cluster 1) are enriched in organelle organization, including Golgi, cytoplasm, and membranes, but depleted in genes coding for regulatory proteins. The piRNAs steadily upregulated (Cluster 3) have potential targets enriched in regulatory genes and inter-cellular signaling. The full lists of the significantly enriched/depleted GO biological process categories for the five clusters are included in the supplementary file 15. We analyzed if specific clusters shown in Figure 5 correlate with subsets of GO categories. The overlaps, for GO BP categories significant at p value <0.01 are listed in Table 6. We found that the early (cluster 0) and the steadily upregulated (cluster 3) piRNAs, for example, were linked to GO subcategories regulating the mediated signal transduction of small GTPase, typically involved in the viral immune response. Many cellular functions such as protein transport, intracellular or membrane protein localization, etc. GO subsets were highly associated with targets of piRNAs from the steady downregulated level (cluster 2), which was also associated with adaptive immune response subset. The piRNA targets from the steadily downregulated or upregulated clusters (2 and three respectively) were interestingly connected to GO subcategories involved in the regulation of cellular signaling or biological quality.

TABLE 6

TABLE 6. Correlation of Gene Ontology (GO) biological processes (BP) subset categories with the clustering of the temporal profiles of X-axis.

Discussion

The present study is the first to report changes in piRNAs expression in response to RSV infection of human primary airway epithelial cells. piRNAs are the largest class of sncRNA molecules expressed in animal cells and they can be derived from multiple sources, including transposons, mRNAs, and ncRNAs, such as lncRNAs, snoRNAs and tRFs. Among these piRNAs, only the function of transposon-derived piRNAs in gonadal cells has been relatively well studied (Weick and Miska, 2014). piRNAs have been shown to play a crucial role to safeguard genome, maintaining the genome complexity and integrity, as they suppress the insertional mutations caused by transposable elements, however, there is growing evidence of their role in controlling gene expression in somatic cells as well. The biogenesis and function of piRNAs is associated with three different proteins, namely PIWI, Argonaute-3 (AGO3), and Aubergine (AUB). PIWI proteins are found both in somatic cells and germ cells, while AGO3 and AUB are observed primarily in germ cells. Precursors of piRNAs are exported from the nucleus into the cytoplasm to be further processed into smaller sequences and reach their partners to form piRNA-PIWI complexes. These complexes can migrate back to the nucleus to block the transcription of target genes. There are several mechanisms underlying piRNA-mediated gene regulation. The better characterized one relates to epigenetic changes, including DNA methylation and histone modifications, whereas piRNA–PIWI complex can bind to the promoter region of a gene and favors recruitment of methylation factors, enabling transcriptional silencing. Some examples of this mode of gene regulation in somatic cells are the silencing of killer immunoglobulin-like receptor, the immune receptor CD1A and the antioxidant enzyme glutathione-S-transferase (Zhang et al., 2016; Zhang et al., 2020). In addition, piRNAs can also regulate gene expression at the post-transcriptional level, as in the case of piRNA-30840 which binds to IL-4 pre-mRNA, leading to its degradation (Zhong et al., 2015), and piR-FTH1, a piRNA found in breast cancer cells, which downregulates ferritin heavy chain mRNA (Balaratnam et al., 2018). One last mechanism by which piRNAs regulate gene expression has been recently identified and it related to the interaction with the cellular RNA methylation machinery, leading to changes in mRNA m6A methylation, a post-transcriptional modification of mRNAs in mammals that plays a fundamental role in RNA stability and translation (Gao et al., 2020).

Although the role of piRNAs was initially confined to gonad development, recent studies have identified specific expression profile of piRNAs from multiple organs and cell types. Currently, piRNA-disease association data are available for several diseases, including various types of cancers, cardiovascular diseases, autoimmune and neurodegenerative disorders (Muhammad et al., 2019). Surprisingly, almost nothing is known regarding piRNA generation/function and PIWI proteins in the context of viral infections. In our recent published study, analysis of RNA content of exosomes isolated from RSV-infected airway epithelial cells identified piRNAs as an important component of the exosomal cargo (Chahar et al., 2018) Next-generation sequencing analysis revealed that 52 piRNAs were commonly present in both mock and RSV exosomes, with 28 upregulated and 24 downregulated piRNAs in RSV exosomes, with 3 and 19 piRNAs that uniquely present in mock and RSV exosomes, suggesting a biological function for this class of sncRNAs during infection. In the current study, we show for the first time that RSV infection can lead to significant changes in piRNAs in infected SAE cells in a time-dependent manner. Our piRNA microarray analysis revealed an enrichment in numbers of up- and down-regulated piRNAs in airway epithelial cells upon RSV infection compared to mock cells, with the highest numbers of differentially expressed piRNAs detected at 24 h p. i. With RSV. Significant changes in the expression profiles of piRNAs was recently reported for Coxsackievirus B3 (CVB3) infection in human HeLa cells using high-throughput sequencing. Similar to RSV, CVB3 infection was associated with time-dependent changes in piRNAs numbers, although the highest number of differentially expressed piRNAs was observed at the early time points post-infection (Yao et al., 2020). In an animal model of ducks infected with avian influenza virus, piRNAs represented the highest class of sncRNAs expressed after infection compared to other snRNA classes, such as microRNAs (miRNAs) and small nucleolar RNAs (snoRNAs) (Samir et al., 2020). There were major organ-specific differences in sncRNA populations at baseline and substantial reprogramming of all sncRNA classes throughout infection.

Target prediction analysis suggested that piR-38587 could possibly target AP5Z1, the adaptin proteins 5 (AP5). Very little is known on AP5 because it is the most recently identified adaptin complex. Although AP5 has not been investigated in the context of viral infection, AP2, another type of AP complex has been shown to function as a target of ocular infection upon human adenovirus species D type 37 (HAdV-D37) infection (Lee et al., 2020) or be involved in the entry process of human enterovirus 71 (HEV71) (Hussain et al., 2011). Similarly, piR-4095 target includes DDX4 which is an RNA helicases enzyme. Different DDX enzymes, as DDX1, DDX3X, and DDX5, have been known to play important roles in the replication cycle of coronaviruses (Taschuk and Cherry, 2020). piR-57640 could target STAT2 which is involved in RSV infection. It is known that RSV is a potent inducer of type I interferon release from infected airway epithelial cells (Garofalo et al., 1996), and the virus specifically downregulates human STAT2 protein expression, thus enabling RSV to evade the host type I interferon response (Elliott et al., 2007). The Gene Ontology functional classification found that the most significant groups of targets of regulated piRNAs are related to long noncoding RNA, DNA- and protein-binding, transcriptional and post-transcriptional regulation.

Although piRNA function in response to an infection in mammalian cells has not been explored yet, murine neural stem cells (NSCs) were recently shown to release piRNA-containing exosomes/microvesicles with antiviral immunity (Ikhlas et al., 2022). When NSCs were exposed to RNA fragments of SARSCoV-2 genome, the derived extracellular vesicles displayed enhanced antiviral activity. The increased antiviral effect of these vesicles was associated with increased expression of piRNA species some of which were predicted to target SARS-CoV-2 genome. Knockout of piRNA-interacting protein PIWIL2 in NSCs led to a reduction in the induced antiviral effect, suggesting that the PIWI-piRNA system is important for these antiviral actions (Ikhlas et al., 2022). A recent study of bacterial pneumonia using a knockout mouse for MIWI2, the mouse homologue of human PIWIL4, revealed an important role of this piwi protein in lung innate immune response. MIWI2 was expressed by a subpopulation of multiciliated airway epithelial cells, and its expression was increased following infection with S. pneumoniae (Wasserman et al., 2017). RNA sequencing found that MIWI2-positive cells expressed a significantly different transcriptome compared with neighboring cells that did not express MIWI2, and MIWI2 −/− was associated with enhanced innate immune responses and faster bacterial clearance (Wasserman et al., 2017).

In conclusion, our study provides insights into changes of the expression profile of piRNAs during RSV infection and highlight the need for more detailed investigations of the piRNA/PIWI protein pathway in the context of viral infection to establish their potential role in regulation of cellular responses. A better understanding of this pathway could potentially identify novel therapeutic targets for future piRNA-mediated strategies to modulate RSV-associated lung disease, as well as potential biomarkers for disease diagnosis and prognosis.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

TC designed and performed the experiments, analyzed the data, wrote the manuscript; AK. contributed to piRNA data analysis and drafting of the manuscript; TL was the technical support of the cell culture. AC contributed to the conception of the manuscript, literature review, financial support, and final manuscript editing.

Funding

This research was funded by NIH grant AI062885, AI125434, AI142570. TC was supported in part by the Jeane B. Kempner Postdoctoral Scholar Award from UTMB.

Acknowledgments

The authors would like to thank Cynthia Tribble for assistance with submission of the manuscript.

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.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2022.931354/full#supplementary-material

References

Balaratnam, S., West, N., and Basu, S. (2018). A piRNA utilizes HILI and HIWI2 mediated pathway to down-regulate ferritin heavy chain 1 mRNA in human somatic cells. Nucleic Acids Res. 46 (20), 10635–10648. 5068339[pii][doi]. doi:10.1093/nar/gky728

PubMed Abstract | CrossRef Full Text | Google Scholar

Chahar, H. S., Corsello, T., Kudlicki, A. S., Komaravelli, N., and Casola, A. (2018). Respiratory syncytial virus infection changes cargo composition of exosome released from airway epithelial cells. Sci. Rep. 8 (1), 387. doi:10.1038/s41598-017-18672-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Domachowske, J. B., Anderson, E. J., and Goldstein, M. (2021). The future of respiratory syncytial virus disease prevention and treatment. Infect. Dis. Ther. 10 (1), 47–60. doi:10.1007/s40121-020-00383-6

CrossRef Full Text | Google Scholar

Elliott, J., Lynch, O. T., Suessmuth, Y., Qian, P., Boyd, C. R., Burrows, J. F., et al. (2007). Respiratory syncytial virus NS1 protein degrades STAT2 by using the Elongin-Cullin E3 ligase. J. Virol. 81 (7), 3428–3436. doi:10.1128/JVI.02303-06

PubMed Abstract | CrossRef Full Text | Google Scholar

Ernst, C., Odom, D. T., and Kutter, C. (2017). The emergence of piRNAs against transposon invasion to preserve mammalian genome integrity. Nat. Commun. 8 (1), 1411. doi:10.1038/s41467-017-01049-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, Q., and Wang, P. J. (2014). Mammalian piRNAs: Biogenesis, function, and mysteries. Spermatogenesis 4, e27889. doi:10.4161/spmg.27889

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, X. Q., Zhang, Y. H., Liu, F., Ponnusamy, M., Zhao, X. M., Zhou, L. Y., et al. (2020). The piRNA CHAPIR regulates cardiac hypertrophy by controlling METTL3-dependent N(6)-methyladenosine methylation of Parp10 mRNA. Nat. Cell Biol. 22 (11), 1319–1331. doi:10.1038/s41556-020-0576-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Garofalo, R., Mei, F., Espejo, R., Ye, G., Haeberle, H., Baron, S., et al. (1996). Respiratory syncytial virus infection of human respiratory epithelial cells up-regulates class I MHC expression through the induction of IFN-beta and IL-1 alpha. J. Immunol. 157 (6), 2506–2513.

PubMed Abstract | Google Scholar

Ghosh, B., Sarkar, A., Mondal, S., Bhattacharya, N., Khatua, S., and Ghosh, Z. (2022). piRNAQuest V.2: an updated resource for searching through the piRNAome of multiple species. RNA Biol. 19 (1), 12–25. doi:10.1080/15476286.2021.2010960

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, B. W., and Zamore, P. D. (2014). piRNAs. Curr. Biol. 24 (16), R730–R733. doi:10.1016/j.cub.2014.07.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Hinrichs, A. S., Karolchik, D., Baertsch, R., Barber, G. P., Bejerano, G., Clawson, H., et al. (2006). The UCSC genome browser database: Update 2006. Nucleic Acids Res. 34, D590–D598. doi:10.1093/nar/gkj144

PubMed Abstract | CrossRef Full Text | Google Scholar

Hussain, K. M., Leong, K. L., Ng, M. M., and Chu, J. J. (2011). The essential role of clathrin-mediated endocytosis in the infectious entry of human enterovirus 71. J. Biol. Chem.. 286 (1), 309–321. doi:10.1074/jbc.M110.168468

PubMed Abstract | CrossRef Full Text | Google Scholar

Ikhlas, S., Usman, A., Kim, D., and Cai, D. (2022). Exosomes/microvesicles target SARS-CoV-2 via innate and RNA-induced immunity with PIWI-piRNA system. Life Sci. Alliance 5 (3), e202101240. doi:10.26508/lsa.202101240

PubMed Abstract | CrossRef Full Text | Google Scholar

Iwasaki, Y. W., Murano, K., Ishizu, H., Shibuya, A., Iyoda, Y., Siomi, M. C., et al. (2016). Piwi modulates chromatin accessibility by regulating multiple factors including histone H1 to repress transposons. Mol. Cell. 63 (3), 408–419. doi:10.1016/j.molcel.2016.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J. S., Mukherjee, S., Lee, J. Y., Saha, A., Chodosh, J., Painter, D. F., et al. (2020). Entry of epidemic keratoconjunctivitis-associated human adenovirus type 37 in human corneal epithelial cells. Invest. Ophthalmol. Vis. Sci.. 61 (10), 50. doi:10.1167/iovs.61.10.50

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H. J., Yu, N., Zhang, X. Y., Jin, W., and Li, H. Z. (2010). Spermatozoal protein profiles in male infertility with asthenozoospermia. Chin. Med. J. 123 (20), 2879–2882.

PubMed Abstract | Google Scholar

Li, H., Ma, Y., Escaffre, O., Ivanciuc, T., Komaravelli, N., Kelley, J. P., et al. (2015). Role of hydrogen sulfide in paramyxovirus infections. J. Virol. 89 (10), 5557–5568. doi:10.1128/JVI.00264-15

PubMed Abstract | CrossRef Full Text | Google Scholar

Mazur, N. I., Martinón-Torres, F., Baraldi, E., Fauroux, B., Greenough, A., Heikkinen, T., et al. (2015). Lower respiratory tract infection caused by respiratory syncytial virus: current management and new therapeutics. Lancet Respir. Med. 3 (11), 888–900. doi:10.1016/s2213-2600(15)00255-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Monga, I., and Banerjee, I. (2019). Computational identification of piRNAs using features based on RNA sequence, structure, thermodynamic and physicochemical properties. Curr. Genomics. 20 (7), 508–518. doi:10.2174/1389202920666191129112705

PubMed Abstract | CrossRef Full Text | Google Scholar

Moyano, M., and Stefani, G. (2015). piRNA involvement in genome stability and human cancer. J. Hematol. Oncol. 8 38. doi:10.1186/s13045-015-0133-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Muhammad, A., Waheed, R., Khan, N. A., Jiang, H., and Song, X. (2019). piRDisease v1.0: a manually curated database for piRNA associated diseases. Database. (Oxford) 2019, baz052. doi:10.1093/database/baz052

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, C., Li, S., Cha, E., and Schindler, C. (2000). Immune response in Stat2 knockout mice. Immunity 13 (6), 795–804. doi:10.1016/s1074-7613(00)00077-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Piuco, R., and Galante, P. A. F. (2021). piRNAdb: A piwi-interacting RNA database. bioRxiv. 2021.2009.2021.461238. doi:10.1101/2021.09.21.461238

CrossRef Full Text | Google Scholar

Quinlan, A. R., and Hall, I. M. (2010). BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26 (6), 841–842. btq033 [pii][doi]. doi:10.1093/bioinformatics/btq033

PubMed Abstract | CrossRef Full Text | Google Scholar

Ransohoff, J. D., Wei, Y., and Khavari, P. A. (2018). The functions and unique features of long intergenic non-coding RNA. Nat. Rev. Mol. Cell Biol. 19 (3), 143–157. doi:10.1038/nrm.2017.104

PubMed Abstract | CrossRef Full Text | Google Scholar

Rima, B., Collins, P., Easton, A., Fouchier, R., Kurath, G., Lamb, R. A., et al. (2017). ICTV Virus Taxonomy Profile: Pneumoviridae. J. Gen. Virol. 98 (12), 2912–2913. doi:10.1099/jgv.0.000959

PubMed Abstract | CrossRef Full Text | Google Scholar

Rayford, K. J., Cooley, A., Rumph, J. T., Arun, A., Rachakonda, G., Villalta, F., et al. (2021). piRNAs as modulators of disease pathogenesis. Int. J. Mol. sci.. 22 (5), 2373. doi:10.3390/ijms22052373

PubMed Abstract | CrossRef Full Text | Google Scholar

Samir, M., Vidal, R. O., Abdallah, F., Capece, V., Seehusen, F., Geffers, R., et al. (2020). Organ-specific small non-coding RNA responses in domestic (Sudani) ducks experimentally infected with highly pathogenic avian influenza virus (H5N1). RNA Biol.. 17 (1), 112 –124. doi:10.1080/15476286.2019.1669879

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, T., Denouelt, A., Tietjen, A. K., Campbell, I., Moran, E., Li, X., et al. (2020). Global Disease Burden Estimates of Respiratory Syncytial Virus-Associated Acute Respiratory Infection in Older Adults in 2015: A Systematic Review and Meta-Analysis. J. Infect. Dis. 222 (Suppl 7), S577–S583. doi:10.1093/infdis/jiz059

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, T., McAllister, D. A., O'Brien, K. L., Simoes, E. A. F., Madhi, S. A., Gessner, B. D., et al. (2017). Global, regional, and national disease burden estimates of acute lower respiratory infections due to respiratory syncytial virus in young children in 2015: A systematic review and modelling study. Lancet. 390 (10098), 946–958. doi:10.1016/s0140-6736(17)30938-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Taschuk, F., and Cherry, S. (2020). DEAD-box helicases: Sensors, regulators, and effectors for antiviral defense. Viruses 12 (2), 181. doi:10.3390/v12020181

PubMed Abstract | CrossRef Full Text | Google Scholar

Théron, E., Dennis, C., Brasset, E., and Vaury, C. (2014). Distinct features of the piRNA pathway in somatic and germ cells: From piRNA cluster transcription to piRNA processing and amplification. Mob. DNA. 5 (1), 28. doi:10.1186/s13100-014-0028-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Ueba, O. (1978). Respiratory syncytial virus. I. Concentration and purification of the infectious virus. Acta Med. Okayama. 32 (4), 265–272.

PubMed Abstract | Google Scholar

Wakisaka, K. T., and Imai, Y. (2019). The dawn of pirna research in various neuronal disorders. Front. Biosci. (Landmark Ed.). 24 1440–1451. doi:10.2741/4789

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Zhang, P., Lu, Y., Li, Y., Zheng, Y., Kan, Y., et al. (2019). piRBase: A comprehensive database of piRNA sequences. Nucleic. Acids Res. 47 (D1), D175–D180. doi:10.1093/nar/gky1043

PubMed Abstract | CrossRef Full Text | Google Scholar

Wasserman, G. A., Szymaniak, A. D., Hinds, A. C., Yamamoto, K., Kamata, H., Smith, N. M., et al. (2017). Expression of Piwi protein MIWI2 defines a distinct population of multiciliated cells. J. Clin. Invest. 127 (10), 3866–3876. doi:10.1172/JCI94639

PubMed Abstract | CrossRef Full Text | Google Scholar

Weick, E. M., and Miska, E. A. (2014). piRNAs: from biogenesis to function. Development 141 (18), 3458–3471. doi:10.1242/dev.094037

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, X., Pan, Y., Fang, Y., Zhang, J., Xie, M., Yang, F., et al. (2020). The Biogenesis and Functions of piRNAs in Human Diseases. Mol. Ther. Nucleic. Acids. 21 108–120. doi:10.3390/ijms22052373

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J. X., Rastetter, R. H., and Wilhelm, D. (2016). Non-coding RNAs: An Introduction. Adv. Exp. Med. Biol. 886, 13–32. doi:10.1007/978-94-017-7417-8_2

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, H., Wang, X., Song, J., Wang, J., Song, Q., and Han, J. (2020). Coxsackievirus B3 infection induces changes in the expression of numerous piRNAs. Arch. Virol.. 165 (1), 105–114. doi:10.1007/s00705-019-04451-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Zamore, P. D. (2010). Somatic piRNA biogenesis. The EMBO J. 29 (19), 3219–3221. doi:10.1038/emboj.2010.232

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, Q., Cai, J., Wan, H., Zhao, S., Tan, Y., Zhang, C., et al. (2021). PIWI-interacting RNAs and PIWI proteins in diabetes and cardiovascular disease: Molecular pathogenesis and role as biomarkers. Clin. Chim. Acta.. 518, 33–37. doi:10.1016/j.cca.2021.03.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, L., Meng, X., Pan, C., Qu, F., Gan, W., Xiang, Z., et al. (2020). piR-31470 epigenetically suppresses the expression of glutathione S-transferase pi 1 in prostate cancer via DNA methylation. Cell. Signal. 67, 109501. doi:10.1016/j.cellsig.2019.109501

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., He, X., Liu, C., Liu, J., Hu, Q., Pan, T., et al. (2016). IL-4 inhibits the biogenesis of an epigenetically suppressive PIWI-interacting RNA to upregulate CD1a molecules on monocytes/dendritic cells. J. Immunol. 196 (4), 1591–1603. doi:10.4049/jimmunol.1500805

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, F., Zhou, N., Wu, K., Guo, Y., Tan, W., Zhang, H., et al. (2015). A SnoRNA-derived piRNA interacts with human interleukin-4 pre-mRNA and induces its decay in nuclear exosomes. Nucleic Acids Res. 43 (21), 10474–10491. doi:10.1093/nar/gkv954

PubMed Abstract | CrossRef Full Text | Google Scholar