Integrative RNA profiling of TBEV-infected neurons and astrocytes reveals potential pathogenic effectors

Abstract

Tick-borne encephalitis virus (TBEV), the most medically relevant tick-transmitted flavivirus in Eurasia, targets the host central nervous system and frequently causes severe encephalitis. The severity of TBEV-induced neuropathogenesis is highly cell-type specific and the exact mechanism responsible for such differences has not been fully described yet. Thus, we performed a comprehensive analysis of alterations in host poly-(A)/miRNA/lncRNA expression upon TBEV infection in vitro in human primary neurons (high cytopathic effect) and astrocytes (low cytopathic effect). Infection with severe but not mild TBEV strain resulted in a high neuronal death rate. In comparison, infection with either of TBEV strains in human astrocytes did not. Differential expression and splicing analyses with an in silico prediction of miRNA/mRNA/lncRNA/vd-sRNA networks found significant changes in inflammatory and immune response pathways, nervous system development and regulation of mitosis in TBEV Hypr-infected neurons. Candidate mechanisms responsible for the aforementioned phenomena include specific regulation of host mRNA levels via differentially expressed miRNAs/lncRNAs or vd-sRNAs mimicking endogenous miRNAs and virus-driven modulation of host pre-mRNA splicing. We suggest that these factors are responsible for the observed differences in the virulence manifestation of both TBEV strains in different cell lines. This work brings the first complex overview of alterations in the transcriptome of human astrocytes and neurons during the infection by two TBEV strains of different virulence. The resulting data could serve as a starting point for further studies dealing with the mechanism of TBEV-host interactions and the related processes of TBEV pathogenesis.

