Integrative transcriptome analysis reveals alternative polyadenylation potentially contributes to GCRV early infection

Introduction Grass carp reovirus (GCRV), a member of the Aquareovirus genus in the Reoviridae family, is considered to be the most pathogenic aquareovirus. Productive viral infection requires extensive interactions between viruses and host cells. However, the molecular mechanisms underlying GCRV early infection remains elusive. Methods In this study we performed transcriptome and DNA methylome analyses with Ctenopharyngodon idellus kidney (CIK) cells infected with GCRV at 0, 4, and 8 h post infection (hpi), respectively. Results We found that at early infection stage the differentially expressed genes related to defense response and immune response in CIK cells are activated. Although DNA methylation pattern of CIK cells 8 hpi is similar to mock-infected cells, we identified a considerable number of genes that selectively utilize alternative polyadenylation sites. Particularly, we found that biological processes of cytoskeleton organization and regulation of microtubule polymerization are statistically enriched in the genes with altered 3’UTRs. Discussion Our results suggest that alternative polyadenylation potentially contributes to GCRV early infection.


Introduction
Aquareoviruses are nonenveloped viruses and classified within the family Reoviridae, a family of double-stranded RNA virus composed of aquareoviruses, mammalian reoviruses (MRV), and the other 13 genera.Aquareoviruses cause infection in aquatic organisms including bony fish, shellfish, and crustacean worldwide (Lupianni et al., 1995).Although most aquareoviruses are isolated from seemingly healthy fish and do not give rise to high mortalities, grass carp reovirus (GCRV) is recognized to be most pathogenic among the isolated aquareoviruses (Rangel et al., 1999).GCRV can cause serious hemorrhagic disease in aquatic organisms.Our previous studies have shown that GCRV can induce cell-cell fusion and produce characteristic cytopathic effect (CPE) consisting of large syncytia within infected cultures (Fang et al., 1989;Ke et al., 1990), and it has been extensively used to understand aquareovirus molecular and structural biology.Seven structural (VP1-VP7) and six nonstructural proteins Tan et al. 10.3389/fmicb.2023.1269164Frontiers in Microbiology 02 frontiersin.org(NS12, NS16, NS26, NS31, NS38, and NS80) of GCRV have been well identified (Guo et al., 2013;Yan et al., 2015).Comparative proteomic analysis of lysine acetylation in fish Ctenopharyngodon idellus kidney (CIK) cells reveals the proteome-wide changes in host cell acetylome with GCRV infection (Guo et al., 2017).MRV can cause chronic infection.It has been revealed that there is a close molecular evolutionary relationship between aquareoviruses and mammalian orthoreoviruses.In addition to morphological similarity, GCRV and MRV share high amino acid conservation.A better knowledge of the interaction during early infection stage between GCRV and host cells will help the understanding molecular pathogenesis of the aquareovirus and other members in the family Reoviridae.
Accumulating evidence has demonstrated that epigenetics is actively involved in host-virus interaction.Epigenetic trait is defined as a stably heritable phenotype resulting from changes in a chromosome without alterations in the DNA sequence (Berger et al., 2009).These chromosomal changes include methylation of cytosine in CpG dinucleotides (often referred to as DNA methylation) and other posttranslational covalent modifications to histones, such as methylation, acetylation, and ubiquitylation.The epigenetic modifications are associated with structural organization of chromatin and transcriptional activities of the affected genes.As intracellular parasites, viruses develop various ways of remodeling epigenetic alterations to facilitate their infection and replication.Through inducing DNA methylation changes in host cells viruses epigenetically manipulate host functions upon virus infection.For instance, Epstein-Barr virus (EBV) infection activates cellular DNA methyltransferases and results in aberrant DNA methylation in host cells (Tsai et al., 2006;Hino et al., 2009).HIV infection can also trigger the differential DNA methylation at cis-regulatory regions of host genomic DNA and inhibit the function of T cells (Pion et al., 2013;Youngblood et al., 2013).Nevertheless, the influence on cellular DNA methylation during GCRV infection remains to be further characterized.
In addition to epigenetic modifications, formation of stress granules is also actively involved in the interaction between viruses and host cells.It has been recognized that the innate immune response of host cells is triggered by upon virus infection to prevent pathogen invasion, partially through stress granules.Some components of stress granules have been identified, such as T-cell-restricted intracellular antigen 1 (TIA-1), TIA-1-related protein (TIAR), Ras GTPase-activating protein-binding proteins (G3BPs)and poly(A)-binding proteins (PABPs).PABPs are a family of RNA recognition motif (RRM)-containing proteins that bind poly(A) tail and regulate translation and stability of mRNA.The previous report has demonstrated that alternative polyadenylation (APA) plays an important role in the antiviral innate immune response (Jia et al., 2017).However, it remains unclear whether APA of host cells is involved in GCRV infection.Thus, in this study we carried out integrative analyses of transcriptome, DNA methylome and APA in GCRV-infected CIK cells for understanding the molecular events in GCRV early infection.

