Porcine endemic diarrhea virus infection regulates long noncoding RNA expression.

Long noncoding RNAs (lncRNAs) have been implicated in various life processes. However, the lncRNA expression and potential functions in porcine endemic diarrhea virus (PEDV) infection and host defense are still poorly understood. In this study, we investigated the lncRNA expression profiles during PEDV infection in intestinal porcine epithelial cell-jejunum 2 (IPEC-J2) cell lines by next-generation sequencing and identified 6188 novel lncRNAs. The functional annotation analysis revealed that these lncRNAs might be associated with many immunity-related genes. We next selected candidate lncRNAs related to immune response pathways and further identified their differential expression in PEDV-infected IPEC-J2 cells and newborn piglets. Our results demonstrated that PEDV infection regulated lncRNA expression patterns in both the IPEC-J2 cell line and piglet ileum. These findings provide the first large-scale survey of lncRNAs associated with PEDV infection, specifically the lncRNAs responsible for the activation of the immune system within the ileum.

Introduction

With the development and improvement of next-generation sequencing (NGS) technology, a number of genomes from different species have been annotated and released to the public. Surprisingly, many studies have revealed that the vast majority of genomes can be transcribed into RNAs but cannot be translated into proteins only less than 2% of mammalian genomes are dedicated to protein coding (Carninci et al., 2005).

Noncoding RNAs (ncRNAs), that are transcribed but not translated into proteins include microRNAs (miRNAs), transfer RNAs (tRNAs), ribosomal RNAs (rRNAs), small interfering RNAs, long noncoding RNAs (lncRNAs), and circular RNAs (circRNAs). The ncRNAs transcribed by eukaryotic cells can be divided into two groups, regulatory and housekeeping RNAs depending on their expression pattern and function. The tRNA and rRNA are considered housekeeping RNAs due to their basic and unchanged functions. All miRNAs, circRNAs and lncRNAs belong to regulatory ncRNAs. In addition, lncRNAs are the main ncRNAs, accounting for more than 80% of the total ncRNAs (Kapranov et al., 2007). Therefore, studies on long noncoding RNAs have been important to the transcriptome. Although an accumulating number of lncRNAs have been discovered and found to be vital in many biological processes, the roles that lncRNAs may play during viral infection still await further exploration.

LncRNAs longer than 200 bp are a newly discovered class of ncRNAs and can be distinguished from short RNAs, such as microRNAs, based on length (Nagano and Fraser, 2011). Some lncRNAs can be longer than 100 kb (Ponting et al., 2009). The majority of discovered lncRNAs have a 5′ cap structure and are polyadenylated. Due to the lack of translated proteins, lncRNAs were considered “transcription noise” or “junk RNA”. However, lncRNAs have already been shown to be a key factor in the regulation of gene expression, especially in cancer research (Gong et al., 2016, Huang et al., 2016, Jiang et al., 2016, Zhang et al., 2015). Recent studies demonstrated that viral infection can alter host lncRNA expression profile. More than 500 lncRNAs were observed with differential expression levels in mice after severe acute respiratory syndrome coronavirus (SARS-CoV) infection (Peng et al., 2010). A microarray analysis showed that the expression levels of more than 4800 lncRNAs were changed during enterovirus 71 (EV71) infection (Yin et al., 2013). Some researchers also reported the distinctive kinetics of lncRNA expression profiles during infectious bursal disease virus and Marek's disease virus infections in chicken (He et al., 2015). It was shown that lncRNAs were also involved in HIV and influenza virus infections (Rice, 2015, Trypsteen et al., 2016, Winterling et al., 2014). These studies suggested that lncRNAs may be involved in regulating the interactions between the host response and viral infection.

