Comprehensive Genomic Characterization Analysis of lncRNAs in Cells With Porcine Delta Coronavirus Infection

Porcine delta coronavirus (PDCoV) is a novel emerging enterocytetropic virus causing diarrhea, vomiting, dehydration, and mortality in suckling piglets. Long non-coding RNAs (lncRNAs) are known to be important regulators during virus infection. Here, we describe a comprehensive transcriptome profile of lncRNA in PDCoV-infected swine testicular (ST) cells. In total, 1,308 annotated and 1,190 novel lncRNA candidate sequences were identified. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis revealed that these lncRNAs might be involved in numerous biological processes. Clustering analysis of differentially expressed lncRNAs showed that 454 annotated and 376 novel lncRNAs were regulated after PDCoV infection. Furthermore, we constructed a lncRNA-protein-coding gene co-expression interaction network. The KEGG analysis of the co-expressed genes showed that these differentially expressed lncRNAs were enriched in pathways related to metabolism and TNF signaling. Our study provided comprehensive information about lncRNAs that would be a useful resource for studying the pathogenesis of and designing antiviral therapy for PDCoV infection.

Introduction

Long non-coding RNAs (lncRNAs), which are transcripts larger than 200 nt in length that lack protein-coding ability, have previously been described in mammalian cells (Kapranov et al., 2007; Mattick and Rinn, 2015). Most of them have a structure similar to mRNA; they have a 5′ methylguanosine cap and are usually spliced and polyadenylated at their 3′ termini. Notably, lncRNA expression shows significant cell and tissue specificity (Mercer et al., 2008; Derrien et al., 2012). Emerging evidence shows that non-coding RNAs have a regulatory role in multiple cellular processes, such as genomic imprinting, chromatin modification, and alternative splicing of RNA (Mercer et al., 2009). Moreover, some diseases such as cancer and neurological disorders are also related to the dysregulated expression of lncRNA (Qureshi et al., 2010; Tsai et al., 2011). Numerous studies have been conducted to ascertain their functional role during viral infection. For example, NRAV can promote influenza virus replication and virulence through negatively regulating the initial transcription of varieties of interferon-stimulated genes (ISGs) (Ouyang et al., 2014). lncRNA-ACOD1, named by its neighboring coding gene aconitate decarboxylase 1, significantly reduces virus multiplication by directly interacting with the metabolic enzyme glutamic-oxaloacetic transaminase (Wang et al., 2017). Neat1, one of the lncRNAs induced by HIV-1 infection, is retained in the nucleus and serves as a scaffold for the nuclear paraspeckle substructure. Importantly, Neat1 deficiency enhances HIV-1 replication (Zhang et al., 2013). Although large amounts of data have proved that several lncRNAs are involved in different kinds of virus infection, the mechanisms by which they act are still largely unknown.

Porcine delta coronavirus (PDCoV), a novel emerging pathogenic enterocytetropic virus, was first discovered from the feces of pigs in Hong Kong in 2012. It is an enveloped, single-stranded positive-sense RNA virus. It belongs to the genus Deltacoronavirus in the family Coronaviridae of the order Nidovirales (Woo et al., 2012). The genome length of PDCoV is approximately 25.4 kb. It is similar in structure to other coronaviruses, with shorter non-coding regions (5′-UTR and 3′-UTR) at both terminals. The 3/4 genome from the 5′ end contains two overlapping open reading frames, ORF1a and ORF1b, encoding the pp1a and pp1b, respectively. The downstream of the genome encodes structural protein spike (S), envelope (E), membrane (M), accessory proteins NS6, structural protein nucleocapsid (N), and accessory proteins NS7 and NS7a. A total of 15 non-structural proteins are encoded by the 5′ terminal of the genome (Fang et al., 2017). PDCoV mainly causes acute, watery diarrhea, vomiting, dehydration, and mortality in suckling piglets, including lesions in the stomach and lungs (Ma et al., 2015). The first outbreak of PDCoV infection was reported in the United States in early 2014 and, to date, it has been detected in Canada, South Korea, China, Thailand, and Vietnam, thus posing huge threat to the swine industry and attracting a great deal of attention (Lee and Lee, 2014; Dong et al., 2016; Lorsirigool et al., 2017; Ajayi et al., 2018; Saeng-Chuto et al., 2019).

During infection, the accessory and non-structural proteins of PDCoV usually perform multiple functions to promote replication in infected cells. A previous study showed that NS6 interaction with RIG-I/MDA5 attenuated the binding activity between RIG-I/MDA5 and double-stranded RNA, resulting in a reduced level of IFN-β production (Fang et al., 2018). Also, the non-structural protein nsp5, a 3C-like protease, is an important molecule in suppressing type I IFN signaling (Zhu et al., 2017b). In addition, NEMO, an essential modulator of NF-κB, can also be cleaved by nsp5, causing inhibition of IFN-β production (Zhu et al., 2017a). Though there are many reports of immune evasion by PDCoV, the precise pathogenic mechanism of PDCoV is largely unclear.

