Temporal transcriptional profiling of host cells infected by a veterinary alphaherpesvirus using Nanopore sequencing

doi:10.21203/rs.3.rs-4601239/v1

Download PDF

Article

Temporal transcriptional profiling of host cells infected by a veterinary alphaherpesvirus using Nanopore sequencing

https://doi.org/10.21203/rs.3.rs-4601239/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

In our research, we performed temporal transcriptomic profiling of host cells infected with Equid alphaherpesvirus 1 by utilizing direct cDNA sequencing based on nanopore MinION technology. The sequencing reads were harnessed for transcript quantification at various time points. Viral infection-induced differential gene expression was identified through the edgeR package. The identified genes were segmented into six groups based on their kinetic characteristics. The initial three clusters encompass immediate-early response genes, typically transcription factors and elements of antiviral signaling pathways. These genes were either upregulated (cluster 1) or downregulated (clusters 2 and 3) during the early infection phase. The remaining three clusters include late response genes. In these categories, it is challenging to determine whether changes in gene expression are functionally linked to the viral infection or merely side effects of the infection. A study of gene associations using the STRINGDB software revealed several gene networks that might be directly impacted by the virus. Lastly, we explored whether gene co-expression could be a result of their collective regulation by upstream transcription factors using the Gene Regulatory Network database.

Biological sciences/Cell biology

Biological sciences/Microbiology

Equid alphaherpesvirus 1 (EHV-1) is a veterinary pathogen affecting horses and belongs to the Varicellovirus genus of alphaherpesviruses1,2. EHV-1 infection in horses is associated with myeloencephalopathy, a neurological condition, and its symptoms include upper respiratory tract disease, abortion, neonatal death, and fatal myeloencephalopathy3. EHV-1 contains a large, linear, double-stranded DNA genome of approximately 150 kbp in size. The genome's short region is flanked by inverted repeats4, and contains 76 protein-coding genes5. Five of these genes (ORF1, 2, 67, 71, and 75) have no homologs in other annotated alphaherpesvirus genomes6. EHV-1 enters cells through a receptor-mediated process7, and this can occur via endocytosis or fusion between the cell membrane and the viral envelope8. Besides productive infection, the virus can also enter a latent state in primary sensory neurons9.

Long-read sequencing (LRS) approaches have been instrumental in analyzing structural and dynamic aspects of transcriptomes. Despite short-read sequencing (SRS) offering high base accuracy and coverage, it falls short in identifying transcription start sites (TSSs), transcription end sites (TESs), splice sites, overlapping transcripts, and polygenic RNA molecules10. LRS techniques developed by Oxford Nanopore Technologies (ONT) and Pacific Biosciences (PacBio)11,12 can sequence full-length cDNA or native RNA (only ONT) sequences, at the cost of a higher sequencing error rate and lower throughput13–15. The ONT approach is especially affected by a high error rate, but this does not pose a significant problem in transcriptome analyses if well-annotated genomes are available. While conserved motifs of gene families create difficulty for SRS in identifying genes and RNA molecules, LRS does not encounter this issue. Moreover, the ONT technique enables the detection of RNA modifications without requiring base conversion16–18.

Several herpesvirus transcriptomes have been studied using different LRS techniques, including the PacBio synthesis-based sequencing platform and the nanopore sequencing from ONT. In this study, we utilized a time-course assay to examine the effect of EHV-1 infection on host cell gene expression. This transcriptome analysis was performed using direct cDNA sequencing, a non-amplification-based technique, on the ONT MinION platform.

Differential expression analysis of host genes