Keywords: A3SS, alternative 3′ splice site; A5SS, alternative 5′ splice site; ACACA, Acetyl-CoA Carboxylase Alpha; AKR1C2, Aldo-Keto Reductase Family 1 Member C2; ANKS1A, Ankyrin Repeat And Sterile Alpha Motif Domain Containing 1A; ANOS1, Anosmin 1; AOX1, Aldehyde Oxidase 1; APOBEC3G, Apolipoprotein B MRNA Editing Enzyme Catalytic Subunit 3G; APOL1/6, Apolipoprotein L1/6; ARID2, AT-Rich Interaction Domain 2; AUTS2, Activator Of Transcription And Developmental Regulator AUTS2; Alternative splicing; Astrocytes; BCL11B, BAF Chromatin Remodeling Complex Subunit BCL11B; BCL9L, BCL9 Transcription Coactivator-like; BDKRB2, Bradykinin Receptor B2; BDNF, Brain Derived Neurotrophic Factor; BEND3, BEN Domain Containing 3; BSA, bovine serum albumin; BST2, Bone Marrow Stromal Cell Antigen 2; CALB1, Calbindin 1; CAMK2A, Calcium/Calmodulin Dependent Protein Kinase II Alpha; CD, complement determinant; CDKN1C, Cyclin Dependent Kinase Inhibitor 1C; CFAP61, Cilia And Flagella Associated Protein 61; CHRNA3, Cholinergic Receptor Nicotinic Alpha 3 Subunit; CHRNB4, Cholinergic Receptor Nicotinic Beta 4 Subunit; CLIC5, Chloride Intracellular Channel 5; CMPK2, Cytidine/Uridine Monophosphate Kinase 2; CNS, central nervous system; CNTN2, Contactin 2; CREG2, Cellular Repressor Of E1A Stimulated Genes 2; CXADR, Coxsackievirus B-Adenovirus Receptor; CYYR1, Cysteine And Tyrosine Rich 1; DACH1, Dachshund Family Transcription Factor 1; DAPI, diamidino-2-phenylindole; DCC, Netrin 1 Receptor; DCX, Doublecortin; DDX60, DExD/H-Box Helicase 60; DDX60L, DExD/H-Box 60 Like; DE, differentially expressed; DENV, Dengue virus; DIRAS2, DIRAS Family GTPase 2; DLX1/5/6, Distal-Less Homeobox 1/5/6; DNMT3B, DNA Methyltransferase 3 Beta; DPYSL2, Dihydropyrimidinase Like 2; EBF1, EBF Transcription Factor 1; EGF, Epidermal Growth Factor; ELAVL2/4, ELAV Like RNA Binding Protein 2/4; EPHB1, EPH Receptor B1; EPSTI1, Epithelial Stromal Interaction 1; ERBB4, Erb-B2 Receptor Tyrosine Kinase 4; ES, exon skipping; ESRRG, Estrogen Related Receptor Gamma; FGFb, Fibroblast Growth Factor 2; FPKM, Fragments Per Kilobase of transcript per Million mapped reads; FUT9, Fucosyltransferase 9; G2E3, G2/M−Phase Specific E3 Ubiquitin Protein Ligase; GABRG2, Gamma-Aminobutyric Acid Type A Receptor Subunit Gamma 2; GAPDH, Glyceraldehyde-3-Phosphate Dehydrogenase; GAS2L3, Growth Arrest Specific 2 Like 3; GAS7, Growth Arrest Specific 7; GATAD2B, GATA Zinc Finger Domain Containing 2B; GFAP, Glial Fibrillary Acidic Protein; GIPC2, GIPC PDZ Domain Containing Family Member 2; GLRA2, Glycine Receptor Alpha 2; GNG2, G Protein Subunit Gamma 2; GO, gene ontology; GOLGA4, Golgin A4; GRIN2A, Glutamate Ionotropic Receptor NMDA Type Subunit 2A; GSEA, gene set enrichment analysis; HERC5/6, HECT And RLD Domain Containing E3 Ubiquitin Protein Ligase 5/6; HEYL, Hes Related Family BHLH Transcription Factor With YRPW Motif Like; HPRT1, Hypoxanthine Phosphoribosyltransferase 1; HS, hot-spot; HSPA6, Heat Shock Protein Family A (Hsp70) Member 6; HUDD (ELAV4), Hu-Antigen D/ELAV Like Neuron-Specific RNA Binding Protein 4; IFI6, Interferon Alpha Inducible Protein 6; IFIH1 (MDA5), Interferon Induced With Helicase C Domain 1/Melanoma Differentiation-Associated Protein 5; IFIT1-3, Interferon Induced Protein With Tetratricopeptide Repeats 1–3; IFITM1/2, Interferon Induced Transmembrane Protein 1/2; IFN, interferon; IGB, Integrated Genome Browser; IL6, Interleukin 6; IR, intron retention; ISG20, Interferon Stimulated Exonuclease Gene 20; ISGF3, Interferon-Stimulated Gene Factor 3 Gamma; ISGs, interferon-stimulated genes; JEV, Japanese encephalitis virus; KCND2, Potassium Voltage-Gated Channel Subfamily D Member 2; KCNK10, Potassium Two Pore Domain Channel Subfamily K Member 10; KCNS2, Potassium Voltage-Gated Channel Modifier Subfamily S Member 2; KIT, KIT Proto-Oncogene, Receptor Tyrosine Kinase; KLHDC8A, Kelch Domain Containing 8A; KLHL13, Kelch Like Family Member 13; KRR1, KRR1 Small Subunit Processome Component Homolog; LCOR, Ligand Dependent Nuclear Receptor Corepressor; LEKR1, Leucine, Glutamate And Lysine Rich 1; LGI1, Leucine Rich Glioma Inactivated 1; LRRTM3, Leucine Rich Repeat Transmembrane Neuronal 3; LSV, local splicing variation; LUZP2, Leucine Zipper Protein 2; MAN1A1, Mannosidase Alpha Class 1A Member 1; MAP2, Microtubule Associated Protein 2; MBNL2, Muscleblind Like Splicing Regulator 2; MCTP1, Multiple C2 And Transmembrane Domain Containing 1; MMP13, Matrix Metallopeptidase 13; MN1, MN1 Proto-Oncogene, Transcriptional Regulator; MOI, multiplicity of infection; MTUS2, Microtubule Associated Scaffold Protein 2; MX2, MX Dynamin Like GTPase 2; MYCN, MYCN Proto-Oncogene, BHLH Transcription Factor; NAV1, Neuron Navigator 1; NCAM1, Neural Cell Adhesion Molecule 1; NDRG4, N-Myc Downstream-Regulated Gene 4 Protein; NEK7, NIMA Related Kinase 7; NFASC, Neurofascin; NKAIN1, Sodium/Potassium Transporting ATPase Interacting 1; NMI, N-Myc And STAT Interactor 2; NRAP, Nebulin Related Anchoring Protein; NRARP, NOTCH Regulated Ankyrin Repeat Protein; NREP, Neuronal Regeneration Related Protein; NRN1, Neuritin 1; NS3, flaviviral non-structural protein 3; NXPH2, Neurexophilin 2; NYNRIN, NYN Domain And Retroviral Integrase Containing; Neurons; Neuropathogenesis; OAS, 2′-5′-Oligoadenylate Synthetase; OASL, 2′-5′-Oligoadenylate Synthetase Like; ONECUT2, ONECUT-2 Homeodomain Transcription Factor; OPCML, Opioid Binding Protein/Cell Adhesion Molecule Like; OTX2, Orthodenticle Homeobox 2; PBS, phosphate buffer saline; PBX1, Pre-B-Cell Leukemia Transcription Factor 1; PCDH18/20, Protocadherin 18/20; PFKFB3, 6-Phosphofructo-2-Kinase/Fructose-2,6-Biphosphatase 3; PIK3C2B, Phosphatidylinositol-4-Phosphate 3-Kinase Catalytic Subunit Type 2 Beta; PIP4P2, Phosphatidylinositol-4,5-Bisphosphate 4-Phosphatase 2; PLCH1, Phospholipase C Eta 1; POU3F4, Brain-Specific Homeobox/POU Domain Protein 4; PPM1L, Protein Phosphatase, Mg2+/Mn2+ Dependent 1L; PPP1R17, Protein Phosphatase 1 Regulatory Subunit 17; PRDM12, PR Domain Zinc Finger Protein 12; PSI, percent selective index; PSRC1, Proline And Serine Rich Coiled-Coil 1; PTPN5, Protein Tyrosine Phosphatase Non-Receptor Type 5; PTPRH, Protein Tyrosine Phosphatase Receptor Type H; RAPGEF5, Rap Guanine Nucleotide Exchange Factor 5; RBFOX1, RNA Binding Fox-1 Homolog 1; RIG-I (DDX58), Retinoic Acid-Inducible Gene 1 Protein; RNF212, Ring Finger Protein 212; RNVU1, RNA, Variant U1 Small Nuclear; RSAD2, Radical S-Adenosyl Methionine Domain Containing 2; RTL8B, Retrotransposon Gag Like 8B; Response to infection; SAMD9, Sterile Alpha Motif Domain Containing 9; SEMA3E, Semaphorin 3E; SH3TC2, SH3 Domain And Tetratricopeptide Repeats 2; SHF, Src Homology 2 Domain Containing F; SHISAL1, Shisa Like 1; SIAH3, Siah E3 Ubiquitin Protein Ligase Family Member 3; SIRPA, Signal Regulatory Protein Alpha; SLITRK5, SLIT And NTRK Like Family Member 5; SNP, single-nucleotide polymorphism; SOGA1, Suppressor Of Glucose, Autophagy Associated 1; SPSB4, SplA/Ryanodine Receptor Domain And SOCS Box Containing 4; ST6GAL1, ST6 Beta-Galactoside Alpha-2,6-Sialyltransferase 1; TBC1D30, TBC1 Domain Family Member 30; TBEV, Tick-borne encephalitis virus; TFAP2A, Transcription Factor AP-2 Alpha; TFAP2B, Transcription Factor AP-2 Beta; THSD7A, Thrombospondin Type 1 Domain Containing 7A; THUMPD2, THUMP Domain-Containing Protein 2/SAM-Dependent Methyltransferase; TIPARP, TCDD Inducible Poly(ADP-Ribose) Polymerase; TM4SF18, Transmembrane 4 L Six Family Member 18; TMC8, Transmembrane Channel Like 6; TMEM229B, Transmembrane Protein 229B; TMTC1, Transmembrane O-Mannosyltransferase Targeting Cadherins 1; TNFSF10, TNF Superfamily Member 10; TRHDE, Thyrotropin Releasing Hormone Degrading Enzyme; TRIM38, Tripartite Motif Containing 38; TSHZ1, Teashirt Zinc Finger Homeobox 1; Tick-borne encephalitis virus; Transcriptomics; USP18, Ubiquitin Specific Peptidase 18/ISG15-Specific-Processing Protease; UTR, untranslated region; UTS2R, Urotensin 2 Receptor; WNV, West Nile virus; XAF1, XIAP Associated Factor 1; XRN1, 5′-3′ Exoribonuclease 1; ZIKV, Zika virus; ZMAT3, Zinc Finger Matrin-Type 3; ZMYM5, Zinc Finger MYM-Type Containing 5; ZNF124, Zinc Finger Protein 124; ZNF730, Zinc Finger Protein 730; gRNA, genomic TBEV RNA; hNSC, human neural stem cells; lncRNA, long non-coding RNA; mRNA, messenger RNA; miRNA; miRNA, micro RNA; ncRNA, non-coding RNA; pc-mRNA, protein-coding mRNA; qRT-PCR, quantitative reverse transcription real-time PCR; snRNP, small nuclear ribonucleoproteins; vd-sRNA, virus-derived small RNA.