Based on the increasing number of reports on host lncRNAs associated with virus infection, we are interested in investigating whether host lncRNAs were involved in PDCoV infection. In this study, genome-wide profiling of lncRNAs in swine testicular (ST) cells infected with PDCoV was performed using RNA-seq. We identified 830 differentially expressed lncRNAs from PDCoV-infected cells. An integrative analysis of lncRNA alterations suggested their putative role in regulating the expression of several key genes in metabolic and TNF signaling pathways during infection. In conclusion, this work supports the role of lncRNAs as important regulators of PDCoV infection.

Materials and Methods

Cells and Viruses

Swine testicular cells and porcine jejunum intestinal epithelial cells (IPEC-J2) were grown in DMEM supplemented with 10% (vol/vol) FBS (Gibco, Carlsbad, CA, United States) at 37°C in a humidified 5% CO₂ atmosphere. The PDCoV-CH-HA3-2017 (MK040455) strain, stored in our laboratory, was propagated in ST cells.

Viral Infection and RNA Extraction

For RNA-seq, ST cells were infected with PDCoV at a multiplicity of infection (MOI) of 10; the medium for PDCoV infection was DMEM containing 0.2 ug/ml Trypsin that had been TPCK-treated (Millipore Sigma, St. Louis, MO, United States) for 11 h. Mock-infected cells were placed in the same volume of DMEM, with the same concentration of TPCK-treated Trypsin. Total RNA was isolated from each group using SuPerfecTRI^TM Total RNA Isolation Reagent (Pufei, Shanghai, China) according to the manufacturer’s instructions. The RNA quality was checked by 1% agarose gel electrophoresis. The purity and concentration of RNA were measured by NanoPhotometer^® spectrophotometer (IMPLEN, München, Germany) and Qubit^® RNA Assay Kit in Qubit^® 2.0 Fluorometer (Life Technologies, Camarillo, CA, United States). RNA integrity was assessed using the RNA Nano6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, United States). For quantitative RT-PCR (RT-qPCR), ST and IPEC-J2 cells were infected or mock-infected with PDCoV at an MOI of 10 and harvested at the indicated time. All experiments were conducted in triplicate.

RNA-Seq and Data Analysis

Sequencing libraries were generated using the rRNA-depleted RNA with a NEBNext^® Ultra^TM Directional RNA Library Prep Kit (New England Biolabs, Ipswich, MA, United States). After determining the quality of the library, RNA-seq was performed using an Illumina HiSeq^TM 4000 (Illumina, San Diego, CA, United States) to generate raw reads. After removing poly-N sequences, adapters, and low-quality reads, clean reads were obtained and the paired-end reads were aligned to Ensemble pig genome (Release 76). lncRNAs were identified using TopHat2 (v2.0.9), and reads that were mapped to the pig genome were assembled using Cufflinks v2.1.1 (Trapnell et al., 2010). Cuffdiff (v2.1.1) was used to calculate the FPKMs of both lncRNAs and coding genes in each sample. Gene FPKMs were computed by summing the FPKMs of transcripts in each gene group, and differentially expressed (DE) transcripts were assigned where there was a statistically significant level of expression (p < 0.05). RNA-seq and data analysis were completed by Novogene.

GO and KEGG Enrichment Analysis

Gene Ontology (GO) enrichment analysis of differentially expressed genes or lncRNA target genes was conducted with respect to biological process, molecular function, and cellular component with the GOseq R package, in which gene length bias was corrected. Kyoto Encyclopedia of Genes and Genomes (KEGG) was used to perform pathway enrichment analysis¹. KOBAS software was used to test the level of statistical significance of enrichment of differentially expressed genes and/or lncRNA target genes in KEGG pathways (Mao et al., 2005).

RT-qPCR, RT-PCR, and Statistical Analysis

To determine the reliability of the RNA-seq data, 15 differentially expressed lncRNAs were randomly selected to test the expression by RT-qPCR. Total RNA was extracted from ST and IPEC-J2 cells using SuPerfecTRI^TM Total RNA Isolation Reagent (Pufei, Shanghai, China). First-strand cDNA was synthesized with a reverse transcriptase kit (Vazyme, Nanjing, China). RT-qPCR was performed with a SYBR Green master mix (Vazyme, Nanjing, China) on the LightCycler 96 (Roche, Basel, Switzerland). The PDCoV M gene was detected by RT-PCR. All the primers are presented in Table 1. Relative expressions were calculated using the 2^–ΔΔCt method with GAPDH as the internal control. Comparisons between groups were made using two-way ANOVA. The data reported are the mean ± SEM. Differences were considered statistically significant when p < 0.05.

TABLE 1

Table 1. Primers used for RT-qPCR validation.

