Identification of novel reproducible combinatorial genetic risk factors for myalgic encephalomyelitis in the DecodeME patient cohort and commonalities with long COVID.

We identified 22,411 double-refined signatures, comprised of 7,555 SNPs mapped to 2,311 genes, that are consistently associated with increased odds of ME in multiple DecodeME and UKB cohorts (Table  2 , Supplementary Tables 6 – 7 ). Ignoring genotypes, the double-refined signatures represent 20,496 combinations of SNPs. 38% of the double-refined signatures have odds ratios greater than 1.1 in the appropriate Refinement dataset, with a maximum odds ratio = 2.55, mean odds ratio = 1.10, and median odds ratios = 1.09 (see Fig.  4 , Supplementary Fig.  4 ). Table 2 Number of signatures, component SNPs, and mapping genes from the initial DecodeME combinatorial analyses after Refinement in UKB and an independent DecodeME cohort # Double-refined signatures # SNPs in double-refined signatures # Genes mapping to double-refined signatures # Candidate core genes 22,411 7,555 2,311 259 Candidate core genes were selected based on significant association with ME and > 20% case prevalence in a DecodeME Refinement dataset. Number of signatures, component SNPs, and mapping genes from the initial DecodeME combinatorial analyses after Refinement in UKB and an independent DecodeME cohort Candidate core genes were selected based on significant association with ME and > 20% case prevalence in a DecodeME Refinement dataset. Fig. 4 Distribution of odds ratios of double-refined signatures in appropriate Refinement datasets. Odds ratio represents the odds that a randomly selected person with the disease signature will be an ME case relative to the odds that a randomly selected person without the disease signature will be an ME case. Y-axis labels are odds ratios plotted on log2 scale for improved readability. Distribution of odds ratios of double-refined signatures in appropriate Refinement datasets. Odds ratio represents the odds that a randomly selected person with the disease signature will be an ME case relative to the odds that a randomly selected person without the disease signature will be an ME case. Y-axis labels are odds ratios plotted on log2 scale for improved readability. Extended Table  1 lists the 7,555 SNPs and 2,311 genes mapped to the double-refined signatures and Extended Table  2 lists locus information for the SNPs and signatures. The SNP combinations predominantly have a ratio of females to males that is closely aligned to the overall sex ratio in the dataset, suggesting that the signatures are not sex-biased (Supplementary Fig.  5 ). Breakdowns of number of signatures and SNPs following each step of the Discovery/Refinement pipeline are included in the Supplementary Materials. Our results highlight the importance of the 6p22.1 band on chromosome 6, which contains the HLA-A , HLA-F , and HLA-G genes, for ME disease biology. 29% of the double-refined signatures contain at least one SNP located within this gene-rich region of low recombination, far more than any other chromosomal band. HLA genes are critical for immune system regulation and elimination of viral infections, and several previous studies have identified a link between HLA genes, including HLA-A variants, and protection against ME and long COVID [ 30 , 50 – 54 ]. To aid in identification of candidate core genes and mechanisms, we filtered the gene list from each analysis by prioritizing genes that had disease odds ratios greater than the median (1.046 in Analysis 1 and 1.050 in Analysis 2), p -values less than 0.05, and mapped to at least 20% of cases in the appropriate Refinement dataset. In total, 259 unique genes (Extended Table  3 ) satisfied these criteria. Pathway enrichment analysis on the candidate core gene set highlighted multiple cellular processes involved in nervous system development and neural signaling, immune response, cellular stress response and calcium signaling (Supplementary Fig.  6 ). Many of these processes have been previously linked to ME, supporting their potential involvement [ 32 , 37 , 55 ]. Additionally, amine-mediated responses may be linked to neuroimmune regulation and autonomic function [ 56 , 57 ], whereas calcium-dependent signaling plays a key role in many biological processes relevant to ME such as neuronal signaling, regulation of Purkinje cell migration and pathobiology of ataxia [ 58 , 59 ]. The candidate core gene set showed overall enrichment of expression in 6 brain regions and spinal cord, which broadly aligned with the pathway enrichment results (see Supplementary Results). Alongside the hypothesis-free pathway enrichment analysis, we functionally characterized the candidate core genes using pre-defined functional groups representing literature-derived biological processes relevant to ME to better understand their functional context. This hypothesis-driven approach showed that more than 50% of the candidate core genes were linked to energy metabolism, followed by neurological processes (13.1%), inflammation or autoimmunity (10%), mitochondrial dysfunction (9.7%) and innate immunity (5.8%) (see Supplementary Results). The candidate core genes include some genes that have previously associated with ME. This includes TLR3 , involved in innate immune response, which has an agonist currently in clinical development for ME treatment [ 60 ], along with CSE1L and DCC which were associated with ME and chronic pain respectively in a GWAS including DecodeME participants of European ancestry [ 30 ]. The remaining 256 genes represent novel disease associations, although 6 genes ( CSMD1 , CYB5RL , MRPL37 , PAX5 , SRGAP1 , RAP1GAP2 ) were previously reported in disease signatures identified in the combinatorial analysis of the fully disjoint UKB ME cohort [ 32 ], with their SNPs co-associated with key genes. We examined associations between gene-defined patient subgroups and self-reported symptoms and comorbidities reported in DecodeME questionnaires. This large-scale phenotype association analysis did not identify associations passing FDR correction. However, exploratory analysis highlighted top-ranked gene-phenotype pairs, described in the Supplementary Results and Extended Table  4 . In order to evaluate the predictive power of the double-refined disease signatures, we used the DecodeME Test dataset (Cohorts D + H) to test whether the count of double-refined signatures is associated with increased odds of ME. Importantly, none of the samples in the Test dataset overlap with samples used for disease signature Discovery or Refinement. This is not intended to be a polygenic risk score or equivalent – it is a simplistic, non-optimized tool to test whether the cumulative set of signatures are associated with increased risk of ME, while minimizing the risk of overfitting. Each feature (i.e., signature in this case) is assumed to additively and independently contribute to a patient’s risk of disease, and all signatures are assumed to have identical effect sizes. Combinatorial signatures are clearly not independent, as they can (and often do) share SNPs that are identical or are sometimes in linkage disequilibrium. Nor do they have identical effect sizes. Nonetheless, this simplistic additive signature count score approach enables unbiased testing of the association of the signatures in aggregate. We observed a bimodal distribution for the number of double-refined signatures possessed by samples in the Test dataset (Fig.  5 a). This reflects signatures that contain at least one of 187 SNPs from the extended major histocompatibility (MHC) region located within the p22.1 cytogenic band of chromosome 6. In total, 6,404 double-refined signatures (29% of total) map to this region of low recombination and are strongly correlated due to shared SNPs and high linkage disequilibrium. As a result, individuals are more likely to possess either many or few signatures mapped to this region than they are to possess intermediate numbers of signatures. Removing signatures that contain at least one 6p22.1 SNP results in a unimodal distribution of double-refined signatures possessed by ME cases (Fig.  5 b). Below we refer to this adjusted count as the signature count score, which removes the potential strong confounding factor introduced by overweighting signal related to 6p22.1 . Fig. 5 Density plots showing distribution of number of double-refined signatures possessed by individuals in the Test dataset cohort. Distributions for cases (pink) and controls (blue) are plotted separately. ( a ) Count of all double-refined ME signatures. ( b ) Count of double-refined ME signatures not mapped to 6p22.1 Density plots showing distribution of number of double-refined signatures possessed by individuals in the Test dataset cohort. Distributions for cases (pink) and controls (blue) are plotted separately. ( a ) Count of all double-refined ME signatures. ( b ) Count of double-refined ME signatures not mapped to 6p22.1 We observed a highly significant correlation between the signature count score and ME in the Test dataset (OR = 1.23 per standard deviation increase, p  = 4 × 10 − 21 ) when including sex and the top 10 genetic principal components as confounders in the logistic regression. The disease odds ratio for individuals with a top 10% signature count score relative to individuals with a bottom 10% signature count score was 1.64 (Fig.  6 ), while the odds ratio for the top 5% vs. bottom 5% was 1.89. Fig. 6 Relative odds of ME across deciles of counts of signatures possessed by individuals in the test dataset. Odds for each decile are standardized to the odds of the bottom 10% decile. Deciles are labeled with their equivalent lower-bound percentile. Error bars represent standard error. Signature count excludes double-refined ME signatures mapped to 6p22.1 to eliminate confounding due to the disproportionate number of strongly correlated signatures mapped to that locus Relative odds of ME across deciles of counts of signatures possessed by individuals in the test dataset. Odds for each decile are standardized to the odds of the bottom 10% decile. Deciles are labeled with their equivalent lower-bound percentile. Error bars represent standard error. Signature count excludes double-refined ME signatures mapped to 6p22.1 to eliminate confounding due to the disproportionate number of strongly correlated signatures mapped to that locus The DecodeME GWAS study identified 8 genome-wide significant loci associated with ME in the DecodeME cohort using 15,579 cases, over 150,000 controls, and more than 8.8 million variants including imputed variants [ 30 ]. For each locus, they identified a lead variant mapped to a gene and classified a total of 37 genes into Tier 1 and Tier 2 based on the strength of genetic evidence and biological relevance [ 30 ]. The study also reported an additional non-genome-wide significant gene linked to chronic pain. The double-refined signatures map to 13 of the 32 Tier 1 or 2 DecodeME GWAS study genes including genes located in 6 of the 8 GWAS loci (Table  3 , Supplementary Table 8 ). The double-refined signatures also contain SNPs that map to the remaining 2 loci but which are not located within protein-coding genes, potentially representing non-coding variants that affect expression of one or more Tier 1 or 2 genes. Of the Tier 1 or 2 GWAS genes, only CSE1L was included among the 259 candidate core genes, along with DCC , which is linked to chronic pain. Table 3 Genes identified by combinatorial analysis that map to GWAS significant and other loci reported by DecodeME study GWAS Significant Locus Gene DecodeME Tier Range of Odds Ratios for Double-Refined Signatures in Refinement Dataset Minimum p -value for Double-Refined Signatures in Refinement Dataset 1q25.1 DARS2 Tier 1 other 1.15–1.79 2 × 10 − 12 6p22.2 BTN2A2 Tier 1 best 1.04–1.78 2 × 10 − 10 ABT1 Tier 1 other 1.03–1.98 9 × 10 − 11 TRIM38 Tier 1 other 1.12–1.43 3 × 10 − 9 HFE Tier 1 other 1.21–1.50 9 × 10 − 9 ZNF322 Tier 1 other 1.13–1.42 5 × 10 − 6 12q24.23 SUDS3 Tier 1 best 1.00–1.81 4 × 10 − 12 VSIG10 Tier 1 other 1.18–1.43 1 × 10 − 7 13q14.3 OLFM4 Tier 2 best 1.13–1.44 2 × 10 − 5 17q22 CA10 Tier 1 best 1.22–1.37 9 × 10 − 6 20q13.13 CSE1L Tier 1 best 1.00–1.90 9 × 10 − 13 DDX27 Tier 1 other 1.01–2.00 9 × 10 − 13 ZNFX1 Tier 1 other 1.01–2.00 9 × 10 − 13 - DCC Chronic Pain 1.11–1.24 1 × 10 − 4 The lead genes ( n  = 5) reported by the DecodeME GWAS are shown in bold. Includes range of odds ratios and minimum p -values for associated double-refined signatures in the corresponding refinement dataset Genes identified by combinatorial analysis that map to GWAS significant and other loci reported by DecodeME study The lead genes ( n  = 5) reported by the DecodeME GWAS are shown in bold. Includes range of odds ratios and minimum p -values for associated double-refined signatures in the corresponding refinement dataset We previously identified 235 long COVID associated genes in three patient cohorts derived from the Sano GOLD study cohort [ 35 ], of which 180 were reproduced in All of Us long COVID cohort [ 36 ]. These included: 43 genes from a ‘Severe’ cohort who reported the most diverse and severe symptoms. 35 genes from a ‘Fatigue Dominant’ cohort who reported mainly fatigue-related symptoms. 165 genes, from a ‘General’ cohort self-reporting ongoing or new symptoms 12 weeks after initial SARS-CoV-2 infection, with symptoms lasting at least 2 months without other explanation. 43 genes from a ‘Severe’ cohort who reported the most diverse and severe symptoms. 35 genes from a ‘Fatigue Dominant’ cohort who reported mainly fatigue-related symptoms. 165 genes, from a ‘General’ cohort self-reporting ongoing or new symptoms 12 weeks after initial SARS-CoV-2 infection, with symptoms lasting at least 2 months without other explanation. 102 (43%) of the 235 ME genes previously associated with long COVID also mapped to the double-refined ME signatures, including 76 of the 180 genes (42%) that reproduced in the All of Us long COVID cohort (Supplementary Table 9 ). Both results represent a significant enrichment of long COVID genes mapping to the ME signatures (Fisher’s exact test p  < 10 − 30 and p  < 10 − 21 respectively). These p -values, however, assume that all genes have similar likelihood of being observed in a combinatorial analysis. As genes vary considerably in length our gene mapping approach could potentially be biased towards identifying longer genes that map to many SNPs. We therefore also conducted a second enrichment analysis to evaluate the overlap between ME and long COVID genes, employing a length-matched permutation test to control for potential biases that favor identification of longer genes (see Supplementary Methods). This analysis showed significant enrichment ( p  < 0.017) across all three gene length bins (lengths lower than 50 kb, 50–250 kb and above 250 kb) after adjusting for multiple-testing correction (see Fig.  7 and Supplementary Table 10 ). This result confirms that the overlap between long COVID and ME genes is not simply an artifact of gene length. Fig. 7 We observed significant overlap between the genes associated with ME by this study and the genes associated with long COVID by a published combinatorial analysis [ 35 ] even after controlling for gene length. Red diamonds represent the number of long COVID genes that were also mapped to the double-refined ME disease signatures. Box plots show the distribution of observed overlap with the ME genes after randomly selecting an equal number of genes to the corresponding long COVID gene lists, repeated over 10,000 randomized permutations. Boxes represent the top 25th and bottom 25th percentile ranks for SNPs assigned to each gene set with solid line denoting median SNP rank and tails denoting the limits of the middle 95th percentile. p -values above each plot denote the fraction of permutations in which the overlap between the randomly selected genes and the set of ME genes is greater than or equal to the observed overlap between the sets of long COVID and ME genes We observed significant overlap between the genes associated with ME by this study and the genes associated with long COVID by a published combinatorial analysis [ 35 ] even after controlling for gene length. Red diamonds represent the number of long COVID genes that were also mapped to the double-refined ME disease signatures. Box plots show the distribution of observed overlap with the ME genes after randomly selecting an equal number of genes to the corresponding long COVID gene lists, repeated over 10,000 randomized permutations. Boxes represent the top 25th and bottom 25th percentile ranks for SNPs assigned to each gene set with solid line denoting median SNP rank and tails denoting the limits of the middle 95th percentile. p -values above each plot denote the fraction of permutations in which the overlap between the randomly selected genes and the set of ME genes is greater than or equal to the observed overlap between the sets of long COVID and ME genes The identified overlap between ME and long COVID is perhaps an underestimate as the originating studies were carried out based on genotype data from two different arrays – Axiom UK Biobank for ME and Illumina Global Screening Array for the Sano GOLD study. The two QCed datasets had very limited SNP overlap (just 10%), and without using imputation, it is only possible to evaluate overlap at a gene level. Since both the DecodeME and All of Us datasets provided fewer than 400,000 SNPs this can only count as a minimum estimate of potential overlap. We also observed a highly significant enrichment of top ranked SNPs in the ML model for ME (i.e., SNPs with high feature importance) assigned to the sets of genes previously associated with long COVID (Fig.  8 ). The top 25th percentile ranks are significantly higher in the sets of SNPs assigned to the gene sets for each long COVID phenotype than the expected top 25th percentile rank for randomly selected sets of the same numbers of genes ( p <  0.001 for all three long COVID phenotypes). SNP rankings in the ML model were particularly strong for the genes previously associated with Severe and Fatigue Dominant long COVID. For example, 25% of SNPs assigned to genes associated with Severe long COVID were included in the top 0.5% of SNPs with the highest feature weights in the ML model for ME. Similarly, 25% of SNPs assigned to genes associated with Fatigue Dominant long COVID were included in the top 2% of SNPs with the highest feature weights in the ML model for ME. These results provide independent validation that a large number of long COVID genes are also biologically important for ME. Alternatively, this result could represent an inherent bias in the methodology that would be expected if longer genes mapping to many SNPs are both more likely to be identified by combinatorial analysis and to map to at least one SNP with strong feature importance in the ML model for ME. However, we did not observe the same extreme enrichment of very top ranked model SNPs mapping to endometriosis and coronary artery disease (CAD) disease signatures that we observed for Fatigue Dominant and Severe long COVID (Supplementary Fig.  10 ). This finding confirms that the very strong enrichment of top ranked ME SNPs that also map to long COVID genes is not an inevitable artifact of any potential gene length biases described above. Genes mapping to disease signatures for endometriosis and CAD also exhibited significantly higher ranks in the ML model for ME relative to randomly selected genes. This pattern is again apparent when comparing the ML model ranks of the top 25th percentile SNPs for each gene set (Supplementary Fig.  7 ). The degree of overlap is, however, substantially lower than that observed for long COVID, confirming that observation of very strong overlaps between top ranked SNPs in the ML model for ME and Fatigue Dominant and Severe long COVID genes are not an inherent artifact of the analytical approach. This result implies that some of the mechanisms of action and gene regulatory networks for these two diseases may also overlap with some of those important to ME. Given the pathologies and broad physiological scope of the mechanisms involved, it is reasonable to expect such overlap, and we have observed evidence for this in combinatorial analyses of these and other diseases (unpublished results). Fig. 8 SNP ranks in the ML model for ME for all SNPs assigned to autosomal protein-coding genes (left) and SNPs assigned to genes associated with Severe, Fatigue Dominant, and General long COVID phenotypes in the Sano GOLD cohort (right three boxplots). Higher SNP ranks (i.e. higher positions on the y-axis) denote high feature importance in the ML model. Boxes represent the top 25th and bottom 25th percentile ranks for SNPs assigned to each gene set with solid line denoting median SNP rank. p -values are the probability of randomly observing a top 25th percentile SNP rank that high based on 10,000 gene subsets of equivalent size randomly sampled from the full list of autosomal protein-coding genes SNP ranks in the ML model for ME for all SNPs assigned to autosomal protein-coding genes (left) and SNPs assigned to genes associated with Severe, Fatigue Dominant, and General long COVID phenotypes in the Sano GOLD cohort (right three boxplots). Higher SNP ranks (i.e. higher positions on the y-axis) denote high feature importance in the ML model. Boxes represent the top 25th and bottom 25th percentile ranks for SNPs assigned to each gene set with solid line denoting median SNP rank. p -values are the probability of randomly observing a top 25th percentile SNP rank that high based on 10,000 gene subsets of equivalent size randomly sampled from the full list of autosomal protein-coding genes This gene overlap provides further evidence of the highly polygenic and heterogeneous nature of ME and long COVID and their partly shared etiology [ 35 , 55 ]. 97 genes (out of 19,146 total annotated autosomal genes) were in the top 20 highest-ranked ME genes in at least one of the 10 folds. These 97 genes included 6 genes that were previously associated with long COVID in the Sano Gold cohort (see Supplementary Table 11 ). Three overlapping genes (i.e., found in both analyses), MAPK9 , NLGN1 , and PTPRD , were associated with Fatigue Dominant long COVID, which is significantly more overlap than expected under the null hypothesis ( p  < 0.001). Each of these genes is associated with a distinct subset of mechanisms previously hypothesized to play an important role in ME and long COVID [ 32 , 35 , 55 ]. MAPK9 , which encodes mitogen-activated protein kinase 9 (c-Jun N-terminal kinase 2, JNK2), is a stress-responsive kinase that regulates immediate-early gene expression through phosphorylation of key transcription factors [ 61 ]. It integrates inflammatory, oxidative, and metabolic signals—processes relevant to the pathophysiology of both ME and long COVID [ 62 ]. Consistent with this, multi-omics analyses in both conditions have identified MAPK9 as differentially regulated or as a component of disease-associated gene interaction networks [ 63 – 65 ], suggesting a potential role in shared disease biology. NLGN1 encodes a postsynaptic adhesion protein critical for glutamatergic synaptic function [ 66 ] and is transcriptionally regulated by CLOCK [ 67 ], linking it to circadian rhythm. NLGN1 variants and altered expression have been associated with ME [ 68 ] and with cognitive disorders such as Alzheimer’s disease [ 69 ]. These findings suggest that NLGN1 may contribute to shared symptoms of ME and long COVID, including cognitive impairment, sleep disruption, and fatigue, although further studies are needed to confirm this. PTPRD , a receptor-type tyrosine phosphatase, is another neuronal adhesion molecule highly expressed in the brain, where it contributes to synaptic specification and cognitive regulation [ 70 ]. Genetic variants in PTPRD have been reported across several human diseases with a neuroimmune component, including Alzheimer’s disease and ME [ 70 , 71 ]. Although direct evidence connecting PTPRD to long COVID is currently lacking, its neuronal functions and genetic association with ME make it a plausible candidate contributing to the fatigue and cognitive symptoms common across post-viral syndromes. Our analysis highlighted over 40 core candidate drug targets and many peripheral genes with existing therapeutic assets in other indications (see Extended Table  5 ), enabling systematic identification using a detailed workflow to evaluate these as potential drug repurposing opportunities in ME based on genetic and other evidence [ 72 ]. Some of these are known and subject of repurposing trials already, but many more are novel. For example, we identified a novel genetic associations for a core candidate gene target TLR3 (odds ratio 1.20, p  = 10 − 4 in the Refinement dataset), which is also linked to medications currently under investigation in ME clinical trials [ 60 , 73 , 74 ]. Toll-like receptor 3 ( TLR3 ) is a key component of the innate immune system that acts as a critical sensor for viral double-stranded RNA in several cell types that are key to host antiviral defense [ 75 ]. Beyond its role in detecting exogenous viral RNA, TLR3 also senses endogenous RNA released by damaged, necrotic, or stressed cells, thereby modulating inflammatory responses [ 76 ]. Dysregulated TLR3 signaling can lead to chronic inflammation and tissue damage, exacerbating conditions such as autoimmune diseases, chronic viral infections, and cancer [ 77 , 78 ]. Rintatolimod is a synthetic double-stranded RNA molecule that acts as a selective agonist of TLR3. On binding to TLR3, rintatolimod activates the MyD88 independent TRIF signaling pathway, leading to the production of interferons and other antiviral proteins without triggering excessive systemic inflammation associated with other dsRNA molecules [ 74 ]. It has been investigated in several Phase II/III clinical trials with ME patients, where it has shown statistically significant improvements in primary endpoint using exercise tolerance and some secondary endpoints when compared to placebo [ 60 , 73 , 74 ]. Despite showing encouraging results and having an acceptable safety profile, rintatolimod still remains an experimental drug in ME in the US, with its use limited to compassionate access programs. The disease signatures associated with TLR3 identified in a patient subgroup (24.1%) in Analysis 1 could help enhance future trial designs to aid higher efficacy and advance toward full regulatory approval. While no statistically significant phenotypes associated with the TLR3 patient subgroup were detected, modest associations were observed (see Supplementary Results), providing an illustrative example of phenotype characterization of patient subgroups. In addition, our findings revealed many other target–drug pairs involving safe, well-tolerated, generic medications. One such example is PDE4B and its modulator apremilast, associated with a subgroup of ME patients. Interestingly, our unpublished results suggest that this medication might be also effective in a subgroup of long COVID patients, which reflects shared disease pathology with ME. Our identification of PDE4B (covering 44 disease signatures with odds ratios 1.05–1.78 and minimum p -value 1 × 10 − 12 in the Refinement dataset) represents a novel genetic link between this gene and ME. We identified signatures associated with PDE4B in around 70% patients in both Analyses 1 and 2. PDE4B encodes phosphodiesterase 4B, an enzyme that hydrolyzes cyclic AMP (cAMP), thereby modulating intracellular signaling related to inflammation, immune activation, and cognition [ 79 , 80 ]. PDE4B is highly expressed in macrophages, lymphocytes, and CNS glial cells [ 81 , 82 ], and its dysregulation has been implicated in multiple inflammatory and neurological conditions [ 83 , 84 ]. We hypothesize that overactivation or upregulation of PDE4B in ME reduces intracellular cAMP levels, disinhibiting NF-κB and TNF-α signaling pathways [ 85 – 87 ]. This may amplify peripheral and central neuroinflammatory cascades, contributing to core ME symptoms such as fatigue, post-exertional malaise, and cognitive dysfunction. Interestingly, our analysis in long COVID patients identified a modest association with PDE4D , another member of the same enzyme family, suggesting a potential shared pathophysiological axis involving dysregulated cAMP signaling (unpublished results). Among available PDE4B modulators, we prioritized apremilast, an orally available small-molecule PDE4 inhibitor approved for psoriasis and other inflammatory disorders, as the most promising candidate for repurposing. Although non-selective across PDE4 isoforms, apremilast’s inhibition of both PDE4B and PDE4D could confer therapeutic benefits in both ME and long COVID. Its safety profile is favorable, with gastrointestinal side effects such as nausea and diarrhea being the most common adverse events. Collectively, these findings highlight apremilast as a promising repurposing candidate for some ME patients based on a novel genetic and mechanistic rationale. Future studies should assess its safety and efficacy in mechanistically stratified clinical trials, using biomarkers to identify the subgroup most likely to benefit.

