Transcriptomic analysis reveals that coxsackievirus B3 Woodruff and GD strains use similar key genes to induce FoxO signaling pathway activation in HeLa cells

Coxsackievirus B3 (CVB3) is a major cause of viral myocarditis in humans. Although there have been studies on CVB3 infection and pathogenesis, the precise disease mechanism is still not clear. In this study, we used RNA-seq technology to compare the transcriptomic profile of virus-infected HeLa cells to that of uninfected cells to identify key genes involved in host-virus interaction. For this, two CVB3 strains, CVB3 Woodruff, an experimental strain, and GD16-69/GD/CHN/2016, a clinical strain, were selected to examine the common mechanisms underlying their infection. Transcriptomic profiles revealed increased expression of the cell cycle genes CCNG2, GADD45B, PIM1, RBM15, KLF10, and RIOK3 and decreased expression of CYBA. The autophagy-related genes ATG12 and YOD1 were found to be upregulated, while the expression of SOD2 and XPO1 increased slightly in infected cells, and only a minor change was observed in GABARAP expression. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed the FoxO signaling pathway to be enriched and showed a close interaction with differentially expressed genes (DEGs) in the protein-protein interaction network. DEGs associated with related pathways such as cell cycle, autophagy, and oxidative stress resistance were also confirmed by qRT-PCR. In summary, the FoxO signaling pathway was activated during infection with both CVB3 strains and was found to have a regulatory role in downstream pathways such as cell cycle, autophagy, oxidative stress resistance, and the antiviral immune response.


Introduction
Coxsackieviruses, a group of non-enveloped single-stranded RNA viruses, belong to the, genus Enterovirus, family Picornaviridae, and cause clinical symptoms ranging from respiratory illness and aseptic meningitis to acute and chronic myocarditis. Coxsackievirus B3 (CVB3) is the most common cause of viral myocarditis in children and young people. Cytopathology induced directly by virus infection and replication and activation of the host immune response are two common causes of myocardial damage. However, the precise mechanism of CVB3-induced viral myocarditis is still not clearly understood [1].
Clinical and experimental studies related to CVB3 infection, the disease mechanism, and the associated signaling pathways have been conducted [1]. One of the most-studied steps of CVB3 infection is invasion, in which the binding of CVB3 to the host is critical. It involves the interaction between the viral capsid protein, the coxsackievirus and adenovirus receptor (CAR), and the co-receptor decayaccelerating factor (DAF), also known as CD55. Once CVB3 enters the host cell, it manipulates the host system for its own reproduction [2]. The virus replication cycle results in the disruption of cellular processes in the host and associated signaling events [3,4].
There are only a few model CVB3 strains for laboratory research, such as CVB3 Nancy and Woodruff, which were developed decades ago in the United States [5]. On the other hand, with each disease outbreak, clinical isolates from patients are evolving with great diversity [6]. China has had CVB3 outbreaks associated with hand, foot, and mouth 1 3 disease (HFMD), and the sequences of viral isolates show regional phylogenetic dynamics, suggesting an association between CVB3 outbreaks and their evolution [7]. In this study, two CVB3 strains with distinct sequences and origins were selected to assess similarities and differences in their infection and pathogenesis mechanisms. The laboratory strain was CVB3 Woodruff (CVB3-W), which was developed in the United States, and the other was the GD16-69/ GD/CHN/2016 (CVB3-GD) strain, which was isolated from a 2016 HFMD outbreak in China and is a representative of the Chinese epidemic strains that have been circulating in recent years. We used RNA-seq technology to analyze the transcriptomic profiles of the host cell after CVB3 infection, with the aim of gaining a broader view of genes involved in the host-virus interaction and better insights into the common disease mechanisms used by both clinical and experimental virus strains to better guide clinical treatment and medicine development.

