Identification of lncRNAs Involved in PCV2 Infection of PK-15 Cells

Porcine circovirus type 2 (PCV2) can cause severe disease in infected pigs, resulting in massive economic loss for the swine industry. Transcriptomic and proteomic approaches have been widely employed to identify the underlying molecular mechanisms of the PCV2 infection. Numerous differentially expressed mRNAs, miRNAs, and proteins, together with their associated signaling pathways, have been identified during PCV2 infection, paving the way for analysis of their biological functions. Long noncoding RNAs (lncRNAs) are important regulators of multiple biological processes. However, little is known regarding their role in the PCV2 infection. Hence, in our study, RNA-seq was performed by infecting PK-15 cells with PCV2. Analysis of the differentially expressed genes (DEGs) suggested that the cytoskeleton, apoptosis, cell division, and protein phosphorylation were significantly disturbed. Then, using stringent parameters, six lncRNAs were identified. Additionally, potential targets of the lncRNAs were predicted using both cis- and trans-prediction methods. Interestingly, we found that the HOXB (Homeobox B) gene cluster was probably the target of the lncRNA LOC106505099. Enrichment analysis of the target genes showed that numerous developmental processes were altered during PCV2 infection. Therefore, our study revealed that lncRNAs might affect porcine embryonic development through the regulation of the HOXB genes.


Introduction
Porcine circovirus (PCV), which belongs to the family Circoviridae, is a nonenveloped virus with a single-stranded circular DNA [1]. Currently, there are three known types of PCVs. PCV1, which was first identified in the pig kidney cell line, is nonpathogenic to pigs [2,3]. PCV3 is a newly identified pathogen with the first report in North America. The afflicted pigs have symptoms of porcine dermatitis and nephropathy syndrome and reproductive problems [4]. Another PCV is porcine circovirus type 2 (PCV2), which is the primary threat to the global swine industry leading to massive economic loss annually [5]. PCV2 is the cause of pig postweaning multisystemic wasting syndrome (PMWS), which is characterized by progressive loss of weight, respiratory difficulty, jaundice, and reproductive failure [6]. Microscopically, PMWS is associated with severe lymphocyte depletion and multiple organ inflammation [7,8]. Thus, the study of PCV2 is of great importance. PCV2 has a genome size of

Cell Culture and PCV2 Infection
PK-15 cells free of PCV2 were purchased from the American Type Culture Collection (ATCC, Manassas, VA, USA), and they were cultured in Dulbecco's modified Eagle medium (HyClone, GE, Beijing, China) supplemented with 10% fetal bovine serum (HyClone, GE, Beijing, China). The cell line was used because it is a widely accepted in vitro infection model, especially in the field of PCV2 virology, from which underlying molecular and cellular mechanisms of PCV2 infection have been learned [26][27][28][29][30][31]. The PCV2 strain 17HNA1 (Accession no.: MT423827) was used in our study. The infection procedures strictly followed a previously published protocol [32] with minor modifications: multiplicity of infection (MOI) was set at 0.1, and the PK-15 cells were harvested at 72 h postinfection. Viral and mock infections were repeated three times. Therefore, the mock infections were named as con-1, con-2, and con-3, while the PCV2 infection repeats were designated as PCV2-1, PCV2-2, and PCV2-3.

RNA Preparation, Quantitative Real-Time PCR (qRT-PCR), and RNA-Seq
Seventy-two hours postinfection, the PCV2-infected and mock-infected cells were washed with ice-cold PBS (Sigma Aldrich, Shanghai, China) three times before being subjected to total RNA extraction using the Mini Best Universal RNA Extraction kit (Takara, Dalian, China). Next, the extracted RNA was tested using the Agilent 2200 Analyzer for quality control (QC), and those samples with an RNA Integrity Number (RIN) greater than 8 passed the QC process. Then, the RNA was subjected to RNA-seq using the Illumina HiSeq X Ten platform (San Diego, CA, USA) according to the pipeline of Novel Bioinformatics Inc. (Shanghai, China). The raw sequencing data were filtered using the fastp [33] to obtain clean reads, which were subsequently mapped to the porcine genome using the Hisat2 [34]. Then, the raw counts were obtained using the FeatureCounts [35]. Next, differentially expressed transcripts were calculated based on the following parameters: |log 2 fold change| > 0.2630 (fold change greater than 1.2 or lesser than 5/6) and adjusted p-value (false discovery rate, FDR) < 0.05 using the DESeq2 tool [36]. The identification of lncRNAs used a strategy described in a previous report [25]. The sequenced reads were submitted to the Gene Expression Omnibus with the accession number: GSE149711.
qRT-PCR was then carried out to confirm the expression levels of the selected transcripts, and it was performed according to a previously described protocol [37]. The transcripts tested were ORF2, NGFR, PIGR, SERPINF1, LOC106505099, LOC106509609, TRIM29, and LOC100512687, the latter two of which served as the negative controls. The ACTB was used as an internal control. All primers used are listed in Table S1. qRT-PCR data were analyzed using the 2 −∆∆Ct method [38].