Construction of the lncRNA-Protein-Coding Gene Co-expression Network

For each lncRNA, the Pearson correlation coefficient of its expression value with that of each protein-coding gene was calculated. Under the conditions of an absolute value of the Pearson correlation coefficient >0.998 and p < 0.00001, the interaction network of the differentially expressed lncRNAs and protein-coding gene co-expression pairs was then constructed using Cytoscape (v3.5.1) (Shannon et al., 2003).

Results

RNA-Seq and lncRNA Screening in PDCoV-Infected Cells

To identify the lncRNAs in PDCoV-infected cells, we sequenced the transcriptomes of the ST cells with or without PDCoV infection using high-throughput RNA sequencing. Robust and reproducible data were obtained from all samples, and more than 1 × 10⁸ clean reads per sample were retained after removing reads containing adapter or poly-N sequences and reads with low quality. Afterward, all clean reads were aligned onto the pig reference genome (Release 76) using TopHat2 and were compared and assembled with Cuffcompare and Cufflinks, respectively, and coverage analysis was performed on those clean reads on different annotated gene types. The distribution of each type of gene was counted according to the expression level. In total, eight categories of RNA were identified, according to database annotation of those transcripts, in which the protein-coding genes were highly represented (66.54% in PS and 69.10% in ST, respectively) (Figure 1A). Next, four software tools, CNCI, CPC, PhyloCSF, and PFAM, were used to calculate the protein-coding potential of assembled-transcripts to screen lncRNAs, then taking the intersection of transcripts with no coding potential in these software products as the novel lncRNA (Figure 1B). In total, 1,308 annotated and 1,190 novel lncRNA candidates were identified (Supplementary Tables S1, S2).

FIGURE 1

Figure 1. RNA-seq and the characteristics of lncRNAs in PDCoV-infected ST cells. (A) RNA categories of quantified genes in RNA-seq. The percentages of each RNA species are shown in infected or mock-infected ST cells. (B) Venn diagrams of the identified lncRNAs shared among four software products. CNCI, CPC, PhyloCSF, and PFAM were used to screen novel lncRNAs. The sum of the numbers in each large circle represents the total number of lncRNA candidates of the indicated software, and the overlapping portions of the circles represented the lncRNAs shared between the software products. The 1190 novel lncRNAs were totally identified in ST cells. (C–E) Comparison of the conservation (C), exon number (D), and ORF length (E) of lncRNAs and protein-coding genes.

It has been reported that lncRNAs, in comparison with protein-coding genes, usually share some common genomic features to their sequences. They are generally shorter in length, have fewer but longer exons, and there is lower evolutionary sequence conservation, with only ∼15% of mouse lncRNAs having homologs in humans. lncRNAs also demonstrate low expression levels (the median is ∼10% of that of protein-coding genes) (Heward and Lindsay, 2014). To further determine the characteristics of the lncRNAs identified in the present study, we compared the transcript length, exon number, and degree of conservation between protein-coding genes and lncRNAs. Conservation analysis of exons, introns, and promoters of lncRNAs and protein-coding genes showed that the exons of protein-coding genes were the most conserved and the exons of lncRNA were far less conserved (Figure 1C). Furthermore, fewer exons and shorter ORFs were found in lncRNAs, which was also consistent with the reported lncRNAs (Figures 1D,E).

Whole Transcriptome-Wide Functional Prediction of lncRNAs in PDCoV-Infected Cells

Long non-coding RNAs sequences are poorly conserved and do not appear to form large homologous families, so it is difficult to infer their common ancestors by sequence similarity (Ponting et al., 2009). Therefore, it is challenging to predict the functions of a type of lncRNA on the basis of its sequence or structure. There have been reports of using genome-wide association analysis between lncRNAs and the co-expressed and/or co-regulated protein-coding genes to characterize the function of the lncRNA (Huarte et al., 2010). To investigate the putative role of lncRNAs, we first analyzed the whole RNA-seq profiles to identify target protein-coding genes whose location or expression was significantly correlated with the candidate lncRNA. For co-located target gene prediction, we searched coding regions 100 k upstream and downstream of lncRNA. In total, 8,812 pairs of lncRNA-protein-coding genes, containing 2,088 lncRNAs and 3,566 protein-coding genes, were identified (Supplementary Table S3). For co-expressed target gene prediction, the expression correlation between lncRNAs and protein-coding genes was evaluated. When the required Pearson correlation coefficient was set above 0.95, 1,048,575 pairs of lncRNA-protein-coding genes, containing 1,730 lncRNAs and 10,581 protein-coding genes, were obtained (Supplementary Table S4). We next performed GO and KEGG pathway analysis for the target genes of lncRNAs. The top 20 GO and KEGG pathways with the highest representation of each term are reported (Figure 2 and Supplementary Tables S5, S6). KEGG enrichment analysis revealed that pathways related to the immune system and metabolism were preferentially targeted. The GO-term analysis was divided into three main categories: cellular component, biological process, and molecular function. Significantly, a large number of biological processes, like protein-DNA complex assembly, DNA packaging and transcription, and the cellular macromolecule metabolic process, were enriched. Furthermore, protein binding and nucleic acid binding and the nucleosome and organelles, belonging to molecular function and cellular component, respectively, were also enriched. GO and KEGG pathway enrichment analysis of target genes revealed that lncRNAs may act in cis or trans to participate in the regulation of expression of multiple important genes in different processes including protein binding, DNA transcription, metabolism, and immune response.

