Multi-cohort cross-omics analysis reveals disease mechanisms and therapeutic targets in HTLV-1-associated myelopathy, a neglected retroviral neuroinflammatory disorder

Abstract

HTLV-1 is an enigmatic retrovirus triggering a debilitating neuroinflammatory disease, HTLV-1-associated myelopathy (HAM), with unknown pathogenesis. Both HTLV-1 infection and HAM predominantly affect women and non-white neglected populations. HAM is lacking disease-modifying treatment, as current treatment is mostly symptomatic and inspired by either HIV-1 or multiple sclerosis therapeutic strategies. We used systems biology analyses of novel and publicly available data comprising (epi)genomics, transcriptomics, metabolomics and proteomics of multi-ancestry cohorts from a total of > 2500 People Living with HTLV-1 from 5 countries (Brazil, Peru, Japan, UK, US). Leveraging an unique admixed Brazilian cohort, genome-wide association study (GWAS) revealed African-specific variants in inflammasome sensor AIM2 with genome-wide significance (p < 5x10^-8). Suggestive loci (p > 5x10^-8) corresponding to metabolic, immune and neuronal genes were validated using published Japanese GWAS. Polygenic risk score and proviral load were independent disease predictors across ancestries. Systems biology analysis revealed neuronal/synaptic signaling, monocyte count, glucose/lipid metabolism, and neurocognition/depression as genetically linked to HAM. In silico drug screening identified estrogen blocker Fulvestrant as the top hit, while also confirming existing (pre)clinical data for HDAC inhibitors and immunosuppressants. Validated GWAS genes were overexpressed in HAM patients' whole blood and CD4 T-cells, as well as in spinal cord astrocytes, oligodendrocytes, and microglia by single-cell RNAseq. We experimentally confirmed decreased ApoA1/lipid/cholesterol levels, higher monocyte levels and lower neurocognitive scores in multi-ancestry cohorts. We found striking biological similarities between retroviral Hbz/Tax overexpression, Hbz interactome and HAM multi-omics findings: enrichment for lipid/cholesterol metabolism, estrogen signaling, neurodegenerative diseases, and viral pathways including EBV, recently identified as the major driver of multiple sclerosis. In conclusion, our data-driven approach uncovers novel disease mechanisms and therapeutic targets, and a validated polygenic risk score allowing targeted surveillance for high-risk individuals. A strong molecular overlap to other neurodegenerative/neuroinflammatory diseases reveals shared neuropathogenic pathways between unrelated viruses.

Keywords: ancestry; multi-omics; neglected disease; neuroinflammation; retrovirus.

Figure 1. GWAS in admixed Brazilians reveals ancestry-specific HAM loci while polygenic risk scores replicate across cohorts

A) The Brazilian cohort displayed mainly European (EUR) and African (AFR) admixture, with smaller contributions from Amerindian (AMR) and Asian (ASN) ancestries, as shown by PCA (left panel) and ADMIXTURE (right panel). B) Manhattan plot from Regenie GWAS showing resulting associations when using HAM status as a phenotype with logistic mixed modeling. The red horizontal line indicates the genome-wide significant threshold, and the blue horizontal line indicates the genome-wide suggestive threshold. C) LocusZoom plot of genome-wide significant variants from Regenie HAM GWAS, showing the significance of tested SNPs with HAM on the left y-axis. The right y-axis is the recombination rate. Each point represents a different SNP tested. The blue line is the level of recombination rate. The color of each point represents the strength of LD with the lead SNP. The shape of each point represents the type of SNP: splice variant (triangle), nonsynonymous (inverted triangle), synonymous (square), untranslated region (square), and none-of-the-above (circle) The lead SNP is labeled (diamond). Below are the genes at their corresponding genomic position. D) ROC curves (left panel) showing significant prediction of AS vs. HAM clinical status by univariate and multivariate logistic regression, comparing PRS and PVL, with and without age and sex correction. No correlation was observed between PRS and PVL (Spearman correlation, right panel)3. E) Violin plot overlaid with individual scores comparing the PRS distribution grouped by disease status: asymptomatics (AS), intermediate syndrome (IS), incident HAM (iHAM) and HAM (p<0.0001, AS vs. HAM, Mann-Whitney test). F) Example of TRACTOR analysis in a n a three-way admixed individual, individual stretches of each ancestry are heterogeneously dispersed throughout all chromosomes. G) Manhattan plots from TRACTOR GWAS in AFR, EUR and AMR ancestries, showing resulting associations when using HAM status as a phenotype with local ancestry-informed analysis. The red horizontal line indicates the genome-wide significant threshold, and the blue horizontal line indicates the genome-wide suggestive threshold. Right panel inset shows validation differential expression of plasma SOD2 protein in an independent Amerindian cohort of PLwHTLV (Lima, Peru, Mann-Whitney, p=0.02).