Graphical abstract

Fig. 1

TBEV infection dynamics assessment in in vitro differentiated human neurons and astrocytes reveals higher neuronal susceptibility. Astrocytes were differentiated for 22 days, neurons for 14–17 days prior to infection with TBEV strains Hypr and Neudoerfl at MOI of 5. (A) Immunofluorescent labelling of selected differentiation markers in astrocytes (GFAP – red and S100B – green) – left and neurons (HuD – red and MAP2 -green) – right. Nuclei were stained with DAPI. Controls are labelled with secondary antibodies only. The scale bar represents 20 µm. Representative images from two independent experiments are shown. (B) TBEV growth curve of strains Hypr and Neudoerfl, virus shedding into cultivation media was quantified by plaque assay using A549 cells. Data summarise three independent experiments and values in graphs are expressed as mean ± SD. Significant difference from control was calculated by unpaired Student’s t-test (p < 0.05 (*), p < 0.01(**)). (C) Relative viability of neurons and astrocytes upon TBEV infection measured by metabolic conversion of alamarBlue and related to the viability of mock-treated cells (1.0 = 100 %, dotted line). Data summarise four biological and two technical replicates experiments, and values in graphs are expressed as mean ± SD. Welch́s t-test showed a significant difference from control with the multiple testing correction by FDR (Benjamini, Krieger, and Yekutieli, q ≤ 0.05), p ≤ 0.05 (*), p ≤ 0.01 (**), p ≤ 0.001 (***). (D) TBEV replication in astrocytes and neurons was assessed by the amount of genomic RNA quantified by qRT-PCR and related to mock with the ΔΔ-ct method. Results are represented as means ± SD of three biological and technical replicates. (E) TBEV C (capsid protein) and GAPDH (loading control) levels in infected astrocytes and neurons at 24 and 72 h p.i. were determined by immunoblotting. Representative data from three independent experiments are shown. N – Neudoerfl strain, H – Hypr strain. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2