FIGURE 2

Figure 2. Whole transcriptome-wide functional prediction of lncRNA identified in ST cells. Representative overrepresented KEGG (top) and GO (bottom) terms of gene-expression clusters. The number of genes within each cluster is shown for significantly enriched terms (p < 0.05). BP, biological process; MF, molecular function; CC, cellular component.

Clustering Analysis Identified Differentially Expressed lncRNAs in PDCoV-Infected Cells

To identify the PDCoV-associated lncRNAs, Cuffdiff software was used to investigate the differentially expressed (DE) lncRNAs in PDCoV-infected cells. The hierarchical clustering heat map in Figure 3A shows the DE lncRNA expression profiling data. Out of the 1,308 annotated and 1,190 novel lncRNAs, we obtained 454 annotated DE lncRNAs (225 up-regulated and 229 down-regulated) and 376 novel DE lncRNAs (252 up-regulated and 124 down-regulated) after PDCoV infection (p < 0.05; Supplementary Table S7). Importantly, we observed 20 lncRNAs whose expression levels were decreased to FPKM = 0 after PDCoV infection, while the FPKMs of another 12 lncRNAs, all novel lncRNAs, were 0 before PDCoV infection (Supplementary Table S7). This suggests that these 32 lncRNAs, though they have very low expression levels, might be strongly associated with the viral infection. Furthermore, to evaluate the reliability of RNA-seq data analysis, 15 lncRNAs were selected for RT-qPCR analysis in PDCoV-infected cells. As shown in Figure 3B, the expression levels of the 15 selected lncRNAs, though exhibiting no significant differences at 4 h post-infection (hpi), were all significantly changed at 11 hpi in ST cells. Also, different expression patterns of lncRNAs were detected in IPEC-J2 cells. As shown in Figure 3C, 11 out of the 15 selected RNA were significantly altered at 11 hpi, and all of them were differently expressed at 24 hpi. For both ST and IPEC-J2 cells, they had a strong expression pattern consistent with the RNA-seq results (Table 2).

FIGURE 3

Figure 3. Identification of differentially expressed lncRNAs during PDCoV infection. (A) Hierarchic clustering analyses of 830 lncRNAs that were differentially expressed in the PDCoV-infected ST cells (PS) compared with the mock-infected ST cells (ST) (p < 0.05). (B,C) Validation of 10 up-regulated (left) and five down-regulated (right) lncRNAs from RNA-Seq in ST (B) and IPEC-J2 cells (C) were performed by RT-qPCR. The PDCoV M gene was confirmed by RT-PCR. Data are shown as mean ± SEM, n = 3. ^∗p < 0.05, ^∗∗p < 0.01.

TABLE 2

Table 2. Expression level of 15 selected lncRNAs in RNA-seq.

Co-location Analysis of DE lncRNAs Revealed Their Potential Regulation of Their Neighboring Protein-Coding Genes

The lncRNA in the genome is not randomly distributed, so locus classification will be an effective first step in analyzing its regulatory functions at the genome level (Luo et al., 2016). In general, lncRNAs function either in cis or in trans to affect the transcription of genes within or far from the same genomic locus (Clark and Blackshaw, 2014). To understand the potential functional association between lncRNAs and cognate genes, we investigated their genomic distribution pattern relative to protein-coding loci and classified all DE lncRNAs to ascertain their potential biological roles. All DE lncRNAs were classified into six categories comprising sense-upstream lncRNA, sense-downstream lncRNA, sense-overlapping lncRNA, antisense-upstream lncRNA, antisense-downstream lncRNA and antisense-overlapping lncRNA. As shown in Figure 4A, 26% of DE lncRNAs were located in the same strand but upstream of protein-coding genes and 24% were located downstream, while antisense-upstream and antisense-overlapping comprised 27 and 1%, respectively, and the remaining 22% were antisense-downstream lncRNAs. Next, in order to define the lncRNA functions more precisely, GO enrichment analysis of the co-located genes of up- and down-regulated lncRNAs were analyzed independently. The results showed that protein-coding genes associated with DE lncRNAs were mainly enriched in terms of molecular function and cellular component, primarily under the category of nucleic acid binding and intracellular membrane-bounded organelle (Figure 4B). Notably, by analyzing the relative expression level, we found that antisense lncRNA and protein-coding genes were specifically co-expressed, in which two pairs showed a positive and three pairs showed a negative correlation in their expression patterns (Figure 4C and Supplementary Table S8). We speculated that these antisense lncRNAs act in cis to modulate the expression of their cognate genes.

