Identication and Characterization of MRNAs, lncRNAs, and MiRNAs Expression Proles and ceRNA Networks in PEDV Infection

Background: Long non-coding RNAs (LncRNAs) are transcripts longer than 200 nucleotides with no protein-coding ability and exert crucial effects on viral infection and host immune responses. Porcine Epidemic Diarrhea Virus (PEDV) is a coronavirus that seriously affects the swine industry. However, our understanding of the function of lncRNA involved in host-PEDV interaction is limited. Results: A total of 1197 mRNA transcripts, 539 lncRNA transcripts, and 208 miRNA transcripts were differentially regulated at 24 h and 48 h post-infection. Moreover, gene ontology (GO) and KEGG pathway enrichment analysis showed that DE mRNAs and DE lncRNAs were mainly involved in biosynthesis, innate immunity, and lipid metabolism. Ten differentially expressed genes were randomly selected and validated by reverse-transcription qRT-PCR. In addition, we constructed a miRNA-mRNA-pathway network followed by a lncRNA-miRNA-mRNA ceRNA network. Conclusions: The present study is the rst to reveal the global expression proles of mRNAs, lncRNAs, and miRNAs during PEDV infection. We comprehensively characterize the ceRNA networks which can provide new insights into the pathogenesis of PEDV.

LncRNAs are a class of RNAs that are more than 200 nucleotides (nt) in length and have no proteincoding capability [17]. Compared with mRNA, lncRNA is more tissue-speci c and less conservative [18].
Several studies have shown that viral infection induces lncRNAs to regulate antiviral innate immunity [19]. LncRNA NEAT1 exerts an antiviral effect by promoting IFN response during Hantavirus infection [20].
On the other hand, another study has shown that lncRNA-CMPK2 acts as a negative regulator of IFN response to promote HCV replication [21]. In recent studies, the competing endogenous RNA (ceRNA) theory has been reported as a new regulatory mechanism between lncRNA and mRNA [22]. LncRNAs act as miRNA sponges to indirectly regulate the degradation or translational inhibition of target mRNA [22].
So far, ceRNA theory has become the mainstream method for studying lncRNA-miRNA-mRNA interaction and has been veri ed by many physiological processes. Lnc-ISG20 binds to miR-326 to enhance the expression of ISG20 and inhibit IAV replication [23]. However, there are no studies on the lncRNA-miRNA-mRNA network during PEDV infection.
In this study, we performed the rst systematic analysis of mRNA, lncRNA, and miRNA expression pro les at two different time points during PEDV infection. We analyzed differentially expressed genes, identi ed the target genes of DE lncRNAs and DE miRNAs, and characterized the mRNA-miRNA-pathway network and immune-related ceRNA network. Our study aims to explore the potential role of ncRNAs during PEDV infection and expand the understanding of the host-virus interactions.

Results
Overview of RNA-sequencing To identify the ncRNAs (lncRNA and miRNA) and protein-coding transcripts associated with viral infections, 18 RNA-seq libraries were sequenced. A total of 995M raw reads were obtained and an average of 983.63 M clean reads per library was produced (Table S1). In addition, at least 89% of transcripts in each library were covered by reads (Table S2). Meanwhile, approximately 71.33M clean reads were obtained from 9 small RNA sequencing libraries (Table S3). In this study, we identi ed 17,500 mRNAs, 5,748 lncRNAs, and 659 miRNAs. The RNA analysis suggested that, except for miRNA, the distribution of mRNA and lncRNA in each treatment group was essentially the same ( Fig. 1a-c, Fig. S1).
The proportion of miRNA expression <1 was greater than 41% in the 48hpi-vs-Mock group. Following that, the properties of mRNA and lncRNA transcripts were assessed. The number of mRNA exons was noticeably greater than that of lncRNA exons. The ndings revealed that mRNAs with more than 7 exons exceeded 60%, while lncRNAs mostly contained only two or three exons (Fig. 1d). The average length of these mRNAs was longer than the average length of lncRNAs, and the majority of the lncRNAs were longer than 2000 bp (Fig. 1e). Furthermore, miRNA length analysis revealed that it was commonly distributed at 22nt (Fig. 1f).

