Comparative transcriptomics analysis on Senecavirus A-infected and non-infected cells

Senecavirus A (SVA) is an emerging virus that causes the vesicular disease in pigs, clinically indistinguishable from other high consequence vesicular diseases. This virus belongs to the genus Senecavirus in the family Picornaviridae. Its genome is a positive-sense, single-stranded RNA, approximately 7,300 nt in length, with a 3′ poly(A) tail but without 5′-end capped structure. SVA can efficiently propagate in different cells, including some non-pig-derived cell lines. A wild-type SVA was previously rescued from its cDNA clone using reverse genetics in our laboratory. In the present study, the BSR-T7/5 cell line was inoculated with the passage-5 SVA. At 12 h post-inoculation, SVA-infected and non-infected cells were independently collected for the analysis on comparative transcriptomics. The results totally showed 628 differentially expressed genes, including 565 upregulated and 63 downregulated ones, suggesting that SVA infection significantly stimulated the transcription initiation in cells. GO and KEGG enrichment analyses demonstrated that SVA exerted multiple effects on immunity-related pathways in cells. Furthermore, the RNA sequencing data were subjected to other in-depth analyses, such as the single-nucleotide polymorphism, transcription factors, and protein–protein interactions. The present study, along with our previous proteomics and metabolomics researches, provides a multi-omics insight into the interaction between SVA and its hosts.

1 Introduction

Senecavirus A (SVA), as an emerging virus, has been demonstrated to be a causative agent for vesicular disease in swine (1–5). SVA-infected pigs develop vesicular lesions mainly on the snout, dewclaw or (and) coronary band. Other signs include lameness, anorexia, lethargy, cutaneous hyperemia, and fever (6, 7). The SVA-induced signs are clinically indistinguishable from those of other vesicular diseases in pigs (8). The outbreak of SVA infection has been recently reported in several countries, including Canada, the United States, Brazil, China, Thailand, Vietnam and Chile. The transmission risk of SVA has attracted a great deal of attention from the pig industry around the world.

SVA is the only member of the genus Senecavirus in the family Picornaviridae (9). The virion is a typical icosahedral particle without envelope. It harbors a positive-sense, single-stranded RNA genome, approximately 7,300 nt in length, composed of 5′ untranslated region (UTR), long encoding region and 3’ UTR. Like those of other picornaviruses, the 5′ terminus of SVA genome does not contain a cap structure. In contrast, a short peptide (VPg) is covalently linked to the 5′ terminus, and plays an essential role in synthesizing the SVA genome. The 5’ UTR bears a type-IV internal ribosome entry site (10), structurally and functionally similar to those of pestiviruses (11), allowing for the initiation of polyprotein translation in a cap-independent manner. The 3’ UTR is approximately 70 nt in length, followed by a variable-length poly(A) tail (12). The encoding region of SVA polyprotein follows the standard “L–VP4–VP2–VP3–VP1–2A–2B–2C–3A–3B–3C–3D” layout. After SVA infection, the viral polyprotein will be translated in cytoplasm, and then gradually cleaved into 12 proteins, namely, L, VP4, VP2, VP3, VP1, 2A, 2B, 2C, 3A, 3B, 3C and 3D (13). The VP4, VP2, VP3 and VP1 as structural proteins participate in viral morphogenesis. The others are nonstructural proteins, albeit uninvolved in the package of virion, playing a crucial role in viral replication (14–16).

RNA sequencing (RNA-seq) is a technique that uses next-generation sequencing to reveal the presence and quantity of RNA molecules in a biological sample, providing a snapshot of gene expression in the sample, also known as transcriptome. A transcriptome is the full range of mRNA molecules expressed by an organism. The RNA-seq technique contributes to identifying a transcriptome in a given population, even in a single cell (17). Comparative transcriptomics facilitates the elucidation of differentiation between two groups (populations, species and so on) in their alternative gene spliced transcripts, post-transcriptional modifications, gene fusion, single-nucleotide polymorphism (SNP) and changes in gene expression over time (18). Large DNA viruses, such as human cytomegalovirus and African swine fever virus, contained very long genomes. Each of these viruses itself has a complicated transcriptome (19, 20) in virus-infected cells. In contrast, some small RNA viruses, such as picornavirus, only have a simple “transcriptome,” i.e., one single RNA genome. Therefore, it is meaningless to uncover a picornaviral “transcriptome” only based on a given picornavirus itself.

SVA can trigger a variety of metabolic and biochemical changes in cells through virus-specific or -nonspecific mechanisms (21–23). For example, SVA 2C protein can target mitochondria and cause release of cytochrome C into cytoplasm, which activates caspase-9 and -3 in a series of signaling cascades to induce the onset of apoptosis (24). In addition, SVA infection is able to affect the level of transcription in hosts. For example, SVA 2C protein can block the transcription of interferon-stimulated gene 56 and interferon-β to weaken the innate immunity in hosts (21). We have demonstrated that SVA infection can lead to significant changes in cellular proteome and metabolome, even at an early stage of infection (25, 26). Virus-caused differentiation of cellular proteome is closely related to the change in cellular transcriptome. Therefore, a comparative transcriptomics analysis was conducted here to uncover a profile of SVA-induced changes in cellular transcriptome at the early stage of infection.

2 Materials and methods

2.1 Cell line and virus