FIGURE 4

Figure 4. Co-location analysis revealed the potential regulation ability of lncRNAs on their neighboring protein-coding genes. (A) Pie chart showing the percentage of each class of lncRNA classified by their location. (B) Over-represented GO terms of up-regulated (left) and down-regulated (right) co-located genes (p < 0.05). BP, biological process; MF, molecular function; CC, cellular component. (C) Heatmap showing the expression correlation of selected antisense lncRNA-protein-coding gene pairs.

Correlation Analysis Provided a Resource to Functionally Identify PDCoV Driver Metabolism- and Immunity-Related lncRNAs

The functional association between regulatory lncRNA and protein-coding gene transcripts can be determined by performing expression correlation analysis coupled with ascertaining their putative role in related physiological processes. To further investigate the potential mechanism of action of the PDCoV-associated lncRNAs, the DE lncRNAs and their predicted target DE protein-coding genes were investigated by delineating lncRNA-protein-coding gene functional interactions. Here we identified 1,048,575 pairs of DE lncRNA-DE protein-coding genes, containing 821 lncRNAs and 8,799 protein-coding genes (p < 0.01). Next, KEGG pathway analysis was repeated once again (Figure 5A), and we found that metabolic and TNF signaling pathways were significantly enriched.

FIGURE 5

Figure 5. Correlation analysis provides a resource to functionally identify PDCoV driver metabolism- and immunity-related lncRNAs. (A) The top 20 overrepresented KEGG pathways of co-expressed genes. (B,C) The network of co-expressed lncRNA-protein-coding gene pairs based on a threshold of a Pearson’s correlation coefficient of 0.998. The network of lncRNAs with metabolic- (B) or immunity-related genes (C) is presented (p < 0.00001). The solid line represents a positive correlation, while the dotted line represents a negative correlation.

The interaction network involving the metabolic and TNF signaling pathways was then constructed. Several key genes in metabolism were positively or negatively regulated by lncRNAs (Figure 5B). Of the significantly enriched genes, ATP5L and ATP5F1, two of the mitochondrial membrane ATP synthase subunits, were regulated by LNC-000625, LNC-001104, ALDBSSCT0000008902, and ALDBSSCT0000006348. In addition, three lncRNAs, LNC-000459, LNC-000258, and ALDBSSCT0000005568, might regulate acyl-coenzyme A thioesterase four expression. These results suggest that these lncRNAs might be involved in the regulation of metabolic processes particularly involving energy and lipid metabolism. Meanwhile, an inducible program of inflammatory gene expression is central to antiviral defense. Many of them, i.e., CCL5, CCL20, CXCL2, CXCL10, MAP3K8, NF-κB1, and interleukin 6 (IL-6), were protein-coding genes known to have roles in the inflammatory response. In the network (Figure 5C), eight lncRNAs have putative regulatory roles in IL-6 expression. Six of them, LNC-000173, LNC-000269, LNC-000242, LNC-000657, ALDBSSCT0000009132, and ALDBSSCT0000001339, might exert positive regulation, while LNC-001173 and ALDBSSCT0000010894 showed the opposite effect. This suggests that these lncRNAs might act as the regulatory module of the circuit that is involved in the inflammatory response.

Discussion

Numerous studies have shown that lncRNAs play a key role during viral infection. The lncRNAs THRIL, NeST, NEAT, and lincRNA-Cox2 can participate in immune responses against viral infection mainly through regulating the expression of TNF-α, IFN-γ, IL8, and inflammatory response, respectively (Carpenter et al., 2013; Gomez et al., 2013; Imamura et al., 2014; Li et al., 2014). PDCoV is an important enteric virus mainly causing diarrhea in suckling pigs. Infection with PDCoV causes changes in the expression levels of several host cell proteins of host innate immune response, but little is known about the critical roles of lncRNAs in these processes.

Here, we performed RNA-seq to identify the lncRNAs involved in PDCoV infection. The results of comparing clean reads to the genome showed that more than 60% of reads are protein-coding genes, and no lncRNA classifications were identified due to the limited lncRNA database annotation in pig. In our results, 1,190 novel lncRNAs were identified. Further analysis showed that the basic characteristics of these novel lncRNAs are consistent with the known ones. Our RNA-seq results further enrich the pig lncRNA database.