LncRNA Target Prediction
LncRNA cis-regulated target genes were identified by searching the protein-coding genes 50 kb up or downstream of each lncRNA. Then, the expression level of each of the identified lncRNAs was compared with all the sequenced reads. Transcripts with Pearson's correlation coefficient >0.99 and adjusted p-value < 0.05 were considered potential trans-targets of lncRNAs.

Enrichment Analysis
The differentially expressed genes (DEGs) identified by RNA-seq and the lncRNA targets were submitted to the Database for Annotation, Visualization and Integrated Discovery (DAVID) [39,40] for the identification of significantly altered Gene Ontology (GO) [41] and Kyoto Encyclopedia of Genes and Genomes (KEGG) [42] terms using the default parameters (EASE score < 0.1).

Statistics
Student's t-tests were used to compare the difference in transcripts between PCV2-infected and mock-infected cells. Data are presented as the mean ± standard errors of the means (SEMs). Type I error was set at 0.05.

Ethical Statement
The experiment was approved by the Animal Welfare Committee of Zhejiang University (No. 11,834 issued on 26 February 2018).

Results and Discussion
Total RNA from the PK-15 cells was extracted and subjected to the RNA-seq. After the QC process and removal of ribosomal RNA sequences, an average of 84% of the sequence reads was mapped to the pig genome (Table S2), demonstrating that RNA-seq provided a reliable result. The positions of the reads were mostly located in the coding sequences or exonic regions ( Figure S1). A total of 38,841 transcripts were identified in the six sequence samples, of which 111 transcripts were differentially expressed; 87 were upregulated, and 24 were downregulated ( Figure 1A and Table S3). Then, a heatmap with clustering analysis using all the DEGs was constructed, in which PCV2-infected cells exhibited significantly different expression patterns compared to the mock-infected cells ( Figure 1B). Within this group of transcripts, six differentially expressed lncRNAs were also observed; 4 were upregulated, and 2 were downregulated ( Figure 1A and Table S4).
Next, qRT-PCR was employed to confirm the expression of some selected genes and lncRNAs. The Cap expression level was significantly higher in the PCV2-infected groups, demonstrating the successful infection of PK-15 cells (Figure 2). Other selected genes and lncRNAs, such as NGFR, PIGR, SERPINF1, LOC106509609, and LOC106505099, showed expression patterns that were similar to the data obtained by RNA-seq ( Figure 2 and Table S3). Accordingly, the negative controls (TRIM29 and LOC100512687) had comparable expression levels in the PCV2-and mock-infected cells. Next, qRT-PCR was employed to confirm the expression of some selected genes and lncRNAs. The Cap expression level was significantly higher in the PCV2-infected groups, demonstrating the successful infection of PK-15 cells ( Figure 2). Other selected genes and lncRNAs, such as NGFR, PIGR, SERPINF1, LOC106509609, and LOC106505099, showed expression patterns that were similar to the data obtained by RNA-seq ( Figure 2 and Table S3). Accordingly, the negative controls (TRIM29 and LOC100512687) had comparable expression levels in the PCV2-and mock-infected cells. Analysis of the underlying molecular mechanisms of PCV2 infection and replication is key to the prevention of the virus. Previously, several high-throughput approaches have been employed to study mRNA, miRNA, and protein changes after PCV2 infection. In these studies, the DEGs analyzed by microarray or RNA-seq were enriched in the immune response, cytokine production, and apoptosis upon viral infection using in vivo or in vitro methods [43][44][45][46][47][48][49]. Later, miRNA was also tested for its role in the disease. Analysis of the targets of miRNAs indicated that disease resistance, MAPK signaling, and transcriptional regulation were involved [47,50]. Direct quantifications of the protein  Next, qRT-PCR was employed to confirm the expression of some selected genes and lncRNAs. The Cap expression level was significantly higher in the PCV2-infected groups, demonstrating the successful infection of PK-15 cells ( Figure 2). Other selected genes and lncRNAs, such as NGFR, PIGR, SERPINF1, LOC106509609, and LOC106505099, showed expression patterns that were similar to the data obtained by RNA-seq ( Figure 2 and Table S3). Accordingly, the negative controls (TRIM29 and LOC100512687) had comparable expression levels in the PCV2-and mock-infected cells. Analysis of the underlying molecular mechanisms of PCV2 infection and replication is key to the prevention of the virus. Previously, several high-throughput approaches have been employed to study mRNA, miRNA, and protein changes after PCV2 infection. In these studies, the DEGs analyzed by microarray or RNA-seq were enriched in the immune response, cytokine production, and apoptosis upon viral infection using in vivo or in vitro methods [43][44][45][46][47][48][49]. Later, miRNA was also tested for its role in the disease. Analysis of the targets of miRNAs indicated that disease resistance, MAPK signaling, and transcriptional regulation were involved [47,50]. Direct quantifications of the protein Analysis of the underlying molecular mechanisms of PCV2 infection and replication is key to the prevention of the virus. Previously, several high-throughput approaches have been employed to study mRNA, miRNA, and protein changes after PCV2 infection. In these studies, the DEGs analyzed by microarray or RNA-seq were enriched in the immune response, cytokine production, and apoptosis upon viral infection using in vivo or in vitro methods [43][44][45][46][47][48][49]. Later, miRNA was also tested for its role in the disease. Analysis of the targets of miRNAs indicated that disease resistance, MAPK signaling, and transcriptional regulation were involved [47,50]. Direct quantifications of the protein amounts after PCV2 infection revealed that there were altered levels of cytoskeletal proteins, stress response proteins, signal transduction proteins, and proteins involved in the ubiquitin-proteasome processes [32,51]. However, only one study had previously been performed to evaluate the effect of lncRNAs on PCV2 infection [25]. Fang et al. used PCV2 to infect the porcine intestinal epithelial cell line IPEC-J2 to mimic porcine enteritis. In their study, some interesting conclusions were drawn; they found that transcriptional regulation, nucleic acid binding, MAPK signaling, and cytokine-related pathways were disturbed upon PCV2 infection. However, they used the p-value instead of the adjusted p-value to identify the altered lncRNAs and mRNAs, which would substantially increase the type I error, leading to false positive findings. For example, in their lncRNA list, no lncRNA identified has an adjusted p-value of less than 0.05 [25]. Therefore, in the current study, we aimed to study whether some differentially expressed lncRNAs could be found using more stringent screening parameters.
Of the 111 differentially expressed transcripts, 105 were assigned to the protein-coding genes. Thus, these DEGs were analyzed using the DAVID tool for the identification of enriched biological processes (BP), cellular components (CC), and molecular functions (MF). As shown in Table 1, the DEGs were mainly involved in cell division, apoptosis, and the cytoskeleton and protein phosphorylation processes. Cell division and apoptosis have been well studied previously [52,53]. The involvement of cytoskeletal processes was also noticed in proteomic studies, indicating that the virus might affect cell division through the mitotic spindle [32,51]. The protein phosphorylation process is associated with numerous critical biological processes, such as the MAPK signaling pathway, which has been shown to be essential in viral infection [30,31]. KEGG analysis indicated that Rap1, Ras signaling pathways, and arginine biosynthesis were significantly enriched. Currently, no research has been conducted regarding the roles of arginine biosynthesis and Rap1 signaling in PCV2 infection. However, Ras is the upstream effector of MAPK signaling, which promotes infected cell proliferation [30,31]. In our enrichment analysis, no immune-related terms were found. This disparity may be because we used an MOI of 0.1, resulting in minimally infected PK-15 cells [26].
With our experimental conditions in combination with the stringent screening parameters, only six lncRNAs were found to be differentially expressed (adjusted p-value < 0.05, Table S4). Thus, these lncRNAs may be sensitive to PCV2 infection. Next, to functionally annotate these lncRNAs, the lncRNA targets were predicted. Currently, there are three main prediction approaches. The first one is cis-regulated target prediction, which is based on the fact that lncRNA can regulate its nearby gene transcription [54]. Thus, the genes adjacent to the lncRNA within an appropriate range are all potentially regulated. The second approach predicts the targets based on the coexpressing pattern. The underlying mechanism is that the interacting lncRNA and mRNA might have a similar expression pattern. The third method is directly computing the interaction potential of lncRNA and its targets using parameters such as free energy change, secondary structural information [55][56][57]. The former two methods can produce fast prediction results; however, the prediction is indirectly based on correlation, which needs to be validated experimentally. The computational approach can calculate the direct interaction of lncRNA and its targets. However, it also needs intensive computational facilities. After weighing the pros and cons, cis-regulated target prediction and coexpression methods were employed in our current study, and the results were then shown in Table 2. Analysis of the neighboring genes showed that 18 genes were potential targets; however, no cis-target was found for LOC102158335. One interesting finding is that lncRNA LOC106505099 was located within the HOXB (Homeobox B) gene family, which is critical for the development of embryos [58][59][60]. Coexpression analysis using stringent parameters predicted that only three genes, DNAAF3, PPP3CB, and CKMT1A, were correlated with LOC102168077, LOC100525935, and LOC102161888, respectively. Thus, the three genes might be regulated by their respective lncRNAs in a trans manner. Next, to predict the potential biological roles of these lncRNAs in the infected PK-15 cells, GO analysis was performed using the potential target genes of the lncRNAs. The enrichment analysis showed that these target genes were mainly located in the nucleus (CC) and were involved in two MF processes: DNA binding and transcription factor regulation (Table 3). GO_BP analysis revealed the involvement of DNA transcription and several developmental processes, e.g., embryonic skeletal development, facial nerve structural organization, and rhombomere 4 development, suggesting the role of PCV2 in organ development. Previous studies on pig embryos showed that PCV2-exposed embryos had a higher mortality rate, especially when the zona pellucida was compromised, resulting in a low pregnancy rate of sows [61][62][63]. These results indicated that PCV2 might affect the reproductive systems through lncRNA-mediated modulation of HOXB genes. Virus infection mediated dysregulation of HOX genes has been reported before. In one study, the hepatitis B virus can induce HOXA10 gene expression in human cells [64]. Another report showed that human papillomavirus gene expression could be modulated by HOXD9 in cervical cancer [65]. The Epstein-Barr virus activates NKL homeobox genes in B-cell lymphoma [66]. Thus, it is reasonable to link PCV2 infection with HOXB gene regulation. Subsequently, we searched the human genome for the presence of possible lncRNA orthologs. Based on syntenic information, LOC106505099 and LOC102161888 have potential human orthologs HOXB-AS1 and SMIM15-AS1, respectively. LOC106505099 and HOXB-AS1 are both located between HOXB2 and HOXB3 and are in the antisense direction, while LOC102161888 and SMIM15-AS1 are adjacent to the SMIM15 gene in a tail-to-tail manner. However, only HOXB-AS1 and LOC106505099 showed some extent of sequence conservation according to a Basic Local Alignment Search Tool (BLAST) search of the pig transcriptome using the HOXB-AS1 sequence (ENST00000435312.5). Functionally, some lncRNAs residing in the HOX gene cluster have been associated with developmental regulation through the modulation of neighboring HOX genes [67,68]. Moreover, two functional studies regarding HOXB-AS1 in two cancer cell lines indicated that the lncRNA might promote proliferation and migration of cancer cells by either acting as a miRNA sponge or by participating in the Wnt/beta-catenin signaling pathways [69,70]. Another research project on mouse ortholog Hoxb3os showed that the lncRNA was downregulated in cystic kidneys, and regulated EIF2, mTOR, and EIF4 signaling pathways [71]. However, the direct link between HOXB-AS1/LOC106505099 and development has yet to be established, and further functional validation is required to unveil the biological roles of HOXB-AS1 or LOC106505099 in porcine embryo development and PCV2 infection.

Conclusions
Using the RNA-seq technique, we identified differentially expressed lncRNAs upon PCV2 infection with stringent screening parameters. Through analysis of the targets of these lncRNAs, we found that LOC106505099 might be involved in the stages of embryonic development through the regulation of HOXB genes. These findings suggest a novel association between PCV2 infection and embryonic development, which would shed new light on the pathological mechanisms of PCV2 infection.