Cell line and virus strains
HeLa cells were propagated in Dulbecco's modified Eagle's medium (DMEM) containing 100 units of penicillin and 100 μg of streptomycin per ml and 10% fetal bovine serum (FBS; GIBCO, Thermo Fisher, Australia). HeLa cells were infected with the coxsackievirus B3 Woodruff strain (CVB3-W) using a standard protocol described in previous publications [41,42]. Coxsackievirus B3 strain GD16-69/GD/ CHN/2016 (CVB3-GD; MT712317) was kindly provided by the WHO WPRO Regional Polio Reference Laboratory, National Health Commission Key Laboratory for Medical Virology, National Institute for Viral Disease Control and Prevention, Chinese Center for Disease Prevention and Control. The virus titer of CVB3-W was 1.9 × 10 9 plaque-forming units (PFU), and that of CVB3-GD was 3.9 × 10 8 PFU.

CVB3 infection of HeLa cells
Strains CVB3-W (MOI = 10) and CVB3-GD (MOI = 7) were propagated in HeLa cells at 37 °C. After viral inoculation, the cells were washed twice with phosphate-buffered saline and maintained in fresh DMEM with 2% FBS. After 2, 4, and 6 h, infected cells and mock-infected cells were harvested for transcriptome analysis and quantitative realtime PCR (qRT-PCR).

Western blot analysis
HeLa cells were lysed using a commercial lysis buffer (Beyotime), and after centrifugation at 4 °C, the sample protein concentration was determined using an Enhanced BCA Protein Assay Kit (Thermo Scientific). Then, 30 μg of protein from each sample was subjected to 10% sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE) and transferred to nitrocellulose membranes (0.02 mm; GE Health). After blocking for 1 h with 5% nonfat milk, the membrane was incubated overnight with primary antibodies against VP1 (Dako) and actin (Cell Signaling Technology) at 4 °C. The membrane was incubated with horseradishperoxidase-conjugated secondary antibodies (Cell Signaling Technology) for 1 h and then washed with PBS containing 0.1% Tween-20. The target proteins were visualized using a chemiluminescence detection kit (PerkinElmer) [4,43].

RNA isolation and next-generation sequencing of mRNA (RNA-Seq)
At selected time points, samples were collected, and total RNA was isolated using TRIzol Reagent (Sigma, USA) according to the manufacturer's instructions. The RNA was quantified in a Bioanalyzer 2100 system using an RNA Nano 6000 Assay Kit (Agilent Technologies, CA, USA). The mean RIN value of the samples was 9.87 (0.83). The library preparation and RNA-Seq were carried out by Novogene Co., Ltd, Beijing, China. A total of 1 µg of RNA per sample was used for library preparation. The clustering of the sequences was performed using a TruSeq PE Cluster Kit v3-cBot-HS (Illumina) on a cBot Cluster Generation System. The libraries were sequenced on an Illumina HiSeq2000 Novaseq platform [43].