Porcine epidemic diarrhea (PED) is a highly devastating enteric disease in pigs, characterized by severe enteritis, vomiting and watery diarrhea. Porcine epidemic diarrhea virus (PEDV), as the causative agent of PED, is a member of the Alphacoronavirus genus within the Coronaviridae family (Pensaert and Bouck, 1978). This virus infects pigs of all ages. The mortality rate is associated with the age of pigs, causing an over 90% death rate in piglets up to one week of age. The mortality decreases to 10% in piglets at the age of day 10 or thereafter. PEDV can only cause mild symptoms in older animals. Virulent PEDV strains have reemerged since 2010 and caused significant economic concern to the global pork industry (Chen et al., 2008, Li et al., 2012b, Mole, 2013, Puranaveja et al., 2009, Stevenson et al., 2013). Although the understanding of the molecular virology, epidemiology and vaccinology of PEDV has greatly improved since the reemergence of these strains, the molecular mechanism and transcriptional dynamics of the host response to PEDV infection remain unclear. There are several reports investigating the mechanisms of immune responses and the pathogenesis of PEDV. However, there are no studies on lncRNAs expression and their potential functions of these molecules in virus-host interactions based on a nucleic acid-based centric view.

In the present study, the lncRNA expression profile in a PEDV-infected porcine ileum epithelium cell line (IPEC-J2 cell line) was investigated by NGS technology. We then selected immune signaling pathway-related lncRNAs and identified their differential expression in the IPEC-J2 cell line and in the newborn piglet ileum induced by PEDV infection using real-time RT-PCR. The lncRNA profiles obtained from IPEC-J2 cells and piglets will definitely contribute to and improve the understanding of PEDV-host interactions and immune responses.

Materials and methods

Cell culture and virus infection

IPEC-J2 cells were cultured in DMEM/F12 (Gibco) medium supplemented with 10% fetal bovine serum (Gemini) without bovine viral diarrhea virus (BVDV) and anti-BVDV antibodies, 100 units/ml penicillin, 5 μg/ml insulin/selenium/transferrin (Life) and 5 ng/ml epidermal growth factors (EGF) (Life) in an incubator (5% CO2, at 37 °C). Vero-E6 cells were cultured in DMEM (Gibco) medium supplemented with 10% fetal bovine serum (Gemini) and 100 units/ml penicillin in a 37 °C incubator supplied with 5% CO2.

The PEDV LJX01/GS/2014 strain was isolated from a farm pig and propagated in our laboratory into a cell-adapted virus. To analyze host responses to PEDV infection, 80% confluent IPEC-J2 cells were inoculated with the PEDV LJX01/GS/2014 strain at an MOI of 0.01 or mock infected. Following incubation at 37 °C for 1 h, IPEC-J2 cells were then washed three times with PBS and supplied with culture medium. At different time points post-infection, both the supernatant and cell pellets were collected for further analysis.

Indirect immunofluorescence assay and virus titration

The titers of PEDV were determined by an indirect immunofluorescence assay (IFA). Briefly, an IFA-based viral titration assay was performed by infecting Vero-E6 cells seeded onto 96-well plates (approximately 5 × 104 cells/per well) with 10-fold serial dilutions and four replicates for each dilution. After further 48-h incubation, Vero-E6 cells were washed with PBS, fixed with 4% paraformaldehyde, and permeabilized with 0.1% Triton X-100. The fixed cells were incubated with anti-PEDV-M protein IgG for 2 h at 37 °C, washed three times with PBS, and then incubated with FITC-labeled goat anti-rabbit IgG (1:100 dilution in PBS; Sigma) for 1 h at 37 °C. After washing three times with PBS, the cells were examined under a fluorescence microscope (TE2000U; Nikon) with a video documentation system. Viral titers were calculated in accordance with the Reed-Muench method and expressed as the number of medium TCID50/ml.

For the titration of PEDV in the ileum from each piglet, 3 g of ileum tissue from each piglet was collected and placed into a 2 ml tube with silica beads and 1 ml PBS. The mixtures were then homogenized and centrifuged at 10,000 × g for 10 min at 4 °C. The supernatant was collected and used for virus titration as described above. The culture medium contains 2 μg/ml of trypsin was used to maintain the virus infection.

Library construction and RNA sequencing