Figure 2. Systems biology analysis of shared genes and pathways between Brazilian and Japanese GWAS reveals novel neuronal and immune links to HAM pathogenesis

A) Out of 2687 genes identified by FUMA using Regenie Brazilian GWAS output, 435 overlapped with the Japanese GWAS at nominal p-values (Venn diagram, left panel). Gene Ontology (GO) analysis of Brazilian GWAS showed a strong enrichment for mitochondria, neuron/synapse and cell membrane (middle panel). Boxes on the right represent the enriched GO biological processes. The size of each box is proportional to the total width of the connections it is a part of. The width of each connection is proportional to the −log10 p-value returned from GO enrichment utilizing the corresponding gene set from the listed GWAS. Enrichment was confirmed when examining only the 435 overlapping genes (right panel), with comprehensive enrichment analysis using publicly available databases (including gene ontology (GO), KEGG/Reactome pathways, transcription factors motifs, human phenotype ontology and drug databases). B) Venn diagram shows a strong overlap at the pathway level, with 40 shared between Brazilian and Japanese GWAS (upper panel). Pathways grouped into three major classes: intracellular signaling (black), metabolic pathways (red), and neuronal pathways (blue). Pathways were highly interconnected when plotted as a network based on overlapping genes (lower panel, boxed inset). C) Enrichment of shared GWAS genes in the CGP (Chemical and Genetic Perturbations) database confirmed and extended the observed metabolic and neuronal associations at the cellular and disease level. D) Spearman correlation analysis between HAM PRS and a complete battery of neuropsychological tests performed in a subgroup of the Brazilian cohort (n=76), significant negative correlations between HAM PRS and working/spatial memory (Animals and R_matrix tests, both p<0.05), a trend was observed for activities of daily living (Lawton test, p=0.057). E) Phenome-wide enrichment analysis in the UK Biobank identified monocyte count as the top hit, alongside associations with broad neurological phenotypes (alcohol use related traits, anhedonia and depression) and neurocognitive skills (Snap-button test). F) Validation of the genetically predicted link of HAM to blood monocyte levels in two independent cohorts (Brazil and Peru). Absolute monocyte counts were increased in HAM vs. AS in both cohorts, but the most significant increase was observed in monocyte percentage (normalized to total leukocytes, p<0.01 for both cohorts). G) Blood monocyte counts at baseline were predictive of clinical outcome during long-time follow-up, with higher monocyte levels predicting worse disease progression (Osame Motor Disability Scale, p=0.004) as well as lower activities of daily living (Lawton, p=0.068).

Figure 3. Validation of GWAS predictions by multi-cohort transcriptomics and immune biomarker screening