Our analytical pipeline used four cohorts of ME cases and controls (Fig.  1 ): A DecodeME ‘Discovery’ dataset used to detect disease signatures via hypothesis-free combinatorial analysis of a cohort of DecodeME cases and UKB controls. A ‘UKB Refinement’ dataset used to filter disease signatures and component SNPs based on reproducibility of disease associations in an independent cohort of UKB cases and controls. A ‘DecodeME Refinement’ dataset used to filter disease signatures and component SNPs based on reproducibility of disease associations in a cohort of independent DecodeME cases and the same controls as the UKB Refinement dataset. A ‘Test’ dataset used to evaluate the predictivity of the final set of ‘double-refined’ signatures (from the UKB Refinement and DecodeME Refinement steps respectively) in an independent cohort of DecodeME cases and UKB controls not included in any other dataset. A DecodeME ‘Discovery’ dataset used to detect disease signatures via hypothesis-free combinatorial analysis of a cohort of DecodeME cases and UKB controls. A ‘UKB Refinement’ dataset used to filter disease signatures and component SNPs based on reproducibility of disease associations in an independent cohort of UKB cases and controls. A ‘DecodeME Refinement’ dataset used to filter disease signatures and component SNPs based on reproducibility of disease associations in a cohort of independent DecodeME cases and the same controls as the UKB Refinement dataset. A ‘Test’ dataset used to evaluate the predictivity of the final set of ‘double-refined’ signatures (from the UKB Refinement and DecodeME Refinement steps respectively) in an independent cohort of DecodeME cases and UKB controls not included in any other dataset. We began by identifying non-overlapping sets of UKB controls to use for Discovery, Refinement, and Testing. Control Cohort B: We first set aside 40,000 samples from the set of potential UKB controls to use for disease signature Discovery (i.e., combinatorial analyses). This number was selected to potentially allow for a minimum 1:2 case: control ratio in a combinatorial analysis of the full DecodeME cohort prior to knowing how many participants had consented to data sharing and passed QC. The set of Discovery controls was selected to be 83.5% female to match available data on the sex-ratio in the DecodeME cohort. Control Cohort C: To create a set of Refinement controls, we first randomly selected half the remaining female UKB potential controls (i.e., excluding Discovery controls). We then randomly selected remaining male UKB potential controls to create a cohort of UKB Refinement controls that is 54% female, matching the overall female: male ratio in UKB. For the UKB Refinement dataset (Step 2), we paired the UKB Refinement controls with ME cases identified using UKB’s pain questionnaire (as described in 32). Control Cohort D: Finally, we created a set of UKB Test controls by taking the remaining potential female UKB potential controls and merging them with randomly selected remaining male UKB potential controls to again produce a controls cohort that is 54% female, again matching the overall female: male ratio in UKB. We conducted two complementary analyses to fully utilize the available samples from DecodeME (Figs.  1 and 2 ), relying on three independent cohorts of DecodeME participants. Case Cohort F was comprised of ME participants from DecodeME Batch 1 (which was the first cohort available for analysis). Case Cohorts G and H were comprised of ME participants from DecodeME Batches 2 and 3 (which were made available after the initial combinatorial analysis had been performed). To mitigate any potential batch effects between Batches 2 and 3, we merged their QCed cases. We then re-split the cases in half to form Case Cohorts G and H, ensuring a consistent ratio of cases from Batch 2 vs. Batch 3 and females vs. males within each batch and case cohort. Case Cohort G was used for Discovery in Analysis 2 and Refinement in Analysis 1, while Case Cohort H was merged with the UKB Test controls (Cohort D). Fig. 1 Study cohorts used as sources of ME cases and controls in the disease signature Discovery, Refinement, and Testing pipeline and the cohort datasets used for Analysis 1 and Analysis 2 Study cohorts used as sources of ME cases and controls in the disease signature Discovery, Refinement, and Testing pipeline and the cohort datasets used for Analysis 1 and Analysis 2 Fig. 2 Disease signature Discovery, Refinement, and testing pipelines and cohort datasets used for Analysis 1 and Analysis 2 Disease signature Discovery, Refinement, and testing pipelines and cohort datasets used for Analysis 1 and Analysis 2 In Analysis 1, we performed combinatorial analysis using a Discovery dataset comprised of ME cases from DecodeME Batch 1 (Cases F) paired with the UKB Discovery controls (Controls B) (see Fig.  2 ). Refinement was first conducted in the UKB Refinement dataset (Cases E paired with Controls C) and then in a Refinement dataset comprised of the first DecodeME Batch 2 + 3 split (Cases G) paired with the UKB Refinement controls (Controls C). In Analysis 2, we inverted the study design from Analysis 1. First, we performed combinatorial analysis using a Discovery dataset comprised of the first DecodeME Batch 2 + 3 split (Cases G) paired with the UKB Discovery controls (Controls B). Refinement was first conducted in the UKB Refinement dataset and then in a Refinement dataset comprised of the DecodeME Batch 1 cases (Cases F) paired with the UKB Refinement controls (Controls C). Pooling the output of Analysis 1 and Analysis 2 provides a set of ‘double-refined’ disease signatures that are associated with ME across three cohorts (i.e., a DecodeME Discovery dataset, a UKB Refinement dataset, and a DecodeME Refinement dataset). Note that after QC, many SNPs that were present in the Analysis 1 Discovery dataset derived from DecodeME Batch 1 were missing from the Test Dataset and Analysis 1 Refinement datasets derived from DecodeME Batches 2 + 3. We therefore imputed these missing SNPs using the approach described in the Supplementary Methods to allow for refinement and validation of signatures containing them. We conducted standard quality control (QC) individually on each batch of participants provided by DecodeME and each merged dataset based on the protocol described in [ 38 ]. Details are provided in the Supplementary Methods. Population substructure within a cohort can result in false positive gene-disease associations when disease prevalence is correlated with patient ancestry, for example if ME is underdiagnosed in non-European genetic ancestry patients. In this scenario, genotypes that differ in frequency between historical human populations either due to selection or genetic drift can be indirectly associated with disease even though they do not directly affect disease biology. Due to the very high proportion of European ancestry participants and low representation of other ancestries, we filtered the DecodeME cohorts to include only participants classified as having European ancestry (see Supplementary Methods). This reduced the potential that the disease signatures associated with increased prevalence of ME in DecodeME are artifacts of population substructure rather than reflecting genetic features directly relevant to disease biology [ 42 ]. Merging DecodeME cases with UKB controls introduces the potential for batch effects to produce large numbers of false positive disease associations [ 43 – 45 ]. That is, differences in signature frequency between cases and controls may reflect technological artifacts in the data rather than differences in disease biology. Although DecodeME used the same genotyping array as UKB, other differences in the genotyping pipeline (e.g., in sample collection, sample preparation and storage, DNA extraction, machine calibration, genotype calling, or bioinformatic pipelines) could result in material differences in genotyping error rates and observed allele frequencies between the two datasets. Because case-control status is perfectly correlated with the dataset source in our data, identifying and removing these batch effects while retaining true genetic associations is challenging. Traditional methods such as those using principal component analysis (PCA) can only be applied when technological artifacts are expected to be the sole reason for differences in allele frequency between batches [ 46 ]. However, in this study, we expect that disease biology, differences in diagnosis and phenotyping, and/or indirect relationships between population substructure and disease prevalence will manifest as differences in allele frequencies between cases and controls. Standardizing the DecodeME case and UKB controls so that they have similar principal component (PC) distributions is therefore likely to remove the signal that we wish to identify in this study. We therefore undertook a series of extensive and highly conservative QC analyses to identify and remove false positive signal associated with differences in allele frequency between DecodeME and UKB and between DecodeME batches that reflect genotyping artifacts. These QC processes are described in detail in the Supplementary Methods. We applied a multi-step analytical pipeline to identify and refine disease signatures for ME using non-overlapping cohorts of ME cases and healthy controls (see Fig.  3 ). We began by using the PrecisionLife combinatorial analytics platform to identify disease signatures significantly enriched in ME cases relative to UKB controls in one of the two Discovery datasets. Fig. 3 Analytical pipeline for discovery and refinement of disease signatures and testing of signature count score for ME. The discovery and refinement datasets are illustrated in Fig.  1 . These processes are described further in the Supplementary Methods section Analytical pipeline for discovery and refinement of disease signatures and testing of signature count score for ME. The discovery and refinement datasets are illustrated in Fig.  1 . These processes are described further in the Supplementary Methods section After combinatorial analysis, we applied a process we call ‘Refinement’ to the set of disease signatures (see Fig.  3 and Supplementary Methods section for more detail). This entailed: Removing signatures and/or component SNP-genotypes that do not show reproducible associations with ME in a second cohort (‘strong components test’). Removing signatures that do not provide insight into disease associations that are independent of other signatures (‘redundant signature filter’). Expanding the set of genotypes reflected in the disease signatures to reflect standard genetic models when appropriate (‘expanded genotypes’). Removing signatures and/or component SNP-genotypes that do not show reproducible associations with ME in a second cohort (‘strong components test’). Removing signatures that do not provide insight into disease associations that are independent of other signatures (‘redundant signature filter’). Expanding the set of genotypes reflected in the disease signatures to reflect standard genetic models when appropriate (‘expanded genotypes’). The above pipeline resulted in a set of ‘double-refined’ disease signatures that were associated with increased odds of ME across three independent cohorts of participants (one Discovery + two Refinement cohorts). To evaluate the aggregate predictive power of these signatures, we tested whether the count of double-refined signatures possessed by a patient is correlated with prevalence of ME in the independent Test dataset, which was comprised of cases and controls not previously used for Discovery or Refinement in either analysis. For each patient, we calculated a signature count score reflecting a simple count of the number of signatures possessed by each patient. As noted in the Results section, we identified thousands of correlated disease signatures containing at least one SNP located within the 6p22.1 locus, which strongly influenced the signature count. We therefore excluded these signatures from the set of signatures evaluated by the signature count score. We then calculated a logistic regression using the glm function in R with CRS as the independent variable, case-control status as the dependent variable, and genetic sex and 10 genetic PCs as confounders to evaluate whether the signature count score is significantly associated with case-control status (see Supplementary Methods). In addition to controlling for the confounding effects of population substructure, inclusion of the PCs as covariates is also expected to distinguish true biological signal from systemic batch effects that were not flagged and removed during QC. We further sorted individuals in the Test dataset by signature count score and calculated the odds of ME (i.e., # cases / # controls) for each signature count score decile. Odds ratios were calculated by dividing the odds for each decile by the odds of ME for individuals falling within the bottom 10% of signature count score values. SNPs identified in the disease signatures were conservatively mapped to protein-coding genes based on the human reference genome (GRCh37). Only SNPs located within the transcription start and stop positions of genes were mapped to the corresponding gene(s). These genes were then used for overlap analyses with other studies. To aid in identifying the candidate core, potentially causal, genes and mechanisms linked to ME, we used the appropriate DecodeME Refinement datasets (C + F and C + G) to filter the mapped genes based on the strength of their associations with ME and the prevalence of associated disease signatures in ME patients. We then constructed networks for each SNP comprised of the set of all double-refined disease signatures that contained that SNP (i.e., the network’s defining SNP). Cases and controls were assigned to a network if they possessed at least one of the constituent disease signatures. This allowed us to calculate a network odds ratio reflecting the odds of ME for samples within and not within each network, and each network’s case prevalence (i.e., the fraction of total ME cases who were assigned to the network). We further calculated a SNP odds ratio reflecting the odds of ME for a network relative to the odds of ME for a matching network comprised of the same constituent signatures but with the defining SNP removed, along with a corresponding Fisher’s exact test p -value. This measure allows us to quantify the association with ME for individual SNPs in a combinatorial setting. Finally, we identified the set of SNPs that define networks with SNP odds ratio p- values less than 0.05 and a case prevalence greater than 20%. Genes mapped to these SNPs, and those with a minimum SNP network odds ratio exceeding the median value across all genes, were identified as the candidate core gene set. Supplementary Table 3 shows the varying SNP and candidate core gene counts at different prevalence thresholds. All genes were annotated using data from over 50 public data sources (see Supplementary Table 4 ) to characterize their biological roles and potential mechanism of action links to ME. We first performed pathway enrichment analysis on the candidate core genes using Gene Set Enrichment Analysis (GSEA) [ 47 ]. Additional functional characterization and tissue enrichment analyses were performed to better understand their functional context and potential relevance to disease biology [ 48 , 49 ]. See Supplementary Methods for detailed description of these analyses. Phenotypic variables were derived from selected responses in DecodeME questionnaires 1 and 2, encompassing symptoms, comorbidities and illness severity. A total of 39 features covering 17 symptoms and 22 comorbidities were generated. See Supplementary Methods for detailed description of phenotype features from DecodeME Questionnaire Data. The constructed phenotypic variables were used to evaluate patient profiles associated with the specific genes identified in the study. This large-scale phenotype association analysis was performed in the Test Dataset, which was comprised of DecodeME Batch 2 + 3 split 2 cases (Cohort H). The enrichment analyses were performed using Fisher’s Exact Test to compare phenotype distributions between gene-associated cases and the remaining cohort. Multiple testing correction was applied using the Benjamini–Hochberg method to control the false discovery rate (FDR). We tested for overlap between the genes identified in our hypothesis-free ME analyses and the genetic associations for long COVID that were previously identified by combinatorial analysis in the Sano GOLD cohort and which also exhibited high rates of reproducibility (92%) in an All of Us cohort [ 35 , 36 ]. These hypothesis-free studies identified genes associated with long COVID in three partially overlapping cohorts from the Sano GOLD dataset [ 35 ]: Patients who have a ‘Fatigue Dominant’ form of long COVID Patients who have a ‘Severe’ form of long COVID Patients who satisfy the ‘General’ criteria for long COVID. Patients who have a ‘Fatigue Dominant’ form of long COVID Patients who have a ‘Severe’ form of long COVID Patients who satisfy the ‘General’ criteria for long COVID. We used a Fisher’s exact test to determine whether we observed a significant enrichment of long COVID genes that also map to the double refined ME disease signatures. We also tested this overlap using an independent method – i.e., not relying directly on the PrecisionLife combinatorial analytics methodology, which may potentially be biased towards outputting larger genes that map to many SNPs on the Axiom UKB genotyping array. To do this, we identified the overlap between the long COVID genes reported in the hypothesis-free analyses above and the output of an independent machine learning (ML) analysis that quantified the relative feature importance of SNPs and genes in models that classify people into ME cases and controls. The accuracy and robustness of the ML model for ME is limited due to the large numbers of included SNPs and their relatively small feature importances. It would not be useful as a discovery tool in its own right, but it does allow identification of genes previously associated with long COVID that are also highly associated with ME. We performed two analyses using the ML model for ME. First, we ranked SNPs by mean feature weight in the ML model for ME and tested whether SNPs mapping to the long COVID genes had significantly stronger ranks (i.e., high feature importances) relative to comparable sets of random genes. To test for potential biases in the methodology, we contrasted the results of the overlap analysis for long COVID genes with the same analysis for genes identified by combinatorial analyses of endometriosis [ 33 ] and coronary artery disease (CAD) ( unpublished results , Supplementary Table 5 ). Second, we identified genes mapping to a small set of SNPs with very high feature importances in at least one of the ten folds of the ML model for ME, then tested whether we observed significantly enriched overlap between these and the long COVID genes. Details of these analyses and the ML model are included in the Supplementary Methods. A PPIE advisory panel of individuals with lived experience of ME and/or long COVID and patient advocates was established in collaboration with the charity Action for ME. Members were recruited through the charity’s networks and reimbursed according to NIHR payment guidelines. The panel met virtually 25 times during the LOCOME project and contributed particularly to both the planning and operational aspects of the study reviewing progress on a monthly basis ensuring clarity in planning future work. In addition, they brought unique perspectives and critical thinking from lived experience and extensive knowledge from the wider patient and academic community, improving analysis and written communications of the results.

DecodeME provided genotype data for 14,767 participants who satisfied the study’s inclusion criteria based on the Canadian Consensus and/or US Institute of Medicine / National Academy of Medicine criteria for ME/CFS using survey responses, who had received a diagnosis of ME, CFS, ME/CFS or CFS/ME from a health professional, and who consented to data sharing with accredited research teams. DecodeME participants were genotyped in four batches [ 30 ], but only the first three were included in this analysis due to the timing of data availability. Although all were genotyped using the Axiom UK Biobank array, non-trivial variations in data generation and quality control (QC) resulted in major differences in SNP counts between the DecodeME batches, and a substantial reduction in the number of participants (cases) and SNPs that could be used in the study. The extensive QC processes that we undertook to remove batch and population structure effects and other data quality issues are described in the Supplementary Methods section. In these batches between 82 and 86% of DecodeME cases and 53–59% of the DecodeME SNPs were able to be used for this analysis. Total counts of cases and SNPs for each batch used in this study are listed in Table  1 below. The number of females and males in the cohorts are listed in Supplementary Table 1 . Table 1 Number of cases and SNPs from the DecodeME and UKB cohorts before and after LOCOME Study QC, broken down by batch ME Cohort # Cases before QC # SNPs before QC # Cases after Study QC (%) # SNPs after Study QC (%) DecodeME Batch 1 3,944 772,374 3,405 (86%) 408,873 (53%) DecodeME Batch 2 4,349 656,807 3,677 (85%) 362,798 (55%) DecodeME Batch 3 4,238 645,542 3,487 (82%) 381,756 (59%) DecodeME Batch 2 + 3 Split 1 -- 3,585 349,734 DecodeME Batch 2 + 3 Split 2 (Test) -- 3,579 350,014 UK Biobank – Pain Questionnaire 2,382 554,063 1,985 (83%) 519,337 (94%) Number of cases and SNPs from the DecodeME and UKB cohorts before and after LOCOME Study QC, broken down by batch UK Biobank – Pain Questionnaire The DecodeME dataset is restricted to participants who satisfy the clinical case criteria for ME, CFS, ME/CFS or CFS/ME [ 29 , 30 ]. We therefore relied on UKB as a source of controls for the genetic association analyses. We used the same criteria as [ 31 ] to identify potential UKB controls, i.e., participants with no evidence in UKB’s Hospital Episodes Statistics (HES), primary care, or self-reported data fields indicating diagnoses of chronic fatigue, post-exertional malaise, post-viral fatigue syndrome or myalgia (Supplementary Table 2 ). We further excluded all UKB participants previously used as controls in the combinatorial analyses described in [ 31 ]. Ensuring we only used non-overlapping controls allowed us to use the previously identified UKB ME dataset for independently assessing the reproducibility of the output of the DecodeME combinatorial analysis.

Despite ME and long COVID’s massive socioeconomic impacts, unmet medical needs, and profound personal and public health consequences [ 11 , 88 ], there are still no treatments that address the root causes of these diseases. There are many scientific and technical reasons for this. The biology of the diseases is complex, involving multiple causes and factors, and a wide range of symptoms. They lack clear diagnostic criteria, an easily measurable disease progression biomarker, and recognized, reliable in vitro/in vivo (cellular or animal) models that accurately represent the diseases’ heterogeneous pathology [ 89 , 90 ]. All of these issues add materially to the uncertainty, time, and risk associated with initiating drug discovery efforts in a new disease indication, and they have substantially affected pharmaceutical companies’ appetite to invest in drug development program for ME and long COVID. This study builds on and confirms findings from a series of studies that provide evidence of significant novel reproducible genetic associations with ME and long COVID, and the partial overlap between them. Even reproducible genetic associations are not however automatically causal or druggable, nor do they all translate equally into clinically useful findings [ 91 , 92 ]. This means that a series of functional experiments and interventional clinical trials will be required for establishing the causality and verifying the druggability of the genes identified. These results may however help to guide, prioritize and accelerate some of those studies, hopefully leading to evidence of clinical utility that encourages new pharmaceutical and diagnostic research efforts in the field. There are several clinically relevant applications of the results, which could potentially help improve diagnosis and/or differential triage of the diseases, build better risk models for the diseases and their key mechanisms/symptoms, and find new therapeutic options personalized to patient subgroups, especially via drug repurposing. Outside of the oncology and rare disease space, genetic tests are not usually definitively diagnostic as they describe a person’s lifetime risk of getting a disease rather than their current disease status. In diseases that have environmental or epidemiological components (such as lifestyle or comorbidities) or which are infection mediated, these other factors (e.g., exposure to a triggering infection) play a significant role alongside genetics. A person’s genetics may however significantly determine their susceptibility and response to infection, their propensity to develop sequelae that persist after the acute phase, which symptoms these may lead to and their severities, and the drugs to which they may respond best (and worst). They may also confer disease resilience in some people. As predictors of personal risk, and in conjunction with clinical symptoms, genetic tests can therefore be very useful clinical tools, especially in complex diseases such as ME or long COVID that have multiple symptoms that could be overlooked or mistaken for other conditions. As shown in this and previous studies, there is a significant overlap of ME disease genes with long COVID, which in some patients may exhibit similar symptoms [ 38 , 39 ]. Long COVID and ME clearly share multiple genetic and mechanistic commonalities based on the results of this study and others [ 13 , 35 , 40 , 55 ] and there likely exist some drugs that will modulate shared targets, and which may therefore prove to be effective for some patients from both disease groups. Many long COVID genes are, however, not biologically associated with ME and vice versa. These results suggest that the diseases are best considered to be partially overlapping but different, consistent with the findings of studies of symptomatic overlap and differentiation [ 93 ] . Similar patterns of significant partial overlap with symptoms and with sets of ME genes have been observed in combinatorial analysis studies of other diseases with a neuroimmune component including fibromyalgia, multiple sclerosis, rheumatoid arthritis, endometriosis, IBD, systemic lupus erythematosus, and hypersomnia (unpublished results). These diseases are typically treated by different clinical specialists who may not be able to differentiate the causes of the patient’s symptoms accurately. When fully developed and validated, the insights into genetic similarities and differences between these diseases could be used to develop a test for the rapid non-invasive differential triage of patients presenting with non-specific neuroimmune symptoms. From a single sample, the disease specific risk and resilience models for multiple conditions could be evaluated to identify which disease and mechanisms are most likely contributing to a patient’s symptoms. In principle, and once validated, such a test could help facilitate quick and accurate referral of patients to the correct disease specialist to achieve a definitive diagnosis. With further validation it may also evaluate a patient’s prognosis for symptoms and severity, and help clinicians select the most appropriate therapies for that patient, both during clinical trial recruitment and when drugs are approved for use in the clinic. Knowing that a patient with a potentially triggering infection such as SARS-CoV-2 or EBV has a high genetic risk of developing long COVID or ME would provide useful information to guide their clinical care. Identifying high risk patients might trigger more proactive methods to avoid reinfection via vaccination boosters [ 94 ], or to limit the initial infection using anti-viral therapy [ 95 , 96 ] and avoidance of post-infectious overexertion via pacing once unwell [ 97 ]. However, building clinically useful risk scores for diseases like long COVID and ME is challenging. The strong correlation between genetics and ME is notable but will always be limited as genes are likely only one contributing factor towards whether an individual develops ME. Such non-genetic factors are expected to weaken the correlation between genetics and disease status. For example, some ‘high-risk’ controls likely have higher genetic predisposition to developing ME but may not have been exposed to sufficient triggers (e.g., viral, bacterial, and/or environmental) to actually develop the condition. Alternatively, such patients may possess ‘actively protective’ disease signatures that confer resilience and wholly or partially mitigate the effects of their ME disease risk signatures (as described below). We expect that combinatorial analysis of additional cohorts of ME patients, especially more diverse non-European cohorts, would provide additional insight and improve our ability to accurately quantify patients’ genetic risk. Any future risk score should also take into account the hypothesis that ME likely represents several distinct mechanistic subtypes, each with its own disease etiology. For example, a patient might have relatively few disease signatures in total, but a significant excess of signatures associated with a single mechanism that is primarily responsible for their disease subtype. Detailed clinical phenotype information is available in DecodeME, and further studies to investigate the genetic and mechanistic underpinnings of specific phenotypes are underway [ 98 ]. Incorporating mechanistic stratification and observed clinical phenotypes into a risk and resilience score framework can potentially provide improved insights not only for quantifying risk, but also for understanding the likely drivers of a patient’s disease and aiding in the selection of precision medicine treatments most likely to be effective for them. Given the lack of tools and appetite for novel drug development, the ME community has widely recognized that traditional drug discovery and development routes are going to be slow at best. Testing safe and well-tolerated compounds directly in humans through proof-of-concept clinical studies may therefore provide a pragmatic alternate route to generate critical evidence that modulation of a novel target in a specific direction can yield a disease modifying therapeutic effect [ 99 ]. Many repurposing studies are already underway including adaptive trial designs and combination therapy trials [ 100 – 103 ]. They are however challenged by the heterogeneity of the disease, which limits the proportion of responders and the degree of efficacy that can be demonstrated in a non-stratified trial population. Access to drugs, especially development candidates and on market compounds, has also proved difficult in some cases [ 104 ]. This study’s results have the potential to extend and enable the use of precision drug repurposing trials. Firstly, the range of candidate core genes identified offers more repurposing targets – a few of these have been highlighted above and there are dozens more potential candidates that could be the subject of future validation studies (Extended Tables). Secondly, genetic stratification of the disease provides an opportunity to enrich potential responders in the trial cohort – i.e., identify an individual patient’s disease drivers (and resilience factors) and therefore predict whether a therapy is likely to be effective (or not) for them. Good inclusion and exclusion criteria are equally important. People with ME and long COVID are often sensitive to medications and adverse effects, and enrolment in trials carries a risk of making their symptoms worse. Better targeted trials would allow selection of only those participants more likely to respond, avoiding those who may be at high risk of worsening of their condition. Precision drug repurposing has the potential to be the most rapid and cost-effective route to delivering effective therapies to specific groups of patients [ 105 ]. Demonstrating such efficacy in a targeted patient cohort would also provide essential clinical validation of the novel target’s druggability and the ability to select likely responders. This would overcome many of the substantial gaps in the translational medicine toolchest for the diseases and remove a key risk that prevents biopharma companies from investing further research and development funds. Generic small molecule drugs (such as apremilast highlighted above) have certain advantages for repurposing – they are easy and cheap to access, have well-known safe formulations, delivery routes, and dosages, and are easy to administer. Initial studies would likely be randomized open-label, proof of concept studies designed to demonstrate disease modifying potential, dose and safety in small genetically targeted cohorts of 30–50 people [ 106 ]. Participants with comorbidity risks that could be affected by the therapy, e.g., excessive hypotension with a blood pressure lowering drug, or other contraindications would be excluded. If successful, follow-up studies could progress to more sophisticated double-blind randomized controlled trials, although for these it will be essential to define the specific phenotypic endpoint clearly and evaluate it consistently [ 107 ]. Core symptoms like fatigue, pain, and cognitive dysfunction are subjective and highly variable day‑to‑day, making them vulnerable to placebo, nocebo, and expectation effects even under blinding. Precision mechanism-based clinical trial recruitment tools could improve the probability of clinical trials’ success, both for repurposing and traditional drug development. They could also inform a precision regulatory strategy that may make their clinical development faster and more likely to succeed with fewer adverse events. The same tools may become complementary diagnostics to guide therapy selection in the clinic in due course. The study used genotype data collected by both UKB and DecodeME in the form of a blood/saliva sample run on a ThermoFisher Axiom UKB genotyping array. This same genotyping platform or an equivalent low-pass whole genome sequence could be reused to form the basis of a new non-invasive test whose use to identify the patients most likely to benefit from a specific intervention during clinical trial recruitment is described above. Potentially (once validated and as part of a new clinical care pathway), such tests could also inform the selection of the therapies most likely to help an individual patient, once their own clinical utility has been established. In a similar vein, such a test could be used alongside the lived experience and reported outcomes of patients on studies such as TreatME, which collated thousands of patient’s responses to over 150 treatments [ 108 ], to further predict a patient’s personal response to a wider range of treatment options. The study’s findings explain significantly more of the disease’s heritable components than have previously been identified, and also why it has historically been so difficult to find or replicate genetic associations for ME using GWAS. The high observed polygenicity implies that the genetics of ME risk factors are quite variable between individuals, i.e., that every ME patient likely has many genetic variants associated with the disease, but that the sets of genetic variants that are responsible for their personal elevated risk of ME differ between most patients. This likely also explains why conventional genetic association studies have struggled to identify replicable disease associations. Each additional causal genetic variant effectively acts as a confounding factor that masks the disease association for any single genetic variant, resulting in very small mean effect sizes that are challenging to statistically validate in available datasets. This does not imply that it is not possible to identify effective treatments that apply to large numbers of ME patients. Many of the genes identified in this study show significant enrichment for druggable mechanisms such as neurological dysregulation, inflammation, cellular stress responses and calcium signaling (Supplementary Fig.  7 ) that have been implicated as playing an important role in ME disease biology for many patients. Focusing on the mechanisms and cellular pathways associated with the highest prevalence reproducible disease signatures can help prioritize testing of potential treatments such as drug repurposing candidates aiming to treat ME (and some long COVID) patients. High polygenicity does however strongly imply that understanding the specific sets of variants at play in an individual is going to be crucial for improving diagnostic tests and the success of clinical trials for new or repurposed medicines via targeted recruitment. It suggests a commonality in disease biology across patients that is consistent with a quasi-‘omnigenic’ model of disease biology. The omnigenic model of disease hypothesizes that many complex traits such as chronic disease are ultimately directly governed by expression of a relatively small set of ‘core’ genes [ 109 – 111 ]. However, identifying those core causal genes via genetic association studies is technically challenging because their expression in the population is directly or indirectly affected by complex gene regulatory networks incorporating a very large number of ‘peripheral’ upstream genes and genetic variants as well as environmental factors [ 112 ]. A genetic variant mapped to a core gene may increase susceptibility to ME by directly altering expression of that gene, but the effect of that variant in some individuals can be wholly or partially offset by the many other variants throughout the genome that also affect gene expression of the core gene indirectly via altered expression of ‘peripheral’ genes. As a hypothetical, simplified example, assume that a disease state results from overexpression of a single core gene beyond a threshold. In this example, a SNP in a cis -regulatory region that directly increases the expression of the core gene will cause some people with the SNP to develop the disease. However, expression of the core gene is also affected by several hormones such as testosterone and estrogen. In some people, levels of circulating hormones decrease the baseline expression of the core gene such that the overall expression remains below the disease threshold despite the increase in expression from the cis -regulatory SNP. These people will not develop disease even though they have the ‘causal’ SNP. Likewise, many people who do not possess the cis- regulatory SNP could nonetheless develop disease due to hormone-mediated increased expression of the core gene. The genetics of the hypothetical disease would be highly complex as any genetic variant or environmental exposure that directly or indirectly alters levels of circulating hormones could potentially cause expression of the core gene to rise above or drop below the disease threshold. Yet even though it is highly polygenic, it would be relatively simple to effectively treat the disease using therapies that reduce expression of the core gene to levels below the disease threshold. We hypothesize that ME exhibits similar omnigenic characteristics, and our findings of 259 candidate core genes with an additional set of 2,052 ‘peripheral’ genes that are also reproducible across 3 ME populations is consistent with this view of the disease. This is not an insurmountable complexity problem for study of the disease. It is possible to stratify patients with complex highly polygenic diseases into subtypes characterized by shared disease etiology, e.g., involving misexpression of key metabolic, neurological, or immune response genes, even if the set of peripheral genetic variants driving core causal gene expression differs between patients [ 113 ]. By representing the subcomponents and non-linear interaction effects of larger gene regulatory networks, combinatorial disease signatures more accurately reflect the omnigenic model of disease biology versus the GWAS and PRS approaches that ignore genetic interactions. Among the core candidate genes identified in our analysis, several are implicated in immune dysregulation and impaired energy metabolism – two key biological mechanisms thought to underlie ME pathophysiology [ 21 , 114 ]. Immune dysregulation associated with ME is characterized by aberrant immune activation and altered cytokine profiles that sustain inflammatory signaling [ 114 ]. In this context, we identified several core candidate genes that have established roles in immune signaling and inflammation, including TLR3 , CD8B , CD22 , and PSMB9 [ 115 – 118 ]. As described in the Results section, TLR3 is an innate immune receptor that detects both viral and damage-associated double-stranded RNA, and its dysregulation can drive chronic inflammation and tissue damage across a range of diseases [ 75 – 78 ]. Experimental modulation of TLR3 influences pro-inflammatory cytokine production and downstream signaling through pathways such as NF-κB and STAT1 [ 116 , 119 ]. Placebo-controlled randomized clinical trials of the TLR3 agonist rintatolimod reported improvements in clinical outcomes in ME patients [ 60 , 74 , 75 ]. CD8⁺ T‑cell dysfunction appears important in ME pathophysiology. High‑dimensional profiling shows that CD8⁺ T‑cell subsets are among the most dysregulated immune populations in ME, with signatures of transcriptional and epigenetic reprogramming toward terminal exhaustion and altered metabolism [ 120 ]. PSMB9 encodes an inducible β‑subunit of the immunoproteasome, which replaces a constitutive proteasome subunit under interferon‑γ and other inflammatory signals. PSMB9 is upregulated as part of a proteostasis defense program under mitochondrial dysfunction, boosting proteasome activity to limit protein aggregation [ 65 , 121 ]. Conversely, experimental depletion of PSMB9 in mice reduces pathological inflammation [ 122 ]. Plasma/immune‑cell proteomics in ME and long COVID cohorts emphasize shared enrichment of pathways of antigen presentation, interferon signaling, mitochondrial dysfunction, and proteostasis, all processes in which the immunoproteasome (including PSMB9) is central. We also identified a genetic signal in CD22 , a negative regulator of TLR3 signaling and B-cell receptor activity that functions to constrain inflammatory responses [ 117 ]. Experimental depletion of CD22 enhances responsiveness to TLR3 ligands in B cells and promotes a pro-inflammatory transcriptional state in macrophages, providing functional evidence for CD22 signaling in immune activation [ 117 , 123 ]. Impaired energy metabolism encompassing mitochondrial dysfunction, altered substrate utilization, and metabolic disturbances associated with fatigue and exercise intolerance is also implicated in ME [ 114 ]. We identified several metabolism-related core candidate genes, including PKM , IGF1R , TBC1D5 , and ABCA1 [ 124 – 127 ]. Notably, PKM encodes pyruvate kinase M1/2 (PKM1/2), which catalyzes the final, rate-limiting step of glycolysis and plays a role in regulating metabolic flux between glycolytic and oxidative pathways [ 124 ]. Experimental studies demonstrate that loss of PKM2 in endothelial cells alters mitochondrial substrate utilization, disrupts tricarboxylic acid (TCA) cycle activity, and induces innate immune signaling [ 128 ]. Moreover, inhibition of PKM2 in macrophages reduces glycolytic flux and attenuates inflammasome activation, further highlighting the link between cellular metabolism and inflammatory responses [ 129 ]. Interestingly, PKM2 has been shown to bind and stabilize the insulin-like growth factor 1 receptor ( IGF1R ), whereas PKM2 knockdown reduces IGF1R stability [ 130 ]. Conversely, IGF1R signaling influences PKM2 activity through AKT-dependent phosphorylation, suggesting mutual regulation between these pathways [ 131 ]. IGF1R is a key regulator of glucose metabolism and systemic energy homeostasis, acting through effects on insulin sensitivity, glucose uptake, and hepatic glucose production [ 125 ]. In line with this role, genetic or pharmacological modulation of IGF1R alters energy expenditure, metabolic adaptation to caloric restriction, and whole-body glucose metabolism in vivo [ 132 , 133 ]. Overall, our findings highlight genetic associations within pathways central to ME biology, particularly immune regulation and energy metabolism, providing mechanistic context for their potential contribution to disease pathophysiology. We can use this deeper understanding of genetic disease risk factors to find people who do not have any form of a disease even though they have a very high level of genetic disease risk factors, have experienced many of the potential disease triggers, and are as old as possible – so they’ve had every chance to contract the disease and be diagnosed. We can call these people ‘protected’ or ‘resilient’, at least for this specific disease (as such protection does not necessarily extend to other unrelated diseases). We have every reason to expect that these protected people should have the disease, so we can speculate that certain biological factors are working behind the scenes to prevent them from showing symptoms. We can find out what makes protected people different by comparing which genes and mechanisms show up significantly more often in them than in diagnosed ME patients. This is an enriched reversal of the usual way we find disease risk genes. In our UKB ME study, we found many protective signatures, mapping to 9 protein-coding genes that we believe merit further study as they may have a protective effect resisting the onset and progression of ME symptoms [ 37 ]. Many of these genes’ functions are consistent with alleviating the ME disease risk mechanisms also identified in this study, e.g., in insulin signaling, stress response and autoimmunity. Some of these have also been associated with mechanisms that are subject of approved drugs. These new ‘actively protective’ or disease resilience genes represent a completely new class of potential drug targets. In the future these may represent opportunities for new or repurposed therapies that could have a prophylactic benefit for many people who have high risk of a specific form of a disease. They could work to reduce the disruption to normal metabolism that causes disease – analogous to statins mitigating cardiovascular risk. Switching on or turning up these beneficial protective mechanisms may well benefit a much wider range of people than just those who have a specific form of the disease. The findings presented represent a preliminary combinatorial analysis of the DecodeME dataset. The studies have been performed on just 3 of the 4 batches of DecodeME data and run in just two small case-control cohort studies as opposed to a complete pooled set. There are several other new and confirmatory studies that we would like to perform. To date neither this study nor the DecodeME study has identified any sex-specific differences between genes associated with ME and we would like to investigate this and associations with HLA region genes in more detail. We have not yet performed actively protective analysis using DecodeME, which we plan to do when the whole dataset has been QCed, and we will also be looking for associations with specific symptoms and/or disease severity. Finally, formal confirmation that the double-refined signatures are broadly associated with disease (e.g., using the signature count approach) will require an independent cohort of people with ME and controls that are genotyped as part of the same study. This will minimize potential artifacts of batch effects that may produce differences in disease signature frequency between cases and controls that are unrelated to disease biology. Mis-phenotyping of the cases and controls used due to unreliable or missed diagnosis reduces the number of features detected and the level of reproducibility that can be observed between cohorts (see Supplementary Methods for a fuller explanation). Similarly, we expect genotyping errors, including errors introduced by SNP imputation, would cause reduced reproducibility of signatures, resulting in improper removal or retention of disease signatures and component SNP-genotypes during the Refinement pipeline, as well as reduced predictivity for the set of disease signatures in the Test dataset. The DecodeME cohort represents the largest and most accurately phenotyped collection of genomic data for ME participants available to date. However, it only contains case samples, and therefore a secondary dataset (UKB) comprised of controls who do not have ME was required to identify disease signatures associated with ME. Although the DecodeME and UKB datasets used for cases and controls used the same Axiom genotyping array [ 134 ], material differences in sample collection and genotyping between the two studies introduced source-specific batch effect artifacts that are inseparable from the case-control phenotype of interest. With the DecodeME team, we conducted extensive and robust QC procedures to remove all SNPs that showed patterns consistent with batch effects. Only 53–59% of the DecodeME SNPs survived this QC and this was further reduced due to lack of overlap between the two sets. However, it is possible that our QC pipeline did not flag all such batch effects, especially those that result in weak disease signal, resulting in potential false positives among the list of identified disease signatures. Conducting secondary screening in UKB mitigates this issue as artifactual signal from lingering batch effects shouldn’t be replicated in the cohort of UKB cases and controls. We believe that any remaining effects are therefore minimal. However, the highly conservative QC approach used to identify and remove potential batch effects and LD haplotypes may have filtered biologically important SNPs from the dataset as well. We opted to remove any potentially mis-genotyped SNPs as we believed that the potential scientific consequences of Type I error (false positives) in this study are more severe than Type II error (false negatives). The Refinement pipeline that was used to identify and remove potential false positive SNPs and signatures is also highly conservative and likely prone to Type II error. For example, a SNP-genotype is only retained in a signature if it is associated with increased odds of ME in the full Refinement dataset and all five of the k-fold subsets. However, many signatures are relatively rare and the 95% confidence intervals associated with the observed odds ratios often overlap due to random sampling effects. As such, it is expected that some true biologically important signatures and SNPs will be removed because random sampling variance makes it appear that they are not associated with increased disease risk in one or more of the comparisons. We again elected to accept this high rate of potential Type II error as a trade-off for minimizing potential Type I error and increased confidence in the final set of reproducible signatures and SNPs. Finally, we applied a conservative gene mapping approach which only assigns SNPs to a gene if they lie within or very close to a gene. This prevents us from identifying any upstream or downstream genetic variants that affect phenotype by altering gene or protein expression. It is likely that these include many of the unannotated SNPs from the double-refined signatures that were not assigned to genes. Alternatively, it is possible that some SNPs located within gene bodies affect phenotype via altered expression of a different gene, or that identified SNPs assigned to one gene are instead in linkage disequilibrium with a true causal SNP mapping to a different gene. The combinatorial analyses presented in this paper relied exclusively on patients with European ancestry, which despite efforts to the contrary, comprise most (93.1%) of the DecodeME study subjects. While not ideal, the exclusion of non-European ancestries from the study design was intended to reduce the likelihood of identifying ‘false-positive’ signatures. These may differ in frequency between cases and controls because they occur at different background frequencies in populations with differing rates of ME diagnosis rather than due to a direct causal link to disease biology. Unfortunately, populations with other ancestries were not sufficiently well-represented in DecodeME to enable separate statistically robust analyses. We note, however, that results from previous combinatorial analyses of long COVID and endometriosis demonstrated broad reproducibility across self-identified Hispanic/Latino and Black/African-American cohorts in All of Us, but sometimes at lower rates than for self-identified White European cohorts [ 33 , 36 ]. In part this high reproducibility results from the disease signatures being comprised of combinations of relatively common variants that are represented more consistently across multiple ancestries than rarer variants. These observations from other diseases suggest that a high proportion of the insights into ME from this study may also be applicable to non-European patients, but also that there may be some ancestry specific effects. This hypothesis remains to be proven. In the omnigenic disease model, effect sizes of individual SNPs are expected to vary between groups with different ancestries due to population-level differences in the frequency of confounding peripheral genetic variants. This pattern has been proposed to explain the poor predictive performance of polygenic risk scores across ancestries [ 135 ]. Combinatorial analyses of non-European ancestry populations may therefore provide insight into novel genes and disease signatures that are associated with ME in all populations but that have relatively greater importance in some patient subgroups due to the complex combinatorial nature of ME. Given the above concerns, development of accurate precision medicine tools like those described above requires confirmation that specific genetic associations should ideally reproduce well across diverse populations of ME and long COVID patients. Combinatorial analyses and other studies of gene function in non-European populations would be particularly helpful. We therefore emphasis the urgent need for the development of well-characterized, publicly available genomic datasets of non-European ME and long COVID cohorts.

The findings of this study substantially increase the number of genetic associations with ME, providing further evidence that ME is a complex multisystemic disease. This is an important step forward in scientific and clinical understanding of the disease and for the patient community who have been overlooked for so long. ME has been shown to be a highly polygenic and heterogenous disease, with a set of over 250 candidate core disease risk genes across a variety of mechanisms. This has significant implications for research strategy and the opportunity to develop more personalized clinical pathways. Stratification of patients by the mechanisms driving their form of the disease will be critical to predicting which patients will benefit from which therapy. It will also provide a crucial tool for developing more accurate diagnostic tests, creating patient selection tools that can make clinical trials more successful, and identifying novel and repurposed drugs. These results confirm the findings from other DecodeME based studies and its role as a key enabler of research into ME. While not without challenges as a dataset, the utility of DecodeME and combinatorial analytics as a research tool for ME and long COVID has been clearly demonstrated, and we hope this will encourage the development of new data collections and studies with more participants, including from diverse ancestries, with more detailed longitudinal clinical histories and deeper multiomic characterization of participants across multiple time points. These are, however, only preliminary findings and there are several on-going analyses to answer more specific questions. We are now beginning to uncover ME targets and tools with the precision required for more coordinated drug research. We hope that this will stimulate research across the community and encourage people with ME and long COVID to continue to contribute their samples, data and time to future studies.

Myalgic encephalomyelitis (also known as ME/CFS or simply ME) is a complex, chronic disease characterized by post-exertional malaise (PEM, sometimes referred to as post-exertional neuroimmune exhaustion PENE [ 1 ] in which symptoms disproportionately worsen, or arise, following minimal physical or mental exertion relative to pre-sickness). The disease also has neurological components (e.g., unrefreshing sleep, pain, neurocognitive impairment, sensory disturbance), immune, gastro-intestinal and/or genitourinary impairment, and impairments to energy metabolism/ion transport. Patients may experience a wide spectrum of other symptoms and comorbidities affecting multiple body systems, including dysautonomia, orthostatic intolerance and postural tachycardia, fibromyalgia, IBS, clinical depression, mast cell disease, and connective tissue differences. ME has a pronounced sex bias, affecting females 3–4 times more frequently than males. Onset of the condition is often linked to an infection (e.g., EBV, enteroviruses, SARS-COV) or other trigger [ 2 ]. Many people have gone on to meet the ME diagnostic criteria after developing long COVID [ 3 ]. Approximately 25% of ME patients have symptoms so severe that they are left house- or bed-bound [ 4 ]. Despite its debilitating impact on the lives of 25–70 million patients worldwide [ 5 ] and its huge societal and economic cost, estimated at over £20B in the UK [ 6 ] and at least $18B - $51B in the US [ 7 ], both scientific research and effective care pathways for ME have been notably lacking [ 8 – 10 ]. This means that the underlying biological mechanisms responsible for ME are still not well understood, posing significant challenges in improving diagnosis and finding effective treatments for patients. Long COVID is a similarly heterogeneous condition that is estimated to have affected 400 million people worldwide and to cost $1 trillion – 1% of global GDP – annually in healthcare costs and lost productivity [ 11 ]. The World Health Organization [ 12 ] defines long COVID as the continuation or development of new symptoms three months after SARS-CoV-2 infection, that last for at least two months with no other explanation. Common long COVID symptoms include fatigue, post-exertional malaise, autonomic dysfunction, cognitive dysfunction and shortness of breath, but over 200 symptoms have been reported that can impact daily life. Genetic studies of ME have yielded few reproducible results to date [ 13 ]. In part this may be related to the broad diagnostic criteria used for case populations, although there may be other explanations. Early studies of CFS (and not specifically ME) found that risk of developing the condition is increased among people who have family members with CFS, suggesting a possible genetic component [ 14 – 16 ]. Estimates of narrow-sense heritability of ME and/or CFS vary greatly however, ranging from 3% to 48% [ 17 – 19 ]. Early genome-wide association studies (GWAS) did not find any replicable genetic loci that were significantly associated with increased or decreased prevalence of ME [ 20 ]. This is likely due to the limited statistical power provided by available datasets as well as the multi-factorial basis of ME [ 21 ]. Of five ME GWAS studies that used data from UK Biobank (UKB) [ 22 ], two found no significant loci [ 23 , 24 ], two found one significant variant each [ 18 , 25 ], and one identified a female-specific and male-specific significant variant ( http://www.nealelab.is/uk-biobank/ ). All five of these studies aimed to identify genetic associations across large numbers of phenotypes in UKB and therefore did not attempt to tailor the case-control criteria specifically for ME. None of the four reported variants, two of which are very rare (minor allele frequency < 0.5%) and another of which is multi-allelic and may represent a genotyping artifact, were able to be replicated across multiple studies even though they all relied on the same base dataset [ 20 ]. Two additional GWAS analyses that relied on extremely small and underpowered datasets (< 50 cases) reported identification of genetic variants associated with ME. In one, none of the reported genetic associations remained significant after applying standard statistical thresholds for GWAS [ 26 ]. The second described 15 variants that were reported to have GWAS-significant disease associations ( p  < 5 × 10 − 8 ) [ 27 ]. The latter is very surprising given the very small dataset size and relatively low heritability of ME, and it is unclear if there were appropriate controls in place for confounding signal of relatedness or population substructure, both of which can result in spurious genetic associations [ 28 ]. Again, none of these purported ME disease associations have yet been replicated in other datasets. DecodeME is a genetic study that recruited 21,620 UK participants over the age of 16 who reported a clinical diagnosis of ME by a healthcare professional [ 29 ]. DecodeME applied rigorous screening criteria to participants’ questionnaire data to identify participants who pass international consensus criteria for ME and who lack any alternative diagnoses that might be responsible for post-exertional malaise, a unique diagnostic feature of ME [ 30 ]. Participants provided self-collected saliva samples that were genotyped using the UKB Axiom array. DecodeME’s recent pre-print of an initial GWAS comparing 15,579 DecodeME cases with European ancestry and 259,909 UKB controls reported 8 loci significantly associated with ME in that cohort [ 30 ]. Annotation of these loci implicated both immunological and neurological processes in ME disease biology. None of these loci were significantly associated with ME-like traits in a replication study of 15,251 cases and 1,878,066 controls assembled across seven independent biobanks. Four of the loci were associated ( p  < 0.05) with post-exertional malaise and fatigue in UKB and the Netherlands Lifelines biobank, although these associations were not statistically significant after false discovery rate (FDR) corrections for multiple testing. The authors of the DecodeME GWAS posit that this lack of replication may reflect the looser case definition criteria used in those replication data sources relative to DecodeME. The PrecisionLife ® combinatorial analytics platform uses a hypothesis-free approach to identify combinations of features (in this case SNP-genotypes), termed ‘disease signatures’, that are significantly enriched in cases relative to controls or vice versa. Unlike GWAS, which assumes that each genetic variant acts independently of all the others, combinatorial analysis captures linear and non-linear interactions between combinations of features such as SNP-genotypes and potentially other exogenous factors [ 31 ]. As such this approach can identify many more significant genetic associations and related biological mechanisms, even in smaller datasets than those required by GWAS [ 32 – 35 ]. The pipeline employed for combinatorial analysis is described in detail elsewhere [ 32 , 34 ]. Briefly, it employs a deterministic heuristic algorithm and a hypothesis-free data analytics framework encompassing integrated network geometry and statistical genetics approaches. These are employed over multiple rounds of increasing combinatorial complexity. During the signature mining phase, the algorithm iteratively combines increasing numbers of features, resulting in multi-SNP disease signatures where each component SNP-genotype contributes to higher disease association. The platform validates these signatures by comparing their properties (e.g., case prevalence, odds ratios) and the properties of their corresponding SNP networks (i.e., the set of all signatures containing a shared SNP – termed the ‘critical’ SNP in the context of that network) against the properties of signatures and SNP networks identified using the same combinatorial analysis methodology across 1,000 fully random permutations of the case-control sample labels. Statistically validating replication of disease associations for combinatorial disease signatures often requires much larger datasets than replication of GWAS results due to the low frequency of samples having the specific combinations of multiple SNP-genotypes. However, as well as direct biological validation of many targets, the disease signatures, novel candidate genes, and potential drug repurposing candidates identified by the PrecisionLife combinatorial analysis approach in complex, heterogenous diseases such as long COVID and endometriosis have been shown to be significantly more likely to also be linked to those diseases across independent multi-ancestry datasets than corresponding GWAS results [ 33 , 36 ]. A previous combinatorial analysis of ME in a UKB cohort identified 84 disease signatures comprised of 199 SNPs and 14 genes that were significantly enriched in ME patients relative to controls [ 32 ]. That study further used those disease signatures to stratify patients into 15 non-exclusive clusters representing potential mechanistic subtypes of ME. These subtypes centered around genes linked to cellular mechanisms hypothesized to underpin ME, including infection response, autoimmune development, mitochondrial dysfunction, neurotransmitter biology, and circadian rhythm. These 84 signatures were then used in a study of actively protective biology that identified 276 ME protective signatures and 9 disease resilience genes that were significantly enriched in UKB controls who remain healthy despite possessing many ME disease risk signatures [ 37 ]. A combinatorial analysis of long COVID patients from the Sano GOLD cohort identified overlap between the mechanisms associated with Severe and Fatigue Dominant long COVID phenotypes (73 genes found) and the mechanisms linked to ME in UKB, including 9 shared genes [ 32 ], a result consistent with observed symptomatic and biological overlap between the two diseases [ 38 – 40 ]. The results of the long COVID combinatorial analysis in the UK Sano GOLD cohort were recently reproduced with up to 92% overlap of genes by the authors of the original study in a cohort of mixed ancestry long COVID patients from the US All of Us dataset [ 36 ]. An independent study likewise found that 12 out of 80 of the ‘critical’ SNPs and 4 genes identified by the long COVID combinatorial analysis from [ 31 ] were also associated with long COVID in at least one of four patient datasets [ 41 ]. This rate of replication is likely underestimated as that study evaluated each SNP individually rather than the networks of disease signatures containing each SNP, which better reflect complex non-linear disease biology. This study formed part of the Innovate UK funded LOCOME project (#10083274), investigating both long COVID and ME. The intention of the project was to compare the use of combinatorial analysis alongside GWAS to improve understanding of the diseases’ biology, inform new diagnostic tests, and identify potential drug repurposing candidates. We applied combinatorial analytics to populations from the DecodeME dataset to identify disease signatures and genes associated with increased risk of ME. The larger sample numbers and rigorous phenotypic screening employed in the construction of the DecodeME patient cohort allowed us to identify many more genetic associations relative to the previous combinatorial analysis, which relied solely on UKB patient data with its less consistent diagnostic criteria from non-disease specific participant surveys. The combinatorial analytics approach allowed identification of many multiply reproducible genetic associations not found by GWAS. These findings were further refined by assessing the reproducibility of the disease signatures and their component SNPs in a UKB cohort as well as in independent case cohorts of participants from DecodeME. This resulted in a set of multiply reproducible signatures that can potentially be used to inform precision medicine patient stratification tools and test new drug repurposing candidates for ME. We also recapitulated and identified broader overlap between genes implicated in ME and those previously associated with long COVID.

Below is the link to the electronic supplementary material. Supplementary Material 1 Supplementary Material 1 Supplementary Material 2 Supplementary Material 2