In this experiment, we scrutinized changes in the expression of host genes post EHV-1 infection over a 48-hour period (mock, 1, 2, 4, 6, 8, 12, 18, 24, and 48 hours post-infection) via direct cDNA sequencing utilizing the ONT MinION platform. From the raw gene counts, we pinpointed 7,308 genes, each having more than ten transcript reads in every one of the three replicates for different time point samples. Principal component analysis (PCA) demonstrates that technical variability between samples from each timepoint is minimal. Moreover, it highlights the biological pattern across timepoints through the markedly segmented distribution of samples between mock (0 hour sample, 0h) and 1 to 6 hours, contrasted by the 8 to 48 hours samples along the first two principal component axes (Fig. 1A). This apparent shift is also well in line with the exponential pattern of viral expansion, indicated by the relative ratio of host to viral read counts during infection (Fig. 1B). It is notable that by 48 hours, a relative rebound is apparent, however, host read counts are still dramatically low and thus might only represent slight technical variation. To investigate the hypothesis that the shift in kinetics might outline distinct host regulatory changes between early and late phases of infection, we have compared significantly differentially expressed genes (DEGs, with log2 fold-change > 1 and adjusted p-value < 0.05) between mock and 1h to 6h samples together (“early contrast”) and between 1h to 6h vs. 8h to 48h samples (“late contrast”). Notably, in the early contrast we find a wide range of significantly altered pathways by Gene Ontology (GO) term enrichment (Fig. 1C). These outline viral entry (“integrin-mediated signaling”, “LDL receptor activity”) and host defense responses, such as “negative regulation of apoptosis”, stress response and proteolysis, alongside several other processes indicating morphological changes and cell recruitment, cytokine stimuli. Interestingly, the late contrast outlines a narrowed-down set of processes dominated by transcription, translation and assembly (Fig. 1D). While this is to be expected, our approach helps to single out key molecular players indicating the host cell’s conversion into an effective viral replication machinery. We have compiled a resource of DEGs between subsequent timepoints, as well as compared to pre-infection, resulting in a list of 2533 unique genes regulated in one or more timepoints (Supplementary Table S1). Host gene expression correlated to the virus-host ratio changes has been determined via Pearson correlation (Supplementary Table S2). Interestingly, the most anti-correlated genes not only include the anti-viral candidate gene SAT1, but several 20S proteasome subunits, implicated in the antigen-presenting capacity of endothelial cells. Unsurprisingly, the most positively correlated genes include ribosome subunits, indicating a dramatic increase in translation capacity. These are, however, joined by an extensive set of long non-coding RNAs as well, which raises an interesting question on their possible role in the modulation of ribosomal biology.

Unsupervised clustering analysis identifies distinct host regulatory patterns

To complement the coarse-level comparison of early and late phases and better determine the directionality of specific alterations, a more refined clustering was carried out, identifying distinct gene expression patterns (Fig. 2). A narrowed-down, strict set of 577 genes exhibiting significant differential expression (FDR < 0.01) between the mock and 1–48 h time points was identified (Supplementary Table S3). We then converted the time series of the expression levels of DE genes to a relative scale, portraying the expression changes between the sampling time points, based on their maximum expression level. Subsequently, genes were clustered according to their scaled expression profiles, facilitating the recognition of genes that displayed co-expression (a similar expression pattern throughout the virus infection). This technique surpasses calculations with absolute expression levels of genes as there is a considerable difference between those. Utilizing this method, we were able to categorize six distinct groups of genes that demonstrated modified expression profiles during the course of infection (Fig. 2, Supplementary Table S4). In addition, Gene Ontology biological processes and GO molecular function annotations were used for functional characterization of separate clusters. The significant outcomes (FDR < 0.05) of the analysis for each cluster and GO terms are detailed in Supplementary Table S5.

Characterization of immediate-early and early response clusters

Cluster 1 was composed of 28 genes that presented low or no expression prior to the viral infection and exhibited a significant, transient up-regulation at 1h, which then significantly decreased at all subsequent time points. This cluster includes genes for transcription factors (ATF3, FOS, and HNF1B), cell-cycle regulators (CDK1), cytokines (SPP1), cell signaling (CTGF, CAP1), tubulin encoding (TUBB), heterochromatin-binding protein (CBX3), stressful growth arrest (GADD45A), translocation through the nuclear pore (KPNA2), cell adhesion (CLDN1), and components of the DNA replication complex (CIMP), among others. Cluster 2 (comprising 33 genes) and Cluster 3 (40 genes) contained genes that exhibited high expression before the viral infection, and their expression was sharply downregulated at 1h and further reduced at all subsequent time points. This downregulation was more pronounced for Cluster 2. Key inflammatory cytokines and acute-phase mediators (CXCL10, SAA3, THBS4) are down-regulated immediately upon infection, alongside TGFB2, low-density lipoprotein receptors, structural constituents and adhesion factors. Interestingly, no over-represented pathways were detected in cluster 1. However, for cluster 2, we identified an over-representation of genes located in the extracellular matrix structural constituent, structural molecule activity, amyloid-beta binding, heparin binding, integrin binding, protein-containing complex binding, cell adhesion molecule binding, and calcium ion binding GO molecular functions. This cluster also showcased a wide variety of GO biological processes, such as cell adhesion, extracellular matrix organization, multi-organ morphogenesis regulation, among others. In cluster 3, over-representation of genes implicated in protein folding in the endoplasmic reticulum, adenylate cyclase-modulating G protein-coupled receptor signaling pathway, steroid metabolic process, lipid metabolic process, platelet degranulation, secretion by cell, and export from cell GO Biological processes were detected. Furthermore, we explored the potential gene-gene associations within these clusters using STRINGDB¹⁹. This examination revealed multiple gene networks in each cluster (Fig. 3A-C). Besides the JUN/FOS and ATF3 immune-response related pathways, tubulin and heat shock protein interactions also feature prominently in the networks, suggesting stress response and viral entry-related activation. Importantly, the network analysis affirmed that co-expression was the primary interaction type among the genes of these clusters (Supplementary Tables S6, S7). To provide further context, building on the predicted virus-host protein interaction network of HSV-1 published in²⁰, we extracted the predicted interactions for significant DEGs in the immediate-early and early clusters (Fig. 3D). Again, the interaction of tubulin, spectrin and several cytoskeletal cross-linking genes, heat shock and related crystallines, alongside immune mediators is demonstrated, with several viral glycoproteins and ribonucleotide reductase.