The BSR-T7/5 cell line, derived from the baby hamster kidney cell, was kindly provided by the China Animal Health and Epidemiology Center. This cell line was cultured at 37°C with 5% CO₂ in Dulbecco’s modified Eagle’s medium (DMEM), supplemented with 4% fetal bovine serum (VivaCell, Shanghai, China), penicillin (100 U/mL), streptomycin (100 μg/mL) and amphotericin B (0.25 μg/mL). The wild-type SVA was rescued previously from a full-length cDNA clone (27), genetically derived from an SVA isolate, CH-LX-01-2016 (28).

2.2 Sample preparation

BSR-T7/5 cells were seeded into six T25 flasks for culture at 37°C. When the cells were 90% confluent, three flasks were randomly selected for incubation with the passage-5 SVA at an MOI (multiplicity of infection) of 2.5. The other flasks, as non-infected controls, were not treated. There were three SVA-infected samples (S1, S2 and S3) and three non-infected controls (C1, C2 and C3). Supernatants were separately removed from the six flasks at 12 h post-inoculation (hpi). Cell monolayers were gently washed with PBS three times, followed by the extraction of total RNAs using the TRIzol reagent (Thermo Fisher, Waltham, MA, United States), as per the manufacturer’s instructions. The concentration, quality and integrity of total RNAs were determined using the NanoDrop spectrophotometer (Thermo Fisher, Waltham, MA, United States). Three μg of RNA was used as input material to prepare RNA sample for each group.

2.3 RNA-seq analysis

The preparation of sequencing libraries was carried out as described previously with modifications (29). The mRNAs were purified from total RNAs using poly-T oligo-attached magnetic beads, further fragmented, and then used as templates to produce cDNAs. The first strand cDNA was synthesized using a system with random hexamer primers and the reverse transcriptase. The second strand cDNA was synthesized via the first strand with dNTP, buffer solution, DNA polymerase I and RNase H. Remaining overhangs were converted into blunt ends through exonuclease/polymerase activities. After adenylation of the 3′ ends of DNA fragments, Illumina paired-end adapter oligonucleotides were ligated to prepare for hybridization. The cDNA fragments of 400 to 500 bp were preferentially size-selected using the AMPure XP system (Beckman Coulter, Beverly, United States). DNA fragments with ligated adaptor molecules on both ends were selectively enriched using Illumina PCR Primer Cocktail in a 15-cycle PCR reaction. Products were purified using the AMPure XP system, and then quantified by the Agilent high sensitivity DNA assay on the Agilent 2,100 Bioanalyzer (Agilent Technologies, CA, United States). The sequencing libraries were subjected to sequencing on the NovaSeq 6,000 platform (Illumina, CA, United States) for obtaining image files.

2.4 Quality control and reads mapping

The image files were transformed by the software of sequencing platform. The original data was generated in a FASTQ format (raw data). Sequencing data contained a number of connectors and low-quality reads. The Cutadapt (v1.15) software was used to filter the sequencing data (30), subsequently obtaining high-quality sequences (clean data) for further analysis. Two reference genomes were those of the golden hamster (Genbank No.: PRJNA77669) and the SVA CH-LX-01-2016 (Genbank No.: KX751945). The filtered reads were separately mapped to both reference genomes using the HISAT2 (v2.0.5) program (31).

2.5 Analysis of differential expression

The analysis of differentially expressed genes (DEGs) was performed as described previously with modifications (32). The HTSeq (v0.9.1) was used to compare the Read Count values on each gene as the original gene expression (33). Gene expression was standardized through the FPKM (Fragments Per Kilobase of transcript per Million mapped reads). DEGs were determined by the DESeq (v1.30.0) with screening parameters as follows: the fold change (FC) > 2 (or <0.5) and the significant p value <0.05 (34). The bi-directional clustering analysis of all DEGs was performed by the Pheatmap (v1.0.8) package. The heatmap was drawn according to the expression level of the same gene in different groups and the expression patterns of different genes in the same group, with the Euclidean method for calculating the distance, and the complete linkage method for clustering.

2.6 Analyses of GO and KEGG enrichments

All genes were mapped to terms in the database of gene ontology (GO). Differentially enriched genes were calculated for each term. The topGO package was designed to perform the GO enrichment analysis on the DEGs. The p value was calculated by the hypergeometric distribution method. The p value <0.05 was determined as the standard of significant enrichment. The GO terms were found with significantly differentially enriched genes, all of which were further classified to determine the main biological functions. The ClusterProfiler (v3.4.4) software was used to carry out the enrichment analysis of DEGs on the KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways. The p value <0.05 was determined as the standard of significant enrichment (35).

2.7 Other analyses on RNA-seq data

2.7.1 Analysis of new transcripts

On the basis of the existing reference genome, the software StringTie (http://ccb.jhu.edu/software/stringtie/) was used to assemble the mapped reads (36). The assembling results were compared with the known transcripts for obtaining unannotated transcripts.

2.7.2 Analysis of alternative splicing events

The rMATS (v3.2.5) software was used to uncover alternative splicing events (37). The main types of alternative splicing events included skipped exon (SE), retained intron (RI), alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), and mutually exclusive exons (MXE).

2.7.3 Analysis of SNP sites

The Varscan program was used to obtain SNP sites (38). The filtering criteria were: (i) SNP site base Q > 20, (ii) the number of reads covering the site >8, (iii) the number of reads supporting the mutation site >2, and (4) the p value of SNP locus <0.01.

2.7.4 Prediction of transcription factors

Transcription factors and their own families were predicted via the comparison with the Animal Transcription Factor Database (AnimalTFDB) (39), a comprehensive database including classification and annotation of genome-wide transcription factors, transcription co-factors and chromatin remodeling factors in numerous animal genomes.

