To investigate the dynamics of host–pathogen interactions based on poly(A) pattern and gene expression profile during PRRSV infection, we used the SAPAS method combined with IVT to generate nine IVT-SAPAS libraries (C0, C6, T6, C12, T12, C24, T24, C36, and T36, where C and T denote control and treatment groups, respectively, and numbers indicate hpi).
A total of 332 million raw reads were generated by Illumina sequencing. After mapping to the monkey genome and internal priming filtering, about 156 million reads were obtained for poly(A) sites analysis (Table 1). Only 39% of reads used University of California at Santa Cruz annotated poly(A) sites, indicating that IVT-SAPAS was effective in detecting novel poly(A) sites, especially in samples with low mRNA expression levels (Fig. 1A). Additionally, 22% and 12% of reads were located in the intergenic region and 1 kb downstream of the gene's 3' UTR. About 2 × 105 variable poly(A) sites—only 13% (Fig. 1B) of which are known—accounted for 58% of reads (Fig. 1A).
Table 1. Summary of IVT-SAPAS data from Illumina sequencing.
Figure 1. Sequencing reads and poly(A) sites in Marc-145 cells. A Genomic location of sequencing reads. B Distribution of poly(A) sites in the genome. C Genes with different numbers of tandem poly(A) sites. D Histogram of the median distance between stop codon and poly(A) sites in genes with single poly(A) sites, and distance between stop codon and closest or longest poly(A) sites in genes with multiple poly(A) sites.
We arranged the annotated 3' UTR APA sites in tandem based on the stop codon and combined these with poly(A) site profiling as previously described (Fu et al. 2011), and found that 6536 genes had tandem 3' UTRs, while 34.33% had at least two tandem 3' UTR poly(A) sites (Fig. 1C). Poly(A) sites switching may regulate a variety of biological functions. The average distribution of distances between poly(A) site and stop codon was 426 bp (median: 270 bp). Median distances between the stop codon and proximal and distal poly(A) sites were 237 and 720 bp, respectively (Fig. 1D). These results reveal the detailed landscape of poly(A) site usage in Marc-145 cells and suggest an important regulatory role for APA in the monkey genome.
The poly(A) signal (PAS) is an important cis-acting element in 3' end processing that is located 10–30 nt upstream of the cleavage site. The six-nucleotide base sequence is recognized by cleavage poly(A) specificity factor (Elkon et al. 2013). To date, 12 types of PAS have been identified, with AATAAA and ATTAAA being the most common (Proudfoot 2011), which was confirmed by our finding that the usage frequency was close to 50%. Our data suggest that poly(A) is dynamic; moreover, 26.6% of the six PAS dimers were not observed (Fig. 2), suggesting that there are other regulatory mechanisms.
The IVT-SAPAS method can be used not only to analyze APA, but also to quantify the expression of a gene based on an analysis of the 3' UTR of its transcripts (Jia et al. 2017). Genes with tandem 3' UTR lengths that differed significantly between libraries were defined as APA switching genes. We identified 19, 17, 52, and 142 such genes at 6, 12, 24, and 36 hpi, respectively (Fig. 3A). In total, there were 185 APA switching genes (false discovery rate [FDR] < 0.05, Rcut ≤ 0.05) after PRRSV infection, of which 141 used shorter 3' UTRs, 32 used longer 3' UTRs, and 12 were dynamically regulated at different time points (Fig. 3B). Genes exhibiting a greater than 1.5-fold difference in expression and had an FDR < 0.05 were considered as differentially expressed. We identified 15, 14, 71 and 361 DEGs at 6, 12, 24, and 36 hpi, respectively (Fig. 3C). In total, there were 393 genes that were differentially regulated between mock- and PRRSV-infected samples, of which 234 were upregulated, 152 were downregulated, and 7 were dynamically regulated at different time points (Fig. 3B).
Figure 3. Analysis of tandem 3' UTR switching and regulation of gene expression in PRRSV-infected Marc-145 cells. A Summary of APA genes between samples at different time points. B Summary of APA genes and DEGs. C Summary of DEGs between samples at different time points. D Wayne chart of APA genes and DEGs. E, F Functional classification of APA genes (E) and DEGs (F).
To clarify the functional relationship between APA sites and differential gene expression in antiviral immunity, we compared APA genes and DEGs and identified 35 genes involved in the simultaneous regulation of the two mechanisms, indicating that the two mechanisms are independent (Fig. 3D). Furthermore, a classification of gene function revealed that there was little difference in gene function between the two gene sets, which covered cellular and metabolic processes, biological regulation, and response to stimulus, among other categories (Fig. 3E, 3F). These results indicate that these two mechanisms are actively involved in the response to viral infection.
To clarify the changes in gene expression over the course of viral infection, we generated a growth curve. Viral titer increased slightly within 12 hpi and then rapidly thereafter (black curves, left Y axis; red curve, right Y axis) unlike mock-infected cells (blue curve, right Y axis), which is consistent with the sequencing data (Fig. 4A). We then analyzed the number of APA genes (Fig. 4B) and DEGs (Fig. 4C) at two sequential time points and found that their total numbers increased during the infection cycle, especially at 12–24 and 24–36 hpi, respectively, when viral replication was in the logarithmic phase, suggesting a robust antiviral response. It is possible that this was due to the induction of host interferon (IFN)-stimulated gene expression after viral infection.
Figure 4. Dynamic regulation in PRRSV-infected Marc-145 cells. A PRRSV growth curve and number of reads. B, C Regulation of APA genes (B) and DEGs (C) at two sequential time points. D, E KEGG analysis of APA genes (D) and DEGs (E).
The KEGG analysis (P ≤ 0.05) showed that APA genes (Fig. 4D) and DEGs (Fig. 4E) were upregulated during the immune response to PRRSV infection. The DEGs were also associated with the adaptive immune response, including B and T cell receptor signaling pathways. Additionally, DEGs were more directly involved in the antiviral innate immune response than APA genes, including apoptosis and activation of IFN-stimulated pathways. These results suggest that the two mechanisms act cooperatively in the host response to PRRSV infection.
To investigate the relationship between APA and DEGs, we identified genes common to the two gene sets at 36 hpi (Table 2). Compared to the negative control group, infected cells preferentially utilized shorter 3' UTRs. Accordingly, most genes involved in APA shortening were upregulated. Proteins associated with innate immunity including JUN and NFKB1 were at the center of the protein interaction map and did not show interactions with other proteins (Fig. 5A). Six mixed APA genes (black columns in Fig. 5B) and DEGs (white columns in Fig. 5B) were validated by qRT-PCR and the results were consistent with the sequencing data.
Table 2. Summary of mixed APA genes and DEGs at 36 hpi.