Integrating human endogenous retroviruses into transcriptome-wide association studies highlights novel risk factors for major psychiatric conditions

Abstract

Human endogenous retroviruses (HERVs) are repetitive elements previously implicated in major psychiatric conditions, but their role in aetiology remains unclear. Here, we perform specialised transcriptome-wide association studies that consider HERV expression quantified to precise genomic locations, using RNA sequencing and genetic data from 792 post-mortem brain samples. In Europeans, we identify 1238 HERVs with expression regulated in cis, of which 26 represent expression signals associated with psychiatric disorders, with ten being conditionally independent from neighbouring expression signals. Of these, five are additionally significant in fine-mapping analyses and thus are considered high confidence risk HERVs. These include two HERV expression signatures specific to schizophrenia risk, one shared between schizophrenia and bipolar disorder, and one specific to major depressive disorder. No robust signatures are identified for autism spectrum conditions or attention deficit hyperactivity disorder in Europeans, or for any psychiatric trait in other ancestries, although this is likely a result of relatively limited statistical power. Ultimately, our study highlights extensive HERV expression and regulation in the adult cortex, including in association with psychiatric disorder risk, therefore providing a rationale for exploring neurological HERV expression in complex neuropsychiatric traits.

Fig. 1. A summary of the retrotranscriptome-wide association study (rTWAS) approach.

A RNA-sequencing and genotype data from individuals of European (EUR, N = 563) or African ancestry (AFR, N = 229) are used to construct (B) single nucleotide polymorphism (SNP) weights. The example depicts a genetic feature more expressed in association with the A-allele from a hypothetical local variant, relative to the alternative a-allele. C GWAS results are then cross-referenced with the SNP weights using a transcriptome-wide association study (TWAS) approach, to identify expression signatures associated with risk. The example illustrates that the A-allele of the hypothetical variant, associated with increased expression of the hypothetical genetic feature, is also associated with trait susceptibility. D Sensitivity analyses including i. conditional analyses and ii. fine-mapping then allow inference of which expression signals are considered E high confidence risk features, as indicated by their ability to independently explain the genetic signal at their respective loci. Created with Biorender.com. This image is published under a CC BY-NC-ND license.

Fig. 2. Retrotranscriptome-wide association studies of major psychiatric disorders.

The Manhattan biplots show the expression signatures significantly associated with (A) schizophrenia, (B) bipolar disorder, and (C) major depressive disorder. We found no HERV expression signatures associated with attention deficit hyperactivity disorder and autism spectrum conditions, so these are omitted. The X-axis indicates genomic location, whereas the Y-axis shows Z score from the TWAS. The horizontal grey lines indicate transcriptome-wide significance, i.e., a threshold adjusted for the number of expressed features using the Bonferroni method (two-sided P value cut-off = 6.10 × 10⁻⁶). Only Bonferroni-significant HERV features are labelled.

Fig. 3. Predicted HERV expression signatures explaining GWAS signals at multiple locations.

A For schizophrenia, we observed instances where the HERV expression signal was the best feature to explain some of the GWAS signal at the locus, e.g., i. MER4_20q13.13 and ii. ERV316A3_5q14.3j, and an instance where more than one expression feature, including a HERV, were associated with risk, e.g., iii. ERV316A3_2q33.1 g. B For bipolar disorder, we also observed the expression of MER4_20q13.13 as a feature explaining the GWAS signal at its locus. C For major depressive disorder, multiple expression signatures correlated with risk on chromosome 1p31 (feature names labelled in blue), but ERVLE_1p31.1c showed independent association with the disorder (feature name labelled in green). Upper part of each image: genomic context. Lower part of each image: a plot in which the X-axis indicates genomic location, and the Y-axis shows -log10(P) of genetic variant associations (from the GWAS, two-sided), before (grey dots) and after (blue dots) conditioning on jointly significant genes in each locus. P-values are not adjusted for multiple testing. Only high confidence risk HERVs are shown.

Fig. 4. Fine-mapping analysis supports high confidence risk HERVs for multiple psychiatric disorders.

The graphs correspond to the HERV expression signals in the fine-mapping analysis that are also significant in the conditional analyses, in relation to (A) schizophrenia, including i. MER4_20q13.13, ii. ERV316A3_5q14.3j, and iii. ERV316A3_2q33.1 g; (B) bipolar disorder, which includes MER4_20q13.13; and (C) major depressive disorder, which includes ERVLE_1p31.1c. Upper part of each image: graph where the Y-axis indicates the TWAS association p value (two-sided), unadjusted for multiple testing, and the X-axis shows genetic features in the linkage disequilibrium block. The size and colour of the points indicate the posterior inclusion probability (PIP), indicating the probability that the expression feature is causal for the association signal at the locus. Lower part: correlation of predicted expression.

Fig. 5. Genomic context of high confidence risk HERVs.

A Expression of HERVs and their nearest canonical genes are shown as median values with interquartile range, with outliers depicted separately (N = 563 biologically independent samples of European ancestry). B Analysis using HOMER indicates that approximately 98% of HERVs from Telescope are in intergenic and intronic regions, whereas the remainder (‘Other’) is located in promoters, untranslated regions, or transcription start or termination sites. C The genomic context of ERV316A3_2q33.1 g and (D) ERV316A3_5q14.3j suggests that these HERVs are likely part of specific isoforms of canonical genes FTCDNL1 and ADGRV1, respectively. On the other hand, (E) MER4_20q13.13 is encoded in the opposite strand of the canonical gene PTGIS, and (F) ERVLE_1p31.1c is intergenic (nearest gene is NERG1), suggesting that they are likely producing novel non-coding RNAs.

Fig. 6. Co-expression analysis identifies HERVs co-expressed with canonical genes and supports their role in a range of biological functions.

A Proportion of HERVs and canonical genes assigned to each co-expression module, including the number of genetic features per module at the top, as detected in the European subset (N = 563 biologically independent samples). B Bubble plot showing top gene ontology (GO) term, per module. The X-axis and colour of the bubbles indicate −log10(P) (two-sided, uncorrected) of the enrichment statistic. The size of the bubbles represents the enrichment ratio. Only Bonferroni-significant GO terms are shown (Bonferroni-adjusted P < 0.05).