In total, we found 454 annotated and 376 novel lncRNAs that were differentially expressed during PDCoV infection. These lncRNAs were classified as sense-upstream lncRNA, sense-downstream lncRNA, sense-overlapping lncRNA, antisense-upstream lncRNA, antisense-downstream lncRNA, and antisense-overlapping lncRNA. Many antisense-overlapping lncRNAs have inverse expression patterns with their sense transcript counterparts. This suggests that these antisense-overlapping lncRNAs may have a negative regulatory effect on them. In contrast, many lncRNAs that do not contain overlapping sequences display expression patterns correlated with their neighboring protein-coding gene transcripts. In the present study, two out of five antisense overlapping lncRNAs were found to have high consistency in their expression (Figure 4). Similarly, the lncRNA Evx1as, which initiates within the first exon of the gene EVX1, has an overlap of eight nucleotides with the EVX1 mRNA and promotes transcription of its neighbor gene by increasing the binding affinity of histone H3 lysine 4 tri-methylation (H3K4me3) and histone H3 lysine 4 acetylation (H3K27ac) at the promoter region. Considering that most lncRNAs might function through their secondary structure rather than the primary one, this suggests that the regulation of antisense transcripts by antisense-overlapping lncRNA may not be mediated through base-complementary pairing.

Correlation analysis of DE lncRNA and protein-coding genes identified a number of DE lncRNA-DE protein-coding gene pairs. The main enriched KEGG pathways of these protein-coding genes were in metabolism and oxidative phosphorylation. In a recent report, 5-day-old neonatal pigs were infected with PDCoV, and transcriptome profile and KEGG pathway enrichment analysis were performed at different stages of infection (Wu et al., 2019). In our study, we found that the lncRNA targeted genes enriched in those pathways that were perturbed during the late stage of infection. In addition, the expression level of transglutaminase 3 (TGM3) and apolipoprotein A-2 (APOA2) in a Wu et al. (2019) study were significantly changed. Similarly, we also found that TGM1 was up-regulated, and APOA1, APOA4, and APOA5 were down-regulated during PDCoV infection (data not shown). Moreover, our data show that many cytokines and chemokines, which elicit an inflammatory response, were differentially expressed in the infected cells compared to mock cells. The inflammation causes injury to the intestinal tissues, resulting in diarrhea or even death. Raised CCL and CXCL10 levels were associated with the severity of virus infection (Betakova et al., 2017; Masood et al., 2018). Here, we identified a number of lncRNAs that may regulate the expression of these inflammatory molecules.

To the best of our knowledge, this is the first study focusing on the expression profile of cellular lncRNAs after PDCoV infection. Our data show the expression landscape of lncRNAs, with special emphasis on the lncRNA-protein modules operating in response to PDCoV infection. Moreover, this study provides a comprehensive genome-wide resource for exploring the molecular and cellular regulatory functions of lncRNAs. This study will also be useful for identifying lncRNAs as potential biomarkers for the diagnosis of PDCoV infection and designing better prophylactic and therapeutic tools against virus infection.

Conclusion

In the present study, the expression profiles of lncRNAs were determined in PDCoV-infected ST cells. In total, 1,190 novel lncRNAs were identified. A total of 830 lncRNAs were differentially expressed between PDCoV-infected or mocked-infected ST cells. KEGG pathway analysis of DE lncRNA co-expressed genes revealed that they might be primarily involved in regulating metabolism and TNF signaling pathways. Our study systematically characterizes lncRNA expression during PDCoV infection and provides a useful resource for identifying and functionally characterizing the cognate gene products of those lncRNAs. This study will also be useful for assigning lncRNAs as potential biomarkers of PDCoV infection and designing better preventive and therapeutic measures against the virus infection, which would be economically beneficial for the pig farming community.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

Author Contributions

JLL, JG, and JZ conceived and designed the experiments. FW, LD, and JL performed the experiments. JLL, YY, YJ, and TY analyzed the data. JLL drafted the manuscript. All authors read and approved the final manuscript.

Funding

This work was supported by the National Key Research & Development Program of China (2016YFD0500102).

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.

Supplementary Material

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

Footnotes

1. ^ http://www.genome.jp/kegg/

References

Ajayi, T., Dara, R., Misener, M., Pasma, T., Moser, L., and Poljak, Z. (2018). Herd-level prevalence and incidence of porcine epidemic diarrhoea virus (PEDV) and porcine deltacoronavirus (PDCoV) in swine herds in Ontario, Canada. Transbound. Emerg. Dis. 65, 1197–1207. doi: 10.1111/tbed.12858

PubMed Abstract | CrossRef Full Text | Google Scholar

Betakova, T., Kostrabova, A., Lachova, V., and Turianova, L. (2017). Cytokines induced during influenza virus infection. Curr. Pharm. Des. 23, 2616–2622. doi: 10.2174/1381612823666170316123736

PubMed Abstract | CrossRef Full Text | Google Scholar

Carpenter, S., Aiello, D., Atianand, M. K., Ricci, E. P., Gandhi, P., Hall, L. L., et al. (2013). A long noncoding RNA mediates both activation and repression of immune response genes. Science 341, 789–792. doi: 10.1126/science.1240925