Data processing
Raw sequencing reads in fastq format were processed using in-house Perl scripts (Novogene Co., Ltd). High-quality clean data were used for downstream analysis. The reads were mapped to the reference human genome sequence GRCh38/hg38 (GRCh38.p12), gene model annotation files were downloaded directly from the genome website ENSEMBL (ftp:// ftp. ensem bl. org/ pub/ relea se-94/ assem bly_ chain/ homo_ sapie ns/), and the alignment was built using Hisat2 v2.0.5. The number of gene-level reads mapped to each gene was determined using FeatureCounts v1.5.0-p3, and the number of uniquely mapped reads was calculated. Reads per kilobase per million mapped (RPKM) reads were determined using FeatureCounts v1.5.0-p3 [44].

Differential expression analysis
Analysis of differential expression between transcripts of the control and virus-infected groups was conducted using the DESeq2 R package (1.20.0). Genes with an adjusted P-value (padj) <0.05 and a fold change >1 detected by DESeq2 were selected as differentially expressed genes (DEGs) [45]. Heat map diagrams and Venn diagrams were constructed using the Pheatmap R package and VennDiagram [46], respectively, to present quantitative differences between DEGs. Integrated protein-protein interaction (PPI) network analysis of DEG was visualized using the STRING website (https:// string-db. org/) [47].

GO and KEGG enrichment analysis of differentially expressed genes (DEGs)
Gene Ontology (GO) enrichment analysis [48] was conducted using the clusterProfiler R package to confirm the biological functions of DEGs. The significance of GO terms was determined using a corrected P-value <0.05. Pathway analysis was performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database [49]. The statistical enrichment of differentially expressed genes in KEGG pathways was assessed using the ClusterProfiler R package.

Quantitative RT-PCR (qRT-PCR)
Total RNA was extracted from cells using TRIzol Reagent (Sigma, USA) for qRT-PCR analysis. The primers used for the analysis of 18 selected host target genes are listed in Supplementary Table S1. First-strand cDNA synthesis was carried out using a PrimeScript RT Reagent Kit (TaKaRa, Japan). Then, qRT-PCR was performed using TB Green ® Premix Ex Taq™ II (TaKaRa, Japan) on a LightCycler ® 480 System (Roche Applied Science) following the manufacturer's recommendations. The qRT-PCR conditions were as follows: 95 °C for 30 s, followed by 40 cycles of 5 s at 95 °C and 30 s at 60 °C. GAPDH was used as an internal control. The expression of each target mRNA relative to the internal control was calculated by the 2 -ΔΔCT method [50].

Detection of CVB3 replication
To detect viral replication, HeLa cells were infected with strain CVB3-W (MOI = 10) or CVB3-GD (MOI = 7) for 2, 4, and 6 h, after which the expression of the CVB3 VP1 protein was examined by western blotting. The expression of VP1 protein of the two strains showed similar trends and increased gradually from 2 h postinfection (hpi) to 6 hpi ( Fig. 1A and B). Thus, the growth rates of the two CVB3 strains were appropriate for comparison of transcriptomic changes in host cells.

Analysis of differentially expressed genes (DEGs)
To investigate the mechanism of cellular damage due to CVB3 infection, we analyzed the transcriptomic changes in HeLa cells infected with the two CVB3 strains at 2, 4, and 6 hpi by RNA sequencing. The transcript expression profiles of CVB3-W-or CVB3-GD-infected HeLa cells and uninfected control cells were compared using threedimensional principal component analysis (PCA), and the results revealed distinct variation between the control samples and the samples infected with two different virus strains in a time-dependent manner ( Fig. 2A). In the differential expression analysis, genes with an adjusted P-value <0.05 and fold change >1 were considered DEGs. Analysis   S1B). We also found that 396 DEGs associated with CVB3-GD infection enriched a different set of pathways. The GO enrichment results showed the significance of pathways such as nuclear export and chromosomal regions responsible for cellular processes, muscle and cardiac ventricular development pathways, actin filament activity, and ribosomal protein binding functions (Supplementary Fig. S1C). KEGG analysis showed that, in CVB3-W infection, pathways involving the cell cycle, FoxO signaling, and viral myocarditis were significantly enriched, in addition to pathways related to the immune response and infection ( Supplementary Fig. S1D).
To identify pathways affected by infection with both viral strains, KEGG pathway enrichment analysis was performed on samples infected with each virus strain at all three time points. As shown in Fig. 2D, the top enriched pathways associated with CVB3-W-infection-induced DEGs were cell cycle, FoxO signaling pathway, RNA transport, endocytosis, mismatch repair, and ubiquitin-mediated proteolysis. In the KEGG pathway enrichment analysis of CVB3-GDinfected samples, DEGs were associated with cell cycle, proteoglycans in cancer, the longevity-regulating pathway, human T-cell leukemia virus 1 infection, RNA transport, and the FoxO signaling pathway (Fig. 2E). Thus, the common pathways associated with the infection of these two virus strains were cell cycle, RNA transport, and the FoxO signaling pathway. Particularly, both strains were able to stimulate the FoxO signaling pathway, which functions as an upstream regulator of enrichment pathways associated with other DEGs [8,9].

Identical DEGs in CBV3-infected HeLa cells
To examine the changes in the expression of identical DEGs induced by CVB3-W and CVB3-GD, all DEGs were compared at 2, 4, and 6 hpi. As shown in the Venn diagram (Fig. 3), 58 common DEGs were found to be induced by two strains at three time points. Of these 58 differentially expressed genes, 43 were upregulated and 15 were downregulated (Supplementary Table S2). GO and KEGG enrichment analyses were also performed to further evaluate the mechanisms associated with infection with both viral strains. The results revealed that identical DEGs were involved in similar cellular pathways, and this was consistent with the individual enrichment analysis results. GO enriched terms such as cellular response to peptide and extracellular stimulus, processes involving cell division and adhesion junctions, and transcriptional regulation were of significance. The top listed KEGG pathways were signaling transduction and viral myocarditis, autophagy, and immune response-related pathways ( Fig. 2D and E, Supplementary Fig. S1). Particularly, many genes among the 58 DEGs were involved in the pathways that coincided with the enrichment results and are downstream of the FoxO signaling pathway. For example, UBE2S, PIM1, and KLF10 are cell cycle regulation factors [10][11][12], and PHAX and RBM27 are RNA-binding proteins related to RNA transport [13,14].

The protein-protein interaction (PPI) network of DEGs
FoxO proteins are a group of transcription factors with essential regulatory roles in the expression of genes associated with cellular processes such as apoptosis, cell-cycle progression, and oxidative stress resistance. To identify the key regulatory pathways and DEGs, we selected genes from the 58 common DEGs, the most enriched KEGG pathways, and the FoxO signaling pathway for protein-protein interaction network analysis using the STRING Protein-Protein Interaction Networks Functional Enrichment Analysis website (https:// string-db. org/). As shown in Fig. 4, the nodes represent DEGs and pathway proteins, and their expression levels are indicated in specific colors. Among them, all four FoxO genes (FoxO1, FoxO3, FoxO4, and FoxO6) in the network were found to interact strongly with each other, indicating the essential role of the FoxO signaling pathway. The enrichment of FoxO signaling pathways in KEGG analysis (Fig. 4B) clearly shows that the pathway regulation correlated with DEGs from samples infected with CVB3-W and CVB3-GD strains in comparison to that in control cells. Both CVB3 strains induced common DEGs in HeLa cells, including those associated with FoxO signaling pathways, upregulated EDN1, RIOK3, PIM1, KLF10, RBM15, OTUD3, YOD1, XPO1, and downregulated CYBA.

Confirmation of transcriptomic data and expression of key DEGs by qRT-PCR
To confirm the expression of selected genes among the 58 common DEGs, we performed qRT-PCR of eight  (Fig. 5B, C, and D). Meanwhile, to assess the significance of CVB3 infection on pathway regulation, the expression of nine FoxO signaling proteins selected from the PPI analysis was examined. In accordance with the RNA-seq and PPI results, expression of FoxO1 (Fig. 5A), CCNG2, ATG12, GADD45B, and SOD2 were found to be upregulated (Fig. 5B, C, and D), whereas that of FoxO6 (Fig. 5A) was gradually downregulated at 2, 4, and 6 hpi.

Discussion
CVB3 infection is able to cause viral myocarditis in individuals of all ages, with clinical symptoms ranging from classical respiratory symptoms and central nervous system infection, which can lead to meningitis and acute flaccid paralysis (AFP), to HFMD-like skin lesions [15,16]. In this study, we selected two CVB3 strains, CVB3-GD, isolated from a Chinese HFMD outbreak, causing mostly HFMDlike symptoms in our cases, and CVB3-W, an experimental strain developed in the United States that has the ability to cause all of the above symptoms, to study the transcriptomic profile of the infected host cells. By comparing differences in their mode of infection, we hope to find common pathogenic mechanisms and appropriate treatments for CVB3 that are more suitable for cases of viral myocarditis. The CVB3 Woodruff strain is evolutionarily similar to the CVB3 Nancy strain, belonging to subtype A, and it is no longer prevalent worldwide, while CVB3-GD is evolutionarily similar to other Chinese epidemic strains, belonging to subtype D (the dominant subtype circulating in China in recent years). In addition to their geographical and evolutionary differences, the specific sequence pattern in the virus-host interface is also notable. In their study, Pan et al. found that the substitution of single amino acid residues in VP3 (E243Q) and VP2 (D138N) of CVB3 could disrupt the binding between CVB3 and the host DAF receptor [17]. In CVB3-W, the corresponding amino acids of the capsid proteins VP3 and VP2 are 234Q and 138N respectively, while those of CVB3-GD are 234D and 138D, respectively. During virus invasion, both strains bind loosely to host DAF receptors with different molecular interactions. Differences in binding between CVB3 strains could result in strain-or cell-specific variation in virus replication, virus spread, and pathology [18]. For example, strong DAF binding of many human enteroviruses does not normally promote lytic infection of cells, and CAR binding usually leads to lytic infection [19]. The interaction of DAF and cardiovirulent CVB3 is important in the pathogenesis of CVB-mediated heart disease [20]. For a better understanding of viral pathogenesis, in this study, we first confirmed by western blotting using a VP1 antibody that the viruses were in a state of gradual proliferation in the host cells (Fig. 1). We then performed RNA-seq for transcriptomics analysis of CVB3-infected HeLa cells to identify genes associated with host-virus interactions, activation of cellular signaling pathways, and the immune response. Since its discovery, extensive research has focused on the pathological mechanisms of CVB3, which have been investigated using HeLa cells. In CVB3-infected HeLa cells, several mechanisms, including cell death, the immune response, and autophagy are at work [21,22]. Our RNA-Seq results reveal that although HeLa is a fairly stable cell line, infection with both CVB3 strains could still stimulate a number of diverse cellular pathways associated with the cell cycle, apoptosis, and RNA transport, and both infections resulted in enrichment of the FoxO signaling pathway (Fig. 2D and E), suggesting that this pathway plays a role in disease progression in infection with both experimental and epidemic CVB3 strains.
FoxO proteins belong to the forkhead protein family, characterized by a conserved "forkhead box" DNAbinding domain, are crucial transcription factors associated with the expression of genes related to apoptosis, cell-cycle progression, and resistance to oxidative stress [23]. These cellular functions play an important role in host-virus interactions and innate immunity. In Sendai virus infection, FoxO1 acts as a negative regulator of retinoic acid-inducible gene-I (RIG-I)-triggered signaling, promotes virus replication, and downregulates type I interferon (IFN) production [24]. Likewise, Japanese encephalitis virus (JEV) infection results in a decrease in FoxO expression and induces apoptosis by inhibiting the STAT3-Foxo-Bcl-6/p21 pathway, which had no effect on JEV replication [25]. In our study, one of the DEGs that were upregulated after infection with both CVB3 strains was EDN1, which encodes endothelin-1, a peptide survival factor in colon cancer that is regulated directly by FoxO1 and NF-κB pathways [26]. In this study, qRT-PCR analysis of samples from HeLa cells infected by both CVB3 strains showed upregulated FoxO1 and downregulated FoxO6, while there was little effect on FoxO3 and FoxO4 expression (Fig. 5A). Western blot analysis confirmed that the FoxO1 protein was downregulated and that the FoxO6 protein was upregulated ( Supplementary Fig. S4). This is likely due to regulation by a negative feedback loop. Indeed, previous studies have shown that FoxO family transcription factors can suppress their own expression by negatively regulating gene expression upon apoptosis [27]. Also, compared with our previous mouse model results, we found a similar regulatory trend of FoxO6 transcriptomic expression (Supplementary Table S2) and FoxO signaling pathway enrichment in mouse myocardial tissue with acute myocarditis induced by CVB3-W [28]. These results show that the FoxO protein is involved in Our RNA-seq data show the significance of varying expression patterns in HeLa cells after CVB3 infection. Combined with KEGG pathway enrichment results, the selected DEGs were found to be closely associated with cell cycle regulation, apoptosis, autophagy, oxidative stress, DNA repair, and immune regulation, which are regulated by the upstream FoxO signaling pathway [23]. We observed that the cell cycle was one of the most significantly enriched KEGG pathways downstream of the FoxO pathway. In fact, qRT-PCR results showed an upregulation of CCNG2 and GADD45B, which encode FoxO cell cycle signaling proteins. CVB3 infection also upregulated four cell-cyclerelated genes: PIM1, RBM15, KLF10, and RIOK3. PIM1 is a serine/threonine-protein kinase that is involved in the cell cycle and cell survival, and its active form can phosphorylate PRAS40, which binds to pFoxO3a to downregulate apoptosis [29]. RNA binding motif protein 15 (RBM15) is crucial in cell survival and cell cycle progression and has an inhibitory effect on apoptosis in CML cells [30]. The multifunctional roles of these genes confirmed the close interaction between the cell cycle and apoptotic regulation in CVB3 infection. Krueppel-like factor 10 (KLF10), also known as Tieg1, regulates the transcription of genes involved in multiple pathways involving proliferation, apoptosis, and cell cycle arrest. A study using a TIEG1 knockout mouse model demonstrated its role in cardiac hypertrophy [12]; therefore, KLF10 may be a potential target gene in CVB3-induced cardiac pathology. Another DEG, encoding the serine/threonine-protein kinase RIOK3, could bridge TBK1 and IRF3, activate the type I IFN pathway, and play a critical role in the defense against both DNA and RNA viruses through the innate immune response [31]. Therefore, it is a potential target gene for CVB3 antiviral treatment.
Autophagy is another signaling event that is regulated by FoxO signaling. Studies have already revealed that replication of CVB3 is promoted by autophagosome and autophagic flux [32]. A study by Giansanti et al. identified autophagy protein mTORC1 signaling as a major regulatory network in CVB3-infected cells, emphasizing the importance of autophagy signaling at the protein level [33]. In our KEGG analysis, autophagy was also one of the most significantly enriched signaling pathways ( Fig. 2D and E). Unlike genes encoding cell cycle proteins, our qRT-PCR results revealed that the autophagy gene ATG12 was only slightly upregulated during infection. Until 6 hpi, this was observed only in CVB3-GD-infected samples, while there was no obvious change in the expression of GABARAP during this infection period. Accordingly, the expression of autophagy-related DEGs and XPO1 did not change significantly during infection. Exportin-1 (XPO1) is a nuclear export protein that can be phosphorylated by STK38 and plays an important role in the export of the autophagy regulator BCEN1 [34]. A recently published study showed XPO1 to be a promising antiviral target for blocking access of the virus to the cell [35]. RNA-seq and qRT-PCR results revealed significantly upregulated YOD1 expression after CVB3 infection. The deubiquitinating enzymatic function of YOD1 is important for clearing damaged lysosomes and promoting autophagosome formation [36]. In the PPI network, the only DEG directly connected with YOD1 was OTUD3, which encodes another deubiquitinating enzyme that is upregulated by CVB3 infection and is important in several cellular processes, including cell cycle regulation. RNA virus infection can stimulate OTUD3 deacetylation by sirtuin 1 (SIRT1), shut down its deubiquitinating function, release mitochondrial antiviral-signaling protein (MAV), and promote the rapid propagation of innate immune signaling for host antiviral defense [37]. Hence, we propose a similar antiviral function during CVB3 infection.
Oxidative stress and DNA repair signaling pathways are also stimulated by CVB3 infection and are regulated by FoxO signaling. As confirmed by qRT-PCR, GADD45B was upregulated, while SOD2 was only slightly upregulated until 6 hpi. The protein encoded by the DEG CYBA is a crucial subunit of NADPH oxidase; it is involved in oxidative stress, and its role in HCV and HBV disease mechanism has been studied [38,39]. In this study, we observed that virus infection could downregulate CYBA expression. A study on HIV showed that HIV Nef could act on NADPH oxidase complex units, thus altering superoxide production in neutrophils and modulating the innate defense mechanism [40]. We predict similar mechanisms in CVB3 infection. Altogether, many of the genes selected in our study are not only components of pathways regulated by FoxO signaling but also have independent roles in other cellular processes or have a direct function in the virus life cycle and may therefore be potential targets for antiviral therapy.

Conclusion
In this study, we used RNA-seq technology to compare the transcriptome profiles of HeLa cells infected with an experimental isolate and a clinical CVB3 isolate. We found that, despite subtle molecular changes due to regional evolution, the FoxO signaling pathway was activated by infection with both of these strains. KEGG enrichment and PPI network analysis of the DEGs indicated a regulatory effect of FoxO proteins on the downstream pathways associated with cell cycle, apoptosis, autophagy, oxidative stress, the DNA repair response, and the antiviral immune response, which could either promote host cell survival or stimulate the antiviral defense response after CVB3 infection. The associations of these key genes with signaling pathways need to be