Total RNA was extracted by TRIzol reagent (Invitrogen) from both PEDV-infected and mock cells at different time points (12 and 36 hpi). Three biological replications of each group were pooled and subjected RNA-seq analysis. After extraction, the mRNA and noncoding RNA were enriched by removing rRNA from the total RNA with the Ribo-Zero™ Gold Kit. The Agilent 2100 Bioanalyzer and the ABI StepOnePlus Real-Time PCR System were used to quantify and qualify the sample library. The mRNA and noncoding RNAs were then enriched and fragmented into short fragments (approximately 200–500 nt). The first-strand cDNA was synthesized with hexamer random primer using the short fragments as templates. The dTTP was substituted by dUTP when the second-strand cDNA was synthesized. Short fragments were then purified and ligated to adapters, and the second-strand cDNA was degraded using UNG (Uracil-N-Glycosylase)(Parkhomchuk et al., 2009). After agarose gel electrophoresis, suitable fragments were selected as templates for PCR amplification. Finally, the library was sequenced with Illumina HiSeq™ 4000 (Fig. 2A).

Fig. 2

Detailed workflow of lncRNA sequencing and basic properties of lncRNAs. IPEC-J2 cells were plated onto six-well plates and infected with the PEDV LJX01/GS/2014 strain at an MOI = 0.01 or mock infected. A) Total RNA was extracted from the cell pellets at 12 and 36 hpi and subjected to a detailed workflow for lncRNA sequencing and primary analysis. The tools and software SOAP, Tophat, Cufflinks, Cuffcompare and CPC were employed in the whole procedure. FPKM stands for fragments per kilobase of transcript per million mapped reads, which represents the relative expression level of a transcript. B) Pie chart showing the different genomic locations of lncRNAs identified. C) Venn diagram showing the differential expression of these putative 6188 ncRNAs between groups. The overlapping expression profile of predicated lncRNA candidates from each sample is depicted in different colors. D) The lncRNAs and coding RNAs were grouped based on length distributions. The blue columns represent lncRNAs while the light green columns represent coding RNAs.

Genome assembly

After sequencing and data cleaning (Li et al., 2009), the reads were mapped to reference genome (Sus scrofa10.2, Ensembl) using an improved version of TopHat2 (Kim et al., 2013), which can align reads across splice junctions without relying on gene annotation. After assembling with Cufflinks (Trapnell et al., 2010), the whole parsimonious set of transcripts was obtained. To detect the novel transcript from the initial assemblies, the assembled transcripts were compared to the reference annotation by utilizing Cuffcompare (Trapnell et al., 2010).

LncRNA predication and differential expression analysis

For the prediction of lncRNAs, protein-coding transcripts, miRNA, tRNA, snoRNA and rRNA were first filtered according to genome annotation information in Ensembl. Next, the lncRNAs were screened according to the following three characteristics transcripts of nonlncRNAs containing more than one exon; transcript lengths must be longer than 200 bp; and the coverage of each transcript must be more than 3 for each transcript. For the ORF prediction, the software HMMER (v3.1b1) was utilized. The e-value was 1e-5. The CPC (Coding Potential Calculator) was employed to predict protein-coding potential (Kong et al., 2007). Differential expression analysis was performed by using the DESeq Bioconductor package (v1.18.0), a model based on the negative binomial distribution. Briefly, the random number of genes we observed as i. The random number of samples we observed as j. The read count is Kij. Therefore, the negative binomial distribution: Kij ⁓ NB (μij, σij2). After being adjusted by Benjamini and Hochberg's approach for controlling the false discovery rate, the P-values of the genes were set to < 0.05 to detect differentially expressed genes.

Animal experiments

The animal experiments were conducted according to the Guide for the Care and Use of Laboratory Animals of Lanzhou Veterinary Research Institute (LVRI), Chinese Academy of Agricultural Sciences, China. Six newborn piglets derived from the same sow were purchased from a PEDV-free farm. After birth, these piglets were immediately transferred to a laboratory isolator without suckling colostrum and fed with supplemental milk. Three piglets were inoculated with 106 TCID50 of the PEDV LJX01/GS/2014 strain orally. The control groups (3 piglets) were inoculated with an equal volume of PBS orally. At 48 h post-challenge, all piglets were euthanized. The ileum was collected and quickly frozen in liquid nitrogen for detection.

RNA extraction and RT-qPCR analysis