Identi cation and analysis of differentially expressed mRNAs
In this study, we identi ed 660 and 669 DE mRNAs in two groups, including 813 up-regulated (61.17%) and 516 down-regulated (38.83%) ( Fig. 2a and b, Table S4). Next, the hierarchical cluster analysis revealed that the expression pro les of mRNAs changed signi cantly during PEDV infection, and 132 common DE mRNAs were found between the two groups through Venny 2.1 (Fig. 2c-e). GO enrichment analysis suggested that these DE mRNAs were signi cantly enriched in the biosynthetic process, metabolic process, cell-cell signaling, and binding ( Fig. 3a and b). Meanwhile, KEGG analysis indicated that several DE mRNAs were enriched in signal transduction pathways, including the Rap1 signal, Wnt signal, cGMP-PKG signal, TGF-beta signal, and NOD-like receptor signaling pathway ( Fig. 3c and d).

Identi cation and analysis of miRNAs
To analyze the characteristics of miRNAs during PEDV infection, miRNA expression levels were compared in different groups. We detected 43.1% mature miRNAs after comparing the clean reads to the GenBank, Rfam database, and reference genome. (Fig. 4a). A total of 616 and 641 miRNAs were detected from 24hpi-vs-Mock and 48hpi-vs-Mock, respectively. Among them, 598 miRNAs were co-expressed, while 18 and 43 miRNAs were speci cally expressed in the two groups. According to fold change > 2 and P-value < 0.05, there were 3 and 206 DE miRNAs, including 100 up-regulated (31.9%) and 109 down-regulated (68.1%) were identi ed ( Fig. 4b-d, Fig. S2a and b, Table S5). A total of 7,352 target mRNAs for DE miRNAs were identi ed, and we identify the potential physiological functions of DE miRNAs through GO and KEGG enrichment analyses. These genes were mainly related to cellular components and molecular functions, including intracellular signal transduction, single-organism component, and binding ( Fig. S3a  and b). KEGG results showed that the enriched pathways involved Endocytosis, FoxO signal, MAPK signal, cAMP signal, and apoptosis ( Fig. S3c and d).
Construction of miRNA-mRNA-pathway network. Among intersecting genes, 12 genes related to cellular immunity and cell growth along with their corresponding 120 DE miRNAs were selected (Table 1) Table 1 (Fig. 6). Table 1 Twelve intersection genes and their corresponding pathways and differentially expressed miRNAs Identi cation and analysis of lncRNAs In this study, CNCI, CPC, Pfam-scan, and PLEK were used to remove potential coding transcripts, and 5,769 lncRNAs were obtained. Through further ltering, we nally identi ed 5,748 lncRNAs, which included 37.2% intergenic lncRNAs and 40.4% intronic lncRNAs ( Fig. 5a and b). Most of the lncRNA transcripts were found on chromosome 20, however, few lncRNAs were found on chromosome Y (  Table S6).
In addition, 1,197 protein-coding genes were found to exhibit signi cant correlations with DE lncRNAs in expression, of which 1,169 were detected as trans-regulated target genes. GO terms analysis was performed, and many terms related to the biosynthetic process and immunity were signi cantly enriched ( Fig. S4a and b). Meanwhile, KEGG analysis showed that DE lncRNAs were mainly involved in Toll-like receptor signal, Jak-STAT signal, NF-kappa B signal, and TGF-beta signal ( Fig. S4c and d).
Functional enrichment analysis revealed that these DE mRNAs were involved in several immune-related pathways, including apoptosis, NOD-like receptor signaling, PI3K-Akt, NF-kappa B signaling, and RIG-I-like receptor signaling pathway.