Viral infection, cytopathic effect observation and plaque assay
The infection assays were carried out as we described previously (Zhang et al., 2019).Briefly, the 80% confluent CIK cells in T-25 flask (Corning Inc., Corning, NY, United States) with a concentration of 2 × 10 6 cells/ml were inoculated with GCRV at a multiplicity of infection (MOI) of 1 in serum-free MEM medium at 28°C for 1 h following the method as previously reported (Guo et al., 2017).For comparison, the mock-infected cells were treated with same amount of medium in the same conditions.Upon adsorption, cells were washed with phosphate-buffered saline (PBS) to remove non-adsorbed virions.The infected cells were maintained in MEM-2 at 28°C and harvested at 0 (mock), 4 and 8 h post infection, respectively.When initial cytopathic effects were observed, the infected cells and mock-infected cells were prepared and harvested for further transcriptome analyses.Three rounds of independent experiments were performed.For MOI determination, plaque assays were done according to our previously described method (Yan et al., 2015;Zhang et al., 2018).
RNA isolation, RNA-seq library construction and deep sequencing CIK Cells were infected by GCRV for 0, 4 and 8 h, respectively.Total RNA was extracted with Trizol reagent (Invitrogen, United States), which was further treated with RNase-free DNase to remove genomic DNA.mRNA was purified with poly(dT) oligoattached magnetic beads and broken down into 200 ~ 400 bp fragments.A strand-specific RNA-seq library was constructed with NEBNext Ultra Directional RNA Library Prep Kit (NEB, New England, United States).Briefly, the fragmented mRNA was reversely transcribed into cDNA with random primers and then the second-strand cDNA was generated.The resulting double-strand DNA fragments were purified with AMPure beads (Beckman Coulter, Brea, CA, United States) and ligated with Illumina adapters.The ligation products were purified by agarose gel electrophoresis to remove adapter dimmers, which were subsequently subjected to HiSeq X sequencing (Illumina, San Diego, CA, United States).The raw sequencing data could be obtained in the EMBL database1 under accession number E-MTAB-13002.

MeDIP-seq library construction and deep sequencing
Genomic DNA of CIK cells was extracted using GenElute™ Mammalian Genomic DNA Miniprep Kit (Sigma, United States).DNA was randomly sheared into fragments of 200 ~ 500 bp and used for library preparation with NEBNext ® Ultra™ II DNA Library Prep Kit for Illumina (NEB), the resulting libraries were purified with 1 × Agencourt AMPure XP beads (Beckman Coulter).The immunoprecipitation procedure was basically performed according to a previous MeDIP protocol (Taiwo et al., 2012).Briefly, the library DNA was denatured at 95°C for 10 min and immediately placed into an ice for 10 min, 1/10 volume of denatured product was set aside as Input.The Protein A + G magnetic beads (Millipore, United States) were incubated with 5-Methylcytosine (5-mC) monoclonal antibody (Epigentek) at 4°C for 2 h and the library was incubated with antibody-bead complexes at 4°C overnight with a slight rotation.The dynabeadantibody-methylated DNA complexes were washed three times, followed by proteinase K (Thermo scientific) treatment for 3 h at 55°C.The immunoprecipitated DNA was extracted by phenol/ chloroform/isoamylalcohol, followed by ethanol precipitation, and resuspended in EB buffer (10 mM Tris-HCl pH 8.0).The enriched methylated DNA and Input DNA were amplified using Q5 High-Fidelity DNA Polymerase (NEB), and subject to Illumina sequencing platforms.The raw sequencing data could be obtained in the EMBL database (see Footnote 1) under accession number E-MTAB-13003. 3

Bioinformatics analysis
The raw reads with low quality and the adapter sequences of RNA-seq and MeDIP-seq data were removed using Cutadapt v4.1 (Kechin et al., 2017).For RNA-seq data, clean reads were mapped to the grass carp reference genome (Wang et al., 2015) using Hisat v2.2.1 (Kim et al., 2019).The Subread toolkits was used to quantify read counts for genes (Liao et al., 2014), and reads per kilobase of transcript per million mapped reads (RPKM) were calculated as expression levels.Differential expression analysis was performed using the edgeR package in R platform v3.6.3 (Robinson et al., 2010).Those genes with an value of p < 0.05 and fold change >1.5 were regarded as differentially expressed genes (DEGs).For MeDIP-seq data, clean reads were mapped to the grass carp reference genome using Bowtie v2.4.5 (Langmead and Salzberg, 2012).PCR duplicate reads were removed with Picard v2.27.4.4DNA methylation peaks were called with MACS2 with deduplicated alignments (Zhang et al., 2008) and the differentially methylated regions (DMRs) were identified with DiffBind and DESeq2 packages (Anders and Huber, 2010).Functional enrichment analysis was performed with DAVID.5