PubMed Abstract | CrossRef Full Text | Google Scholar

Clark, B. S., and Blackshaw, S. (2014). Long non-coding RNA-dependent transcriptional regulation in neuronal development and disease. Front. Genet. 5:164. doi: 10.3389/fgene.2014.00164

PubMed Abstract | CrossRef Full Text | Google Scholar

Derrien, T., Johnson, R., Bussotti, G., Tanzer, A., Djebali, S., Tilgner, H., et al. (2012). The GENCODE v7 catalog of human long noncoding RNAs: analysis of their gene structure, evolution, and expression. Genom. Res. 22, 1775–1789. doi: 10.1101/gr.132159.111

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, N., Fang, L., Yang, H., Liu, H., Du, T., Fang, P., et al. (2016). Isolation, genomic characterization, and pathogenicity of a Chinese porcine deltacoronavirus strain CHN-HN-2014. Vet. Microbiol. 196, 98–106. doi: 10.1016/j.vetmic.2016.10.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, P., Fang, L., Hong, Y., Liu, X., Dong, N., Ma, P., et al. (2017). Discovery of a novel accessory protein NS7a encoded by porcine deltacoronavirus. J. Gen. Virol. 98, 173–178. doi: 10.1099/jgv.0.000690

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, P., Fang, L., Ren, J., Hong, Y., Liu, X., Zhao, Y., et al. (2018). Porcine deltacoronavirus accessory protein NS6 antagonizes interferon beta production by interfering with the binding of RIG-I/MDA5 to double-stranded RNA. J. Virol. 92:e00712-18. doi: 10.1128/JVI.00712-18

PubMed Abstract | CrossRef Full Text | Google Scholar

Gomez, J. A., Wapinski, O. L., Yang, Y. W., Bureau, J. F., Gopinath, S., Monack, D. M., et al. (2013). The NeST long ncRNA controls microbial susceptibility and epigenetic activation of the interferon-gamma locus. Cell 152, 743–754. doi: 10.1016/j.cell.2013.01.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Heward, J. A., and Lindsay, M. A. (2014). Long non-coding RNAs in the regulation of the immune response. Trends Immunol. 35, 408–419. doi: 10.1016/j.it.2014.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Huarte, M., Guttman, M., Feldser, D., Garber, M., Koziol, M. J., Kenzelmann-Broz, D., et al. (2010). A large intergenic noncoding RNA induced by p53 mediates global gene repression in the p53 response. Cell 142, 409–419. doi: 10.1016/j.cell.2010.06.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Imamura, K., Imamachi, N., Akizuki, G., Kumakura, M., Kawaguchi, A., Nagata, K., et al. (2014). Long noncoding RNA NEAT1-dependent SFPQ relocation from promoter region to paraspeckle mediates IL8 expression upon immune stimuli. Mol. Cell 53, 393–406. doi: 10.1016/j.molcel.2014.01.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Kapranov, P., Cheng, J., Dike, S., Nix, D. A., Duttagupta, R., Willingham, A. T., et al. (2007). RNA maps reveal new RNA classes and a possible function for pervasive transcription. Science 316, 1484–1488. doi: 10.1126/science.1138341

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S., and Lee, C. (2014). Complete genome characterization of korean porcine deltacoronavirus strain KOR/KNU14-04/2014. Genom. Announc. 2:e01191-14. doi: 10.1128/genomeA.01191-14

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., Chao, T. C., Chang, K. Y., Lin, N., Patil, V. S., Shimizu, C., et al. (2014). The long noncoding RNA THRIL regulates TNFalpha expression through its interaction with hnRNPL. Proc. Natl. Acad. Sci. U.S.A. 111, 1002–1007. doi: 10.1073/pnas.1313768111

PubMed Abstract | CrossRef Full Text | Google Scholar

Lorsirigool, A., Saeng-Chuto, K., Madapong, A., Temeeyasen, G., Tripipat, T., Kaewprommal, P., et al. (2017). The genetic diversity and complete genome analysis of two novel porcine Deltacoronavirus isolates in Thailand in 2015. Virus Genes 53, 240–248. doi: 10.1007/s11262-016-1413-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, S., Lu, J. Y., Liu, L., Yin, Y., Chen, C., Han, X., et al. (2016). Divergent lncRNAs regulate gene expression and lineage differentiation in pluripotent cells. Cell Stem Cell 18, 637–652. doi: 10.1016/j.stem.2016.01.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, Y., Zhang, Y., Liang, X., Lou, F., Oglesbee, M., Krakowka, S., et al. (2015). Origin, evolution, and virulence of porcine deltacoronaviruses in the United States. mBio 6:e00064. doi: 10.1128/mBio.00064-15

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, X., Cai, T., Olyarchuk, J. G., and Wei, L. (2005). Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics 21, 3787–3793. doi: 10.1093/bioinformatics/bti430