Discussion
PEDV caused huge economic losses to the livestock industry, and commercial vaccines cannot provide complete protection for pigs [7]. Although research on PEDV has made great progress in recent decades, the pathogenic mechanism of PEDV has not been fully studied. Recent research has demonstrated that lncRNAs can act as miRNA sponges, regulating the expression of host genes and hence in uencing virus replication. However, the characterization of the ceRNA mechanism in host cells during PEDV infection has not been elucidated. In this study, we detected changes in the expression pro les of mRNAs, lncRNAs, and miRNAs in PEDV-infected MARC-145 cells. The highest number of DE lncRNAs was discovered at 48 hpi, and there were more up-regulated genes were identi ed. These characteristics reveal that the host response to PEDV infection involves activating a large number of genes. Comparative analysis of the genomic characteristics of mRNAs and lncRNAs reveals that lncRNAs have shorter transcripts, fewer exons, and lower expression levels, which is consistent with the results of transcriptome analysis in different species [24][25][26]. In addition, the expression changes of randomly selected differentially expressed genes were con rmed by qPCR analysis, indicating the accuracy of the sequencing results.
Recent studies suggest a close relationship between cholesterol metabolism and innate immunity [27][28][29]. In addition, the ATP binding cassette (ABC) transporters play a key role in many metabolic processes and contribute to resisting the penetration of xenobiotic [30]. By analyzing these DE mRNAs, we found that the expression levels of ABCG1, ABCC3, LPIN1, DHCR7, HMGCR, and INSIG1 genes were signi cantly up-regulated during PEDV infection. ABCG1 and ABCC3 are members of the superfamily of ABC transporters. ABCG1 can catalyze the e ux of cholesterol and plays a role in cellular lipid homeostasis [31]. The lipin protein family contains three members, among which LPIN1 is a key gene that acts as a magnesium-dependent phosphatidic acid phosphatase to regulate fatty acid metabolism [32]. DHCR7 catalyzes the nal step of cholesterol biosynthesis, which is of great signi cance in maintaining cholesterol homeostasis [33]. Xiao et al. reported that reducing the expression of DHCR7 or adding 7dehydrocholesterol enhanced macrophage-mediated immune responses via increased IFN-I production [34]. INSIG1 is an endoplasmic reticulum membrane protein that is mainly involved in lipogenesis, which is combined with RNF145 to ubiquitinate and degrade HMGCR, thereby regulating cholesterol synthesis [35]. These genes related to cholesterol metabolism were signi cantly up-regulated, suggesting that the host may resist PEDV infection by regulating lipid metabolism.
We also studied the changes in the expression of genes associated with immunity. The sequencing results revealed that toll-like receptors (TLR2 and TLR10), some in ammatory factors (IL24, IL1RN, and IL37), and several antiviral genes such as IFIT1, ISG15, and TNK1 were differentially expressed. These upregulated DE mRNAs were mainly focused on the immune-related pathways, such as Toll-like receptor signaling, p38 MAPK signaling, IFN signaling, and TNF signaling pathway. IFIT1 is an interferon-induced antiviral RNA binding protein that reduces the viral mRNA production by speci cally binding singlestranded RNA with a 5'-triphosphate group [36]. In addition, IFIT1 can also bind 5'capped 2'-O unmethylated mRNA to inhibit mRNA translation and thus inhibit virus replication [37]. TNK1, a member of the Ack kinase family, is an essential regulator of intracellular signal transduction pathways, mediating important physiological activities such as cell proliferation and development [38]. Activated TNK1 inhibits virus replication by enhancing JAK-STAT signaling to induce the expression of ISGs [39]. The data presented here may indicate that the host resists viral infection by up-regulating the genes described above in 48h. As we all know, PEDV can use a variety of strategies to antagonize innate immune responses to ensure self-proliferation. However, no signi cant changes in mRNA levels of typical immune signaling factors such as MAVS, TRAF3, and IRF3 were detected during the early stages of PEDV infection. Therefore, we hypothesize that this may be one of the ways PEDV evades the host innate immune response. Further in vivo and in vitro studies are required to con rm these ndings. Chen et al. compared the lncRNA pro les in PEDV infected IPEC-J2 cells at different time points and obtained the DE lncRNAs related to immunity [41]. There are several differences between our study and other studies. Firstly, we conducted a synthetic analysis of the mRNA, lncRNA, and miRNA expression pro les of PEDV-infected cells and performed KEGG and GO enrichment analysis. Secondly, we predicted the DE mRNAs and DE lncRNAs interacting with DE miRNAs through the database, and the ceRNA network was built. Subsequently, the enrichment analysis of DE mRNAs and DE lncRNAs showed the enrichment of cell growth, steroid synthesis, fatty acid synthesis, and other lipid metabolism pathways, which is different from previous reports.
In our study, 660 and 669 signi cantly DE lncRNAs were obtained at two different time points in PEDVinfected cells. We found that the position of lncRNA in the genome is non-random, which indicates a connection between the lncRNA function and genomic position. Ørom et al. found that a group of human lncRNAs can perform enhancer-like functions and increase the expression of neighboring protein-coding genes through cis-acting [42]. Furthermore, another report clari ed that lncRNA-mediated cis-regulation of adjacent gene transcription is common in mammalian gene regulation [43]. Therefore, knowledge of the cis-regulation of lncRNA is essential to reveal its biological effects. In our study, several lncRNAs (TCONS_00031360, TCONS_00030751, TCONS_00005777, TCONS_00021191, etc.) were involved in the cis-regulation of autophagy-related genes during PEDV infection. IFN also plays a vital role in innate immunity defense, and several studies have reported that numerous viral proteins of PEDV have been identi ed as IFN antagonists [7]. Our data also indicated that some lncRNAs targeted immunity-related genes such as DDIT3, TXNIP, and FOSL1. DDIT3 is a stress-inducible nuclear protein that induces apoptosis in response to ER stress. Moreover, it inhibits the antiviral immune response by promoting the degradation of MAVS [44]. Although lncRNA-mediated cis-regulation is ubiquitous, some lncRNAs can also function in trans. For example, the expression of several lncRNAs (TCONS_00039521, TCONS_00019044, TCONS_00024478, and TCONS_00021764) was signi cantly increased, and their target genes were ADA2, FADS1, MEA1, ADGRE1, MAN1C1, PNPLA3, and RND2, which were located in the trans-position.
Several studies have demonstrated the extensive involvement of miRNAs in viral pathogenesis as well as in multiple immune responses. In our study, we found that the NF-κB signaling pathway was shown to be regulated by miR-16-5p and miR-24-3p. It has been reported that during IBDV infection, up-regulated miR-16-5p can enhance IBDV-induced apoptosis by inhibiting Bcl-2 and promoting IBDV replication [45]. A similar promotion of virus replication by miR-24-3p was observed during PRRSV infection [46]. Upregulated miR-24-3p inhibited the expression of Heme oxygenase-1 (HO-1) through targeting its 3'UTR, facilitating PRRSV replication. It is universally known that IRF3 is a member of the IRF family that can induce the expression of type I IFN to regulate the antiviral innate immune response. Interestingly, our sequencing data showed that the expression levels of miR-296-3p, miR-486-3p, miR-940, and miR-1306-5p, whose target gene is IRF3, were signi cantly down-regulated. Therefore, we hypothesize that the host innate immune system was effectively activated to down-regulate miRNA expression in response to viral infection. Furthermore, several DE miRNAs targeting immune-related genes were identi ed, but with different expression trends. This nding also indicates that these miRNAs are differentially regulated and perform diverse functions during viral infection. In addition, we identi ed several DE miRNA target genes enriched in lipid-associated pathways, such as fatty acid biosynthesis and metabolism, lipid binding, and PPAR signaling pathway.
More and more research has shown that lncRNAs can act as miRNA sponges, and the construction of ceRNA networks is a new strategy to investigate gene function. Therefore, we selected 3 DE miRNAs to build a ceRNA network to explore the exact roles and mechanisms of DE lncRNAs. Integration analysis of the lncRNA-associated ceRNA network revealed the functions of lncRNAs and their roles in cell proliferation and innate immune response. We identi ed many lncRNA-miRNA-mRNA regulatory axes that may be involved in the regulation of PEDV replication, for example lnc_00037724-miR-455-5p-MEFV, lnc_00036554-miR-27b-3p-ID3, and lnc_00010748-miR-145-5p-TNFAIP3. It has been reported that MEFV regulates innate immunity by degrading multiple in ammasome components, including CASP1, NLRP1, and NLRP3 [47]. Other studies have documented that ID3 is involved in several biological processes, including cell proliferation, senescence, and apoptosis [48]. TNFAIP3 is a ubiquitin editing enzyme with immunosuppressive properties, preventing NF-κB activation and TNF-mediated apoptosis [49,50]. Furthermore, several lncRNA-miRNA-mRNA axes with similar expression patterns were identi ed, which were enriched in different biological processes. Although the results predicted by the software need to be further investigated, our ndings provide new evidence for the involvement of lncRNAs in host antiviral immunity.