Total RNA was extracted from PEDV-infected and mock-infected IPEC-J2 cells and animal specimens using TRIzol reagent (Invitrogen) and treated with DNase I (TaKaRa) to remove potential genomic DNA contamination. The synthesis of cDNA was performed with 2000 ng of total RNA and 50 U of Moloney murine leukemia virus (M-MLV) reverse transcriptase (Promega), 200 μM dNTPs (TaKaRa), and 10 μl of 5× M-MLV reverse transcriptase buffer in a total volume of 50 μl. The primers for IFN and OASL detection and for validation of differentially expressed lncRNAs related to immune responses (Supplementary Table 1) were designed with Primer 5 software and synthesized from Genewiz Biotech Ltd. Co. (Suzhou, China). Real-time qPCR was performed using the ABI 7500 system with Fast Start Universal SYBR Green Master Mix (Roche), and the primer sets are listed in Supplementary Table 1. Briefly, the reactions were incubated at 95 °C for 10 min, followed by 40 cycles at 95 °C for 15 s and 60 °C for 30 s. All reactions were run in triplicate. The ∆∆Ct method was used to measure the expression levels of target lncRNAs and genes.

Statistical analysis

The differences between matched groups were examined for statistical significance using Student's t-test. An unadjusted P value of less than 0.01 was considered to be highly significant.

Results

PEDV infection activated the immune system

To determine the cellular lncRNA expression profile upon PEDV infection, we first carried out virus replication kinetics to optimize the sample collection time points for lncRNA sequencing. To this end, IPEC-J2 cells were infected with PEDV at an MOI = 0.01 or mock-infected and harvested at 12, 24, 36, 48, 60 and 72 h post infection (hpi). The culture supernatant was collected for virus titration while the cell pellets were used for gene expression analysis. The PEDV initiated its replication soon after infection and maintained a logarithmic phage replication trend between 12 and 48 hpi, peaked at 48 hpi and remained in the plateau phase until 72 hpi ( Fig. 1A). Total RNA was extracted from the cell pellets and subjected to real-time PCR analysis to determine the expression levels of IFN-α, IFN-β and 2′-5′-oligoadenylate synthetase-like protein (OASL), which represents early host immune responses to viral infection. The results showed that the expression level of all three target genes was significantly upregulated, peaking at 36 hpi, and declining to normal levels thereafter (Fig. 1B). All of these data indicated that the host immune system responded to PEDV infection at an early stage of infection. Based on these results, total RNA from 12 hpi and 36 hpi were selected for library construction and lncRNA sequencing.

Fig. 1

PEDV infection activates type 1 IFNs and ISG expression at the early stage. IPEC-J2 cells plated onto six-well plates were infected with the PEDV LJX01/GS/2014 strain at MOI = 0.01 or mock infected. A) The culture supernatant was collected at the indicated time points and titrated for virus growth curve. B) The cell pellets were used for RNA extraction. The cDNA was prepared with hexamer random primers and subjected to real-time PCR analysis for the expression of type 1 IFNs, IFN-α and IFN-β, and OASL at 12, 24, 36, 48, 60 and 72 hpi.

Expression of cellular lncRNAs in PEDV infected IPEC-J2 cells

After RNA extraction, the mRNA and noncoding RNAs were enriched and fragmented into short fragments (approximately 200–500 nt). Then, cDNA was synthesized and detected by agarose gel electrophoresis. The suitable fragments were selected as templates for the PCR amplification. Finally, the library was sequenced with Illumina HiSeq™ 4000 ( Fig. 2A).