Finally, to investigate whether the co-expression of genes in clusters 1–3 might be due to shared regulatory elements, we gathered all high-confidence associations between transcription factors and genes from the Gene Regulatory Network database (GRNdb; human embryo E-MTAB-3929, Supplementary Table S8). Our analysis discovered numerous transcription factors (TFs) shared among the genes of these clusters. Furthermore, cluster 1, which displayed a transient up-regulation, comprises three transcription factors (ATF3, FOS, and HNF1B). These transcription factors, as potential pivotal regulators, may control a significant number of downstream genes involved in antiviral defense.

Characterization of early-late and late response clusters

Cluster 4, containing 141 genes, demonstrated low expression in mock and showed an increasing expression throughout the course of the infection. This cluster correlates moderately with the virus-host ratio, however it is also dominated by ribosomal genes and translation/elongation facilitators. Cluster 5 comprised 129 genes that had low expression in the mock and showed a later response to viral infection that peaked at 48 hours. Accordingly, genes are highly positively correlated with the virus-host ratio and are dominated by mostly ribosomal genes and long non-coding RNAs, complementing cluster 4 and the genes found among “late contrast” DEGs. Even though the expression kinetics differed slightly between clusters 4 and 5, they shared several overlapping pathways, such as ribosomal RNA assembly, mRNA processes, transcription initiation, viral transcription, translation initiation, protein localization, protein transport, and others that account for the complete replication and packaging of the herpesvirus within the host cell. Lastly, Cluster 6, with 200 genes, had medium expression in the mock that slightly elevated until 6h, after which they were extensively down-regulated until the end of the infection. Importantly, this large cluster is dominated by cellular respiration, splicing and key life-supporting metabolic pathways, and follows the kinetic shift in viral expansion starting between 6- and 8-hours post-infection. We found pathways related to the host cell's energy and metabolic functions, responsible for cell energy supply, catabolic processes, and various signaling pathways. To a lesser extent, pathways in cluster 6 that overlapped with some pathways in clusters 4 and 5 (RNA binding, protein transporter activity) were also identified.

In this investigation, we employed a long-read sequencing method (ONT MinION, dcDNA-Seq) to explore the host cell's response to EHV-1 infection over a span of 48 hours. We utilized annotated sequencing reads to compute the abundance of transcripts at the specified time points. First, we inspected the major differences in gene expression between the early and late stages of infection by forming two coarse contrasts, with a separation point at 6 hours. This decision was informed by the principal component-derived distribution of samples. We then further categorized the identified genes into six clusters according to their scaled expression profiles. Each cluster includes genes demonstrating similar expression patterns throughout the entire period of viral infection we examined. We note that in the later stage of infection, it's challenging to discern if the changes in gene expressions are functional or if they simply signify transcriptional noise in cells undergoing disintegration. The most intriguing clusters are clusters 1–3, which contain genes that are either transiently upregulated (cluster 1) or downregulated (clusters 2 and 3) in the very early stages of infection.

Cluster 1 contains the transcription factors ATF3, FOS, and HNF1B, which may govern downstream genes involved in antiviral defense. The ATF/CREB transcription factor family, comprising nearly twenty members, embodies a collection of leucine zipper proteins²¹. The Activating Transcription Factor 3 (ATF3), encoded by the ATF3 gene, is a marker for cellular stress and regenerative response and has been seen to be induced upon infection by the Japanese encephalitis virus, regulating cellular antiviral pathways in the absence of type I interferons²². The c-fos nuclear phosphoprotein, encoded by the FOS gene, forms a complex with the JUN/AP-1 transcription factor, unlinked by covalent bonds. This immediate-early gene is initially expressed in response to various external stimuli, such as cytokines, growth factors, hormones, and stress²³, playing a crucial role in controlling cell proliferation, differentiation, and survival, as well as immune responses and more. The hepatocyte nuclear factor 1-beta protein, encoded by the HNF1B gene, has a broad range of functions, including suppressing the innate immune response²⁴. The exact role of this gene in defending against the virus remains unclear. On the other hand, marked upregulation of immediate-early genes is observed in secreted phosphoprotein 1 (osteopontin, SPP1), which strongly correlates with interferon gamma (IFNG) expression in a variety of inflammatory conditions and infections^25–27. However, during the course of infection, these usually well-correlated Th1 cytokines show alternating patterns in IE upregulation of SPP1 contrasted by the late upregulation of IFNG (Supplementary Figure S1).