A) Venn Diagram showing overlap between Brazilian GWAS genes identified by FUMA, transcriptomic analysis of differentially expressed genes in HAM in CD4 cells (Japanese cohort, “CD4 JP”) and in whole blood (UK cohort, “WB UK”), and KEGG-defined pathways (“DNA inflamm”: cytosolic sensing of pathogen DNA and inflammasome activation). B) From the cytosolic DNA sensing and inflammasome activation KEGG-defined genes, 15 genes overlapped with Brazilian GWAS genes, of which 8 were consistently upregulated in both “CD4 JP” and “WB UK “ data sets (all p<0.05 GEO2R analysis). C) Digital transcriptomic analysis of the US cohort showed higher pathway scores for IL-1 signaling and cytosolic DNA sensing in 2 unique incident HAM (iHAM) cases (One-sample t-test, p<0.05 for all), as compared to uninfected healthy controls (HC) or asymptomatic PLwHTLV (AS), which remained elevated after HAM diagnosis (individual iHAM cases are connected with dotted lines before and after HAM diagnosis). D) Transcription factor "upstream” analysis of 428 “HAM CD4” genes (as defined in A) across several tissues confirmed enrichment in brain and spinal cord. E) Circos plot combining GWAS, transcriptome data and interactome data showing AIM2, KIT, CD7 and NK cell-related genes (KLRs, NCRs) as major hubs in the HAM protein-protein interaction network (STRING database, middle part, all FDR-corrected p<0.05). Significant protein-protein interaction with previously known HAM biomarkers STAT1^, and CXCL10^- confirms the pathobiological significance of GWAS findings. F) Tabula Sapiens whole-body single-cell atlas (CellxGene) was used to map HAM genes prioritized by both GWAS and transcriptomics. Among immune cells (upper panels), we found sparse expression of AIM2 in B-cells > neutrophil > plasma cells >activated T-cells, while the “HAM inflammasome” gene module (as defined in B) was ubiquitously expressed with Tregs>macrophages>monocytes>NK cells. Exploring a spinal cord single-cell dataset (CellXgene), we found no detectable AIM2 but the “HAM inflammasome” gene module (as defined in B) was modestly expressed in microglial cells>endothelial cells>macrophages>oligodendrocyte precursors>astrocytes (lower panels). Exploring the “broad HAM” gene module (95 genes shared between Brazilian and Japanese GWAS, all overexpressed in HAM CD4 cells), we observed modest expression among whole-body immune cells (upper right panel: highest in macrophages>monocytes>microglial cells>B-cells), but strong and widespread expression across all cell types in the spinal cord (lower right panel: oligodendrocyte precursor cells (OPC) >astrocytes >oligodendrocytes >microglia >mural cells> vascular smooth muscle cells> glutamergic neurons> macrophages> cerebral granule cells>GABAergic neurons >capillary endothelial cells >neurons). Inset (lower right panel) shows the quantification of all 95 genes in spinal cord cell types as density plots. G) Spearman correlation analysis of HAM PRS and circulating plasma cytokine levels for IL-6 (p=0.038) and plasma IFN-γ (p=0.061) in 66 patients for which paired samples were available.

Figure 4. Genetic links to lipid metabolism and proviral load are validated by in vivo metabolomics and in vitro retroviral Tax/Hbz overexpression and interactome

A) Using log-transformed PVL as a phenotype, 841, 27, 11, 402, and 162 GO terms were enriched for AFR, AMR, ASN, EUR, and Regenie datasets, respectively. AFR and EUR PVL GWAS showed the strongest enrichment, highlighting small molecule and lipid metabolism as well as oxidation-reduction processes. B) ¹H-NMR metabolomic analysis in a subset of the Brazilian GWAS cohort (n=110) showed lower levels of lipid metabolites in HAM compared to AS, including ApoA1 (p<0.001), ApoB (p<0.05), total (p<0.01) and HDL-cholesterol (p<0.05). ApoA1 levels were already decreased in patients with prodromal intermediate syndrome (Intermed, *p<0.05), suggesting an early event in HAM pathogenesis, Consistent with PVL GWAS enrichment analysis, proviral load was inversely correlated with total and clinical LDL-cholesterol levels (p<0.05 for both). C) Pathway enrichment analysis of all genes upregulated by both Tax and Hbz overexpression in Jurkat cells. D) Validation of a shared molecular imprint of viral Tax and Hbz using digital transcriptomic analysis in the US cohort. Total viral RNA (Tax+Hbz) levels were positively correlated to cell cycle (p=0.0012) and metabolism pathway scores (p=0.0026), as well as several transcripts corresponding to GWAS genes (CTLA4, CADM1, BIRC5, al p<0.05). E) Enrichment analysis of the HBZ human interactome revealed striking biological similarities to HAM GWAS findings: cell cycle, lipid metabolism and oxidative phosphorylation pathways were confirmed, as well as several viruses that can cause neurodegeneration (HIV-1), encephalitis (HSV-1), and Epstein-Barr virus (EBV), recently identified as the major driver of neuroinflammation in multiple sclerosis. EBV infection was the most central pathway (yellow nodes) in the network, with strong molecular overlap to other neurodegenerative diseases (spinocerebellar ataxia, ALS, Alzheimer, Parkinson and prion disease). F) Proposed HAM disease mechanism(s) of spinal cord infiltration and inflammation triggered by the combined effects of AIM2 inflammasome activation, monocyte and microglial activation, proinflammatory cytokine secretion (IL1γ, IL-6, TNF, IFN-γ), accompanied by decreased antioxidant capacity mitochondrial (SOD2) and reprogramming of (systemic and/or local) lipid metabolism.