2.7.5 Analysis of differential exon usage

The DEXSeq package was used to analyze the RNA-seq data for identifying the differential exon usage, as described previously (40).

2.7.6 Interaction analysis in protein network

The STRING database (https://string-db.org/) was used to unveil putative protein–protein interactions (PPI) (41), contributing to clarifying the relationship among genes of interest.

2.8 Validation of gene expression by RT-qPCR

Four representative genes, namely, SVA genome, Nfkbia, Phlda2 and Txnip, were selected for validating the profile of gene expression. The GAPDH (glyceraldehyde-3-phosphate dehydrogenase) gene was used as an internal reference control. Gene-specific primers were listed in Supplementary file 18 for RT-qPCR validation. The RT-qPCR analysis was performed with three technical repeats, using the AceQ qPCR SYBR Green Master Mix (Vazyme, Nanjing, China) based on the LightCycler 480® Real-time PCR System (Roche, Rotkreuz, Switzerland), as per the manufacturer’s instructions. The RT-qPCR results were analyzed through the 2^−ΔΔCt method for relatively quantifying the four genes of interest (42). The GraphPad Prism (v8.0) was used for statistical analysis by two-tailed Student′s t-test with Welch′s correction. Data were shown as means ± standard deviations of three independent experiments.

3 Results

3.1 Sequencing for de novo transcriptome assembly

The BSR-T7/5 cell monolayers showed no obvious cytopathic effect (CPE) at 12 hpi (Figure 1). Cells were separately collected from SVA-infected and non-infected groups to extract total RNAs for the construction of high-quality cDNA libraries. The primary library-related data were listed in Supplementary file 1. The cDNA libraries were subjected to sequencing to obtain image files, subsequently transformed into raw data for statistical classification, as shown in Table 1. Because the raw data contained a number of connectors and low-quality reads, the Cutadapt (v1.15) software was used to filter the raw data for obtaining high-quality clean sequences, as listed in Supplementary file 2. There were three key parameters, namely, base mass, base content, and average mass of reads. Their distributions were independently shown in Supplementary 3, 4 and 5.

Figure 1

Figure 1. SVA-infected and non-infected cell monolayers at 12 hpi. S1, S2, and S3: SVA-infected groups; C1, C2 and C3: non-infected controls.

Table 1

Table 1. Statistic data of RNA-seq for all six groups.

3.2 Transcriptomic mapping

The filtered reads were mapped to both reference genomes, those of the golden hamster and the SVA CH-LX-01-2016, using the HISAT2 (v2.0.5) program. The results of RNA-seq mapping were listed in Table 2. The global profile of reads was subjected to the further statistical analysis on the distribution of reads mapped to both genomes, in which genetic elements included coding sequence, intron, intergenic spacer, and UTR. The mapping results were listed in Table 3, and shown in Supplementary file 6. Supplementary file 7 exhibited the coverage distributions of reads mapped to genes. To sum up, a high-quality dataset of RNA-Seq was harvested here, meeting a standard for the further bioinformatics analysis.

Table 2

Table 2. Statistical data of RNA-seq mapping.

Table 3

Table 3. Distribution of read-mapped regions.

3.3 Profile of gene expression

A total of 20,374 genes were identified in all six groups (Supplementary file 8), but these genes contained more than 3,000 components with FPKM value = 0. FPKM was a simple method for normalizing the read count data, based on gene length and the total number of mapped reads. The FPKM-normalized expression level was divided into different intervals (Supplementary file 9) for the six groups, as shown in Supplementary file 10. The number of genes, either co-identified in different groups or recognized in a single group, was shown in Figure 2A. The density distribution of FPKM values was displayed in Figures 2B,C, as two different forms. To validate whether the sequencing depth of RNA-seq was sufficient for the analysis of gene expression, the saturation analysis was performed for all six groups, as shown in Supplementary file 11. The correlation analysis of gene expression, based on calculation of Pearson’s correlation coefficients, was carried out among the six groups (Figure 2D). The closer to 1.0 the correlation coefficient was, the higher the similarity of expression pattern was among the six groups. Principal component analysis made it possible to project a high-dimensional dataset onto two or three dimensions, as shown in Figure 2E. The closer the distance was, the higher the similarity was among groups.

Figure 2

Figure 2. Profile of gene expression in all groups. The numbers of genes that are either co-identified in different groups or recognized in a single group (A). Violin plot of FPKM distributions in all groups (B). Distributions of FPKM densities in all groups (C). The correlation analysis of gene expression via the calculation of Pearson’s correlation coefficients (D). Principal component analysis on all groups (E). PC1: principal component 1. PC2: principal component 2.

3.4 Analysis of differential expression

DEGs were determined by the DESeq (v1.30.0) with screening parameters, |log₂FC| > 1 and the significant p value <0.05. A total of 628 DEGs, including 565 upregulated (Supplementary file 12) and 63 downregulated (Supplementary file 13) components, were identified here. The basemean was described as the “mean of normalized counts of all samples.” The basemean values of DEGs, corresponding to the group C and S, were exhibited in Figures 3A,B, respectively. The distributions of sequence length, p value and log₂FC were exhibited in Figures 3C–H. A single asterisk in Figures 3G,H indicated the exclusion of positive or negative infinity (“Inf” in Supplementary files 12, 13), respectively. The distribution and degree of differential expression were graphically shown in a volcano plot and a heatmap, respectively.

Figure 3

Figure 3. Profile and analysis of DEGs. Distribution of basemean values in group C (A) and S (B). Distributions of sequence lengths of upregulated (C) and downregulated (D) DEGs. Distributions of p values of upregulated (E) and downregulated (F) DEGs. Distributions of log₂FC values of upregulated (G) and downregulated (H) DEGs. *Excluding the SVA-related data. Volcano plot of p value versus FC for all identified genes but excluding outliers (I). The threshold values are set as |log₂FC| > 1 and p value <0.05. Heatmap based on bi-directional clustering analysis of all DEGs (J). Clustering analysis on expression patterns of DEGs (K). Grey lines indicate expression patterns. Each blue line represents the average value in each cluster.

The volcano plot (Figure 3I), drawn by the GraphPad Prism software, revealed the p value versus the FC for all identified genes. The threshold values were set as |log₂FC| > 1 and p value <0.05. The upregulated, downregulated and stably-expressed genes were indicated by red, green and grey circles, respectively. The R language Pheatmap (1.0.8) software package was used for the bi-directional clustering analysis to draw the heatmap (Figure 3J), which provided a visual depiction for hierarchical clustering of all 628 DEGs from the six groups. The red and green labels represented the upregulated and downregulated DEGs, respectively. The intensity of color reflected the degree of differentiation in gene expression. All 628 DEGs, based on their differences in expression patterns, were classified into nine different clusters (Figure 3K), in which grey lines indicated expression patterns, and each blue line represented the average value in each cluster.

3.5 GO enrichment analysis

The topGO package was used for performing the GO enrichment analysis on DEGs. The GO terms were found with significantly differentially enriched genes. The numbers of category BP, CC and MF were 3,093, 357 and 644, respectively. The complete GO data were listed in Supplementary file 14 in detail. Figure 4A displayed the top-10 statistically significant GO terms for each GO category. The false discovery rate (FDR), ranging from 0 to 1, was associated with the degree of GO enrichment. The lower the FDR was, the more significant the enrichment degree was. The GO terms with the top-20 lowest FDRs were shown in a bubble plot (Figure 4B). Each GO category was organized further as a directed acyclic graph (Figures 4C–E), in which parental terms described more general functional categories than their next-generation terms. GO terms with the top-10 lowest FDRs were framed with rectangles, and the others were indicated by ellipses. The more statistically significant a GO term was, the darker its color was.

Figure 4

Figure 4. GO and KEGG enrichment analyses of DEGs. Top-10 statistically significant GO terms of three categories (A). Bubble plot of top-20 statistically significant GO terms (B). Directed acyclic graphs of top-10 statistically significant GO categories (C–E). GO terms with the top-10 lowest FDRs are framed with rectangles, and the others are indicated by ellipses. The more statistically significant a GO term is, the darker its color is. Top-30 statistically significant KEGG pathways, classified into four categories (F). Bubble plot of top-20 statistically significant KEGG pathways (G).

3.6 KEGG enrichment analysis

The analysis of KEGG pathway enrichment was performed to uncover DEG-related pathways. The result showed that DEGs were enriched in a total of 275 KEGG pathways (Supplementary file 15). Figure 4F displayed the top-30 statistically significant KEGG pathways (p value <0.05), classified into four categories, namely, environmental information processing, human diseases, metabolism, and organismal systems. According to the result of KEGG enrichment, the degree of enrichment was evaluated through the rich factor, FDR, and the number of DEGs enriched in a given pathway. The higher the rich factor was, the more significant the enrichment degree was. The lower the FDR was, the more significant the enrichment degree was. The KEGG pathways with the top-20 lowest FDRs were shown in a bubble plot (Figure 4G).

3.7 Other analyses on RNA-seq data

The StringTie was used to assemble the mapped reads. The assembling results were compared with the known transcripts to obtain unannotated transcripts. The transcripts from class codes j, i, u, and x, regarded as new transcripts, were functionally annotated and listed in Supplementary file 16. The proportion of each class code was shown in Figure 5A. Five types of alternative splicing events were analyzed by the rMATS (v3.2.5) software. The SE and RI exhibited the most and the least alternative splicing events, respectively (Figure 5B). The SNP sites were analyzed by the Varscan program. The numbers of heterozygous and homozygous variants were shown in Figure 5C. Transcription factors and their own families were predicted through the comparison with those in the AnimalTFDB. Figure 5D showed the number of transcription factors in each family. Out of the identified families, 20 were demonstrated to contain upregulated, downregulated or both components (Figure 5E). The DEXSeq package was used to analyze the RNA-seq data for identifying the differential exon usage, as shown in Supplementary file 17. Figure 5F revealed a representative gene with differential exon usage. DEGs were comprehensively analyzed in the STRING database to unveil potential PPIs (Score > 0.95) for constructing a PPI network, which including 26 nodes and 17 edges (Figure 5G). Red and green nodes indicated upregulated and downregulated genes, respectively.

Figure 5

Figure 5. In-depth analysis of RNA-seq data. Pie chart of new transcripts (A). All new transcripts are classified into four categories, x, u, j and i. Analysis of alternative splicing events (B). X and Y axes indicate the types of alternative splicing events, and the number of new transcripts, respectively. S.D.: significant differentiation. The numbers of SNP sites in all groups (C). Profile of transcription factor families (D). The Y axis indicates the number of transcription factors in each family. Transcription factor families with significantly differential transcription factors (E). A representative gene with differential exon usage (F). The differential exon usage is marked with a purple rectangle. PPI network with 26 nodes and 17 edges (G). The score of PPI is set to be more than 0.95. Red and green nodes indicate upregulated and downregulated genes, respectively.

3.8 Validation of gene expression by RT-qPCR

Three upregulated and one downregulated DEGs were selected for validating the profile of gene expression through RT-qPCR. The three upregulated genes included SVA genome, Nfkbia and Phlda2 (Supplementary file 12); the downregulated gene was Txnip (Supplementary file 13). The RT-qPCR detection demonstrated that the expression trend of DEGs was consistent with the result of RNA-seq analysis (Figure 6). Due to the group C without SVA inoculation, the expression differentiation of SVA genome was extremely significant between both groups (Figure 6, Left upper). There was no need for the statistical analysis on it.

Figure 6

Figure 6. RT-qPCR validation of gene expression. The 2^-ΔΔCt method is used for relative quantification. The GAPDH gene is an internal control for normalization. Data are shown as means ± standard deviations of three independent experiments. Statistical significance is determined by two-tailed Student′s t-test with Welch′s correction. *p<0.05.

4 Discussion

The Picornaviridae is a well-characterized family within the plus-strand RNA viruses. SVA is a typical picornavirus. Its genome is only a positive-sense, single-stranded mRNA, harboring a 3′ poly(A) tail but no 5′ capped structure. In other words, an SVA virion has a single mRNA, which however is not the viral transcript. Picornaviruses, albeit structurally simple, possibly have significant effects on physiological functions in their hosts. After entrance of virion into a cell, a picornaviral genome will be released into cytosol. This genome either relies on the host translation machinery to initiate the translation of polyprotein, or serves as a template to synthesize an antigenome, which will be used as a template for synthesizing another genome. The nascent genome can be used as a template for the next round of translation or replication, and alternatively is packaged into a virion (43). Therefore, although SVA as such has no concept of viral transcriptome, its infection may exert a significant impact on the cellular transcriptome. This prompted us to conduct the present study for uncovering the transcriptomic change in SVA-infected cells.

A replication-competent SVA was previously rescued from its cDNA clone in our laboratory (27). The passage-5 SVA was used here as a model virus. Despite SVA inoculation with MOI of 2.5, three cell monolayers showed no visible CPE at 12 hpi (Figure 1). Because we demonstrated previously that SVA infection led to significant cellular changes both in proteomic and in metabolomic profiles at 12 hpi, it could be postulated that the cellular transcriptome would be also affected at 12 hpi to some extent. The RNA-seq recognized totally 20,374 genes in the six groups, but containing more than 3,000 genes with FPKM value = 0. The correlation of gene expression is an important indicator to demonstrate the reliability of experiment, and the reasonability of samples. A certain correlation coefficient, if between 0.8 and 1.0, would indicate the extremely strong correlation between two groups. The current correlation analysis displayed the extremely strong intra-group correlation, but the weak inter-group correlation (Figure 2D), implying the RNA-seq data that were reliable.

RNA-seq data were subjected to the further analysis on the differentiation expression. The result totally recognized 565 upregulated and 63 downregulated DEGs here. Out of these DEGs, five representative genes (three upregulated and two downregulated genes) were selected out for RT-qPCR analysis to validate preliminarily the profile of DEGs. Except the downregulated DEG, Tcta gene (data not shown), the other four showed their expression trends consistent with the result of RNA-seq analysis (Figure 6). In our previous study on comparative proteomics between SVA-infected and non-infected cells, we identified totally 305 upregulated and 56 downregulated DEPs (differentially expressed proteins) (25). Regardless of the present or the previous study, the number of upregulated components was much higher than that of downregulated ones. Such a result was consistent with our postulation that DEGs shared a similar regulation trend with DEPs between SVA-infected and non-infected groups. Out of the identified DEGs in group S, the SVA genome was most statistically significant in the expression level (Supplementary file 12). Due to the group C without SVA inoculation, both GO and KEGG enrichment analyses excluded the data of SVA genome, and others with positive or negative infinity (Supplementary files 12, 13). The KEGG enrichment analysis showed that several DEGs were significantly enriched in many immunity-related pathways, such as TNF signaling pathway, IL-17 signaling pathway, Toll-like receptor signaling pathway, and B cell receptor signaling pathway (Figure 4F). The GO enrichment analysis also revealed a few statistically significant terms associated with immune responses, e.g., the response to interleukin-1 and the cellular response to tumor necrosis factor (Figure 4A). Like the conclusion drawn in a previous report (44), the current results also suggest that SVA infection may be able to induce significantly immune responses, especially the innate immune response, in hosts at an early stage of infection.

Further, RNA-seq data were subjected to the in-depth analyses, concerning SNPs, transcription factors, PPI and so on. The analysis of SNP events indicated that there was no significant differentiation in the number of SNP events between group C and S (Figure 5C), implying that SVA infection had no ability of inducing the occurrence of SNP events in the host genome. It is worth noting that out of the 20 statistically significant families of transcription factors, most of them only contain upregulated components (Figure 5E). This result implies that SVA infection can notably stimulate multiple transcriptional pathways, resulting in upregulated DEGs far more than downregulated DEGs. The STRING database was used here to unravel putative PPIs in SVA-infected cells. The result revealed no formation of complicated interaction network among 26 putative DEPs (Figure 5G). Although the 26 DEPs include no SVA-related protein, the possibility that SVA-related proteins interact with cellular proteins cannot be ruled out, because the information of SVA proteins has not been deposited in the STRING database.

SVA emerged in many countries and regions over the past 20 years. It has been still considered as an emerging virus. Natural selection has been a primary evolutionary force affecting SVA codon usage bias (45). Multi-omics analysis provides an integrated approach to facilitate in-depth studies on the virology, especially on the interaction of viruses with their hosts. Based on our previous researches on proteomics and metabolomics, it was demonstrated that SVA infection could lead to significant changes in cellular intrinsic components even at an early stage of infection (25, 26). In order to comparatively analyze transcritpomic profiles between SVA-infected and non-infected cells, we conducted the present study. To sum up, the current results revealed that most of the DEGs were upregulated genes, indicating that SVA infection positively stimulated the transcription initiation in cells. GO and KEGG enrichment analyses demonstrated that SVA could markedly affect immunity-related pathways in cells, whereas the mechanism remained to be elucidated.

Data availability statement

The data presented in the study are deposited in the NCBI repository, accession number PRJNA1100277.

Author contributions

YL: Writing – original draft, Methodology, Formal analysis. HC: Writing – original draft, Methodology, Investigation, Data curation. YJ: Writing – review & editing, Formal analysis. ZL: Writing – review & editing, Investigation. JW: Writing – review & editing, Software. FL: Writing – review & editing, Writing – original draft, Supervision, Conceptualization.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported by the National Natural Science Foundation of China (Grant No. 32273000), the Qingdao Demonstration Project for People-benefit from Science and Techniques (Grant No. 24–2-8-xdny-4-nsh), and the Postgraduate Innovation Program of Qingdao Agricultural University (QNYCX24043).

Acknowledgments

We gratefully thank the Shanghai Bioprofile Biotechnology Co., Ltd. (Shanghai, China) for providing technical assistance in RNA-seq. We also thank Xianggan Cui for his help in technical support.

Conflict of interest

YJ was employed by the Qingdao Zhongren-OLand Bioengineering Co., Ltd.

The remaining 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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

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

References

1. Pasma, T, Davidson, S, and Shaw, SL. Idiopathic vesicular disease in swine in Manitoba. Can Vet J. (2008) 49:84–5.

PubMed Abstract | Google Scholar

2. Singh, K, Corner, S, Clark, S, Scherba, G, and Fredrickson, R. Seneca Valley virus and vesicular lesions in a pig with idiopathic vesicular disease. J. Vet. Sci. Technol. (2012) 3:1–3. doi: 10.4172/2157-7579.1000123

Crossref Full Text | Google Scholar

3. Leme, RA, Zotti, E, Alcantara, BK, Oliveira, MV, Freitas, LA, Alfieri, AF, et al. Senecavirus a: an emerging vesicular infection in Brazilian pig herds. Transbound Emerg Dis. (2015) 62:603–11. doi: 10.1111/tbed.12430

PubMed Abstract | Crossref Full Text | Google Scholar

4. Sun, D, Vannucci, F, Knutson, TP, Corzo, C, and Marthaler, DG. Emergence and whole-genome sequence of Senecavirus a in Colombia. Transbound Emerg Dis. (2017) 64:1346–9. doi: 10.1111/tbed.12669

PubMed Abstract | Crossref Full Text | Google Scholar

5. Wu, Q, Zhao, X, Bai, Y, Sun, B, Xie, Q, and Ma, J. The first identification and complete genome of Senecavirus a affecting pig with idiopathic vesicular disease in China. Transbound Emerg Dis. (2017) 64:1633–40. doi: 10.1111/tbed.12557

Crossref Full Text | Google Scholar

6. Leme, RA, Alfieri, AF, and Alfieri, AA. Update on Senecavirus infection in pigs. Viruses. (2017) 9:170. doi: 10.3390/v9070170

PubMed Abstract | Crossref Full Text | Google Scholar

7. Leme, RA, Oliveira, TE, Alcantara, BK, Headley, SA, Alfieri, AF, Yang, M, et al. Clinical manifestations of Senecavirus a infection in neonatal pigs, Brazil, 2015. Emerg Infect Dis. (2016) 22:1238–41. doi: 10.3201/eid2207.151583

PubMed Abstract | Crossref Full Text | Google Scholar

8. Fernandes, MHV, de Lima, M, Joshi, LR, and Diel, DG. A virulent and pathogenic infectious clone of Senecavirus a. J Gen Virol. (2021) 102:001643. doi: 10.1099/jgv.0.001643

PubMed Abstract | Crossref Full Text | Google Scholar

9. Oliveira, TES, Leme, RA, Agnol, AMD, Gerez, JR, Pelaquim, IF, Miyabe, FM, et al. Seneca Valley virus induces immunodepressionin suckling piglets by selective apoptosis of B lymphocytes. Microb Pathog. (2021) 158:105022. doi: 10.1016/j.micpath.2021.105022

PubMed Abstract | Crossref Full Text | Google Scholar

10. Hellen, CU, and de Breyne, S. A distinct group of hepacivirus/pestivirus-like internal ribosomal entry sites in members of diverse picornavirus genera: evidence for modular exchange of functional noncoding RNA elements by recombination. J Virol. (2007) 81:5850–63. doi: 10.1128/JVI.02403-06

PubMed Abstract | Crossref Full Text | Google Scholar

11. Willcocks, MM, Locker, N, Gomwalk, Z, Royall, E, Bakhshesh, M, Belsham, GJ, et al. Structural features of the Seneca Valley virus internal ribosome entry site (IRES) element: a picornavirus with a pestivirus-like IRES. J Virol. (2011) 85:4452–61. doi: 10.1128/JVI.01107-10

PubMed Abstract | Crossref Full Text | Google Scholar

12. Zhao, D, Li, Y, Li, Z, Zhu, L, Sang, Y, Zhang, H, et al. Only fourteen 3′-end poly(a)s sufficient for rescuing Senecavirus a from its cDNA clone, but inadequate to meet requirement of viral replication. Virus Res. (2023) 328:199076. doi: 10.1016/j.virusres.2023.199076

PubMed Abstract | Crossref Full Text | Google Scholar

13. Liu, F, Wang, Q, Huang, Y, Wang, N, and Shan, H. A 5-year review of Senecavirus a in China since its emergence in 2015. Front Vet Sci. (2020) 7:567792. doi: 10.3389/fvets.2020.567792

PubMed Abstract | Crossref Full Text | Google Scholar

14. Choudhury, SM, Ma, X, Zeng, Z, Luo, Z, Li, Y, Nian, X, et al. Senecavirus a 3D interacts with NLRP3 to induce IL-1β production by activating NF-κB and Ion Channel signals. Microbiol Spectr. (2022) 10:e0209721. doi: 10.1128/spectrum.02097-21

PubMed Abstract | Crossref Full Text | Google Scholar

15. Liu, H, Li, K, Chen, W, Yang, F, Cao, W, Zhang, K, et al. Senecavirus a 2B protein suppresses type I interferon production by inducing the degradation of MAVS. Mol Immunol. (2022) 142:11–21. doi: 10.1016/j.molimm.2021.12.015

PubMed Abstract | Crossref Full Text | Google Scholar

16. Zhao, K, Guo, XR, Liu, SF, Liu, XN, Han, Y, Wang, LL, et al. 2B and 3C proteins of Senecavirus a antagonize the antiviral activity of DDX21 via the caspase-dependent degradation of DDX21. Front Immunol. (2022) 13:951984. doi: 10.3389/fimmu.2022.951984

PubMed Abstract | Crossref Full Text | Google Scholar

17. Duhan, L, Kumari, D, Naime, M, Parmar, VS, Chhillar, AK, Dangi, M, et al. Single-cell transcriptomics: background, technologies, applications, and challenges. Mol Biol Rep. (2024) 51:600. doi: 10.1007/s11033-024-09553-y

PubMed Abstract | Crossref Full Text | Google Scholar

18. Maher, CA, Kumar-Sinha, C, Cao, X, Kalyana-Sundaram, S, Han, B, Jing, X, et al. Transcriptome sequencing to detect gene fusions in cancer. Nature. (2009) 458:97–101. doi: 10.1038/nature07638

PubMed Abstract | Crossref Full Text | Google Scholar

19. Zeng, J, Cao, D, Yang, S, Jaijyan, DK, Liu, X, Wu, S, et al. Insights into the transcriptome of human cytomegalovirus: a comprehensive review. Viruses. (2023) 15:1703. doi: 10.3390/v15081703

PubMed Abstract | Crossref Full Text | Google Scholar

20. Deb, R, Sengar, GS, Sonowal, J, Pegu, SR, Das, PJ, Singh, I, et al. Transcriptome signatures of host tissue infected with African swine fever virus reveal differential expression of associated oncogenes. Arch Virol. (2024) 169:54. doi: 10.1007/s00705-023-05959-4

Crossref Full Text | Google Scholar

21. Wen, W, Yin, M, Zhang, H, Liu, T, Chen, H, Qian, P, et al. Seneca Valley virus 2C and 3C inhibit type I interferon production by inducing the degradation of RIG-I. Virology. (2019) 535:122–9. doi: 10.1016/j.virol.2019.06.017

PubMed Abstract | Crossref Full Text | Google Scholar

22. Hou, L, Dong, J, Zhu, S, Yuan, F, Wei, L, Wang, J, et al. Seneca valley virus activates autophagy through the PERK and ATF6 UPR pathways. Virology. (2019) 537:254–63. doi: 10.1016/j.virol.2019.08.029

PubMed Abstract | Crossref Full Text | Google Scholar

23. Li, L, Bai, J, Fan, H, Yan, J, Li, S, and Jiang, P. E2 ubiquitin-conjugating enzyme UBE2L6 promotes Senecavirus a proliferation by stabilizing the viral RNA polymerase. PLoS Path. (2020) 16:e1008970. doi: 10.1371/journal.ppat.1008970

PubMed Abstract | Crossref Full Text | Google Scholar

24. Liu, T, Li, X, Wu, M, Qin, L, Chen, H, and Qian, P. Seneca Valley virus 2C and 3C(pro) induce apoptosis via mitochondrion-mediated intrinsic pathway. Front Microbiol. (2019) 10:1202. doi: 10.3389/fmicb.2019.01202

PubMed Abstract | Crossref Full Text | Google Scholar

25. Liu, F, Ni, B, and Wei, R. Comparative proteomic profiling: cellular metabolisms are mainly affected in Senecavirus A-inoculated cells at an early stage of infection. Viruses. (2021) 13:1036. doi: 10.3390/v13061036

PubMed Abstract | Crossref Full Text | Google Scholar

26. Liu, F, Ni, B, and Wei, R. Senecavirus A- and non-infected cells at early stage of infection: comparative Metabolomic profiles. Front Cell Infect Microbiol. (2022) 11:736506. doi: 10.3389/fcimb.2021.736506

PubMed Abstract | Crossref Full Text | Google Scholar

27. Liu, F, Huang, Y, Wang, Q, Li, J, and Shan, H. Rescue of Senecavirus a to uncover mutation profiles of its progenies during 80 serial passages in vitro. Vet Microbiol. (2021) 253:108969. doi: 10.1016/j.vetmic.2020.108969

PubMed Abstract | Crossref Full Text | Google Scholar

28. Zhao, X, Wu, Q, Bai, Y, Chen, G, Zhou, L, Wu, Z, et al. Phylogenetic and genome analysis of seven senecavirus a isolates in China. Transbound Emerg Dis. (2017) 64:2075–82. doi: 10.1111/tbed.12619

PubMed Abstract | Crossref Full Text | Google Scholar

29. Ding, L, Li, Q, Chakrabarti, J, Munoz, A, Faure-Kumar, E, Ocadiz-Ruiz, R, et al. MiR130b from Schlafen4(+) MDSCs stimulates epithelial proliferation and correlates with preneoplastic changes prior to gastric cancer. Gut. (2020) 69:1750–61. doi: 10.1136/gutjnl-2019-318817

PubMed Abstract | Crossref Full Text | Google Scholar

30. Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. (2011) 17:10–2. doi: 10.14806/ej.17.1.200

Crossref Full Text | Google Scholar

31. Kim, D, Paggi, JM, Park, C, Bennett, C, and Salzberg, SL. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat Biotechnol. (2019) 37:907–15. doi: 10.1038/s41587-019-0201-4

PubMed Abstract | Crossref Full Text | Google Scholar

32. Qiao, J, Li, H, and Li, Y. Dietary Clostridium butyricum supplementation modifies significantly the liver transcriptomic profile in weaned piglets. J Anim Physiol Anim Nutr. (2020) 104:1410–23. doi: 10.1111/jpn.13326

PubMed Abstract | Crossref Full Text | Google Scholar

33. Anders, S, Pyl, PT, and Huber, W. HTSeq--a Python framework to work with high-throughput sequencing data. Bioinformatics (2015) 31:166–169. doi: 10.1093/bioinformatics/btu638, HTSeq—a Python framework to work with high-throughput sequencing data

PubMed Abstract | Crossref Full Text | Google Scholar

34. Li, D, Zand, MS, Dye, TD, Goniewicz, ML, Rahman, I, and Xie, Z. An evaluation of RNA-seq differential analysis methods. PLoS One. (2022) 17:e0264246. doi: 10.1371/journal.pone.0264246

PubMed Abstract | Crossref Full Text | Google Scholar

35. Yu, G, Wang, LG, Han, Y, and He, QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | Crossref Full Text | Google Scholar

36. Kovaka, S, Zimin, AV, Pertea, GM, Razaghi, R, Salzberg, SL, and Pertea, M. Transcriptome assembly from long-read RNA-seq alignments with StringTie2. Genome Biol. (2019) 20:278. doi: 10.1186/s13059-019-1910-1

PubMed Abstract | Crossref Full Text | Google Scholar

37. Shen, S, Park, JW, Lu, ZX, Lin, L, Henry, MD, Wu, YN, et al. rMATS: robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. Proc Natl Acad Sci USA. (2014) 111:E5593–601. doi: 10.1073/pnas.1419161111

PubMed Abstract | Crossref Full Text | Google Scholar

38. Koboldt, DC, Chen, K, Wylie, T, Larson, DE, McLellan, MD, Mardis, ER, et al. VarScan: variant detection in massively parallel sequencing of individual and pooled samples. Bioinformatics. (2009) 25:2283–5. doi: 10.1093/bioinformatics/btp373

PubMed Abstract | Crossref Full Text | Google Scholar

39. Zhang, HM, Chen, H, Liu, W, Liu, H, Gong, J, Wang, H, et al. AnimalTFDB: a comprehensive animal transcription factor database. Nucleic Acids Res. (2012) 40:D144–9. doi: 10.1093/nar/gkr965

PubMed Abstract | Crossref Full Text | Google Scholar

40. Anders, S, Reyes, A, and Huber, W. Detecting differential usage of exons from RNA-seq data. Genome Res. (2012) 22:2008–17. doi: 10.1101/gr.133744.111

PubMed Abstract | Crossref Full Text | Google Scholar

41. Szklarczyk, D, Morris, JH, Cook, H, Kuhn, M, Wyder, S, Simonovic, M, et al. The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. (2017) 45:D362–8. doi: 10.1093/nar/gkw937

PubMed Abstract | Crossref Full Text | Google Scholar

42. Livak, KJ, and Schmittgen, TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods. (2001) 25:402–8. doi: 10.1006/meth.2001.1262

Crossref Full Text | Google Scholar

43. Kinobe, R, Wiyatno, A, and Artika, IM. Insight into the enterovirus A71: a review. Rev Med Virology. (2022) 32:e2361. doi: 10.1002/rmv.2361

Crossref Full Text | Google Scholar

44. Wang, J, Mou, C, Wang, M, Pan, S, and Chen, Z. Transcriptome analysis of senecavirus A-infected cells: type I interferon is a critical anti-viral factor. Microb Pathog. (2020) 147:104432. doi: 10.1016/j.micpath.2020.104432

PubMed Abstract | Crossref Full Text | Google Scholar

45. Zeng, W, Yan, Q, Du, P, Yuan, Z, Sun, Y, Liu, X, et al. Evolutionary dynamics and adaptive analysis of Seneca Valley virus. Infect Genet Evol. (2023) 113:105488. doi: 10.1016/j.meegid.2023.105488

PubMed Abstract | Crossref Full Text | Google Scholar