The inflammatory cytokine CXCL10 is found in cluster 2 and is the only CXC family cytokine expressed in our model system. It is instrumental to initiating host response through its receptor CXCR3 and has been implicated in a variety of viral infections^28,29. Interestingly, here we find CXCL10 among the top down-regulated genes at the first hour of infection, suggesting the silencing of immune cell recruitment during early response, and the expression stays similarly suppressed throughout the time-course (Supplementary Figure S1).

TGFB2 shows IE downregulation in cluster 2, thus the widespread and context-dependent effects of TGFB signaling³⁰ appear to be suppressed upon infection. Importantly, the inhibition of TGFB signaling affected not only the lytic infection efficacy, but also the latency period of HSV-1 in mouse models³¹ CTGF, a main target of TGFB signaling also displays steady downregulation following a small peak in expression at 2h. Two of the main SMAD pathway elements, SMAD2 and SMAD4 are less robustly regulated during infection, with SMAD2 showing minor upregulation at late timepoints (Supplementary Figure S1).

Another representative of the immediate-early downregulation pattern is the low-density lipoprotein receptor (LDLR) gene. While LDLR has been known as a receptor for some RNA viruses³², it has recently been shown that LDLR expression effectively attenuates Gammaherpesviral infection in a cell-type specific manner³³, while it is downregulated in infected cells, a pattern closely replicated in our study as well (Supplementary Figure S1). In addition, the very low-density lipoprotein receptor (VLDLR) is known to serve as an evolutionarily conserved binding site for the E1 and E2 viral glycoproteins of a range of alphaherpesviruses³⁴, and shows a linearly decreasing trend during infection (Supplementary Figure S1). Interestingly, a further key inflammatory player that shows dramatic and lasting IE dowregulation is the apolipoprotein Serum amyloid A3 (SAA3), an acute phase reactant long known to be highly induced upon infection³⁵ and tested as a candidate marker of equine alphaherpesvirus infection³⁶. However, as our results show, this gene is also effectively silenced in infected cells (Supplementary Figure S1). Integrins are well known as herpesvirus receptors³⁷, and we observe the IE downregulation especially in ITGAV, known to bind the gH glycoprotein³⁸, but also ITGAX and ITGA3.

Late clusters, on the other hand, show the overrepresentation of ribosomal complexes alongside protein synthesis and assembly, while genes anti-correlated with the virus-host ratio represent declining mitochondrial function and putative anti-viral agents. One such agent, spermidine-spermine N1-acetyltransferase (SAT1), has been shown to deplete polyamines, thereby effectively reducing the replication of RNA viruses^39,40; however, they are less well characterized in Herpesviridae⁴¹.

Taken together, our high resolution time-course analysis of host response to EHV-1 infection offers a detailed resource on the significant differential gene expression, co-expressed gene clusters and inferred host and viral protein interactions, that may provide future targets and intervention opportunities in herpesvirus research.

Cells and viruses

In this study, we used equid alphaherpesvirus 1 (EHV-1) strain MdBio (EHV-1-MdBio)⁴². EHV-1 was propagated in confluent rabbit kidney (RK-13) epithelial cell line (ECACC: 00021715). RK-13 cells were cultivated in DMEM (Sigma), supplemented with 10% fetal calf serum (Gibco) and 80 μg of gentamycin per ml (Gibco) at 37 °C in the presence of 5% CO₂. Virus stock solution was prepared by infecting the freshly prepared cells with 0.1 multiplicity of infection [MOI = plaque-forming units (pfu)/cell]. Infected cells were harvested when complete cytopathic effect was observed. As a next step, three successive cycles of freezing and thawing of infected cells were performed to release of viruses from the infected cells. For the time-course sequencing experiments, the cells were infected with MOI=4 of EHV-1-MdBio in three technical replicates. For the synchronization of infections, cells were incubated for 1 h at 4 °C, followed by the removal of the virus suspension and washing the cells with phosphate-buffered saline. As a next step, fresh culture medium was added to the infected cells, which were incubated for 1, 2, 4, 6, 8, 12, 18 24 and 48 hours post-infection. After the incubation, the culture medium was removed, and the infected cells were frozen at -80 °C until further use.

RNA isolation