Conclusions
In summary, we constructed the lncRNA-associated ceRNA network related to cellular immunity by comprehensive bioinformatics analysis. To our knowledge, this is the rst comprehensive analysis report about the expression pro les of mRNAs, lncRNAs, and miRNAs during PEDV infection. Therefore, the data acquired in this study provides a bene cial reference for understanding the interaction between the host and PEDV.

Cells culture and virus infection
MARC-145 cells were maintained in Dulbecco's modi ed Eagle's essential medium (DMEM; Hyclone, USA) supplemented with 10% fetal bovine serum (FBS; Gibco, USA) and cultured at 37 °C with 5% CO 2 . PEDV strain CH/SXYL/2016 (GenBank: MF462814.1) used in this study was isolated as previously described [51]. The cells were infected with PEDV at an MOI of 1 and incubated in DMEM containing 2μg/mL Trypsin for 1h at 37 °C with 5% CO 2 . After incubation, the cells were washed with phosphate-buffered saline (PBS) and cultured in DMEM supplemented with 2% FBS.

RNA extraction
MARC-145 cells were divided into three groups: cells infected with PEDV for 24 hours (S24hpi_1, 2, 3), cells infected with PEDV for 48 hours (S48hpi_1, 2, 3), and mock-infected cells as controls (Mock_1, 2, 3). According to the manufacturer's instructions, the total RNA of PEDV-infected cells and uninfected cells was extracted using the TRIzol ® reagent (Invitrogen, USA). RNA quality and concentration were tested by NanoDrop 2000 (Thermo Scienti c, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies, USA). The samples with RNA Integrity Numbers (RIN) ≥ 7 were subjected to the subsequent analysis.