Analysis of APA with RNA-seq data
The APAs were identified by using DaPars algorithm (Masamha et al., 2014) based on RNA-seq data.Briefly, the observed sequence coverage was represented as a linear combination of annotated 3'UTRs.For each transcript with annotated proximal adenylation site (PAS), a regression model was used to infer the end point of alternative novel PAS within this 3' UTR at single nucleotide resolution, by minimizing the deviation between the observed read coverage and the expected read coverage based on a two-PAS model in both mock-infected and GCRV-infected samples simultaneously.A percentage dPAS usage index (PDUI) was utilized to define shortening (negative index) or lengthening (positive index) of 3'UTR and thus capable of quantifying the degree of difference in 3'UTR usage between mock-infected and GCRV-infected CIK cells.The greater PDUI means that the more distal PAS of a given transcript is used and vice versa.

Grass carp reovirus infection-induced cytopathic effects at early stage
To characterize the interaction of GCRV and host cells for integrative analyses of transcriptome in GCRV-infected CIK cells, we firstly carefully examined the cytopathic effects induced by GCRV infection at early stage.In both mock-infected cells and GCRV-infected cells at 4 h, we did not observe obvious CPE.As infection progressed, we detected an initial characteristic CPE on the monolayers of CIK cells at 8 h post infection (hpi) by comparing to mock-infected cell (Figure 1), which suggests that efficient infection was obtained, and the harvested infected cell lysates were suitable for follow-up transcriptome related assays.

Transcription program associated with grass carp reovirus early infection
To detect the molecular events at the GCRV early infection, we performed RNA-seq analysis of CIK cells 0, 4 h and 8 hpi.Totally we generated 197.5 millions raw sequencing reads in the groups of mock, 4 and 8 hpi.Among these reads 94.7% are mappable and are used for downstream analysis.
Totally we identified 15,255 expressed genes in three groups.We then used gene set enrichment analysis (GSEA) to compare the transcriptome data between CIK cells 8 hpi and the MOCK-infected cells.We found that several gene sets were significantly enriched in cells 8 hpi comparing with the MOCK, such as defense response to virus, immune response, and cholesterol metabolic process (Figure 2A).Comparing with MOCK, we identified 675 differentially expressed genes in CIK cells 8 hpi (Figure 2B).Gene ontology (GO) analysis indicates that the biological processes of defense response to virus, cholesterol metabolic process (Figure 2C) and mitogen-activated protein kinase (MAPK) signaling pathway (Figure 2D) are significantly enriched in these differentially expressed genes.

DNA methylation pattern in grass carp reovirus early infected Ctenopharyngodon idellus kidney cells
It has been reported that DNA methylation contributes to resistance to GCRV infection (Shang et al., 2017).Here, we asked whether DNA methylation is involved in GCRV early infection.To address this issue, we performed MeDIP-seq analysis with MOCK and CIK cells 8 hpi.We totally generated 232 million sequence reads, and 96% are mappable.The data sets of biological replicates are highly correlated (Supplementary Figure).Among 9,426 identified methylation sites 37.7% are located in intergenic regions, 32% in exons, 21% in introns and only 9.3% in promoter regions (Figure 3A).We examined DNA methylation signal around transcription starting site (TSS) and found the obviously enriched methylation signal at TSS regions both in MOCK and 8 hpi groups (Figure 3B).It is well recognized that DNA methylation is negatively associated with gene expression.We then examined the correlation between methylated genomic regions and transcription levels.We observed that the methylated regions at promoters, exons and introns are weakly and negatively correlated with transcription (Figure 3C).Compared with the MOCK, we found the DNA methylation of CIK cells 8 hpi is very similar to MOCK (Figure 3D), suggesting that DNA methylation pattern is less functionally involved in early GCRV infection.