The cells underwent RNA extraction following the spin column-based procedure detailed in the NucleoSpin RNA kit from Macherey-Nagel. Initially, the cells were bathed in a lysis buffer with chaotropic ions, a step which served to neutralize RNases. This was followed by the adhesion of DNA and RNA molecules to the silica membrane. All samples were subsequently treated with DNase I to eradicate any lingering genomic (g)DNA. The total RNA was then drawn out into nuclease-free water. Any potentially persistent gDNA was removed with the TURBO DNA-free™ Kit from Invitrogen. The concentration of the samples was ascertained using the Qubit 4.0 fluorometer and Qubit Broad Range RNA Assay Kit, both from Invitrogen. The Agilent TapeStation 4150 was employed for quality control, and only those samples with RIN scores of 9.2 or higher were selected for cDNA production for subsequent procedures.

Purification of polyadenylated RNA

The poly(A)⁺ RNA fraction was isolated from the total RNA samples using the Oligotex mRNA Mini Kit from Qiagen. Briefly, the samples' volumes were adjusted to 250 µL using RNase-free water. The samples were then treated with 15 µL of Oligotex suspension and 250 µL of OBB buffer, both supplied in the Qiagen kit. The samples were subsequently heated to 70°C for a duration of 3 minutes, before being cooled to 25°C for 10 minutes. Following this, the mixtures were spun down at 14,000×g for 2 minutes, after which the supernatants were disposed of. A quantity of 400 µL of OW2 wash buffer, also from the Oligotex kit, was added to the samples before they were transferred onto the spin columns from the Qiagen kit. They were then spun at 14,000×g for 1 minute. This wash step was repeated once more. Ultimately, the polyadenylated RNA fraction was collected from the membrane by introducing 50 µl of heated elution buffer (EB, Qiagen kit). The RNA was collected in 60 µl EB, and a secondary elution step was also performed to ensure maximum yield.

ONT MinION – dcDNA sequencing

Both poly(A)+-enriched and poly(A)+-enriched samples processed with Terminator were employed to construct libraries for direct cDNA sequencing on the ONT MinION device. The ONT Direct cDNA Sequencing Kit (SQK-DCS109) was used according to the instructions provided in the kit’s manual. Briefly, the RNA samples were combined with the VN primer (VNP; from the ONT kit) and 10 mM dNTPs, then heated to 65°C for 5 minutes. Following this, 5x RT Buffer, RNaseOUT (from Thermo Fisher Scientific), and Strand-Switching Primer (SSP; from the ONT Kit) were incorporated, with the mixtures then heated to 42°C for 2 minutes. The first cDNA strand was produced by adding Maxima H Minus Reverse Transcriptase enzyme (from Thermo Fisher Scientific) to the samples. RT and strand-switching reactions took place at 42°C for 90 minutes, before the enzyme was heat inactivated at 85°C for 5 minutes. The RNase Cocktail Enzyme Mix (from Thermo Fisher Scientific) was used to remove RNA from the RNA-cDNA hybrids, with this step carried out at 37°C for 10 minutes.

For the synthesis of the second strand, LongAmp Taq Master Mix [from New England Biolabs (NEB)] and PR2 Primer (PR2P) were utilized. The specific details of the PCR reactions can be found in the provided reference ⁴². The fragmented DNAs underwent end-repair and dA-tailing using the NEBNext End repair /dA-tailing Module (NEB) at 20°C for 5 minutes and then 65°C for 5 minutes. This was succeeded by the ligation of the sequencing adapter at room temperature for 10 minutes, with the NEB Blunt /TA Ligase Master Mix (NEB) used for this process. The ONT dcDNA libraries were marked using barcodes as detailed in the provided reference ⁴² and using the ONT Native Barcoding (12) Kit following the manufacturer's instructions. The prepared and tethered cDNA libraries (200 fmol/flow cell) were purified and introduced into ONT R9.4.1 SpotON Flow Cells. In total, five flow cells were used for the dcDNA sequencing.

In order to circumvent possible "barcode hopping," the earlier time point samples were sequenced separately from those taken at later time points. AMPure XP Beads were utilized after each additional enzymatic process. The samples were then eluted in UltraPure™ nuclease-free water (from Invitrogen), and their concentration was determined with a Qubit 4.0 fluorometer, in conjunction with the Qubit dsDNA HS Assay kit.

Pre-processing and data analysis

The MinION data underwent base calling using the Guppy base caller version 3.4.1, with the --qscore_filtering feature enabled. Reads achieving a Q-score of more than 7 were mapped onto a combined host genome, which comprised the Oryctolagus cuniculus (rabbit) GCF_000003625.3 and Equid alphaherpesvirus 1 (EHV-1) NC_001491.2 reference genomes. This was accomplished using the minimap2 software, with the “-ax splice -Y -C5” options. We chose to exclude MAPQ=0, secondary, and supplementary alignments from all subsequent analysis. Reads aligned to the host genome were associated with host genes based on the GCF_000003625.3_OryCun2.0_genomic.gff genome coordinates. We only used reads that corresponded with the exon structure of the host reference genes (employing a +/- 5 base pair window for matching exon start and end positions) for the computation of raw gene counts. A summary of mapped reads, virus/host ratio, read length/quality for each time point and experiment can be found in Supplementary Table S9.