PubMed Abstract | CrossRef Full Text | Google Scholar

Masood, K. I., Jamil, B., Rahim, M., Islam, M., Farhan, M., and Hasan, Z. (2018). Role of TNF alpha, IL-6 and CXCL10 in Dengue disease severity. Iran J. Microbiol. 10, 202–207.

PubMed Abstract | Google Scholar

Mattick, J. S., and Rinn, J. L. (2015). Discovery and annotation of long noncoding RNAs. Nat. Struct. Mol. Biol. 22, 5–7. doi: 10.1038/nsmb.2942

PubMed Abstract | CrossRef Full Text | Google Scholar

Mercer, T. R., Dinger, M. E., and Mattick, J. S. (2009). Long non-coding RNAs: insights into functions. Nat. Rev. Genet. 10, 155–159. doi: 10.1038/nrg2521

PubMed Abstract | CrossRef Full Text | Google Scholar

Mercer, T. R., Dinger, M. E., Sunkin, S. M., Mehler, M. F., and Mattick, J. S. (2008). Specific expression of long noncoding RNAs in the mouse brain. Proc. Natl. Acad. Sci. U.S.A. 105, 716–721. doi: 10.1073/pnas.0706729105

PubMed Abstract | CrossRef Full Text | Google Scholar

Ouyang, J., Zhu, X., Chen, Y., Wei, H., Chen, Q., Chi, X., et al. (2014). NRAV, a long noncoding RNA, modulates antiviral responses through suppression of interferon-stimulated gene transcription. Cell Host Microb. 16, 616–626. doi: 10.1016/j.chom.2014.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Ponting, C. P., Oliver, P. L., and Reik, W. (2009). Evolution and functions of long noncoding RNAs. Cell 136, 629–641. doi: 10.1016/j.cell.2009.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Qureshi, I. A., Mattick, J. S., and Mehler, M. F. (2010). Long non-coding RNAs in nervous system function and disease. Brain Res. 1338, 20–35. doi: 10.1016/j.brainres.2010.03.110

PubMed Abstract | CrossRef Full Text | Google Scholar

Saeng-Chuto, K., Jermsutjarit, P., Stott, C. J., Vui, D. T., Tantituvanont, A., and Nilubol, D. (2019). Retrospective study, full-length genome characterization and evaluation of viral infectivity and pathogenicity of chimeric porcine deltacoronavirus detected in Vietnam. Transbound. Emerg. Dis. doi: 10.1111/tbed.13339 [Epub ahead of print].

CrossRef Full Text | PubMed Abstract | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genom. Res. 13, 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., van Baren, M. J., et al. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. doi: 10.1038/nbt.1621

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsai, M. C., Spitale, R. C., and Chang, H. Y. (2011). Long intergenic noncoding RNAs: new links in cancer progression. Cancer Res. 71, 3–7. doi: 10.1158/0008-5472.CAN-10-2483

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, P., Xu, J., Wang, Y., and Cao, X. (2017). An interferon-independent lncRNA promotes viral replication by modulating cellular metabolism. Science 358, 1051–1055. doi: 10.1126/science.aao0409

PubMed Abstract | CrossRef Full Text | Google Scholar

Woo, P. C., Lau, S. K., Lam, C. S., Lau, C. C., Tsang, A. K., Lau, J. H., et al. (2012). Discovery of seven novel Mammalian and avian Coronaviruses in the genus Deltacoronavirus supports bat Coronaviruses as the gene source of Alphacoronavirus and Betacoronavirus and avian Coronaviruses as the gene source of Gammacoronavirus and Deltacoronavirus. J. Virol. 86, 3995–4008. doi: 10.1128/JVI.06540-11

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, J. L., Mai, K. J., Li, D., Wu, R. T., Wu, Z. X., Tang, X. Y., et al. (2019). Expression profile analysis of 5-day-old neonatal piglets infected with porcine Deltacoronavirus. BMC Vet. Res. 15:117. doi: 10.1186/s12917-019-1848-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Q., Chen, C. Y., Yedavalli, V. S., and Jeang, K. T. (2013). NEAT1 long noncoding RNA and paraspeckle bodies modulate HIV-1 posttranscriptional expression. mBio 4:e596-12. doi: 10.1128/mBio.00596-12

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, X., Fang, L., Wang, D., Yang, Y., Chen, J., Ye, X., et al. (2017a). Porcine Deltacoronavirus nsp5 inhibits interferon-beta production through the cleavage of NEMO. Virology 502, 33–38. doi: 10.1016/j.virol.2016.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, X., Wang, D., Zhou, J., Pan, T., Chen, J., Yang, Y., et al. (2017b). Porcine Deltacoronavirus nsp5 antagonizes type I interferon signaling by cleaving STAT2. J. Virol. 91:e00003-17. doi: 10.1128/JVI.00003-17

PubMed Abstract | CrossRef Full Text | Google Scholar