Schematic representation of an experimental design performed in this study. (A) hNSCs were differentiated into astrocytes and neurons, infected with TBEV strains Hypr and Neudoerfl at the MOI of 5 and total RNA was isolated at 24 and 72 h p.i. with RNA Blue reagent. The quality of RNA was verified using Bioanalyzer 2100 and (B) Diagram describing the generation of miRNA/mRNA/lncRNA/vd-sRNA networks. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3

Differential expression analysis of poly-(A) enriched RNAs in TBEV-infected neurons and astrocytes reveals strain-specific and cell-specific responses. (A) Graphic summary of poly-(A) classes identified as differentially expressed upon TBEV infection. (B) Graphic summary of DE genes (protein-coding mRNAs and non-coding RNAs) in the respective dataset. (C) Heatmap visualizing the overall expression pattern of protein-coding mRNAs in analysed samples. (D) Heatmap visualizing the overall expression pattern of non-coding RNAs in analysed samples. (E) Verification of RNA-seq data by qRT-PCR. Relative expressions of HSPA-6, OASL, RNVU1, and RSAD2 (viperin) genes were processed via the ΔΔ-ct method with HPRT1 as a reference gene. All samples were analysed in biological and technical triplicates. Spearman correlation was evaluated with GraphPad Prism software. (F) Gene set enrichment analysis (DAVID tool; GO Biological processes; p < 0.001) of significantly up-regulated genes in analysed samples performed by DAVID tool. (G) Gene set enrichment analysis (DAVID tool; GO Biological processes; p < 0.001) of significantly down-regulated genes in analysed samples performed by DAVID tool.