Differential expression analysis of host genes

EdgeR_3.24.3 29 along with R version 3.5.1 was utilized for the differential expression (DE) analysis. We eliminated host genes that had less than ten reads across any of the three technical replicates, retaining 7308 genes above our chosen threshold (Supplementary Table S10). The three technical replicates were then categorized in EdgeR (utilizing the DGEList function) according to their respective time points (mock, 1h, 2h, 4h, 6h, 8h, 12h, 18h, 24h, 48h). We carried out data normalization (with the calcNormFactors function) using the method="TMM" option, and the robust=True option in the subsequent analysis (with the estimateDisp, glmQLFit, and glmQLFTest functions). We sought to identify key genes involved in the earliest response to viral infection in the host by testing for DE between the mock and the 1h time points. To identify genes that exhibited significantly altered expression levels (utilizing the decideTests function), we set a false discovery rate (FDR) threshold of 0.01, adjusting p-values using the Benjamini & Hochberg procedure (with the p.adjust function and method="BH" option, Hochberg and Benjamini, 1990). In addition, PCA was run using the prcomp R function on the normalized count matrix. The first two principal components were used for a coarse DEG analysis between the mean gene expression of early (mock vs. 1h-6h timepoints) and late (between mean values of 1h-6h vs. 8h-48h samples). Functional enrichment of the resulting DEG lists was done in Cytoscape using the ClueGO plugin (Reactome Pathway database, significant hits with adjusted p-value < 0.05) and manually edited for readability. The normalized pseudo-counts of DE genes, used for clustering and visualization demonstrated superb reproducibility across the three technical replicates (an average R=0.994 in all cross correlations between the three replicates and various time points), as reinforced by the first two PCA coordinates. Supplementary table S1 gathers significant DEGs across all timepoints and in reference to mock samples, and may be easily browsed and filtered according to the “contrast” column. Furthermore, Pearson correlation between normalized gene counts in the virus-host ratio across time points was determined by the corr.test function of the “psych” R package, and significantly correlated genes collected in Supplementary Table S2.

Examination of host genes exhibiting differential expression during viral infection

In order to categorize genes based on changes in their relative expression (rather than absolute expression levels) throughout the entire viral infection, we adjusted the expression of DE genes in relation to their peak expression levels. We executed cluster analysis utilizing the amap_0.8-16 R package Kmeans function with the Euclidean distance method to pinpoint genes with similar expression kinetics during viral infection. This resulted in the identification of six clusters of genes with unique expression profiles (Supplementary Table S4). We visualized the gene expression heatmap of these clusters using the pheatmap R package (version 1.0.12). The median of the relative gene expressions for each time point within each of the identified gene clusters was visualized with the ggplot2 R package (version 2.3.3.2) using the geom_smooth function.

