Transcriptome-wide 5-methylcytosine modification profiling of long non-coding RNAs in A549 cells infected with H1N1 influenza A virus
BMC Genomics volume 24, Article number: 316 (2023)
In recent years, accumulating evidences have revealed that influenza A virus (IAV) infections induce significant differential expression of host long noncoding RNAs (lncRNAs), some of which play important roles in the regulation of virus-host interactions and determining the virus pathogenesis. However, whether these lncRNAs bear post-translational modifications and how their differential expression is regulated remain largely unknown. In this study, the transcriptome-wide 5-methylcytosine (m5C) modification of lncRNAs in A549 cells infected with an H1N1 influenza A virus was analyzed and compared with uninfected cells by Methylated RNA immunoprecipitation sequencing (MeRIP-Seq).
Our data identified 1317 upregulated m5C peaks and 1667 downregulated peaks in the H1N1 infected group. Gene ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses showed that the differentially modified lncRNAs were associated with protein modification, organelle localization, nuclear export and other biological processes. Furthermore, conjoint analysis of the differentially modified (DM) and differentially expressed (DE) lncRNAs identified 143 ‘hyper-up’, 81 ‘hypo-up’, 6 ‘hypo-down’ and 4 ‘hyper-down’ lncRNAs. GO and KEGG analyses revealed that these DM and DE lncRNAs were predominantly associated with pathogen recognition and disease pathogenesis pathways, indicating that m5C modifications could play an important role in the regulation of host response to IAV replication by modulating the expression and/or stability of lncRNAs.
This study presented the first m5C modification profile of lncRNAs in A549 cells infected with IAV and demonstrated a significant alteration of m5C modifications on host lncRNAs upon IAV infection. These data could give a reference to future researches on the roles of m5C methylation in virus infection.
RNAs are subject to a range of covalent modifications at the single nucleotide level, such as pseudouridine, N6-methyladenosine (m6A) , N1-methyladenosine (m1A) , 5-methylcytosine (m5C) . So far, over 100 distinct RNA modifications have been described. As a widespread modification of RNA, 5-methylcytosine (m5C) has received considerable attention in recent years. Increasing numbers of studies have demonstrated that RNA m5C modification plays important roles in multiple biological processes including RNA processing [4, 5], RNA stability [6, 7], RNA transport , and mRNA translation . With the development of high-throughput sequencing, it is now possible to identify and quantify m5C modifications in low-abundance RNA species such as non-coding RNAs (ncRNAs), and transcriptome-wide identifications of m5C modification in different types of RNAs have been reported [10,11,12].
Evidence is emerging that RNA modifications are involved in the regulation of virus infection. Epitranscriptomic marks such as m6A, m5C, N4-acetylcytidine (ac4C) and 2ʹO-methylated nucleosides (Nm) have been reported to promote viral replication by upregulating viral mRNA stability or translation, or by preventing the recognition of viral RNA through modulation of host RNA-specific innate immunity factors including RIG-I and MDA5 . As a relatively common epitranscriptomic mRNA modification, m5C has been found to present at higher levels in retroviral transcripts than in cellular mRNAs, and the modification can regulate RNA splicing and promote the translation of viral mRNAs . A more recent study has demonstrated that RNA m5C methylation can control antiviral innate immunity through modulating the m5C methylome of noncoding RNAs (ncRNAs) and their expression, which regulate type I interferons .
As an important pathogen for seasonal respiratory illness, influenza A viruses (IAV) have caused human pandemic outbreaks such as those occurred in 1918, 1957 and 1968, and they are still continuing to threaten public health [16, 17]. IAV primarily infects respiratory epithelial cells and causes pulmonary diseases. If uncontrolled, the infection can cause loss of lung function and even mortality [18, 19]. Understanding the process of IAV infection is critical to target IAV-induced pathogenesis and develop effective anti-viral approaches. Recent studies have revealed that host long noncoding RNAs (lncRNAs) are key regulators of host-virus interactions during viral infection [20, 21]. Thousands of lncRNAs have been identified to be differentially expressed during IAV infection [22,23,24]. Some of them affect IAV infection by regulating the host innate immune responses [25,26,27,28,29,30], modulating cellular metabolism  or directly interacting with viral proteins [32, 33].
Epitranscriptomic regulation has been shown as an important mechanism to control lncRNAs expression and tissue specificity . Although m5C is a common epitranscriptomic modification found in RNAs, knowledge surrounding the prevalence and transcriptome-wide distribution of m5C in lncRNA is still very limited, and the roles of m5C modification during IAV infection have not yet been explored. It remains unknown whether m5C modification plays a role in regulating lncRNA expression during IAV infection.
In this study, the m5C methylation of lncRNAs in H1N1-infected A549 cells were globally mapped by Methylated RNA immunoprecipitation sequencing (MeRIP-seq), using the uninfected cells as controls. Marked alterations in the amount and distribution of m5C peaks in lncRNAs were detected between the two groups, suggesting that m5C modifications could play important regulatory roles during IAV replication.
Differentially expressed lncRNAs in H1N1 infected and uninfected A549 cells
RNAseq was performed to determine whether the lncRNA expression profile changed upon IAV infection. In total, 12,497 lncRNAs were detected in the IAV infected and uninfected A549 cells, including 1779 novel lncRNAs which were predicted as noncoding transcripts by all the four coding prediction softwares used in this study (Fig. 1A). Among the identified lncRNAs, 5435 expressed in both the infected and uninfected cells, 403 were uniquely detected in the uninfected group and 6659 lncRNAs were only detected in the H1N1-infected group, suggesting that H1N1 infection remarkably changed the expression profile of lncRNAs in the host cell (Fig. 1B). In both the infected and uninfected cells, the expressed lncRNAs were widely distributed across all chromosomes. The H1N1-infected group had a higher number of lncRNAs detected in all chromosomes than the uninfected group, demonstrating that the changes in lncRNA expression profile had taken place in all chromosomes (Fig. 1C).
Hierarchical clustering was conducted to analyze the lncRNA expression profile between the H1N1-infected and uninfected groups. Obviously, the expression levels had significantly altered upon H1N1 infection (Fig. 1D). The differential expression profile between the two groups was shown in Volcano plot in Fig. 1E. In total, there are 1566 lncRNAs that were significantly differentially expressed (DE), including 1327 up-regulated and 239 down-regulated lncRNAs (|FC|>2, P-value < 0.05) (Additional file 1).
LncRNAs are classified into intergenic and genic types which are further divided into six subtypes. Classification analysis showed that IAV infection did not dramatically change the proportional distribution of the different subtypes, and more than half of the DE lncRNAs were genic type (68%) and only 32% belonged to intergenic type (Supplementary Fig. 1 in Additional file 4). When the up-regulated and down-regulated lncRNAs were separately analyzed, it turned out that overlapping genic subtype accounted for about half of the DE transcripts in both groups. The ratios of nested genic subtype and same strand intergenic subtype were higher in the down-regulated lncRNAs than that in the up-regulated ones, while the other three subtypes, especially the divergent intergenic subtype, accounted for a larger proportion in the up-regulated lncRNAs than the down-regulated lncRNAs (Fig. 1F).
To verify the reliability of the RNA-Seq results, we selected 20 DE lncRNAs (10 upregulated and 10 downregulated) and detected their expression levels by RT-qPCR. In the A549 cells infected by 0.1 MOI of IAV at 36 hpi, the qPCR results were consistent with the RNA-Seq results (Fig. 1G-H). To evaluate whether the differential expression of these lncRNAs was the early response upon the virus infection, A549 cells infected by 1 MOI of IAV at 8 hpi were also examined by RT-qPCR (Fig. 1G-H), which identified the differential expression of 15 lncRNAs. The rest 5 lncRNAs showed consistent trend of upregulation or downregulation of expression at 8 hpi with that of 36 hpi, but their expression levels in IAV-infected cells were not significantly changed from the levels in uninfected cells, indicating that the differential expression of these lncRNAs at 36 hpi might be follow-up responses to other cellular factors that altered in early infection.
Transcriptome-wide m5C methylation of lncRNAs in IAV infected cells
To investigate the changes of m5C methylation in lncRNAs upon IAV infection, Methylated RNA immunoprecipitation sequencing (MeRIP-seq) was performed on H1N1-infected and uninfected A549 cells. We identified 2343 m5C peaks on lncRNA transcripts in A549 cells infected with IAV, and 4022 m5C peaks in uninfected cells (Fig. 2A). Of the total 5796 identified m5C peaks, 1774 were specific to IAV-infected cells, comparing with 3453 were specific to the uninfected control, with 569 sites commonly modified in both groups. These methylated peaks were mapped in 1809 and 2688 annotated lncRNAs in the IAV-infected and control groups, respectively. Among them, 606 m5C modified lncRNAs were detected in both groups (Fig. 2A). To evaluate the reliability of our data, LncRNA RPPH1, VTRNA1-1 and SCARNA2, which have been reported harboring m5C modification sites [15, 34], were examined and the m5C peaks detected in our data were visualized by IGV (Supplementary Fig. 2A in Additional file 4). It has been reported that NSUN2 is the main methyltransferase responsible for the m5C methylation of ncRNAs. To validate the specificity of the methylated RNA immunoprecipitation assay used in this study, MeRIP-seq was performed on NSUN2 knockout cells. With 879 m5C peaks identified in uninfected cells and 953 peaks in IAV-infected cells, transcriptome-wide mapping of m5C revealed a significant reduction in m5C modification (Supplementary Fig. 2B in Additional file 4). Notably, the number of m5C modification sites in uninfected cells dropped more dramatically than the number in IAV-infected cells. The m5C peaks on LncRNA RPPH1, VTRNA1-1 and SCARNA2 identified before did not appear in the NSUN2 knockout cells (Supplementary Fig. 2A in Additional file 4). These results demonstrated that the m5C peaks identified by MeRIP-seq are generally reliable.
HOMER software was used to search the consensus motifs in the m5C peak regions. Two consensus motifs were identified in the m5C peak sites, which were RYCGRH and CCGYUB (Fig. 2B). According to the gene annotation in the reference genome, the distribution of these m5C peaks in the whole lncRNA transcriptome was studied. The metagene analysis demonstrated that m5C peaks in lncRNAs were obviously enriched around the transcription start site (TSS), and this enrichment was slightly reduced after infection (Fig. 2C).
Classification analysis of the identified m5C lncRNAs revealed that the ratios of all the three subtypes of intergenic lncRNAs increased upon H1N1 infection, and the total ratio raised from 15.86 to 35.45%. Accordingly, the proportions for the subtypes of genic lncRNAs decreased in the virus infected cells, and the ratio of m5C modified genic lncRNAs decreased from 84.14 to 64.55% (Fig. 2D).
Analysis of the chromosome distribution showed that the m5C modified lncRNAs appeared more frequently on chromosomes 1, 2, 17, and 19 than the others, while few modified lncRNAs was located in the Y chromosome (Fig. 2E). Notably, we found that most (about 80%) of the methylated lncRNAs had only one m5C peak, whereas about 10% of them contained two peaks and 10% had three or more peaks (Fig. 2F). IAV infection did not significantly change the number of m5C modification in each lncRNA and the chromosome distribution of the modified lncRNAs.
Differentially methylated lncRNAs in IAV infection
When the distribution of m5C methylation sites on the whole genome was analyzed by R-circlize package, it was shown that the distribution and density of m5C peaks on each chromosome were very different between the infected and uninfected cells (Fig. 3A). Clustering analysis was performed to compare the m5C modifications in the two groups. In the heatmap, the lncRNAs differentially modified in infected cells could be clearly distinguished from those in the control cells (Fig. 3B), demonstrating the significant alteration of m5C modification upon IAV infection.
Compared to the uninfected group, there were 2984 differentially modified (DM) m5C sites identified on lncRNAs in IAV-infected cells, including 1317 hypermethylated and 1667 hypomethylated sites (Additional file 2). The 2984 differentially methylated m5C sites were located across 1956 lncRNAs, of which 914 were hypermethylated and 1042 were hypomethylated (Fig. 3C). The top ten lncRNAs, in which the m5C methylation was upregulated or downregulated with the highest fold change values, were respectively listed in Tables 1 and 2. Representative examples of hypermethylated and hypomethylated peaks were visualized by IGV and shown in Fig. 3D. Approximately 83.1% (1623/1956) of the DM lncRNAs contained one DM site, 14.2% (278/1956) had two DM sites and 2.7% had three or more altered methylation peaks (Fig. 3E). Analysis of the distribution of altered m5C peaks on chromosomes showed that these alterations were distributed on all chromosomes, with the highest number in chromosome 2 (Fig. 3F).
GO Enrichment and KEGG pathway analyses of lncRNAs harboring DM m5C sites
To explore the potential functions of lncRNA m5C methylation in IAV infection, the nearest protein-coding genes paired with the differentially modified lncRNAs were searched out and used for Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses. GO analysis based on biological processes (BP) revealed that the target-genes of up-methylated lncRNAs (H1N1-infected vs. uninfected cells) were significantly enriched in protein modification by small protein removal, protein deubiquitination, and calcium-dependent cell-cell adhesion via plasma membrane cell adhesion molecules, while the down-methylated genes were largely enriched in non-membrane-bounded organelle assembly, establishment of organelle localization and nuclear export. The molecular functions (MF) output showed that the genes targeted by the up-methylated lncRNAs were notably involved in cysteine-type deubiquitinase, HMG box domain binding and alcohol dehydrogenase (NAD(P)+) activity. In contrast, the down-methylated lncRNAs were predominantly associated with genes playing roles in transcription coregulator activity, transcription coactivator activity and cadherin binding. The cellular components (CC) data showed that up-methylated lncRNAs mainly targeted genes in SAGA-type complex, MLL1 complex and chromocenter, while the down-methylated lncRNAs were primarily associated with the genes in actin cytoskeleton, cell-substrate junction and focal adhesion (Fig. 4A).
The KEGG pathway analysis revealed that some up-methylated lncRNAs were significantly associated with genes involved in metabolism (e.g. Fatty acid metabolism, Tyrosine metabolism, Pyruvate metabolism, Retinol metabolism, Cholesterol metabolism). The down-methylated lncRNAs were notably related to genes in Central carbon metabolism and Adherens junction. A number of target genes of hypomethylated lncRNAs were involved in some infection related pathways (e.g. Epstein-Barr virus infection, Salmonella infection) (Fig. 4B).
Conjoint analysis of DM and DE lncRNAs
Conjoint analysis was then performed to explore the relationship between m5C epitranscriptomic modification and lncRNA expression level. Conjoint analysis of the DM and DE lncRNAs identified 143 hyper-up genes with both m5C modification and expression levels up-regulated, 4 hyper-down lncRNAs with up-regulated modification and decreased expression levels, 81 hypo-up lncRNAs with down-regulated modification and increased expression levels, and 6 hypo-down lncRNAs with down-regulated modification and decreased expression levels (fold change > 2, P-value < 0.05) (Fig. 5A and Additional file 3). An example of a hyper-up lncRNA and a hypo-up lncRNA visualized by IGV were presented in Fig. 5B.
GO and KEGG analyses of the genes targeted by the hyper-up and hypo-up lncRNAs were performed to investigate the potential roles of these lncRNAs in IAV infection. For the hyper-up lncRNAs, BP analysis showed that they are mainly participated in chromatin assembly, protein-DNA assembly and nucleosome assembly. MF data revealed the involvement of these lncRNAs in IgG binding, protein tyrosine kinase binding and protein heterodimerization activity. For CC, the target genes were mainly enriched in protein-DNA complex, DNA packaging complex and nucleosome. As for hypo-up lncRNAs, viral process and regulation of interferon production pathways were highlighted by BP analysis. For MF, target genes were enriched in cysteine-type endopeptidase activity involved in apoptotic signaling pathway, ATP-dependent activity and ribonucleoprotein complex binding. The CC analysis found that genes targeted by hypo-up lncRNAs were mainly enriched in ER ubiquitin ligase complex, nuclear envelope and small ribosomal subunit (Fig. 5C). KEGG pathway analyses highlighted the association of hyper-up lncRNAs with Neutrophil extracellular trap formation, and the link of hypo-up lncRNAs with Nucleocytoplasmic transport and Hepatitis C pathway. It is worth noting that most of the genes targeted by these lncRNAs were involved in pathogen recognition and disease process, and both hyper-up and hypo-up lncRNAs targeted genes in NOD-like receptor signaling pathway and Coronavirus disease - COVID-19 pathway (Fig. 5D). These data implied the potential regulation roles of the DM and DE lncRNAs in IAV infection and pathogenesis.
Researches have demonstrated that lncRNAs can act as microRNA (miRNA) sponges to regulate the expression level of protein-coding genes. To explore the possible regulatory roles of the DM and DE lncRNAs in IAV infection, five lncRNAs including three hyper-up lncRNAs (ENST00000587826.1, ENST00000656493.1, ENST00000570843.1) and two hypo-up lncRNAs (ENST00000573866.2, ENST00000552784.1) were picked out by retrieving lncbook and miRwalk databases, and their competing endogenous RNA (ceRNA) network were constructed. As shown in Fig. 6, the 5 lncRNAs were found to associate with 9 miRNAs, which could regulate the transcription of dozens of mRNAs that playing distinct roles in cell biology.
It has been reported that IAV infection alters the expression of numerous host lncRNAs, suggesting that this class of noncoding RNAs may play important regulation roles in the virus-host interaction [22,23,24]. Previous reports have revealed that some lncRNAs, such as TSPOAP1-AS1 , lnc-MxA  and lnc-Lsm3b , are upregulated in response to IAV infection and they promote virus replication by suppressing the production of type I interferon. LncRNA PAAN  and IPAN  have been shown to be induced to promote viral RNA synthesis by association with viral PA and PB1 proteins, respectively. Lnc45 was found to be highly upregulated by different IAV subtypes, and it suppressed influenza virus replication by inhibiting viral polymerase activity and retaining NP and PA in the nucleus of infected cells through its stem-loop arm . The expression of LncRNA ACOD1 is elevated during IAV infection to facilitate virus replication by modulating cellular metabolism . LncRNA NRAV could promote IAV replication by suppressing the initial transcription of several key ISGs, including IFIT2, IFIT3, OASL, IFITM3 and MxA, but it is dramatically downregulated during IAV infection . Several IAV-upregulated lncRNAs have been found to inhibit IAV replication. For example, lnc-ISG20  and ISR  function as interferon-stimulated genes (ISGs), and LncRNA-155 , NEAT1 , SAAL , RDUR , AVAN  and IVRPIE  can suppress IAV replication by promoting innate immune responses to viral infection. However, it remains unknown how the expression levels of these lncRNAs are regulated and whether RNA modifications play roles in the lncRNA regulation in response to the virus infection.
In this study, we examined the expression of lncRNAs in H1N1-infected cells, and detected more lncRNAs on each chromosome in the infected group than the uninfected group, indicating an overall upregulation of lncRNAs expression after H1N1 infection. The phenomenon that the majority of DE lncRNAs have elevated expression levels upon IAV infection had been observed in H1N1 and H3N2 infected cells [24, 26]. Here, by differential expression analysis, we identified 1566 DE lncRNAs including 1327 upregulated and 239 downregulated ones in H1N1-infected A549 cells. It is worth exploring how the abundance of these lncRNAs are regulated and what roles they play in IAV infection.
In RNA metabolism, m5C has been reported to act as a modulator for the stability , nuclear export  and protein translation  of cellular mRNAs. Studies have also found that m5C modification of enhancer RNA  and lncRNA  can increase the stability of these noncoding RNAs. m5C and m6A modifications of LncRNA NKILA have been found to facilitate cholangiocarcinoma growth and metastasis through the miR-582-3p-YAP1 axis . Recently, it was reported that m5C methylation in lncRNAs played an important role in the regulation of type I interferons and antiviral innate immunity .
By MeRIP-Seq, we studied how the m5C modification in lncRNAs was affected by H1N1 IAV infection. The results obtained in A549 cells showed that the m5C modification sites and distribution on lncRNAs were globally altered after IAV infection. A total of 1317 upregulated and 1667 downregulated m5C peaks were detected upon H1N1 infection. GO analysis showed that the genes targeted by the m5C differentially modified (m5C-DM) lncRNAs in IAV infection were notably enriched in biological processes of calcium-dependent cell-cell adhesion via plasma membrane cell adhesion molecules, protein deubiquitination, nuclear export and organelle localization. KEGG pathway analyses suggested that these m5C-DM lncRNAs were mainly associated with metabolism (such as fatty acid metabolism, and tyrosine metabolism) and some infectious pathogens related pathways which were known to be important for virus infections. These data indicated that m5C modification on lncRNAs could play important roles in the host responses to IAV infection.
By conjoint analysis of the epitranscriptomic profile and expression profile, we identified 143 ‘hyper-up’, 4 ‘hyper-down’, 81 ‘hypo-up’, and 6 ‘hypo-down’ lncRNAs. GO analysis highlighted the enrichment of target genes of DM and DE lncRNAs in nucleosome assembly, and regulation of type I interferon production pathway which is consistent with the recent report of the involvement of m5C methylation in the regulation of lncRNAs expression and type I interferon production . KEGG pathway analyses revealed that these DM and DE lncRNAs were predominantly associated with pathogen recognition and pathogenesis pathways. These results suggest that m5C modifications may be involved in multiple pathways to regulate IAV replication by modulating the expression and/or stability of lncRNAs.
Different influenza viruses may influence the host cells in altered pathways or in different degrees. The same virus may also behave differently in varied cells. Our data in this study demonstrated the association of the expression and m5C modification with IAV replication. It remains elusive which of these DE and DM lncRNAs identified are common for IAV infection and which of them specifically respond to H1N1 infection in A549 cells. Furthermore, the samples subjected to MeRIP-seq analyses were harvested at 36 hpi. Part of the DE and DM lncRNAs could be indirectly regulated by the cellular factors altered in early infection. Validation of the altered expression and modification of lncRNAs in early infection could be helpful for screening out the key players in IAV infection. Further studies are also needed to answer some key questions on how IAV induces the differential modification level and distribution of m5C, what are the roles of the modifications in the regulation of lncRNA expression and function, and what are the functions of m5C modification of lncRNAs in IAV replication and pathogenesis.
Recent studies have revealed that the viral genomic and messenger RNAs of retroviruses including HIV-1 and Murine Leukemia Virus are heavily modified with m5C, and the modification plays positive roles in ribosomal recruitment and RNA splicing to benefit viral gene expression and virus replication [14, 48]. Post-transcriptional modification analyses performed by Mass spectrometry have found that viral RNAs from positive-sense RNA viruses such as Zika virus, Dengue virus, hepatitis C virus and poliovirus bear m5C modification as well as HIV-1 . For IAV, m6A modifications on its mRNAs and vRNAs have been mapped. The modification has been shown to be able to increase viral RNA expression, and IAV HA m6A mutants show reduced pathogenicity in mice . It will be important to map m5C sites on IAV transcripts and determine whether the modifications are involved in virus replication and pathogenesis.
In conclusion, we performed a transcriptome-wide 5-methylcytosine modification analysis of lncRNAs in A549 cells infected with influenza A virus, and demonstrated the significant alteration of m5C modification upon IAV infection, suggesting a link between lncRNA m5C modification and IAV infection. The data obtained in this study could provide new insights into our understanding of virus-host interactions in influenza virus infection, and prompt further studies to explore the potential of lncRNAs as diagnostic markers and therapeutic targets for the virus infection caused diseases.
Cell culture and virus infection
Human alveolar basal epithelial adenocarcinoma cell line (A549) was purchased from Shanghai Institute of Biochemistry and Cell Biology, Chinese Academy of Sciences. Influenza virus isolate A/WSN/33 (H1N1) was donated by Dr. Long Liu. A549 cells were cultured in Ham’s F-12 K (BasalMedia) containing 10% fetal bovine serum (FBS, Hyclone), 100 U/ml penicillin and 100 µg/ml streptomycin (Sangon Biotech) at 37℃ in 5% CO2. A549 cells were seeded in 100 mm cell culture dishes. When the confluence reached 80–90%, the cells were infected with influenza A virus strain A/WSN/33 at an MOI of 0.1, incubated at 37℃ for 2 h, and then the media was replaced with F-12 K supplemented with 1% bovine serum albumin (BSA, Solarbio), antibiotics and TPCK-treated Trypsin (Sigma). Cell samples were harvested at 36 hpi for subsequent total RNA extraction. All samples are analyzed in triplicate.
RNA extraction and fragmentation
Three repeats of IAV-infected and uninfected samples were obtained. RNA was extracted from cells using RNAprep pure Cell / Bacteria Kit (TIANGEN) following the manufacturer’s instructions. Denaturing agarose gel electrophoresis was used to measure RNA integrity and gDNA contamination (Supplementary Fig. 3 in Additional file 4). The concentration of total RNA was determined by NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific). The quality of RNA was assessed by the ratio of OD260/OD280, and the samples with the value between 1.8 and 2.1 were marked as qualified (Supplementary Table 1 in Additional file 4).
RNA Library construction and sequencing
Transcriptome high throughput sequencing was performed by DIATRE Biotech (Shanghai, China). Briefly, rRNAs were removed from the total RNA with NEBNext® rRNA Depletion Kit (New England Biolabs, Inc., Massachusetts, USA). MeRIP-Seq and RNA-Seq were performed by DIATREBiotech Inc. (Shanghai, China). Immunoprecipitation of m5C RNA was performed with the GenSeqTM m5C RNA IP Kit (GenSeq Inc) by following the manufacturer’s instructions. Briefly, fragmented RNA was incubated with anti-m5C polyclonal antibody (ABclonal) in IPP buffer for 1 h at 4℃. The mixture was then immunoprecipitated by incubation with protein-A beads (Thermo Fisher) at 4℃ for an additional 2 h. Then, bound RNA was eluted from the beads with Proteinase K for 30 min at 55℃, and purified by RNA clean&concentrantorTM-5 (ZYMO Research). RNA-seq libraries for the input and IP RNA samples were generated using NEBNext® Ultra II Directional RNA Library Prep Kit (New England Biolabs). The library quality was evaluated with BioAnalyzer 2100 system (Agilent Technologies) (Supplementary Tables 2 and Supplementary Fig. 4 in Additional file 4). Library sequencing was performed on an illumina NovaSeq 6000 instrument with 150 bp paired-end reads. The sample preparation processes are illustrated in the flowchart in Fig. 7.
Quality control and transcriptome assembly
Quality control of the paired-end reads was performed with Q30, which was followed by trimming of the 3’ adaptor and removal of low-quality reads using Cutadapt software (v1.9.3) . After that, clean reads of all libraries were aligned to the reference genome (Homo sapiens.GRCh38) by Hisat2 software (v2.0.4) . Then, input data reads were aligned to the reference genome and used for transcriptome assembly and quantification by Stringtie (v2.2.1) . Novel transcripts were annotated using GffCompare software (v0.11.2) .
Identification and quantification of lncRNAs
Within the annotated transcripts, lncRNAs were filtered by the classcodes of x (exonic overlap with reference on the opposite strand), u (intergenic transcripts), i (transcripts entirely within intron), j (at least one splice junction is shared with a reference transcript), and o (other same strand overlap with reference exons). Transcripts with length > 200 bp and exon number > 2 were maintained. The coding potential of these transcripts was predicted by CNCI (v2) , CPC2 , CPAT (v3.0.0)  and PLEK (v1.2)  software, and transcripts with coding potential were removed. After these steps, identified lncRNAs were used for quantification analysis by FeatureCounts (v2.0.1) . Deseq2 software was then applied to search differentially expressed lncRNAs (|FC|>2, pvalue < 0.05) . The distribution of lncRNA on gene was analyzed by FEELnc .
MeRIP-seq data analysis
After transcriptome assembly, annotated lncRNA transcripts were used for m5C peak calling and lncRNA quantification. Methylated sites of lncRNAs (peaks) were identified by MACS2 software . Differentially methylated sites were determined using DiffBind software : m5C peaks with a log2FC > 1 (P-value < 0.05) in infected cells were considered to be up-regulated, and those with log2FC < -1 were down-regulated.
The distribution of lncRNA m5C peaks on chromosome was analyzed by R package Circlize . The normalized read counts were used for clustering analysis by R package Pheatmap . Guitar package was used for analyzing the distribution of m5C peaks on lncRNAs . The read alignments on genome were visualized using the interactive analysis tool Integrative Genomics Viewer (IGV) . HOMER software was used to search the motifs in the m5C peak regions. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)  analyses were performed based on the cis-target genes of lncRNAs to explore the functions of these lncRNAs by R package ClusterProfiler . The interactions between lncRNAs and miRNAs were retrieved from lncbook (https://ngdc.cncb.ac.cn/lncbook/home) , and the association between miRNA and mRNA was predicted by miRWalk (http://mirwalk.umm.uni-heidelberg.de/) . The lncRNA-miRNA-mRNA network was constructed by Cytoscape based on the interactions of lncRNAs, miRNAs and mRNAs.
Quantitative real-time PCR
Total RNAs from H1N1-infected and uninfected A549 cells were used to synthesize cDNA using the HiFiScript gDNA Removal cDNA Synthesis Kit (Cwbio, China). Quantitative Real-time PCR (qPCR) was carried out using ChamQ Universal SYBR qPCR Master Mix (Vazyme, China) according to the manufacturer’s instructions. GAPDH was used as a normalization control. The relative expression level of each lncRNA was calculated with the 2−ΔΔCt method. All samples were analyzed in triplicate. The primer sequences used in the qPCR are listed in Table 3.
The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive in BIG Data Center, Beijing Institute of Genomics (BIG), Chinese Academy of Sciences, under accession number HRA003296 which is publicly accessible at https://ngdc.cncb.ac.cn/gsa-human. For editors and reviewers, please access the sequence data from https://ngdc.cncb.ac.cn/gsa-human/s/1d1Ht4xw.
All other data supporting the findings of this study are available from the corresponding author on reasonable request.
Influenza A virus
Human alveolar basal epithelial adenocarcinoma cell line
Long non-coding RNAs
Methylated RNA immunoprecipitation sequencing
Kyoto Encyclopedia of Genes and Genomes
Fetal bovine serum
Bovine serum albumin
Phosphate buffered saline
Meyer KD, Jaffrey SR. The dynamic epitranscriptome: N6-methyladenosine and gene expression control. Nat Rev Mol Cell Biol. 2014;15:313–26.
Li X, Xiong X, Wang K, Wang L, Shu X, Ma S, et al. Transcriptome-wide mapping reveals reversible and dynamic N(1)-methyladenosine methylome. Nat Chem Biol. 2016;12:311–6.
Hussain S. The emerging roles of Cytosine-5 methylation in mRNAs. Trends Genet. 2021;37:498–500.
Hussain S, Sajini AA, Blanco S, Dietmann S, Lombard P, Sugimoto Y, et al. NSun2-mediated cytosine-5 methylation of vault noncoding RNA determines its processing into regulatory small RNAs. Cell Rep. 2013;4:255–61.
Yuan S, Tang H, Xing J, Fan X, Cai X, Li Q, et al. Methylation by NSun2 represses the levels and function of microRNA 125b. Mol Cell Biol. 2014;34:3630–41.
Tuorto F, Liebers R, Musch T, Schaefer M, Hofmann S, Kellner S, et al. RNA cytosine methylation by Dnmt2 and NSun2 promotes tRNA stability and protein synthesis. Nat Struct Mol Biol. 2012;19:900–5.
Zhang X, Liu Z, Yi J, Tang H, Xing J, Yu M, et al. The tRNA methyltransferase NSun2 stabilizes p16INK4 mRNA by methylating the 3’-untranslated region of p16. Nat Commun. 2012;3:712.
Yang X, Yang Y, Sun BF, Chen YS, Xu JW, Lai WY, et al. 5-methylcytosine promotes mRNA export - NSUN2 as the methyltransferase and ALYREF as an m(5)C reader. Cell Res. 2017;27:606–25.
Xing J, Yi J, Cai X, Tang H, Liu Z, Zhang X, et al. NSun2 promotes cell growth via elevating cyclin-dependent kinase 1 translation. Mol Cell Biol. 2015;35:4043–52.
David R, Burgess A, Parker B, Li J, Pulsford K, Sibbritt T, et al. Transcriptome-wide mapping of RNA 5-Methylcytosine in Arabidopsis mRNAs and noncoding RNAs. Plant Cell. 2017;29:445–60.
Huang T, Chen W, Liu J, Gu N, Zhang R. Genome-wide identification of mRNA 5-methylcytosine in mammals. Nat Struct Mol Biol. 2019;26:380–8.
Squires JE, Patel HR, Nousch M, Sibbritt T, Humphreys DT, Parker BJ, et al. Widespread occurrence of 5-methylcytosine in human coding and non-coding RNA. Nucleic Acids Res. 2012;40:5023–33.
Tsai K, Cullen BR. Epigenetic and epitranscriptomic regulation of viral replication. Nat Rev Microbiol. 2020;18:559–70.
Courtney DG, Tsai K, Bogerd HP, Kennedy EM, Law BA, Emery A, et al. Epitranscriptomic addition of m(5)C to HIV-1 transcripts regulates viral gene expression. Cell Host Microbe. 2019;26:217–227e6.
Zhang Y, Zhang L-S, Dai Q, Chen P, Lu M, Kairis EL, et al. 5-methylcytosine (m(5)C) RNA modification controls the innate immune response to virus infection by regulating type I interferons. Proc Natl Acad Sci U S A. 2022;119:e2123338119.
Krammer F, Smith GJD, Fouchier RAM, Peiris M, Kedzierska K, Doherty PC, et al. Influenza Nat Rev Dis Primers. 2018;4:3.
Medina RA, García-Sastre A, Influenza. A viruses: new research developments. Nat Rev Microbiol. 2011;9:590–603.
Sanders CJ, Vogel P, McClaren JL, Bajracharya R, Doherty PC, Thomas PG. Compromised respiratory function in lethal influenza infection is characterized by the depletion of type I alveolar epithelial cells beyond threshold levels. Am J Physiol Lung Cell Mol Physiol. 2013;304:L481–8.
Taubenberger JK, Morens DM. The pathology of influenza virus infections. Annu Rev Pathol. 2008;3:499–522.
Basavappa M, Cherry S, Henao-Mejia J. Long noncoding RNAs and the regulation of innate immunity and host-virus interactions. J Leukoc Biol. 2019;106:83–93.
Meng XY, Luo Y, Anwar MN, Sun Y, Gao Y, Zhang H, et al. Long non-coding RNAs: emerging and versatile regulators in host-virus interactions. Front Immunol. 2017;8:1663.
Peng X, Gralinski L, Armour CD, Ferris MT, Thomas MJ, Proll S et al. Unique signatures of long noncoding RNA expression in response to virus infection and altered innate immune signaling. mBio. 2010;1.
Wang J, Cen S. Roles of lncRNAs in influenza virus infection. Emerg Microbes Infect. 2020;9:1407–14.
Zhang Y, Yu T, Ding Y, Li Y, Lei J, Hu B, et al. Analysis of expression profiles of long noncoding RNAs and mRNAs in A549 cells infected with H3N2 swine influenza virus by RNA sequencing. Virol Sin. 2020;35:171–80.
Chai W, Li J, Shangguan Q, Liu Q, Li X, Qi D et al. Lnc-ISG20 inhibits Influenza A Virus replication by enhancing ISG20 expression. J Virol. 2018;92.
Maarouf M, Chen B, Chen Y, Wang X, Rai KR, Zhao Z, et al. Identification of lncRNA-155 encoded by MIR155HG as a novel regulator of innate immunity against influenza a virus infection. Cell Microbiol. 2019;21:e13036.
More S, Zhu Z, Lin K, Huang C, Pushparaj S, Liang Y, et al. Long non-coding RNA PSMB8-AS1 regulates influenza virus replication. RNA Biol. 2019;16:340–53.
Ouyang J, Zhu X, Chen Y, Wei H, Chen Q, Chi X, et al. NRAV, a long noncoding RNA, modulates antiviral responses through suppression of interferon-stimulated gene transcription. Cell Host Microbe. 2014;16:616–26.
Pan Q, Zhao Z, Liao Y, Chiu SH, Wang S, Chen B et al. Identification of an Interferon-Stimulated long noncoding RNA (LncRNA ISR) involved in regulation of Influenza A Virus Replication. Int J Mol Sci. 2019;20.
Wang Q, Zhang D, Feng W, Guo Y, Sun X, Zhang M, et al. Long noncoding RNA TSPOAP1 antisense RNA 1 negatively modulates type I IFN signaling to facilitate influenza a virus replication. J Med Virol. 2019;94:557–66.
Jacob R, Zander S, Gutschner T. The Dark side of the Epitranscriptome: chemical modifications in long non-coding RNAs. Int J Mol Sci. 2017;18.
Wang J, Wang Y, Zhou R, Zhao J, Zhang Y, Yi D et al. Host long noncoding RNA lncRNA-PAAN regulates the replication of Influenza A Virus. Viruses. 2018;10.
Wang J, Zhang Y, Li Q, Zhao J, Yi D, Ding J, et al. Influenza virus exploits an Interferon-Independent lncRNA to preserve viral RNA synthesis through stabilizing viral RNA polymerase PB1. Cell Rep. 2019;27:3295–3304e4.
Schumann U, Zhang H-N, Sibbritt T, Pan A, Horvath A, Gross S, et al. Multiple links between 5-methylcytosine content of mRNA and translation. BMC Biol. 2020;18:40.
Li X, Guo G, Lu M, Chai W, Li Y, Tong X et al. Long noncoding RNA Lnc-MxA inhibits Beta Interferon transcription by forming RNA-DNA triplexes at its promoter. J Virol. 2019;93.
Jiang M, Zhang S, Yang Z, Lin H, Zhu J, Liu L, et al. Self-recognition of an inducible host lncRNA by RIG-I feedback restricts Innate Immune Response. Cell. 2018;173:906–919e13.
Zhang L, Zheng X, Li J, Wang G, Hu Z, Chen Y, et al. Long noncoding RNA#45 exerts broad inhibitory effect on influenza a virus replication via its stem ring arms. Virulence. 2021;12:2443–60.
Wang P, Xu J, Wang Y, Cao X. An interferon-independent lncRNA promotes viral replication by modulating cellular metabolism. Science. 2017;358:1051–5.
Imamura K, Imamachi N, Akizuki G, Kumakura M, Kawaguchi A, Nagata K, et al. Long noncoding RNA NEAT1-dependent SFPQ relocation from promoter region to paraspeckle mediates IL8 expression upon immune stimuli. Mol Cell. 2014;53:393–406.
Liu Q, Yang H, Zhao L, Huang N, Ping J. A novel lncRNA SAAL suppresses IAV replication by promoting innate responses. Microorganisms. 2022;10:2336.
Chen Y, Hu J, Liu S, Chen B, Xiao M, Li Y, et al. RDUR, a lncRNA, promotes innate antiviral responses and provides Feedback Control of NF-κB activation. Front Immunol. 2021;12:672165.
Lai C, Liu L, Liu Q, Wang K, Cheng S, Zhao L, et al. Long noncoding RNA AVAN promotes antiviral innate immunity by interacting with TRIM25 and enhancing the transcription of FOXO3a. Cell Death Differ. 2021;28:2900–15.
Zhao L, Xia M, Wang K, Lai C, Fan H, Gu H, et al. A long non-coding RNA IVRPIE promotes host antiviral Immune responses through regulating Interferon β1 and ISG expression. Front Microbiol. 2020;11:260.
Yang Y, Wang L, Han X, Yang WL, Zhang M, Ma HL, et al. RNA 5-Methylcytosine facilitates the maternal-to-zygotic transition by preventing maternal mRNA decay. Mol Cell. 2019;75:1188–1202e11.
Aguilo F, Li S, Balasubramaniyan N, Sancho A, Benko S, Zhang F, et al. Deposition of 5-Methylcytosine on enhancer RNAs enables the Coactivator function of PGC-1α. Cell Rep. 2016;14:479–92.
Sun Z, Xue S, Zhang M, Xu H, Hu X, Chen S, et al. Aberrant NSUN2-mediated m(5)C modification of H19 lncRNA is associated with poor differentiation of hepatocellular carcinoma. Oncogene. 2020;39:6906–19.
Zheng H, Zhu M, Li W, Zhou Z, Wan X. m5C and m6A modification of long noncoding NKILA accelerates cholangiocarcinoma progression via the mir-582-3p-YAP1 axis. Liver Int. 2022;42:1144–57.
Courtney DG, Chalem A, Bogerd HP, Law BA, Kennedy EM, Holley CL et al. Extensive epitranscriptomic methylation of A and C residues on murine leukemia virus transcripts enhances viral gene expression. mBio. 2019;10.
McIntyre W, Netzband R, Bonenfant G, Biegel JM, Miller C, Fuchs G, et al. Positive-sense RNA viruses reveal the complexity and dynamics of the cellular and viral epitranscriptomes during infection. Nucleic Acids Res. 2018;46:5776–91.
Courtney DG, Kennedy EM, Dumm RE, Bogerd HP, Tsai K, Heaton NS, et al. Epitranscriptomic enhancement of Influenza A Virus Gene expression and replication. Cell Host Microbe. 2017;22:377–386e5.
Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. 2011. 2011;17:3.
Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357–60.
Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33:290–5.
Pertea G, Pertea M. GFF Utilities: GffRead and GffCompare. F1000Res. 2020;9.
Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, et al. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;41:e166.
Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L et al. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007;35 Web Server issue:W345–9.
Wang L, Park HJ, Dasari S, Wang S, Kocher JP, Li W. CPAT: coding-potential Assessment Tool using an alignment-free logistic regression model. Nucleic Acids Res. 2013;41:e74.
Li A, Zhang J, Zhou Z. PLEK: a tool for predicting long non-coding RNAs and messenger RNAs based on an improved k-mer scheme. BMC Bioinformatics. 2014;15:311.
Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30:923–30.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Wucher V, Legeai F, Hédan B, Rizk G, Lagoutte L, Leeb T, et al. FEELnc: a tool for long non-coding RNA annotation and its application to the dog transcriptome. Nucleic Acids Res. 2017;45:e57.
Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008;9:R137.
Stark RBG. DiffBind: differential binding analysis of ChIP-Seq peak data. 2011;http://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf.
Gu Z, Gu L, Eils R, Schlesner M, Brors B. Circlize implements and enhances circular visualization in R. Bioinformatics. 2014;30:2811–2.
R K. pheatmap: pretty heatmaps. https://cran.r-project.org/web/packages/pheatmap/pheatmap.pdf. 2019.
Cui X, Wei Z, Zhang L, Liu H, Sun L, Zhang SW, et al. Guitar: an R/Bioconductor Package for Gene Annotation guided transcriptomic analysis of RNA-Related genomic features. Biomed Res Int. 2016;2016:8367534.
Thorvaldsdóttir H, Robinson JT, Mesirov JP. Integrative Genomics viewer (IGV): high-performance genomics data visualization and exploration. Brief Bioinform. 2013;14:178–92.
Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2023;51:D587–92.
Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innov (Camb). 2021;2:100141.
Ma L, Cao J, Liu L, Du Q, Li Z, Zou D, et al. LncBook: a curated knowledgebase of human long non-coding RNAs. Nucleic Acids Res. 2019;47:D128–d134.
Sticht C, De La Torre C, Parveen A, Gretz N. miRWalk: an online resource for prediction of microRNA binding sites. PLoS ONE. 2018;13:e0206239.
We thank the Teaching and Research Core Facility at College of Life Science, NWAFU for their support in this work.
This work was supported by NWAFU and the Principle Investigator Program at Hubei University of Medicine (HBMUPI202102).
Ethics approval and consent to participate
Consent for publication
All authors declare that they have no conflict of interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Jiang, S., Hu, J., Bai, Y. et al. Transcriptome-wide 5-methylcytosine modification profiling of long non-coding RNAs in A549 cells infected with H1N1 influenza A virus. BMC Genomics 24, 316 (2023). https://doi.org/10.1186/s12864-023-09432-z