Library construction and sequencing
A total of 1μg RNA from each sample was used as input for RNA-seq library preparation. According to the manufacturer's instructions, ribosomal RNA (rRNA) was removed using an Epicenter Ribo-Zero rRNA Removal Kit (Illumina, USA), and small RNA libraries were constructed using TruSeq Small RNA Sample Prep Kit (Illumina, USA). Meanwhile, sequencing libraries of mRNA and lncRNA were generated using NEBNext® Ultra™ II Directional RNA Library Prep Kit for Illumina® (NEB, USA). Thereafter, the library products were assessed for purity using an Agilent 2100 Bioanalyzer and then sequenced using the Illumina HiSeq X Ten platform by OE Biotech Co.

Analysis of differentially expressed mRNAs
First, low-quality reads were removed from the raw data using Trimomatic (v0.36) [52] to obtain clean reads. Then, the clean reads were mapped to the Chlorocebus sabaeus genome using HISAT2 (v2.2.1) [53], and the mapped reads of each sample were assembled using StringTie (v2.1.1) [54]. Cu inks (v2.2.1) [55] was used to determine the FPKM [56] of each gene, and HTSeq-count (v0.9.1) [57] was used to obtain the read counts of each gene. Analysis of DE mRNAs was performed using the DESeq R package. Fold change > 2 and P-value < 0.05 were set as the thresholds for signi cantly differential expression.
miRNA identi cation and analysis Cutadapt (v1.14) was used to remove adapter dimers, repetitive sequences, and sequences shorter than 18 nucleotides in the raw data. The clean data was aligned by ACGT101-miR (v4.2) to eliminate RNA families (rRNA, tRNA, snRNA, snoRNA) and repeat sequences. Then, the clean reads were mapped to the miRBase database (v22.0) to identify known miRNAs and novel miRNAs. The miRNA expression was calculated and normalized as transcripts per million (TPM), and DE miRNAs were selected according to fold change > 2 and P-value < 0.05.