Employing the identified subset of genes, we carried out an over-representation analysis for each cluster using PANTHER (version 16.0, using the 2020-12-01 data set release; https://doi.org/10.1093/nar/gkaa1106) software tool. The reference was the 7308 genes that were above the selected expression threshold. We compiled the results of the over-representation analysis (FDR<0.05) using the Gene Ontology (GO) biological processes and GO molecular functions annotation datasets.

We carried out a STRINGDB network analysis (version 11.5) to uncover possible gene networks in clusters 1-3, where the over-representation analysis revealed little or no enrichment in known pathways. Within these specific gene clusters, we also pinpointed the shared regulatory elements of our DE genes based on the Gene Regulatory Network database (GRNdb) data (Fang et al. 2021). For this analysis, we used the Human Embryo (E-MTAB-3929) transcription factor (TF) database, which consists of 1529 single-cell experiments from samples across five different embryonic stages, resulting in 170256 TF-gene associations between 18355 target genes and 695 TFs. We only included the high-confidence TF-gene pairs (annotated in the original SCENIC cisTarget database or inferred by orthology (Aibar et al., 2017) in our investigation.

Ethics declaration: Not applicable

Author Contributions: Conceptualization, D.T., Z.M. and Z.B.; Methodology, D.T., Z.M. B.K., P.O. and Z.B.; Validation, D.T., T.K. Á.D., Z.C. and Z.B.; Formal analysis, G.T., G.G., D.T., B.K., Á.D., P.O. and Z.B.; Investigation, D.T., Z.C., Á.D.; Writing the original draft, D.T., Z.M. and Z.B.; Writing, review and editing, Z.B.; Visualization, G.T., G.G., B.K., P.O. and D.T; Supervision, Z.M. and Z.B.; Funding acquisition, D.T., and Z.B. All authors have read and agreed to the published version of the manuscript.

Funding: This study was supported by the National Research, Development and Innovation Office grant: K 128247 to ZB and FK 128252 to DT. The APC charge was covered by the University of Szeged Open Access Fund: 6997. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Conflicts of Interest: The authors declare no conflict of interest.

Data availability: The sequencing datasets generated in this study are available at the European Nucleotide Archive under the accession: PRJEB52190.

Granoff, A. & Webster, R. G. Encyclopedia of Virology. (Elsevier, 1999).
Oladunni, F. S. et al. Equid Herpesvirus 1 Targets the Sensitization and Induction Steps To Inhibit the Type I Interferon Response in Equine Endothelial Cells. J. Virol. 93, (2019).
Carroll, C. L. & Westbury, H. A. Isolation of equine herpesvirus 1 from the brain of a horse affected with paresis. Aust. Vet. J. 62, 345–346 (1985).
Roizmann, B. et al. The family Herpesviridae: an update. The Herpesvirus Study Group of the International Committee on Taxonomy of Viruses. Arch. Virol. 123, 425–449 (1992).
Telford, E. A., Watson, M. S., McBride, K. & Davison, A. J. The DNA sequence of equine herpesvirus-1. Virology 189, 304–316 (1992).
Reed, S. M. & Toribio, R. E. Equine herpesvirus 1 and 4. Vet. Clin. North Am. Equine Pract. 20, 631–642 (2004).
Osterrieder, N. Construction and characterization of an equine herpesvirus 1 glycoprotein C negative mutant. Virus Res. 59, 165–177 (1999).
Frampton, A. R. J. et al. Equine herpesvirus 1 enters cells by two different pathways, and infection requires the activation of the cellular kinase ROCK1. J. Virol. 81, 10879–10889 (2007).
Warda, F. F. et al. Application of equine herpesvirus-1 vaccine inactivated by both formaldehyde and binary ethylenimine in equine. Vet. World 14, 1815–1821 (2021).
Boldogkői, Z., Moldován, N., Balázs, Z., Snyder, M. & Tombácz, D. Long-Read Sequencing - A Powerful Tool in Viral Transcriptome Research. Trends Microbiol. 27, 578–592 (2019).
Workman, R. E. et al. Nanopore native RNA sequencing of a human poly(A) transcriptome. Nat. Methods 16, 1297–1305 (2019).
Marx, V. Method of the year: long-read sequencing. Nat. Methods 20, 6–11 (2023).
Laver, T. et al. Assessing the performance of the Oxford Nanopore Technologies MinION. Biomol. Detect. Quantif. 3, 1–8 (2015).
Irimia, M. et al. A highly conserved program of neuronal microexons is misregulated in autistic brains. Cell 159, 1511–1523 (2014).
Tombácz, D. et al. Full-Length Isoform Sequencing Reveals Novel Transcripts and Substantial Transcriptional Overlaps in a Herpesvirus. PloS One 11, e0162868 (2016).
Balázs, Z. et al. Template-switching artifacts resemble alternative polyadenylation. BMC Genomics 20, 824 (2019).
Luo, G. X. & Taylor, J. Template switching by reverse transcriptase during DNA synthesis. J. Virol. 64, 4321–4328 (1990).
Cocquet, J., Chong, A., Zhang, G. & Veitia, R. A. Reverse transcriptase template switching and false alternative transcripts. Genomics 88, 127–131 (2006).
Szklarczyk, D. et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 47, D607–D613 (2019).
Lian, X. et al. Prediction and analysis of human-herpes simplex virus type 1 protein-protein interactions by integrating multiple methods. Quant. Biol. 8, 312–324 (2020).
Hai, T. The ATF Transcription Factors in Cellular Adaptive Responses. in Gene Expression and Regulation (ed. Ma, J.) 329–340 (Springer New York, New York, NY, 2006). doi:10.1007/978-0-387-40049-5_20.
Sood, V. et al. ATF3 negatively regulates cellular antiviral signaling and autophagy in the absence of type I interferons. Sci. Rep. 7, 8789 (2017).
Healy, S., Khan, P. & Davie, J. R. Immediate early response genes and cell transformation. Pharmacol. Ther. 137, 64–77 (2013).
He, J. et al. Hepatocyte nuclear factor 1A suppresses innate immune response by inducing degradation of TBK1 to inhibit steatohepatitis. Genes Dis. 10, 1596–1612 (2023).
Kovács, M. Á. et al. Elevated Osteopontin and Interferon Gamma Serum Levels and Increased Neutrophil-to-Lymphocyte Ratio Are Associated With the Severity of Symptoms in Schizophrenia. Front. Psychiatry 10, 996 (2019).
Lund, S. A., Giachelli, C. M. & Scatena, M. The role of osteopontin in inflammatory processes. J. Cell Commun. Signal. 3, 311–322 (2009).
Park, J. W. et al. Osteopontin depletion decreases inflammation and gastric epithelial proliferation during Helicobacter pylori infection in mice. Lab. Investig. J. Tech. Methods Pathol. 95, 660–671 (2015).
Wuest, T. R. & Carr, D. J. J. Dysregulation of CXCR3 signaling due to CXCL10 deficiency impairs the antiviral response to herpes simplex virus 1 infection. J. Immunol. Baltim. Md 1950 181, 7985–7993 (2008).
Elemam, N. M., Talaat, I. M. & Maghazachi, A. A. CXCL10 Chemokine: A Critical Player in RNA and DNA Viral Infections. Viruses 14, (2022).
Massagué, J. TGFβ signalling in context. Nat. Rev. Mol. Cell Biol. 13, 616–630 (2012).
Allen, S. J. et al. Adaptive and innate transforming growth factor beta signaling impact herpes simplex virus 1 latency and reactivation. J. Virol. 85, 11448–11456 (2011).
Finkelshtein, D., Werman, A., Novick, D., Barak, S. & Rubinstein, M. LDL receptor and its family members serve as the cellular receptors for vesicular stomatitis virus. Proc. Natl. Acad. Sci. U. S. A. 110, 7306–7311 (2013).
Aurubin, C. A., Knaack, D. A., Sahoo, D. & Tarakanova, V. L. Low-Density Lipoprotein Receptor Suppresses the Endogenous Cholesterol Synthesis Pathway To Oppose Gammaherpesvirus Replication in Primary Macrophages. J. Virol. 95, e0064921 (2021).
Clark, L. E. et al. VLDLR and ApoER2 are receptors for multiple alphaviruses. Nature 602, 475–480 (2022).
Mitchell, T. I., Jeffrey, J. J., Palmiter, R. D. & Brinckerhoff, C. E. The acute phase reactant serum amyloid A (SAA3) is a novel substrate for degradation by the metalloproteinases collagenase and stromelysin. Biochim. Biophys. Acta 1156, 245–254 (1993).
Pusterla, N. et al. Investigation of The Usefulness of Serum Amyloid A in Characterizing Selected Disease Forms of Equine Herpesvirus-1 Infection. J. Equine Vet. Sci. 104, 103699 (2021).
Campadelli-Fiume, G., Collins-McMillen, D., Gianni, T. & Yurochko, A. D. Integrins as Herpesvirus Receptors and Mediators of the Host Signalosome. Annu. Rev. Virol. 3, 215–236 (2016).
Cheshenko, N. et al. Herpes simplex virus type 2 glycoprotein H interacts with integrin αvβ3 to facilitate viral entry and calcium signaling in human genital tract epithelial cells. J. Virol. 88, 10026–10038 (2014).
Pozzi, B. et al. Dengue virus targets RBM10 deregulating host cell splicing and innate immune response. Nucleic Acids Res. 48, 6824–6838 (2020).
Mounce, B. C. et al. Interferon-Induced Spermidine-Spermine Acetyltransferase and Polyamine Depletion Restrict Zika and Chikungunya Viruses. Cell Host Microbe 20, 167–177 (2016).
Firpo, M. R. & Mounce, B. C. Diverse Functions of Polyamines in Virus Infection. Biomolecules 10, (2020).
Tombácz, D. et al. Hybrid sequencing discloses unique aspects of the transcriptomic architecture in equid alphaherpesvirus 1. Heliyon 9, e17716 (2023).

No competing interests reported.

Supplementaryinformation.zip

Download PDF

Editorial decision: Revision requested
30 Jul, 2024
Reviews received at journal
10 Jul, 2024
Reviews received at journal
09 Jul, 2024
Reviewers agreed at journal
09 Jul, 2024
Reviewers agreed at journal
26 Jun, 2024
Reviewers invited by journal
23 Jun, 2024
Editor assigned by journal
23 Jun, 2024
Editor invited by journal
22 Jun, 2024
Submission checks completed at journal
21 Jun, 2024
First submitted to journal
18 Jun, 2024

You are reading this latest preprint version

Temporal transcriptional profiling of host cells infected by a veterinary alphaherpesvirus using Nanopore sequencing

Status:

Version 1

Abstract

Figures

Introduction

Results

Differential expression analysis of host genes

Unsupervised clustering analysis identifies distinct host regulatory patterns

Characterization of immediate-early and early response clusters

Characterization of early-late and late response clusters

Discussion

Materials and methods

Cells and viruses

RNA isolation

Purification of polyadenylated RNA

ONT MinION – dcDNA sequencing

Pre-processing and data analysis

Differential expression analysis of host genes

Examination of host genes exhibiting differential expression during viral infection

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1