Alternative polyadenylation profile in grass carp reovirus early infected Ctenopharyngodon idellus kidney cells
Since DNA methylation is less involved in GCRV early infection, we next investigated other mechanisms.Alternative polyadenylation (APA) modulates gene expression and has been reported to be involved in antiviral response.We then examined the APA patterns between MOCK and 8 hpi group.Comparing with the MOCK, we identified 404 genes with the APA-derived altered 3'UTRs, including 201 genes with lengthened 3'UTRs and 203 genes  with shortened 3'UTRs (Figure 4A).When examining the 3'UTR alterations and transcription levels, we observed that the overall transcription levels of the genes with shortened 3'UTRs is higher than the those with lengthened 3'UTRs (Figure 4B).Through GO analysis with the genes containing altered 3'UTRs we found the biological processes of cytoskeleton organization and regulation of microtubule polymerization are statistically enriched (Figure 4C).In particular, we observed that Camsap1b, a gene involved in microtubule formation and stability, preferentially utilized the proximal poly(A) sites in GRCV-infected CIK cells when comparing MOCK (Figure 4D).These observations suggest that alternative poly(A) usage is potentially involved in the early infection of CIK cells.

Discussion
Viral infections involve intensive interactions between viruses and host cells.As obligate intracellular parasites, viruses misappropriate host cellular machinery to allow their replication; while host cells also orchestrate the transcriptional programs to repress viral infection.For example, in our previous studies we found that aquareovirus NS38 (the GCRV nonstructural protein expressed in host cells as early as 3 h post infection) interacts with host translation initiation factor eIF3A for virus replication (Shao et al., 2013;Zhang et al., 2019).Meanwhile, host innate immune response would be activated after virus infection.Consist with the reported studies (Chen et al., 2012;Shi et al., 2014;Wan and Su, 2015;Dang et al., 2016;Xu et al., 2016;Chen et al., 2018), we observed that the host genes related to defense response to virus and immune response are differentially expressed in CIK cells 8 hpi (Figures 2A,C).In addition to these immune-related genes, we found the genes involved in cholesterol metabolic process and cholesterol biosynthetic process are also activated (Figure 2C), supporting our previous report that cellular membrane cholesterol is required for GCRV productive infection (Zhang et al., 2018).Activation of MAPK signaling pathway has been reported to be required for cell entry of avian reovirus (Huang et al., 2011).Interestingly, in this study we found that this pathway is most significantly enriched among all identified cellular signaling pathways (Figure 2D), suggesting MAPK signaling pathway is involved in GCRV infection.
DNA methylation has been reported to control the resistance and susceptibility to GCRV infection in CIK cells (Shang et al., 2017).In this study we examined the DNA methylome of CIK cells 8 hpi and found that the DNA methylation pattern of infected cells is very similar to the MOCK (Figure 3D).The statistically enriched biological processes of differentially methylated genes do not include defense response to virus or immune response (data not shown).These findings indicate that DNA methylation is less functionally involved in early GCRV infection.
Alternative polyadenylation functionally contributes to antiviral immune response (Jia et al., 2017).Some poly(A) binding proteins are the components of stress granules, the membrane-less ribonucleoprotein (RNP)-based cellular compartments in the cytoplasm triggering antiviral immune response.Moreover, alternative polyadenylation is involved in chronically infected disease (Su et al., .Previously, we have performed extensive alternative polyadenylation analysis to understand its functional relevance in tumorigenesis (Lai et al., 2015;Tan et al., 2018Tan et al., , 2021)).In this study we identified a considerable number of genes that selectively utilize alternative poly(A) sites in GCRV-infected CIK cells (Figures 4A,B).Among the genes with altered 3'UTRs we identified the biological processes of cytoskeleton organization, regulation of microtubule polymerization (Figures 4C,D).Interestingly, our recent study reported that microtubules are required for productive GCRV infection (Zhang et al., 2020), which is similar to MRV infection (Mainou et al., 2013).These observations suggest that alternative polyadenylation is potentially involved in GCRV early infection.
Taken together, our study provides evidence of molecular events during early infection of dsRNA viruses for understanding their pathogenesis.

FIGURE 1
FIGURE 1 Grass carp reovirus (GCRV)-induced cytopathic effect in CIK cells at different time points.CIK cells were mock-infected (Left panel) or infected with GCRV(Right panel), and the images were taken at 0, 4 and 8 hpi, respectively.Red arrow shows the representative CPE.Sale bar: 200 μm.

FIGURE 2
FIGURE 2Transcriptional programs in GCRV-infected CIK.(A) Enrichment plot for defense-and immune-related genes.(B) Volcano plot of statistically significant differentially expressed genes in early infection.(C, D).GO (C) and KEGG pathway (D) analyses of differentially expressed genes.

FIGURE 3 DNA
FIGURE 3 DNA methylome in early GCRV infection.(A) Genomic distribution of methylated regions in CIK cells.(B) DNA methylated signal around TSS regions.(C) Correlation between transcriptional signal and methylated regions.(D) DNA methylation heatmaps of MOCK and two biological replicates of 8 hpi group. 2