Primary sequencing data produced by Illumina Hiseq™ 4000 were first filtered into clean reads by removing adapters and low-quality reads. Approximately fifty-four million valid reads per sample were generated by 150 bp paired-end sequencing after filtering by SOAP analysis (Li et al., 2009). The filtered reads were then mapped to the pig genome (GCA_000003025.4, Ensembl). A total of 6188 novel lncRNAs and 9004 protein-coding transcripts were identified from all samples. The coding RNA was the most abundant RNA identified in our samples, accounting for approximately 59.3% of all sequenced RNAs. The rest RNA molecules were lncRNAs. The genomic locations of these lncRNAs were determined. According to our analysis, 69% of lncRNAs were located in intergenic regions, 4% lncRNAs were antisense transcripts and 7% lncRNAs were generated from intronic regions. The remaining 20% of lncRNAs were overlap-lncRNAs (Fig. 2B). Within the treatment groups at different time points, 5492 lncRNAs were stably expressed in each group and 696 lncRNAs were expressed in different groups (Fig. 2C). The differential expression analysis was performed by using the DESeq Bioconductor package (v1.18.0), which identified 153 differentially expressed lncRNAs at 12 hpi and 218 differentially expressed lncRNAs at 36 hpi. We further analyzed their lengths and found that approximately 18% of lncRNAs were either 200–300 bp in size or longer than 2000 bp. The longest lncRNA was > 16 kb. In contrast, 58% of coding transcripts were longer than 2000 bp (Fig. 2D).

PEDV infection significantly influenced the expression of host genes and lncRNAs

Previous studies have proven that there are strong expression correlations between lncRNAs and neighboring protein-coding genes than for randomly selected protein-coding gene pairs (Cabili et al., 2011, Guttman et al., 2011, He et al., 2015, Pauli et al., 2012, Peng et al., 2010). To test whether lncRNAs were preferentially located in the vicinity of protein-coding genes with specific functions, KEGG and Gene Ontology (GO) enrichment analyses were performed for the neighboring genes within 2 kb of differentially expressed lncRNAs. KEGG enrichment analysis revealed that PEDV infection influenced seven pathways related to immune and antiviral responses, including JAK/STAT, T cell receptor, NF-κB, IgA production, mitogen-activated protein kinase (MAPK), cytokine to cytokine receptor and transforming growth factor-β (TGF-β) signaling pathways ( Fig. 3A). GO analysis showed that seventeen pathways were influenced by PEDV infection (Fig. 3B).

Fig. 3

Host genes and lncRNA signatures in response to PEDV infection. Computational lncRNA functional annotation associated with immune responses was subjected to KEGG A) and GO B) analyses. The x-axis indicates enrichment score [-log10 (P value)]. KEGG analysis was also performed on mRNA data at 12 C) and 36 hpi D). The x-axis indicates enrichment score [log2 (P value)].

We next analyzed the mRNA expression influenced by PEDV infection. The result demonstrated that a total of 9792 genes were expressed during PEDV infection. KEGG analysis revealed that the host cell might utilize JAK/STAT, antigen processing and presentation, and the ubiquitin-proteasome pathways to fight against virus infection at both 12 and 36 hpi (Fig. 3C and D). For the PEDV infection group, the major DE genes were involved in the metabolic pathway, ubiquitin-proteasome pathway and viral reproduction between 12 and 36 hpi (Supplementary Fig. 1). The JAK/STAT pathways were identified from both analyses of the genes in the vicinity of lncRNAs and mRNAs, indicating their vital roles during PEDV infection. Antigen processing and presentation, cytokine-cytokine receptor interaction, and metabolic pathways were also identified in another study (Zhang et al., 2018), suggesting these molecules may also be important to host responses against PEDV infection.

Validation of lncRNA differential expression

Among all DE lncRNAs, 34 lncRNAs that predicted to be involved in JAK/STAT, T cell receptor, NF-κB, IgA production, MAPK, cytokine to cytokine receptor and TGF-β signaling pathways were selected for further analysis. Among these molecules, 18 lncRNAs were differentially expressed at both 12 and 36 hpi, and another 16 lncRNAs were differentially expressed at a single time point. The expression level of their neighboring coding transcripts was also displayed in the heatmap ( Fig. 4A). The FPKM values of both lncRNAs and mRNAs in different groups were compared with each other and transformed to log2. With the combination of two heatmaps, several mRNAs were found to be more active than others. Additionally, these mRNA neighboring regions would generate more lncRNAs. These results further confirmed our previous findings that there were strong correlations between lncRNA expression and their respective neighboring mRNAs.

Fig. 4

