Characterizing changes in transcriptome and kinome responses in testicular cells during infection by Ebola virus

Abstract

Ebola virus (EBOV) is able to persist and actively replicate in the reproductive tract of male disease survivors months or years after recovery from Ebola virus disease (EVD)¹. Persistent EBOV infections are usually asymptomatic and can be transmitted sexually, but the host and viral factors that mediate these infections have not been characterized^2,3. We investigated the interaction between host and viral factors during EBOV infection of the blood testis barrier (BTB), with a focus on Sertoli cells as a potential reservoir for viral persistence. We assessed viral replication kinetics and host responses of mouse testicular Leydig cells and Sertoli cells infected with EBOV Makona (i.e. infectious EBOV) and collected samples up to 28 days post-infection. Viral replication was apparent in both cell lines, but intracellular early viral loads were much higher in Leydig cells compared to Sertoli cells. We used RNAseq analysis to characterize transcriptomic responses of Leydig cells and Sertoli cells to EBOV infection over time. Further investigation of early interactions between host cells and EBOV was performed using virus-like particles (EBOV trVLP) and assays of phosphorylation-based cell signaling. Our findings indicate that virus-treated Sertoli cells responded more rapidly and robustly than Leydig cells, and with a particular emphasis on detection of, and response to, external stimuli. We discuss how the roles played by Sertoli cells in immune privilege and spermatogenesis may affect their initial and continued response to EBOV infection in a manner that could facilitate asymptomatic persistence.

Fig. 1. EBOV replication kinetics in Leydig cells and Sertoli cells.

Quantitation by RT-PCR of VP40 gene copies in cell pellet and supernatant of Leydig cells (A) and Sertoli cells (B) treated with infectious EBOV. Primer sequences for strand-specific RNA reverse transcription and quantitation of the VP40 gene are listed in Table 1. The cutoff for minimum threshold of detection was at 36 cycles. Each data point represents the average of three biological replicates, for each of which two technical replicates were performed during quantitation. For each cell type, significant changes in viral load between neighboring time points were determined using the Student’s t test. Asterisks indicate a significant change in viral load, where *, **, ***, and **** indicates that P ≤ 0.05, P ≤ 0.01, P ≤ 0.001, and P ≤ 0.0001, respectively.

Fig. 2. Cluster analysis of differentially expressed genes in EBOV-infected Leydig cells.