Fig. 4

Differential expression analysis of miRNAs in TBEV-infected neurons and astrocytes revealed strain-specific and cell-specific responses. (A) Graphic summary of DE miRNAs in analysed samples. (↑) up-regulated miRNAs, (↓) down-regulated miRNAs. (B) Venn diagram describing the cell specificity of miRNA expression upon TBEV infection. (↑) up-regulated miRNAs, (↓) down-regulated miRNAs. (C) Venn diagram describing the strain specificity of miRNA expression upon TBEV infection. (H) – TBEV Hypr strain, (N) – TBEV Neudoerfl strain. (D) Heatmap visualizing the overall expression pattern of mi RNAs in analysed samples. (E) Verification of RNA-seq data by qRT-PCR. Relative expression of hsa-miR-1248 and hsa-miR-145-5p was calculated using the ΔΔ-ct method with hsa-miR-103a as a reference gene and Sp6 as an internal control. All samples were analysed in biological and technical triplicates. Spearman correlation was evaluated with GraphPad Prism software. (F) Frequency distribution of 21–23 nt vd-sRNAs mapped to TBEV Hypr genome in N72_H (black) and A72_H (red) samples. HS – hot-spot. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5

TBEV Hypr infection in astrocytes results in a stronger innate immune response than in neurons. (A) Heatmap summarizing results from gene set enrichment analysis (DAVID tool; GO Biological processes; p < 0.05) of DE genes 1) uniquely expressed in either A72_H or N72_H, and 2) genes expressed with >1.5 log₂ difference in one of the samples. (B) Relative expression of OASL and RSAD2 (viperin) genes in TBEV infected cells. Calculations were performed via the ΔΔ-ct method with HPRT1 as a reference gene. All samples were analysed in biologcal and technical triplicates. (C) Expression pattern of IFN-β and IFN-λ in TBEV-infected cells based on RNA-seq data. (D) Overview and comparison of the expression pattern of selected immune response-related genes from this study and works of Fares et al. [12] and Pokorna et al. [13] (E) Number of mapped reads for vd-sRNAs (16–50 nt) in the respective sample. Data summarise three independent experiments, and values in graphs are expressed as mean ± SD.

Fig. 6

Comparison of TBEV Hypr/Neudoerfl infection in neurons reveals disruption of proper neural development. (A) Heatmap summarizing results from gene set enrichment analysis (DAVID tool; GO Biological processes; p < 0.05) of DE genes 1) uniquely expressed in either N72_H or N72_N, and 2) genes expressed with >1.5 log₂ difference in one of the samples. (B) Number of mapped reads for vd-sRNAs (16–50 nt) in the respective sample. Data summarise three independent experiments and values in graphs are expressed as mean ± SD. (C) Frequency distribution of 21–23 nt vd-sRNAs mapped to TBEV Hypr/Neudoerfl genome in N72_H (black) and N72_N (red) samples. HS – hot-spot. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7

TBEV infection induces changes in the splicing of host pre-mRNAs. (A) Graphical summary of identified LSVs (local splicing variants) in RNA-seq data using MAJIQ tool. Binary LSVs include exon skipping (ES), intron retention (IR), alternative 3′/5′ splice site (A3SS/A5SS) and involve only two exons or two splice sites in the same exon. Complex LSVs combine several binary events originating from or targeting the same site. (B) Classification of binary and complex LSVs identified in TBEV-infected astrocytes and neurons (combined). (C) Heatmap summarizing results from gene set enrichment analysis (DAVID tool; GO Biological processes; p < 0.05) of genes in which > 1 LSV was identified (for samples A72_H, N72_H, and N72_N). (D) MAJIQ-generated Sashimi plots of top-ranked genes in N72_H (CFAP61) and A72_H (MTUS2) samples. Numbers above splicing event arcs represent the number of mapped junction reads. For the sake of readability, cropped areas with identified LSVs are shown. (E) Analysis of read depth distribution across CFAP61 and MTUS2 genes (areas with identified LSVs) in N72_H/M and A72_H/M samples, respectively, using the integrated genome browser.

Supplementary figure 1

Supplementary figure 2

Supplementary figure 3

Supplementary figure 4

Supplementary figure 5