PEDV infection regulates lncRNA expression. The lncRNA data related to immune signaling pathways were selected and visualized on a heatmap A). The FPKM value of lncRNA in different groups was compared to each other and transformed to log2. Another heatmap of the expression level of mRNAs in different groups was also constructed according to the same method. Different colors represent the differential expression levels in the different compared groups. The white bar of 22370 shows that they were not detected at 36 hpi. IPEC-J2 cells were infected with PEDV at MOI = 0.01 and collected at 0, 12 and 36 hpi for lncRNA validation by real-time RT-qPCR assay. The validated results from RT-qPCR analysis were normalized by the data at 0 hpi and shown in B) 12 hpi and C) 36 hpi. Different columns represent different immune response-related pathways potentially affect by lncRNAs.

In the RT-qPCR results, we identified twenty lncRNAs displaying temporal expression changes in response to PEDV infection. To further verify PEDV infection to the expression pattern of lncRNAs, the results from PEDV-infected samples at 0 hpi were used as the baseline to normalize the data obtained at 12 and 36 hpi. We found that most of the lncRNAs were upregulated at 12 and 36 hpi while lncRNAs 161 and 162 were downregulated at 12 hpi (Fig. 4B, C). Although lncRNA 159, 160, 161, 162, 164 and 165 shared the same neighboring coding transcript, their expression levels varied from each other. Their divergence was even more significant at 36 hpi. These lncRNAs might play different roles during PEDV infection. The expression levels of lncRNA 9606 and 44432 were upregulated by over 100-fold at both 12 and 36 hpi. Together with lncRNA 88504, 88529 and 251608, these two lncRNAs belonged to the intestinal IgA production and cytokine-to-cytokine receptor signaling pathways. These two pathways were significantly important to intestinal mucosal immunity. Combining the PEDV growth kinetics (Fig. 1B), these findings indicate that PEDV infection regulates host lncRNA expression, especially the lncRNA candidates involved in antiviral responses.

Differential expression levels of lncRNAs in newborn piglets

After validation of the differential expression level of lncRNAs in vitro, we next measured their expression in the ileum of newborn piglets challenged by PEDV. We inoculated 3 newborn piglets with 106 TCID50 of the PEDV LJX01/GS/2014 strain orally and 3 newborn piglets with 1 ml PBS. At 12 hpi, the piglets infected with PEDV displayed vomiting and watery diarrhea. After 24 hpi, the piglets refused to eat and showed dehydration. At 48 hpi, the piglets could not stand and were on the verge of death. The three mock-infected piglets were healthy and had good appetites. We then euthanized the infected and mock-infected piglets and collected the ileum from each piglet for virus titration and lncRNA validation. The results demonstrated that the piglets were successfully infected with PEDV while the mock-treated piglets were not infected ( Fig. 5A). All of the lncRNAs selected were highly upregulated in the ileum of PEDV-infected newborn piglets compared to those in PBS inoculated newborn piglets (Fig. 5B). Both in vitro and in vivo results showed that lncRNA-9606 was most significantly upregulated, with a 300- to 1300-fold increase (Fig. 4 and Fig. 5). This lncRNA was predicted to be involved in the intestinal immune network for IgA production. IgA was indicated to be the most important antibody against PEDV infection in the small intestine. All other lncRNAs were upregulated approximately 60-fold. LncRNAs predicted to be involved in the T cell receptor signaling pathway showed higher levels than the lncRNAs predicted to be involved in other pathways related to immune responses. These results suggest that the lncRNAs detected in this study may play important roles in host antiviral responses. Additionally, the T cell receptor signaling pathway and IgA production may be the key pathways involved in the activation of the immune system in the intestine.

Fig. 5

PEDV infection regulates lncRNA expression in newborn piglets. Newborn piglets were inoculated with 106 TCID50 of the PEDV LJX01/GS/2014 strain or mock infected and sacrificed at 48 h after challenge. A) The viral titers in the ileums from both infected- and mock-infected piglets. B) The expression level of selected lncRNAs in the ileums. The ileums of each piglet were collected for the detection of immune response-related lncRNAs by real-time RT-qPCR assay. Different columns represent different immune response-related pathways affect by lncRNAs.

Discussion