A Counts of significantly upregulated (log2 FC ≥ 2; P ≤ 0.05) and downregulated (log2 FC ≤ −2; P ≤ 0.05) genes for Leydig cells. Values were calculated from three replicates for each set of variables, using normalized feature counts determined by RNAseq analysis. For each time point, separate queries of upregulated and downregulated genes were performed against the Gene Ontology Biological Processes (https://geneontology.org; accessed October 5, 2023) and Reactome Gene Sets (https://reactome.org; accessed October 5, 2023) to identify associated ontological terms^–. Significance of associated terms was determined by calculating accumulative hypergeometric p-values (cutoff: 0.01) and enrichment factors (minimum 1.5). Samples were clustered based by significant terms using Kappa-statistical similarities (0.3 kappa threshold) among upregulated (B) and downregulated (C) member genes. A heatmap was generated and colored by p-values, where white cells indicate the lack of enrichment for that term in the corresponding gene list. The term with the best p-value for each term cluster was used as a representative on the heatmap. Protein-protein interaction networks were generated based on upregulated (D) and downregulated (E) genes, and the MCODE algorithm was applied to these networks to identify neighborhoods where proteins are densely connected. Enrichment analysis was applied to each MCODE network to extract “biological meanings” from the network component, where top three best p-value terms were retained and represented based by one term per cluster as a label. A breakdown of all MCODE interpretations for each time point is available as supplementary data. Nodes represent individual proteins, node pie sectors indicate which samples differentially expressed each protein, and edges represent interactions between proteins. Analysis was performed using the Metascape online tool (https://metascape.org; accessed October 5, 2023), which incorporates network visualization by Cytoscape^,. L01, Leydig + infectious EBOV 1 dpi; L02, Leydig + infectious EBOV 2 dpi; L04, Leydig + infectious EBOV 4 dpi; L07, Leydig + infectious EBOV 7 dpi; L14, Leydig + infectious EBOV 14 dpi.

Fig. 3. Cluster analysis of differentially expressed genes in EBOV-infected Sertoli cells compared to mock-infected cells.

A Counts of significantly upregulated (log2 FC ≥ 2; P ≤ 0.05) and downregulated (log2 FC ≤ −2; P ≤ 0.05) genes for Sertoli cells. Values were calculated from three replicates for each set of variables, using normalized feature counts determined by RNAseq analysis. For each time point, separate queries of upregulated and downregulated genes were performed against the Gene Ontology Biological Processes (https://geneontology.org; accessed October 5, 2023) and Reactome Gene Sets (https://reactome.org; accessed October 5, 2023) to identify associated ontological terms. Significance of associated terms was determined by calculating accumulative hypergeometric p-values (cutoff: 0.01) and enrichment factors (minimum 1.5)^–. Samples were clustered based by significant terms using Kappa-statistical similarities (0.3 kappa threshold) among upregulated (B) and downregulated (C) member genes. A heatmap was generated and colored by p-values, where white cells indicate the lack of enrichment for that term in the corresponding gene list. The term with the best p-value for each term cluster was used as a representative on the heatmap. Protein-protein interaction networks were generated based on upregulated (D) and downregulated (E) genes, and the MCODE algorithm was applied to these networks to identify neighborhoods where proteins are densely connected. Enrichment analysis was applied to each MCODE network to extract “biological meanings” from the network component, where top three best p-value terms were retained and represented based by one term per cluster as a label. A breakdown of all MCODE interpretations for each time point is available as supplementary data. Nodes represent individual proteins, node pie sectors indicate which samples differentially expressed each protein, and edges represent interactions between proteins. Analysis was performed using the Metascape online tool (https://metascape.org; accessed October 5, 2023), which incorporates network visualization by Cytoscape^,. S01, Sertoli + infectious EBOV 1 dpi; S02, Sertoli + infectious EBOV 2 dpi; S04, Sertoli + infectious EBOV 4 dpi; S07, Sertoli + infectious EBOV 7 dpi; S14, Sertoli + infectious EBOV 14 dpi; S28, Sertoli + infectious EBOV 28 dpi.

Fig. 4. Leydig cells and Sertoli cells are permissive to trVLP entry and transcription.

A Quantification of VP40 gene copies in Leydig cells treated with EBOV trVLP. B Luminescence of the firefly luciferase reporter in mock-treated and EBOV trVLP-treated Leydig cells. C Quantitation of VP40 gene copies in Sertoli cells treated with EBOV trVLP. D Luminescence of the firefly luciferase reporter in mock-treated and EBOV trVLP-treated Sertoli cells. For each time point, EBOV VP40 RNA and luciferase luminescence was quantified based on three biological replicates and two technical replicates. Primer sequences for strand-specific RNA reverse transcription and quantitation of viral VP40 are listed in Table 1. Significant differences in firefly luciferase luminescence between mock-treated and EBOV trVLP-treated cells was determined using 2-way ANOVA. Significant changes in quantity of EBOV VP40 between neighboring time points were determined using the Student’s t test. Asterisks indicate a significant difference, where *, **, ***, and **** indicates that P ≤ 0.05, P ≤ 0.01, P ≤ 0.001, and P ≤ 0.0001, respectively.

Fig. 5. Assessment of kinome response differences between Leydig cells and Sertoli cells treated with EBOV trVLP.

For each peptide, the average of three technical replicate signal intensity values for EBOV trVLP-treated samples was compared against that of time-matched controls; changes in phosphorylation were identified by paired t tests. Mathematical analysis of data was performed using the Platform for Integrated, Intelligent Kinome Analysis 2 (PIIKA2) software. The Pearson correlation distance metric and the McQuitty linkage method were used for hierarchical clustering to compare phosphorylation at all peptides (A), and then by comparing only peptides that reported consistent results across biological replicates (B). Protein-protein interaction networks were generated based on increased (C) and decreased (D) phosphorylation of peptides. Ontological terms were assigned based on Reactome Gene Sets (https://reactome.org; accessed October 5, 2023) and the MCODE algorithm was applied to these networks to identify neighborhoods where proteins are densely connected. Enrichment analysis was applied to each MCODE network to extract “biological meanings” from the network component, where top three best p-value terms were retained and represented by one term per cluster as a label. Nodes represent individual proteins, node pie sectors indicate which samples differentially expressed each protein, and edges represent interactions between proteins. Protein-protein interaction analysis was performed using the Metascape online tool (https://metascape.org; accessed October 5, 2023), which incorporates network visualization by Cytoscape^,. L01, Leydig cells 1 h post treatment; L06, Leydig cells 6 h post treatment; S01, Sertoli cells 1 h post treatment; S06, Sertoli cells 6 h post treatment.

Fig. 6. Characterization of longitudinal cell membrane integrity in Leydig and Sertoli cells following EBOV trVLP treatment.

Cells were seeded 48 h prior to treatment, and the figures in question display membrane resistance (A, C) and cell capacitance (B, D) from 24 h to 96 h after seeding. Data points represent the average value of nine wells per treatment per cell line, as three of the original twelve columns were removed from the experiment prior to treatment due to failed quality control readings. The point at which treatments were added is represented by a downward-pointing arrow (↓), and the point at which media was replaced is indicated by an upward-pointing arrow (↑). Readings were paused while the ECIS plate was manipulated to apply treatments replace media. Resistance readings are displayed at 8.00 × 10³ Hz and capacitance readings are displayed at 6.40 × 10⁴ Hz.