lncRNA identi cation and analysis
The raw data were rst ltered by Trimomatic (v0.36) to remove low-quality reads and adapter sequences, which were then mapped to the Chlorocebus sabaeus genome using HISAT2 (v2.2.1). Afterward, transcripts that overlapped with known mRNAs and shorter than 200 bp were discarded. Next, we utilized CNCI (v1.0), CPC (0.9-r2), Pfam-scan (v30), and PLEK (v1.2) to predict the coding ability of transcripts. Transcripts with coding potential predicted by the four tools listed above were discarded, and those without coding potential were classi ed as lncRNAs. The target genes of lncRNAs were predicted using cis-and trans-regulation analyses [58]. Cu inks (v2.2.1) was used to calculate the expression level of each transcript by calculating FPKM. DE lncRNAs were selected with fold change > 2 and P-value < 0.05 by DESeq R package.
The DAVID online tool was used to perform Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of differentially expressed genes and target genes. GO terms and KEGG pathways with P-value < 0.05 were considered to be signi cantly enriched.
Construction of lncRNA-miRNA-mRNA network According to the ceRNA hypothesis, lncRNAs act as miRNA sponges to relieve the inhibition of the expression of target mRNAs [22]. TargetScan (v7.2), miRTarBase (v8.0), and miRanda (v 3.3a) were employed to predict the relationship between miRNAs and mRNAs as well as between lncRNAs and miRNAs. Cytoscape (v3.5.1) was used to create and visualize the lncRNA-miRNA-mRNA networks.

Quantitative Real Time-PCR (qRT-PCR) veri cation
To validate the data obtained from RNA-seq, qRT-PCR was performed for the randomly selected RNAs. To detect the differentially expressed mRNAs and lncRNAs, total RNA from PEDV-infected and mock-infected cells was extracted using TRIzol reagent (Invitrogen, USA), and reversely transcribed using the Fast Quant RT Kit (Tiangen, China). Then, qRT-PCR assays were performed with SYBR Green qPCR Master Mix (Bio-Rad, USA). The β-actin gene was selected as the internal control for the normalization of mRNA and lncRNA. The relative quanti cation of mRNAs and lncRNAs was determined using the 2 -ΔΔct method. The sequence of the primers used in this study is listed in Table S7.

Statistical analysis.
GraphPad Prism 6.0 software was used for statistical analyses. Differences were determined to be statistically signi cant at values of P < 0.05 by Student's t-test for two groups and one-way ANOVA analysis for more than two groups. (Data are shown as means ± standard deviation (SD). n = 3).       The miRNA-mRNA-pathway network. The ellipse, round rectangle, and diamond symbols represent mRNA, miRNA, and the pathway, respectively. Up-regulated and down-regulated genes are represented by red and green, respectively.

Figure 7
An overview of the lncRNA-miRNA-mRNA ceRNA network. The ellipse, triangle, and round rectangle symbols represent mRNA, lncRNA, and miRNA respectively.