In this study, NGS was used to investigate the whole transcriptome of porcine ileum epithelium cells over two time points corresponding to critical phases of PEDV infection. More than 6100 candidate lncRNAs were identified. These lncRNAs were characterized by lower expression levels and shorter transcript lengths. These properties were also observed in other reported lncRNAs within the mammalian genome (Cabili et al., 2011, He et al., 2015, Peng et al., 2010, Ulitsky et al., 2011, Young et al., 2012, Boltana et al., 2016, Li et al., 2012a, Okazaki et al., 2002, Ulitsky et al., 2011). In addition, GO term enrichment analysis and KEGG pathway analysis indicated that various immune signaling pathways and antiviral responses may be regulated by cellular lncRNAs in response to PEDV infection. Real-time PCR evaluated the differential expression of twenty lncRNAs at different time points in vitro and in vivo. These findings suggest that PEDV infection could alter the expression of cellular lncRNAs; thus, these lncRNAs may have potential to regulate host immune signaling pathways and antiviral responses.

Evolutionary studies have demonstrated that humans and other animals may have the same number of protein-coding genes. However, the size of the genome varies from one species to another (Kapusta and Feschotte, 2014). The number of ncRNAs is the key factor that determines the size of the genome. The genome of primitive species such as worms and yeast contains a much smaller number and a lower ratio of lncRNAs to protein-coding genes than advanced species. Thus, the more complicated the life activity is, the more lncRNAs their genome contains and these molecules are involved in nearly every process of the life cycle.

The immune system is gradually established and perfected in biological and evolutionary processes. Invertebrates only have nonspecific phagocytic function and inflammatory reaction abilities. In invertebrate lymphocytes, the degree of differentiation and division is still very poor. In mammals, there are specific immune organs, immune cells and immune molecules for identifying and eliminating antigens. The increasing number of lncRNAs may be partly related to the immune system in higher organisms. In our study, many lncRNAs were coexpressed with genes that are involved in host immune signaling and antiviral responses, including JAK/STAT, NF-κB, MAPK, cytokine to cytokine receptor, TGF-β, T cell receptor signaling pathway and IgA production signaling pathways, which supports the above hypothesis. The former three pathways are related to the adapted immune response, and the latter four pathways play important roles in the mucosal immune response.

The whole intestine is covered with mucous membranes, which are responsible for interfacing with the outside environment and modulating the immune response of an organism to microbes approaching the mucosal surfaces. Thus, the mucosal immune response is the most powerful defense system against PEDV. IgA is secreted by the mucosal immune system for neutralizing antigens (Blutt and Conner, 2013). The activation of cytokines, TGF-β and T cells is essential for the production of IgA. TGF-β1 plays a pivotal role in intestinal cellular and immune regulation (Min and Kim, 2003, Ohtsuka and Sanderson, 2000). In the presence of TGF-β, peritoneal B cells switch to IgA more potently than other peritoneal B cell subpopulations (Roy et al., 2013). With the help of TGF-β, the costimulation from the CD4+ T helper cell subset via the CD40-CD40L interaction contributes to the switch from an IgM to an IgA isotype (Cazac and Roes, 2000, Tsuji et al., 2009). Various cytokines, such as IL-13, IL-17A and IL-21, also contribute greatly to this process (Cao et al., 2015, Dann et al., 2015, Wu et al., 2016). LncRNAs 159, 160, 164, 165, coexpressed with genes involved in the T cell receptor signaling pathway, were all upregulated at 12 hpi and 36 hpi in IPEC-J2 cells and in the ileum. LncRNA-9606, which are coexpressed with a gene related to IgA production was significantly up-regulated in both PEDV-infected IPEC-J2 cells and piglets. This result strongly proves that PEDV infection may activate the specific expression of several lncRNAs and that these lncRNAs may potentially activate the immune response.

It is of great interest to find that lncRNA 9606 was upregulated 400-fold in the ileum and at least 300-fold in the IPEC-J2 cells during PEDV infection, which showed the highest upregulation among all identified lncRNAs in this study. Previous studies have shown that there are strong expression correlations between lncRNAs and neighboring protein-coding genes than for randomly selected protein-coding gene pairs (Cabili et al., 2011, Guttman et al., 2011, He et al., 2015, Pauli et al., 2012, Peng et al., 2010). We searched the neighboring coding transcript of lncRNA 9606 and found that the T cell receptor alpha chain V region was located within a 2 kb region from lncRNA 9606. Additionally, the KEGG analysis indicated that lncRNA 9606 was included in the intestinal immune network for the IgA production pathway. Therefore, we speculated that lncRNA 9606 may activate the TCR to sense the antigen and then promote IgA production. We performed an alignment search and the results indicated that lncRNA 9606 shared a query coverage of over 60% but lower than 75% between various genes, such as adenylate kinase 3-like 1 (AK3L1), interferon α receptor 1 (IFNAR1), solute carrier family 46 member 3 (SLC46A3), and toll-like receptor 2 (TLR2). IFNAR1 and TLR2 play important roles in host defense against viral infection. LncRNA 9606 may also potentially function in these processes. Notably, to guarantee the reliability of lncRNA detection, all primer sets were designed according to their unique regions to avoid the interference of other similar genes.

LncRNAs regulate gene expression by multiple mechanisms. After transcription, lncRNAs can be associated with chromatin-modifying enzymes and transcription factors or form ribonucleoprotein complexes (Mattick et al., 2009, Mercer et al., 2009, Nagano and Fraser, 2011, Satpathy and Chang, 2015). Although proteins are thought to be the executors of life activities, lncRNAs may perform much better than proteins. These molecules can regulate gene expression much faster than proteins because the production of RNAs needs only transcription. A correlation of expression between lncRNAs and their neighboring protein-coding genes has been widely described in previous studies (Cabili et al., 2011, Guttman et al., 2011, He et al., 2015, Pauli et al., 2012, Peng et al., 2010) and observed in our study. The association of lncRNAs with chromatin-modifying complexes has been recognized as the reason for the coexpression (Khalil et al., 2009, Peng et al., 2010). By using this method, lncRNAs may regulate the expression of their neighboring genes. Once the lncRNAs are knocked down, they can change the neighboring gene expression patterns (Guttman et al., 2011). This mechanism is indeed a fast and reliable method for the regulation of gene expression.

Among the identified lncRNAs, the differential expression of lncRNAs was not only affected by PEDV infection. There were also many differentially expressed lncRNAs identified between mock-infected cells of different culture times (Fig. 3B). This effect may be due to the time and spatial-temporal expression of lncRNAs. The expression of lncRNAs can be affected by many factors, such as cell type, tissue type, development stage, culture time and stimulants (Cabili et al., 2011, Gaiti et al., 2015, Guo et al., 2016, Hart and Goff, 2016, Pauli et al., 2012, Peng et al., 2010, Zhao et al., 2016). The reasons for this effect can be explained as above. Every life activity needs the participation of lncRNAs.

Previous studies focusing on virus-host interactions were mostly based on a protein-centric view. Since the discovery of miRNA and lncRNA, there has been increasing evidence suggesting that RNA could also directly or indirectly interact with viruses. Considering the large number and rapid transcription of RNA, these molecules may play a greater role in the regulating host-virus response. In this study, we first employed RNA-seq technology to identify lncRNAs differentially expressed during PEDV infection. Several lncRNAs potentially involved in the host response to PEDV infection were selected and identified by RT-qPCR. The JAK/STAT signaling pathway, cytokine-cytokine receptor interaction, MAPK signaling pathways, and NF-κB signaling pathway were found to be regulated during PEDV infection in our study and the proteome study (Lin et al., 2017). Our study further discovered a number of DE lncRNAs that are potentially involved in the TGF-beta signaling pathway, T cell receptor signaling pathway and intestinal immune network for IgA production. These data have provided novel targets for studies on host-virus interactions. Therefore, an RNA-centric view of virus-host interactions should also be considered for a better understanding of virus infection. To our knowledge, this study is the first to illustrate the expression profile of lncRNAs in newborn piglets against PEDV infection.

In conclusion, we characterized and identified the lncRNA expression profiles in PEDV-infected IPEC-J2 cells and newborn piglets and discovered a total of 6188 novel porcine lncRNAs, which would provide a basis for understanding the possible roles of lncRNA during virus-host interactions.