Shared Genetic Liability across Systems of Psychiatric and Physical Illness

Domain: chemistry
Published: 2026-02-21
Editorial outcome: minor_revision

Abstract
Epidemiological literature describes widespread psychiatric–physical comorbidity. However, the extent to which observed bivariate risk is transdiagnostic remains uncertain. We apply Genomic SEM and introduce Genomic E-SEM, an exploratory factor-analytic extension of Genomic SEM, here applied to 73 physical outcomes (~1.9 million cases) across eight medical domains. We then modeled latent physical illness factors alongside previously identified psychiatric factors. Thought/psychotic and compulsive factors exhibited few associations with physical illness factors. Conversely, internalizing, neurodevelopmental, and substance-use factors showed substantial genome-wide genetic correlations with all physical illness systems. Multivariate GWAS of a transdiagnostic physical illness factor, defined by 21 physical disorders, identified 27 genomic risk loci and exhibited substantial predictive validity in subsequent PheWAS. This transdiagnostic factor again exhibited substantial genome-wide genetic correlations with internalizing, neurodevelopmental, and substance-use factors. Our results reveal pervasive risk sharing between specific groups of psychiatric and physical conditions and suggest the need for additional nosological frameworks.

Introduction
Individuals with one psychiatric disorder are at such an elevated risk for a second psychiatric diagnosis that comorbidity is the norm rather than the exception. Critically, this problem extends to physical disorders with multimorbidity (≥2 chronic conditions) observed in ~38% of the global population. The separation of psychiatric from physical diseases reflects a traditional dichotomy distinguishing disorders with characterizable physical changes (e.g., organic) from those without (e.g., functional). However, prior research utilizing large-scale national registries indicates that pervasive comorbidity extends beyond this theoretical distinction, with a psychiatric diagnosis also associated with an increased likelihood of developing a physical illness. In isolation, psychiatric conditions already contribute substantially to excess mortality, diminished health-related quality of life, reduced functional status, and increased primary care burdens. When comorbid psychiatric and physical illnesses are present, this further exacerbates negative clinical outcomes. For example, early mortality amongst diabetes cases is 38% higher when comorbid with a diagnosis of major depression. Given these significant clinical implications, understanding the etiological factors underlying these comorbidities is crucial.

Pervasive comorbidity indicates that there are likely transdiagnostic risk pathways shared beyond single psychiatric–physical disorder pairings. For example, literature investigating patterns of physical comorbidity among patients with depression has found an increased risk for coronary artery disease. However, it may be that the risk pathways shared between coronary artery disease and depression are shared broadly across other internalizing disorders (e.g., anxiety) that phenotypically and genetically overlap. However, a hindrance to multivariate applications of modeling broader systems of comorbidity patterns (e.g., structural equation modeling) has historically been the practical limitations of obtaining large samples that encompass a variety of psychiatric and physical outcomes. The current study leverages cutting-edge multivariate genomic methods applied to a wide range of physical and psychiatric datasets to circumvent these pragmatic barriers and evaluate whether genetic overlap across pairs of psychiatric and physical disorders cuts across diagnostic boundaries.

Genomic structural equation modeling (SEM) extends bivariate frameworks for estimating genome-wide genetic correlations, such as LD-score regression (LDSC), into a multivariate approach designed to model genome-wide genetic correlations across multiple traits. As genome-wide association study (GWAS) summary statistics need not come from the same participant sample, Genomic SEM offers the unique opportunity to examine genetic overlap across a broad range of rare and potentially mutually exclusive outcomes. Here, we apply Genomic SEM and introduce and validate a new extension, genomic exploratory SEM (E-SEM), to explore the relationships across psychiatric (compulsive, psychotic/thought, neurodevelopmental, internalizing, and substance use) and newly identified latent genomic factors of physical illness. Genomic E-SEM is developed for genomic exploratory factor analyses and enables examining the multivariate genomic architecture of complex traits for which there are no a priori model specifications or when theoretically derived model specifications do not provide a reasonable fit to the observed data. The combination of these exploratory and confirmatory genomic modeling approaches allows us to include a broad range of psychiatric and physical health phenotypes ascertained from different participant samples in the same statistical model, conduct exploratory and confirmatory factor analyses to identify genomic factors indexing shared risk across groupings of traits, and delineate transdiagnostic risk pathways from more circumscribed shared risk pathways that can reflect overlap unique to a pair of disorders. These results provide a critical context for interpreting the consistent observations of psychiatric–physical comorbidity in the extant literature.

Results
Factor structures across physical illness domains
Through a comprehensive analysis of physical illness phenotypes spanning eight domains (neurological, respiratory, circulatory, digestive, endocrine/metabolic, genitourinary, musculoskeletal, and cancer), we identified 73 physical disorders across 1,889,850 case observations. We first estimated the genome-wide genetic correlations (rg's) using LDSC, which specifically captures the average direction of risk sharing across common variants in the genome. The genetic correlations were used to estimate a single-factor model in each domain in Genomic SEM. When this simple structure did not provide adequate model fit, we utilized Genomic E-SEM to conduct exploratory factor analyses (EFAs). We validate key features of Genomic E-SEM through a series of simulations and demonstrate its flexibility to identify complex factor structures across a wide range of physical outcomes. We further describe the features and validation procedures in the “Method” below. Using a pre-registered modeling pipeline (https://doi.org/10.17605/OSF.IO/YBA58), we identified a common factor for genitourinary disorders and complex, multifactorial structures for the other seven domains (see Supplementary Figs. 1–8). These models all provided a good fit to the data, evidenced by confirmatory fit indices (CFIs)>0.9 and standardized root mean square residuals (SRMR)<0.1. In multifactorial solutions within systems, each latent factor represents subgroups of physical illnesses that exhibit higher levels of genome-wide genetic correlation compared to the other conditions within that domain. For example, in the respiratory system, we identified a factor defined by traits with an atopic component (asthma, nasal polyps) and another factor characterized by high genetic correlation in traits with central pulmonary pathology (pneumonia, chronic airway obstruction). We report model fit indices for the systems-level models in Supplementary Data 1.

Relationships between physical and psychiatric illness latent factors
We took forth the factor structures identified in each physical illness system and modeled them alongside five previously established latent psychiatric factors representing the genetic correlation between compulsive, psychotic/thought, neurodevelopmental, internalizing, and substance use factors defined by subclusters of 13 psychiatric disorders (see Fig. 1). To determine the extent to which factor correlations are operating through common pathways or pathways driven by specific pairs of traits we utilized the Q heterogeneity index. A significant Q statistic indicates that the inter-factor correlation fails to recapture the pairwise genome-wide genetic correlations across the conditions that define the factors. For example, a Q statistic would likely be significant if two physical illnesses on one factor had directionally opposing relationships with the disorders defining a psychiatric factor. Q thereby differentiates between broad risk-sharing patterns reflected in the factor correlation and risk-sharing more appropriately modeled at the more circumscribed disorder level.

Across all eight physical illness domains, we identified 46 significant inter-factor correlations with latent psychiatric factors with a Bonferroni corrected threshold of p<4.5×10−4 (0.05/110 inter-factor correlations) that were not Q significant at the same threshold. Figure 2 visually depicts the inter-factor genetic correlations across the physical and mental illness domains, with the corresponding numeric values presented in Supplementary Data 2. Overall, our results revealed distinct patterns of relationships depending on the psychiatric domain.

For the compulsive factor (defined by anorexia nervosa, obsessive–compulsive disorder, and Tourette’s syndrome), we observed a protective effect for Digestive Factor 1, defined by acute pancreatitis, appendicitis, cholelithiasis/cholecystitis, and peptic ulcers (inter-factor rg: −0.26, SE=0.07). For the psychotic/thought disorders factor (defined by bipolar disorder and schizophrenia), we observed significant genome-wide genetic correlations with physical illness factors in the digestive, musculoskeletal, and respiratory systems (mean rg: 0.15; range: 0.10 to 0.20). We further observed significant Q results with factors in digestive, neurological, endocrine/metabolic, and genitourinary systems, indicating divergent associations across bipolar disorder (BD), schizophrenia (SCZ), and the physical illnesses defining these factors. Bivariate rg's between schizophrenia, bipolar disorder, and physical disorders revealed that BD is likely driving heterogeneous factor relationships indicated by Q, with BD exhibiting greater mean rg's than SCZ.

In contrast, we observed pervasive and substantial genome-wide genetic correlations between neurodevelopmental, internalizing, substance use disorders, and a range of physical illness systems. Within the neurodevelopmental domain (defined by autism spectrum disorder, attention deficit hyperactivity disorder, posttraumatic stress disorder, major depression, and Tourette’s syndrome), we observed moderate to high inter-factor correlations with factors across all physical illness systems (Mean rg: 0.51; range: 0.25–0.76). We also observe heterogeneous associations with factors in the circulatory, digestive, musculoskeletal, respiratory, neurological, and endocrine/metabolic systems. Inspection of the pairwise genome-wide genetic correlations between physical outcomes and the disorders that define the neurodevelopmental factor indicates that these Q results were driven by larger rg's with attention deficit hyperactivity disorder (ADHD) and posttraumatic stress disorder (PTSD).

For internalizing disorders (defined by posttraumatic stress disorder, major depression, and anxiety disorders), we further observed significant inter-factor correlations with specific factors across all physical illness systems (Mean rg: 0.32; range: 0.13–0.55). In addition, we find heterogeneous (Q) relationships between the internalizing factor and factors in the circulatory, digestive, and respiratory systems. Genetic correlations between internalizing and physical illnesses indicate that major depression (MD) and PTSD were driving heterogeneous inter-factor relationships.

In line with neurodevelopmental and internalizing disorders, we observed substantial positive relationships between the substance use disorders (defined by alcohol use disorder, cannabis use disorder, and opioid use disorder) and factors across all physical illness systems (Mean rg: 0.39; range: 0.18–0.91), with one significant Q result in the cancer system. Unlike neurodevelopmental and internalizing disorders, the relative absence of Q findings across most illness systems for the substance use factor indicates that the shared genetic signal across this group of psychiatric disorders is shared broadly with the genetics underlying physical illness.

Providing broad characterization for the specific psychiatric disorders driving heterogeneous factor relationships, we conducted bivariate analyses between psychiatric and physical outcomes (Supplementary Fig. 9). Two findings of note here are that: (i) compared to all other psychiatric disorders in psychiatric domains with multiple Q significant relationships, PTSD, ADHD, and MD had the highest average rg with physical illnesses and (ii) perhaps most surprising was that the average rg between these disorders and physical illnesses was higher than the average pairwise rg between physical disorders and other physical disorders (Fig. 3). This pattern of results collectively indicates that the genetic signal for PTSD, ADHD, and MD indexes an expansive set of risk pathways that transcend psychiatric and physical illness diagnostic boundaries that are likely driving the observed instances of heterogeneous factor relationships (as indexed by Q).

Factor analysis across all physical traits
In the primary analysis, we utilized physical illness domains based on prior epidemiological work, grouping traits according to the body system where the disease manifests. However, this categorization system assumes that physical outcomes manifesting in the same body system share a common etiology. We evaluated this assumption by conducting a follow-up analysis using a data-driven, system-agnostic model across all 73 curated physical traits. We began by fitting a single-factor model across all physical illness measures and then removed any traits with a standardized factor loading below 0.6. This process yielded a common factor across 21 physical traits spanning six of eight domains of illness, which provided an adequate fit to the observed genome-wide genetic correlations (CFI=0.91; SRMR=0.11; see Supplementary Data 3 and Fig. 4A).

Multivariate GWAS and PheWAS of broadly transdiagnostic physical illness factor
We applied Genomic SEM to conduct a multivariate GWAS on the physical illness factor, which estimates the effects of single-nucleotide polymorphisms (SNPs) on the shared liability across the 21 physical disorders defining the factor. We removed genetic loci that were genome-wide significant (p<5×10−8) for Q_SNP, a heterogeneity metric that is conceptually similar to the Q metric described above, that functions to remove SNPs that are unlikely to operate via the factor (e.g., if the SNP is only associated with one of the 21 physical disorders). After removing these Q_SNP loci, we identified 27 genome-wide significant risk loci (Supplementary Data 4 and Fig. 4B).

We then conducted a phenome-wide association study (PheWAS) in the Penn Medicine BioBank (PMBB; N=29,883), testing for associations between multivariate GWAS results for the transdiagnostic physical illness factor (expected sample size N_hat=106,118.7) and a wide range of disease phecodes (see Fig. 4C). Providing external validity for our risk loci, we found substantial effects across multiple physical illness systems (e.g., chronic airway obstruction, type 2 diabetes, obesity, hypertension, and ischemic heart disease). Additional significant associations emerged for sleep-related breathing disorders (sleep apnea, shortness of breath), emphysema, GERD/diseases of the esophagus, chronic pain, and viral hepatitis C. PheWAS results revealed no significant associations with psychiatric outcomes. This result may be due to more limited power to predict outcomes, such as psychiatric disorders, not used to define the 21-disorder common physical factor in Genomic SEM. Further, the ratio of cases to controls was much smaller for many of the psychiatric disorders in the PMBB, resulting in a much smaller effective sample size, relative to many of the physical outcomes. We provide a summary of PheWAS results and the number of cases for all associations with a p<0.05 in Supplementary Data 5.

To investigate this power-based interpretation of null psychiatric associations, we went on to examine better powered 34-indicator common medical factor (N_hat=129,650.2), expanded to include a combination of symptoms and disease proxies (Supplementary Fig. 10a). We find that this updated medical factor significantly predicted psychiatric illnesses (see Supplementary Fig. 10b). However, the lack of associations for the 21-disorder model should still be kept in mind when interpreting genetic correlations reported below.

We then modeled the 21-disorder common physical factor alongside the five psychiatric factors utilized in the primary analysis, using the same Bonferroni corrected threshold (p<4.5×10−4) to identify significant inter-factor correlations. Consistent with the large, pervasive genetic correlations observed in the systems-level analysis (see Fig. 5; Supplementary Data 6), we identified substantial inter-factor correlations with substance use (rg=0.63, SE=0.04), internalizing (rg=0.61, SE=0.03), and neurodevelopmental factors (rg=0.68, SE=0.05). In addition, we identified significant Q results with the psychotic/thought disorders factor.

Methods
Trait selection
All research presented in this study complies with the institutional review board at the University of Colorado Boulder. A prior large-scale epidemiological investigation of psychiatric and physical comorbidity utilized diagnostic bins encompassing various physical conditions available in Danish national registries. Building upon this, we curated an initial list of physical traits by intersecting the diagnostic bins from this prior study with the diagnostic categories described in the UK Biobank (UKB), which are based on the World Health Organization’s International Classification of Diseases and Related Health Problems (ICD) codes (versions 9 and 10). This resulted in nine diagnostic categories (neurological, respiratory, circulatory, digestive, endocrine/metabolic, genitourinary, musculoskeletal, cancer, and hematopoietic) being analyzed in the present study. However, the hematopoietic domain did not contain sufficient traits passing the quality control (QC) filters outlined below. As a result, we excluded hematopoietic disorders from further modeling. To focus on more common physical conditions that frequently coexist in the population, we restricted this list to traits with over 2000 cases in the UKB. The initial list included broad illness definitions (e.g., heart and valve disorders) for each physical trait and their subtypes (e.g., aortic aneurysm). Leveraging Genomic SEM’s capability to integrate traits from multiple samples, we updated this list by incorporating the most recent GWAS summary statistics for each trait from the GWAS Catalog (a comprehensive database of genomic studies-https://www.ebi.ac.uk/gwas/) for either (a) the broad disorder category obtained in the UKB or (b) a specific disorder within the broader category, depending on which had more cases. We utilize the number of cases rather than the total sample size due to the disproportionate number of controls for some traits. Because the UKB may not contain sufficient case numbers for theoretically relevant but statistically rare disorders, we surveyed the extant literature to assess for disorders theoretically relevant in their respective domain or exhibiting previous overlap with psychiatric disorders. We included coronary artery disease, venous thromboembolism, Crohn’s disease, Celiac disease, urinary incontinence, Duparten’s disease, psoriatic arthritis, melanoma, Alzheimer’s disease, Parkinson’s disease, amyotrophic lateral sclerosis, insomnia, and childhood obesity. Next, we pruned symptom measures (e.g., nonspecific chest pain, joint pain) overlapping with other disorders in their system and ambiguously specified disorders (e.g., other urethra/urinary tract disorders, synovium/tendon/bursa disorders-other) that likely contained symptoms represented by more specifically measured phenotypes. This allowed for more precisely defined factors at systems-level analyses and avoided model fit inflation from multicollinearity. Supplementary Data 7 contains a list of all disorders included in our modeling pipeline.

We utilized a previous psychiatric factor structure—containing five correlated psychiatric dimensions-with specific indicators (major depression, bipolar disorder, and posttraumatic stress disorder) updated for newer publicly available GWAS. Encompassing compulsive, psychotic/thought, neurodevelopmental, internalizing, and substance use disorders, the psychiatric factors represent groupings of genome-wide genetic liability across 13 major psychiatric disorders, ascertained with internationally recognized diagnostic criteria and participant self-report of symptom presentation. Supplementary Data 8 list the disorders defining each psychiatric factor, with specific ascertainment methods presented in the corresponding paper of each discovery GWAS.

LD-score regression
We applied a standard set of QC filters to all GWAS summary statistics using the munge function available within Genomic SEM. Munging GWAS summary statistics involves reducing the number of SNPs to HapMap3 SNPs and then pruning by minor allele frequency (MAF) less than 0.01 and an imputation score (INFO) below 0.9 when this information is available in the summary statistics. “Munged” summary statistics were then input for the multivariable version of LDSC implemented in Genomic SEM. LDSC regression models utilized 1000 Genomes Phase 3 European linkage disequilibrium (LD) weights for estimation. These weights exclude the major histocompatibility complex (MHC) as an intricate LD structure in this region that may bias estimates. Following the guidelines from original LDSC developers for producing InterData genetic covariance estimates, we removed physical traits with SNP-based heritability (h2_SNP) Z-statistics < 4.

LDSC produces two meaningful sets of output: the genetic covariance and sampling covariance matrix. The genetic covariance matrix contains SNP-based heritability estimates on the diagonal and genetic covariances on the off-diagonal. These estimates were converted to the liability scale, using the corresponding population prevalence from the original publication and the sum of effective sample sizes across contributing cohorts. The sampling covariance matrix contains squared standard errors on the diagonal (e.g., the sampling variances) and sampling covariances on the off-diagonal, which reflect sampling dependencies that emerge when participant samples overlap across the included traits. The sampling covariance matrix is estimated directly from the GWAS data using a block jackknife resampling procedure. This matrix allows for appropriate parameter estimates and standard errors for traits with varying and unknown degrees of participant sample overlap and varying levels of statistical power across the included GWAS. SNP-based heritability Z-statistics ranged from 4.14 to 26.40 across included physical disorders.

Genomic SEM
For each physical illness domain, we began by supplying the genomic matrices from LDSC to Genomic SEM to fit a common factor in confirmatory factor analysis (CFA) to identify if this simple factor structure adequately recaptures the observed genetic covariance between the traits within each physical illness domain. All models were fit in Genomic SEM using the diagonally weighted least squares estimator, which produces unbiased parameter estimates and standard errors for models with GWAS summary statistics. If the common factor provided an adequate model fit (as determined by CFI>0.9 and SRMR<0.1), then this parsimonious, single-factor model was carried forward for that physical illness domain. If a common factor did not provide adequate fit, we conducted an exploratory factor analysis (EFA) using Genomic E-SEM (described below) to examine multifactorial solutions, a novel extension described and validated via extensive simulations in a separate section below. We determined the number of factors to model in each EFA using the median number of factors selected by the optimal coordinates, acceleration factor, and Kaiser rule when applied to 1000 genetic covariance matrices simulated from the observed genetic covariance and sampling covariance matrices. This novel simulation-based approach to determining the number of genomic factors to model is introduced and validated below. Next, we modeled traits in the EFA with standardized factor loadings exceeding our preregistered cutoff of 0.4 in subsequent confirmatory models.

We assessed the associations between each physical illness and psychiatric disorders by modeling each physical illness domain model alongside five previously established psychiatric factors. The inter-factor correlations within this model represent the degree to which shared signals captured by the physical and psychiatric factors are across factors. However, these factor correlations may not always sufficiently describe patterns of pairwise genome-wide risk sharing of the disorders across factors. We evaluated factor-level heterogeneity by utilizing the Q metric for each physical and psychiatric factor correlation. Q is a heterogeneity metric that becomes more significant when the factor correlation cannot recapture the individual associations across disorders. For example, if two disorders loading on one factor have directionally discordant associations with the disorders loading on another, a factor correlation insufficiently describes this divergent pattern of relationships. Q is a χ2 distributed test statistic with degrees of freedom equal to one less than the cross-factor indicator correlations. We also calculate pairwise associations to follow up on Q results. We first calculated the mean bivariate genetic correlation (rg) between each of the 13 psychiatric disorders and all other psychiatric and physical conditions. Next, we calculated the mean rg amongst all psychiatric and physical disorders.

Genomic E-SEM
Overview. Genomic exploratory structural equation modeling (E-SEM) is the first software developed for genomic exploratory factor analysis, providing a robust framework for identifying the genomic factor structure among a set of variables with limited or no a priori structure. This method enables the identification of latent variables from genomic data without predefining the entire factor structure. Moreover, Genomic E-SEM allows for incorporating subsets of traits with a pre-defined factor structure into exploratory analysis. Genomic E-SEM implements a weighted least squares (WLS) estimator, prioritizing the reduction of misfits in cells of the genetic covariance matrix that are more precisely estimated (e.g., have a smaller SE). Although a WLS estimator will produce reliable model parameter estimates with model identification, the “naïve” standard errors (SEs) and fit statistics generated in stage 2 estimation will be invalid from the full sampling covariance matrix (V_S), not utilized when estimating either estimate. Therefore, robust correction techniques are necessary to obtain valid SEs and model fit statistics. We calculate the correct sampling covariance matrix (V_θ) in stage 2 estimation through a sandwich correction:

V_θ = Δ^0 Γ−1 Δ^−1 Δ^0 Γ−1 V_S Γ−1 Δ^ Δ^0 Γ−1 Δ^−1

where Δ^ is a matrix of model derivatives at each parameter estimate and Γ is the “naïve” weights matrix estimated in stage 2 before correction.

Simulations. We validated key parameters of Genomic E-SEM using simulated LDSC objects for 10% heriData phenotypes with a univariate LDSC intercept of 1.05 to model slight levels of population stratification. These simulations validated the sandwich-corrected SEs and determined the power to discern different population structures in an EFA conducted with Genomic E-SEM. For each model, we simulated participant samples of 50,000 or 500,000. We validated sandwich-corrected SEs on simulated population sizes of 50,000 and the power to discern different population structures on both sets of simulated population sizes. The first population-generating structure was a two-factor model, where each factor had three indicators with a standardized loading of 0.71. We ran three iterations of this model using inter-factor correlations of 0.3, 0.5, and 0.7, with 100 simulations per iteration (e.g., 300 simulations). The second was a single-factor model using standardized factor loadings across all indicators of either 0.32, 0.5, or 0.71 with 100 simulations per iteration (e.g., 300 simulations).

The power to detect the ground-truth structures across the population-generating models is presented in Supplementary Fig. 11, with the validation of sandwich-corrected SEs presented in Supplementary Fig. 12. Across simulations, we find valid sandwich-corrected SEs, with ratios of sandwich-corrected SEs to empirical SEs (calculated as the standard deviation of the parameter estimates across the 100 simulations) close to 1. In addition, we find adequate power to discern population structures, with power scaling in the expected direction. That is, we observed greater power to detect the ground-truth factor structure as factor loadings increased in the case of a single-factor generating model and as inter-factor correlations increased in the case of a two-factor generating model. At a simulated population size of 50,000, power was ≥78% when the standardized inter-factor correlation was 0.5 in a two-factor generating model and when the standardized factor loading was 0.71 (~50% of the variance explained by the factor) in a common-factor generating model. When the simulated population size reached 500,000, power was ≥99% across all generating structures.

Incorporating sampling variability into factor retention heuristics
As the GWAS used as input to LDSC and Genomic SEM can vary with respect to statistical power, it is important to account for these precision differences for all steps in the modeling pipeline. A first step in EFA is determining the number of exploratory factors to model. Here, we introduce an approach to incorporate sampling variability (i.e., power differences) in applied factor retention heuristics using simulated genetic covariance matrices (s*):

s* = s + Lz

where s is the vectorized form of the genetic covariance matrix, and L is the lower triangle Cholesky factor of the sampling covariance matrix (V). We assume:

z ~ N(0, I)

so that, marginally:

s* ~ N(s, V)

Three different factor retention heuristics (optimal coordinates, acceleration factor, Kaiser rule) can then be applied to each simulated matrix s*, with the median of each heuristic’s distribution used as the final retention criterion. By selecting the median factor count from simulated matrices that incorporated differences in power across GWAS of included traits, we capture the central tendency of sampling variability in our retention criterion rather than relying on a single, potentially unrepresentative estimate.

We validated the utility of incorporating sampling variation in determining the number of factors to model by first generating 100 covariance matrices for computational efficiency from an eight-factor generating model containing realistic variability in parameter estimates, variably powered indicators, and a dense implied correlation structure to obscure factor boundaries. For each of the 100 simulated covariance matrices, we then calculated factor retention heuristics using 100 perturbed covariance structures, with the median heuristic per dataset used to capture central tendency in the perturbed covariance matrices. We find that the applied factor retention heuristics on simulated perturbations of the genetic covariance matrix significantly reduce mean absolute error across all three heuristics (p’s<0.05; see Supplementary Fig. 13).

Factor analysis across all physical illness traits
The physical illness domains were chosen based on prior epidemiological work on psychiatric and physical comorbidity. These domains largely reflect physical illness grouped by the body system where the cardinal symptom or disease manifests, with the implicit assumption that these disorders will, in turn, have a shared etiology. Given prior associations with psychiatric illness observed in the phenotypic space using these illness domains, we retained these broad system-based categories for the primary analysis to produce a set of results that are comparable to extant phenotypic literature. In addition, we also conducted a follow-up analysis to explore the factor structure across all curated physical illness traits (e.g., the largest publicly available GWAS; h2_SNP Z-statistic>4) and subsequent relationships with the five psychiatric factors utilized in the primary analysis.

Our data-driven follow-up analysis involved multiple stages. First, we employed the acceleration factor, optimal coordinates, and Kaiser rule to determine the number of factors. Since each index suggested a different number of factors, we chose the most parsimonious factor structure—a common factor solution indicated by the acceleration factor. Next, we conducted a confirmatory factor analysis (CFA) on the 73 curated physical illness traits using Genomic SEM and applied a standardized loading cutoff of 0.6 to decide which traits to include in further modeling. We then fit two subsequent CFAs: (i) to assess the model fit for our data-driven common physical illness factor, and (ii) to evaluate the relationship between the common physical illness factor and the five psychiatric factors used in the primary analysis.

Multivariate GWAS and biological annotation of broadly transdiagnostic physical illness factor
We applied a standard set of QC filters to the individual-level GWAS using the sumstats function in Genomic SEM, which involved filtering SNPs with a minor allele frequency > 1% present in the 1000 Genomes Phase 3 reference panel, aligning to the same reference allele, and ensuring that univariate GWAS results were all partially standardized with respect to the variance in the outcome (the physical illness) but not the predictor (the SNP). The user GWAS function in Genomic SEM was then used to perform a multivariate GWAS on the physical illness factor (expected sample size N_hat=106,118.7) to assess the effects of 6,867,052 SNPs on the shared signal across the 21 physical disorders defining the factor. We utilized FUMA v1.5.2 to identify genomic risk loci with a significance threshold of p<5×10−8. A genetic locus was defined as genome-wide significant SNPs not in linkage disequilibrium (LD) with other significant SNPs at an R²>0.1 based on LD structure from the 1000 Genomes Phase 3 European ancestry reference panel. In addition, genomic loci within 250 kilobases (kb) of one another were merged into a single locus, as these proximal SNPs may not have been identified as in LD due to a mismatch between the physical illness GWAS samples and the reference panel. As a quality control step, we also estimated heterogeneity metrics for each SNP (Q_SNP), which evaluate if SNP effects do not conform with the identified factor structure (e.g., if a SNP is only associated with one of the physical illnesses that defined the common factor). We pruned SNPs with a Q_SNP p-value<5×10−8 and any SNPs within a one-megabase window upstream or downstream of a Q_SNP-significant SNP.

We then applied MAGMA v1.08 to conduct gene property analyses, identifying tissue-specific gene expression, and gene set analyses, identifying previously established sets of genes with shared functional properties or associations with external traits. Gene expression data were obtained from GTEx v8 for 54 body tissues, and gene set data were obtained from MSigDB v7. No tissue-specific gene expression or gene sets passed our multiple-testing corrected significance threshold.

Phewas of physical illness factor in PMBB
First, a polygenic score weight file was generated from the summary statistics by using the “auto” setting of LDpred2 and a linkage equilibrium reference matrix for HapMap3+ variants. Principal component-normalized scores were then calculated for participants of the Penn Medicine BioBank (PMBB) using pgsc_calc. The PMBB is a biobank that links to the electronic health records of consenting participants who receive care within the Penn Medicine health system. Finally, a phenome-wide association study (PheWAS) was performed in individuals genetically similar to a European reference population in the Penn Medicine BioBank (PMBB) using age at enrollment, sex, and the first 10 principal components of genetic ancestry as covariates, testing for associations with disease phecodes.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this study, Lawrence and colleagues aimed to shed light on the extensive comorbidity between psychiatric and physical illnesses through analysis of GWAS summary datasets. For this purpose, they applied their own analytical tool, the multivariate method Genomic SEM, and a novel extension, exploratory Genomic SEM. They first constructed genomic factors that explain shared variance among 76 physical disorders within 8 domains. Next, they estimated genetic correlations between these physical factors and six previously defined psychiatric factors. The authors report genetic correlations between multiple psychical and psychiatric factors, in particular internalizing, neurodevelopmental, and substance use factors. The authors suggest that the findings call into question the bifurcation of psychiatric and physical illness. They also identify a physical factor that explain shared variance across 34 physical illnesses and conduct a GWAS on this factor to identify and characterize the associated loci. Overall, the study is interesting, it represents a considerable body of work and the manuscript is well written. However, I have some comments to the manuscript in its current form.\nThe first comment relates to the selection of the traits, which form the basis of the study. There are several notable absences of major physical disorders with psychiatric comorbidity, such as hypertension, hyperthyroidism, inflammatory bowel diseases and dementia, which should have been included. At the same time, many vague unspecific phenotypes like Respiratory abnormalities, Dysphagia, Chemotherapy, Symptoms involving digestive system, are included and it is unclear why. Preferably, the authors should focus on actual diseases and not symptom groups or other proxies for illness. Moreover, the domains are unequally accounted for by a representative set of disorders, which affect the interpretation of the results. For example, the “neurological factor” explains shared variance in Carpal tunnel syndrome, Epilepsy, Migraine, Headache, Sleep disorders and Nerve root and plexus disorders, but it is unclear how this factor relate to other neurological disorders. Why are not major neurological disorders like Alzheimer’s disease, ALS or multiple sclerosis included, which also have publicly available GWAS and display psychiatric comorbidity? Likewise, for the cancer category, only six phenotypes are included, with four of them relating to the digestive system, as well as one phenotype named “cancer suspected or other” and “chemotherapy” as the final phenotype. Why are not other major cancers like lung cancer included? On the other hand, multiple phenotypes are included within the genitourinary domain. Do these disorders show more pervasive psychiatric comorbidity? I would also recommend splitting diabetes into type 1 and 2, as they have different genetic architectures.\nThe authors write that the GWAS datasets included “represent the largest publicly available GWAS summary statistics for each trait”. However, this is not correct. There are multiple publicly available GWAS for many of the included traits that are larger than those they used. Importantly, this includes several of the psychiatric and substance use disorders, like depression (PMID: 34045744), OUD (PMID: 36171425), PTSD (PMID: 38637617), etc, for some of which the included GWAS are very small. GWAS sample size impacts statistical power and the ability to detect and estimate genetic correlations, with smaller GWAS having lower accuracy.\nBefore constructing physical factors using Genomic SEM, the authors predefined eight physical domains and assigned phenotypes to each domain. However, I think it would be more interesting if they applied their new exploratory multivariate method first across all physical disorders and then defined genomic factors that explained variance across these disorders, irrespective of domain. Their predefined domains prevent the identification of factors that explain risk across organ-specific domains, like auto-immune disorders. Consequently, they cannot investigate correlated risk between auto-immune and psychiatric factors, despite known comorbidity with psychiatric disorders. Furthermore, their current definition of each domain was apparently based on the presumed organ where the illness manifests, similar to the organization of the ICD categories. However, they were not consistent in this. For example, they included stroke in the circulatory domain, contrasting the ICD definition. If they argue that stroke is circulatory in origin that is fine, but then e.g. the auto-immune diseases should perhaps be in a separate category as well reflecting presumed etiology not affected organ?\nGiven all these limitations above, it is unclear to what extent the findings move the field forward and it is difficult to interpret the findings, and compare them across physical domains that are unequally represented by their major disorders and have unequally powered datasets. Furthermore, there is almost a complete lack of references to the existing literature on the shared genetic relationship between psychiatric and physical disorders, which is extensive. The authors should reference the existing literature and their findings should be interpreted in the context of this existing knowledge.\nThe authors repeatedly write that Genomic SEM and LD score regression estimates “genetic overlap”, and use this term throughout the manuscript. However, this is unprecise. Genomic SEM models and estimate correlation between traits and factors, but ignores overlap at the local level or when there is abundant overlap with a mixture of antagonistic and agonistic directional allelic effects. There is a number of other statistical tools that can evaluate shared genetic associations using GWAS datasets, e.g. mendelian randomization, LAVA, PRS, which have revealed shared genetic associations between physical and psychiatric disorders, that the authors do not refer to. Importantly, the terms used also affect the interpretation of the results. For example, the authors write that “psychotic/thought disorders have null genetic relationships with factors across the eight physical illness domains, except for positive genetic correlations with subsets of respiratory and gastrointestinal outcomes.” However, this conclusion is not supported by their findings. Rather, they demonstrate that there are no indications of correlations, but they cannot discard overlap in genes, loci, variants, i.e. a genetic relationship between these factors, which is likely present to some degree.\nThe study introduces a novel extension of Genomic SEM, exploratory Genomic SEM, but it is not discussed at all in the Discussion. It would be good to include a methodological discussion about its usability, limitations, etc.\nThe authors write that their findings call into question the bifurcation of psychiatric and physical conditions, calling for a reevaluation of diagnostic classification system. However, I don’t believe there are sufficient grounds to claim this based on the current analysis. For example, could not the high correlations between ADHD and physical disorders be accounted for by shared associations with low education that also increase risk of these physical disorders? Are the authors suggesting that e.g. ADHD and cardiovascular illness should be classified in the same diagnostic category? With regards to psychotic disorders, for which they do not find much correlations despite extensive somatic comorbidity, are these then different from neurodevelopmental disorders in terms of the bifurcation with physical illness?\nThe motivation for the post-GWAS analysis of the common factor for physical illness is unclear, as this is not connected to the investigation of shared genetic risk with psychiatric disorders. For example, they could look up existing associations of the identified loci with psychiatric disorders using the GWAS catalog.\nSince Genomic SEM builds on the structure of pairwise genetic correlations between traits, it would be natural to display these correlations, for example in a supplementary table or a heat map. Furthermore, the results that account for Fig. 3 are not provided anywhere in the manuscript and are thus difficult to interpret.\nThe authors name their depression phenotype MDD, but the sample is not based on clinically ascertained cases. Instead, this phenotype should be called major depression or depression like the original GWAS (PMID: 30718901). Please also note that there are publicly available and larger MDD GWAS datasets that the authors could use instead, at least if they want to use the largest data available.\nIt would be good if the specific psychiatric factors were specified in the introduction to help readers unfamiliar with previous Genomic SEM work.\nThere are a few instances of double spaces in the text. e.g. page 11, and misspellings, e.g. genomic ssem on page 6\nFig. 6. It would be good to display the categories on y-axis instead of in a box", "reviewer_2": "Dr. Lawrence et al., explored the shared genetic liability across psychiatric and physical illness. They evaluated the genetic overlap between latent genomic factors from 76 physical outcomes, encompassing 1.9 million cases across eight domains. They found that internalizing, neurodevelopmental, and substance use factors share genetic risk with factors across all physical illness domains.\nThis study is of interest but the authors must address the issues outlined below.\n1. Which diagnostic criteria for the psychiatric disorders has been used in the present study? According to the ICD or DSM diagnostic criteria, common psychiatric disorders including anxiety disorder, major depressive disorder, obsessive compulsive disorder, schizophrenia, bipolar disorder, alcohol and drug use disorders, etc.\n2. Why did the authors focus on the compulsive, psychotic/thought, neurodevelopmental, internalizing, and substance-use disorders, which may not be the typical types of the psychiatric disorders? How did the authors assess the above-mentioned psychiatric characteristics?\n3. The authors used the Genomic SEM and Genomic E-SEM software for analyses. Genomic E-SEM is the first software developed for genomic exploratory factor analyses that enables examining the multivariate genomic architecture of complex traits. How about the statistic power for the genomic E-SEM when the compulsive, psychotic/thought, neurodevelopmental, internalizing, and substance-use disorders, are not independent diagnostic units?\n4. The authors conducted a confirmatory factor analysis (CFA) on the 76 curated physical illness traits using Genomic SEM. Have the authors tried to classified 76 curated physical illness traits as several system diseases of the digestive system, cardiovascular system, and respiratory system, etc., before that they explored the shared genetic liability across 13 psychiatric and 76 physical illness?\n5. When the authors explored the shared genetic liability across 13 psychiatric and 76 physical illness, how about the multiple testing corrections?"}, "author_response": "We thank the reviewer for pointing out the absence of key disorders and the presence of symptom overlap in the systems-level models that reduce the interpretability of results and potentially inflate model fit from the presence of symptom overlap amongst physical illness indicators. We now include hypertension, inflammatory bowel diseases, Alzheimer’s disease, and ALS, and split diabetes into type 1 and 2. We further include other key disorders in their respective systems (e.g., coronary artery disease, Parkinson’s disease, melanoma, Addison’s disease). Supplementary Table 6 presents the list of disorders included in our modeling pipeline. We note that excluding hyperthyroidism from modeling is due to underpowered h2 (Z-statistic < 4).\nSNP\nIn addition, we now write on p. 10:\n“Because the UKB may not contain sufficient case numbers for theoretically relevant but statistically rare disorders, we surveyed the extant literature to assess for disorders theoretically relevant in their respective domain or exhibiting previous overlap with psychiatric disorders. We included coronary artery disease, venous thromboembolism, Crohn’s disease, Celiac disease, urinary incontinence, Duparten’s disease, psoriatic arthritis, melanoma, Alzheimer’s disease, Parkinson’s disease, amyotrophic lateral sclerosis, insomnia, and childhood obesity. Next, we pruned symptom measures (e.g., nonspecific chest pain, joint pain) overlapping with other disorders in their system and ambiguously specified disorders (e.g., other urethra/urinary tract disorders, synovium/tendon/bursa disorders-other) that likely contained symptoms represented by more specifically measured phenotypes. This allowed for more precisely defined factors at systems-level analyses and avoided model fit inflation from multicollinearity.”\nRegarding the reviewer’s comment about the absence of key disorders in the systems-level model, we highlight that the inclusion or exclusion of specific disorders in confirmatory models reflects the genome-wide overlap that exists amongst disorders in each system. Since SEM requires a shared component with the other indicators in the model to be present, disorders were not taken forward into the confirmatory modeling stage of our analysis if they exhibited low levels of genome-wide overlap (standardized factor loading < 0.4) within their physical classification system. This may be the case for certain diseases that have a very unique pathology. We also note that the number of disorders in each system before modeling is a function of ascertainment in large-scale biobanks and the power to detect meaningful genetic signals as opposed to their theoretical relevance to psychiatric and physical comorbidity. We believe the value of this inclusive approach to physical trait selection is revealed post-hoc, where we observed high genetic correlations across clusters of psychiatric outcomes and a wide range of physical traits, including those not typically described in the comorbidity literature.\n2. The authors write that the GWAS datasets included “represent the largest publicly available GWAS summary statistics for each trait”. However, this is not correct. There are multiple publicly available GWAS for many of the included traits that are larger than those they used. Importantly, this includes several of the psychiatric and substance use disorders, like depression (PMID: 34045744), OUD (PMID: 36171425), PTSD (PMID: 38637617), etc, for some of which the included GWAS are very small. GWAS sample size impacts statistical power and the ability to detect and estimate genetic correlations, with smaller GWAS having lower accuracy.\nWe thank the reviewer for pointing out newer psychiatric GWAS for the psychiatric model. We do not include the OUD GWAS as it reflects an analysis of the MVP sample, which we understand to only be available through an access request via dbGAP. We exclude restricted access samples such as these with the goal of producing finds that can be more readily tested or brought forward for future applications. We now include updated GWAS for major depression, bipolar disorder, and posttraumatic stress disorder, with the corresponding information for each psychiatric GWAS presented in Supplementary Table 7. We now write on p. 10:\n“We utilized a previous psychiatric factor structure15 - containing five correlated psychiatric dimensions - with specific indicators (major depression, bipolar disorder, and posttraumatic stress disorder) updated for newer publicly available GWAS. Encompassing compulsive, psychotic/thought, neurodevelopmental, internalizing, and substance use disorders, the psychiatric factors represent groupings of genome-wide genetic liability across 13 major psychiatric disorders, ascertained with internationally recognized diagnostic criteria and participant self-report of symptom presentation. Supplementary Table 7 lists the disorders defining each psychiatric factor, with specific ascertainment methods presented in the corresponding paper of each discovery GWAS.”\n3. Before constructing physical factors using Genomic SEM, the authors predefined eight physical domains and assigned phenotypes to each domain. However, I think it would be more interesting if they applied their new exploratory multivariate method first across all physical disorders and then defined genomic factors that explained variance across these disorders, irrespective of domain. Their predefined domains prevent the identification of factors that explain risk across organ-specific domains, like auto-immune disorders. Consequently, they cannot investigate correlated risk between auto-immune and psychiatric factors, despite known comorbidity with psychiatric disorders. Furthermore, their current definition of each domain was apparently based on the presumed organ where the illness manifests, similar to the organization of the ICD categories. However, they were not consistent in this. For example, they included stroke in the circulatory domain, contrasting the ICD definition. If they argue that stroke is circulatory in origin that is fine, but then e.g. the auto-immune diseases should perhaps be in a separate category as well reflecting presumed etiology not affected organ?\nWe agree that the predefined domains represent, at times, arbitrary means of classifying disorders, in internationally utilized diagnostic systems and large-scale biobanks. However, the initial domain-level grouping of traits was non-arbitrary, reflecting a consistent application of the groupings appearing in the UK Biobank (UKB). These can differ from ICD 11 categorizations as the UKB has large amounts of data coded with the ICD 9-10. We retain their use to make the presented findings comparable to existing literature and to evaluate the extent to which disorders classified with this approach do indeed exhibit genome-wide levels of risk-sharing.\nRegarding modeling physical illnesses irrespective of their domain, we approach this concern by modeling all physical disorders at a systems-agnostic level, identifying a common factor representing genome-wide liability from disorders across all eight domains. A newly included PheWAS provides external validity for the cross-disorder risk captured by this factor, predicting disease phecodes across multiple illness systems. Further discussion of the PheWAS is in the response to point 8 below.\n4a. Given all these limitations above, it is unclear to what extent the findings move the field forward and it is difficult to interpret the findings and compare them across physical domains that are unequally represented by their major disorders and have unequally powered datasets.\nWe thank the reviewer for pointing out limitations in our previous approach and believe that the included changes have significantly increased the interpretability of our results. Among these changes, we note the inclusion of key disorders within organ-based domains, the removal of disease proxies and ambiguously defined phenotypes, and the out-of-sample predictive capabilities of the common-factor physical illness GWAS across multiple illness domains. In addition, we highlight the utility of Q in adding much-needed context for the specificity of the Factor pervasive risk-sharing observed between psychiatric and physical illnesses across both genetic and phenotypic investigations. These results, alongside methods introduced and validated in this same manuscript, collectively make several novel contributions to the field.\nRegarding unequally powered GWAS data, Genomic SEM and E-SEM take these differences in precision into account (via the sampling covariance matrix [V]) when producing parameter estimates and their corresponding standard errors. The resubmission has also been updated to include validation and introduction of methods for determining the number of factors to model in EFA that similarly incorporates differences in power across the GWAS using the V matrix. When determining the number of factors to model in a covariance matrix, traditional factor retention heuristics (e.g., Kaiser rule, optimal coordinates, and acceleration factor) treat the covariance matrix (S) as fixed and ignore sampling variability across included GWAS. Here, we propose and validate a method of propagating uncertainty by generating perturbed replicates of S and applying eigen-based factor retention heuristics. We now write on ps. 13-14:\n“Incorporating Sampling Variability into Factor Retention Heuristics\nAs the GWAS used as input to LDSC and Genomic SEM can vary with respect to statistical power, it is important to account for these precision differences for all steps in the modeling pipeline. A first step in EFA is determining the number of exploratory factors to model. Here, we introduce an approach to incorporate sampling variability (i.e., power differences) in applied factor retention heuristics using simulated genetic covariance matrices (s*):\n𝑠∗ = 𝑠+𝐿𝑧\nwhere 𝑠 is the vectorized form of the genetic covariance matrix and 𝐿 is the lower triangle Cholesky factor of the sampling covariance matrix (𝑉). We assume:\n𝑧 ∼ 𝒩(0,𝐼)\nso that, marginally:\n𝑠∗ ∼ 𝒩(𝑠,𝑉)\nThree different factor retention heuristics (optimal coordinates, acceleration factor, Kaiser rule) can then be applied to each simulated matrix s*, with the median of each heuristic’s distribution used as the final retention criterion. By selecting the median factor count from simulated matrices that incorporate differences in power across GWAS of included traits, we capture the central tendency of sampling variability in our retention criterion rather than relying on a single, potentially unrepresentative estimate.\nWe validated the utility of incorporating sampling variation in determining the number of factors to model by first generating 100 covariance matrices for computational efficiency from an eight-factor generating model containing realistic variability in parameter estimates, variably powered indicators, and a dense implied correlation structure to obscure factor boundaries. For each of the 100 simulated covariance matrices, we then calculated factor retention heuristics using 100 perturbed covariance structures, with the median heuristic per dataset used to capture central tendency in the perturbed covariance matrices. We find that the applied factor retention heuristics on simulated perturbations of the genetic covariance matrix significantly reduce mean absolute error across all three heuristics (p’s < .05; see Supplementary Fig. 12).”\nWe believe that incorporating sampling variability across exploratory and confirmatory analyses allows for valid conclusions when modeling traits with datasets that may be unequally powered but are pre-selected to have sufficient power for producing interpretable estimates. We agree that these differences in power are an important interpretive consideration and have updated the Method section to discuss how these power differences affect model estimates (p. 11):\n“The sampling covariance matrix is estimated directly from the GWAS data using a block jackknife resampling procedure. This matrix allows for appropriate parameter estimates and standard errors for traits with varying and unknown degrees of participant sample overlap and varying levels of statistical power across the included GWAS. SNP-based heritability Z-statistics ranged from 4.14 to 26.40, across included physical disorders.\nFor each physical illness domain, we began by supplying the genomic matrices from LDSC to Genomic SEM to fit a common factor in confirmatory factor analysis (CFA) to identify if this simple factor structure adequately recaptures the observed genetic covariance between the traits within each physical illness domain. All models were fit in Genomic SEM using the diagonally weighted least squares estimator (DWLS), which produces unbiased standard errors when GWAS vary in sample size and prioritizes reducing model misfit for genetic covariance estimates with greater precision. This could include, for example, prioritizing producing a lower factor loading for the best-powered GWAS if it has low genetic correlations with the other traits.”\nIn addition, we now write in the Discussion (p. 8):\n“The absence of significant results should also be interpreted in the context of power differences for the GWAS that defined the psychiatric or physical factors.”\n4b. Furthermore, there is almost a complete lack of references to the existing literature on the shared genetic relationship between psychiatric and physical disorders, which is extensive. The authors should reference the existing literature, and their findings should be interpreted in the context of this existing knowledge.\nWe appreciate the reviewer pointing out the lack of reference to previous molecular genetic investigations of psychiatric and physical comorbidity and have included several additional citations. We now include new citations and write in five instances on ps. 6-7:\n“Genome-wide investigations have corroborated these findings. The extant literature has also linked common-variant liability between psychotic/thought disorders to metabolic and gastrointestinal liability30,31.”\n“Conversely, large-scale genomic analyses have found modest negative genome-wide relationships between compulsive disorders and body-mass index, with null relationships for cardiometabolic traits34.”\n“Our results support prior literature that finds extensive pleiotropy between substance use GWAS and metabolic, inflammatory, and pulmonary traits30.”\n“The internalizing factor also exhibited genome-wide overlap with at least one factor in all eight physical illness domains, extending extant literature that finds pervasive genome-wide risk between internalizing disorders and cardio-metabolic37 and respiratory traits38.”\n“We observed a similar broad and disorder-specific risk-sharing pattern for the neurodevelopmental factor, with disorders defining this factor previously found to exhibit substantial genome-wide risk-sharing with physical illnesses30.”\n5. The authors repeatedly write that Genomic SEM and LD score regression estimates “genetic overlap” and use this term throughout the manuscript. However, this is unprecise. Genomic SEM models and estimate correlation between traits and factors, but ignores overlap at the local level or when there is abundant overlap with a mixture of antagonistic and agonistic directional allelic effects. There is a number of other statistical tools that can evaluate shared genetic associations using GWAS datasets, e.g. mendelian randomization, LAVA, PRS, which have revealed shared genetic associations between physical and psychiatric disorders, that the authors do not refer to. Importantly, the terms used also affect the interpretation of the results. For example, the authors write that “psychotic/thought disorders have null genetic relationships with factors across the eight physical illness domains, except for positive genetic correlations with subsets of respiratory and gastrointestinal outcomes.” However, this conclusion is not supported by their findings. Rather, they demonstrate that there are no indications of correlations, but they cannot discard overlap in genes, loci, variants, i.e. a genetic relationship between these factors, which is likely present to some degree.\nWe agree that the term genetic overlap is overly vague with respect to level of analysis (genome-wide) and the fact that genetic correlation captures average risk sharing that may sum to 0 when there are mixtures of loci with convergent and divergent effects across the traits. We have updated the language across the paper and now use the term genome-wide genetic correlation to more accurately describe what we are modeling. We also have removed any claims that genetic relationships are entirely null based on our findings; as the reviewer notes, we cannot rule out shared signal at more granular levels of biological analysis. We now write in the Introduction (p. 2):\n“Genomic Structural Equation Modeling (SEM16) extends bivariate frameworks for estimating genome-wide genetic correlations, such as LD-score regression (LDSC17), into a multivariate approach designed to model genome-wide genetic correlations across multiple traits.”\nWe have also updated the Results section to explicitly state (p. 3):\n“We first estimated the genome-wide genetic correlations from LDSC, which specifically captures the average direction of risk sharing across common variants the genome.”\nFinally, we have updated the text to note that value of applying local methods like LAVA in future work (p. 8):\n“The current findings specifically index genetic risk sharing using genome-wide genetic correlations estimated from common variants (e.g., SNPs with a minor allele frequency > 1%). Future work should reevaluate our observed patterns of risk-sharing in rare genetic variation and at more localized levels of analysis. The importance of analysis at more granular levels is highlighted by recent applications of local analysis of [co]variant association (LAVA42) that find significant local, but not genome-wide correlations, for psychiatric disorders and insulin-resistant conditions (e.g., type 2 diabetes43) and neurodegenerative diseases (e.g., Alzheimer’s44).”\n6. The study introduces a novel extension of Genomic SEM, exploratory Genomic SEM, but it is not discussed at all in the Discussion. It would be good to include a methodological discussion about its usability, limitations, etc.\nWe agree that a methodological discussion of Genomic E-SEM will aid its future applications. We now write on p. 6:\n“Across eight systems of medical illness, Genomic E-SEM characterized the latent structure across high-dimensional genetic landscapes. Combined with robust applications of factor retention heuristics, Genomic E-SEM uncovered patterns of global genetic signal shared between an array of physical illnesses, with subsequent confirmatory models exceeding conventional fit criteria. Genomic E-SEM offers a user-friendly, scalable pipeline with valid standard errors, allowing for more reliable estimates than naïve exploratory approaches. Nevertheless, it shares the typical exploratory caveats – specification of the number of latent factors and dependence on investigator judgment for factor loading cut-offs, underscoring the necessity for future replication. Our identification of well-fitting latent structures across all eight systems in a single analytic framework reflects Genomic E-SEM’s practical utility for mapping the genomic architecture of complex multimorbidity.”\n7. The authors write that their findings call into question the bifurcation of psychiatric and physical conditions, calling for a reevaluation of diagnostic classification system. However, I don’t believe there are sufficient grounds to claim this based on the current analysis. For example, could not the high correlations between ADHD and physical disorders be accounted for by shared associations with low education that also increase risk of these physical disorders? Are the authors suggesting that e.g. ADHD and cardiovascular illness should be classified in the same diagnostic category? With regards to psychotic disorders, for which they do not find much correlations despite extensive somatic comorbidity, are these then different from neurodevelopmental disorders in terms of the bifurcation with physical illness?\nWe agree that the previous language made particular claims without external validation and did not accurately express our argument for reevaluating diagnostic systems.\nWith the inclusion of PheWAS results that uncover substantial cross-disorder predictive validity in an external sample, we believe that the updated language is appropriate for the observed convergence across results. We now write on ps. 7-8:\n“Our findings challenge this dichotomy. Within physical illness domains, we identify shared genome-wide risk through a common factor that cuts across six illness systems and exhibits substantial cross-disorder predictive validity in external samples. We observed pervasive genome-wide genetic correlation between neurodevelopmental, internalizing, and substance use factors and physical illness factors in systems and systems-agnostic levels of analysis - with the average genetic correlations between ADHD, MD, PTSD, and physical disorders higher than the mean pairwise relationships within the physical domain.”\nRegarding our arguments for updated classification systems, we now write on p. 8:\n“These findings suggest a need for additional classification systems conceptualizing illnesses as outcomes from systems of interconnected risk, likely unrelated to particular phenotypic presentations.”\nWith the inclusion of key disorders within each physical domain and the exclusion of ambiguously defined phenotypes and disease proxies, we now find increased associations with the psychotic/thought disorders factor, providing support and refining our understanding of its risk patterns. We now write on ps. 6-7:\n“Our results support prior literature and extend our knowledge of relationship specificity across these disease systems. At the factor level, we find that psychotic/thought disorders (BD, SCZ) exhibit modest genome-wide genetic correlations with factors in digestive, musculoskeletal, and respiratory domains, indicating that liability shared between these disorders further cuts across diagnostic boundaries in specific domains of physical illness. Certain risks, however, appear to be disproportionally driven by single psychotic/thought disorders, suggesting that signal from BD is driving risk-sharing with genitourinary disorders and specific factors in neurological (sleep disorders), endocrine/metabolic (e.g., hypothyroidism, diabetes type 1), and digestive (e.g., irritable bowel syndrome) systems.”\nAnd thank the reviewer for the suggested inclusion and exclusion of traits.\n8. The motivation for the post-GWAS analysis of the common factor for physical illness is unclear, as this is not connected to the investigation of shared genetic risk with psychiatric disorders. For example, they could look up existing associations of the identified loci with psychiatric disorders using the GWAS catalog.\nWe appreciate the reviewer’s comments and now extend our results by performing a PheWAS with the common-factor physical illness GWAS in the Penn Medicine Biobank (PMBB). Providing external validity that the observed risk loci exert pleiotropic effects across diagnoses, we find substantial predictive validity across multiple physical illness domains. We now write on p. 5:\n“We then conducted a phenome-wide association study (PheWAS) in the Penn Medicine BioBank (PMBB), testing for associations between multivariate GWAS results for the transdiagnostic physical illness factor and a wide range of disease phecodes (see Fig. 6). Providing external validity for our risk loci, we found substantial effects across multiple physical illness systems (e.g., chronic airway obstruction, type 2 diabetes, obesity, hypertension, and ischemic heart disease). Additional significant associations emerged for sleep-related breathing disorders (sleep apnea, shortness of breath), emphysema, GERD/diseases of the esophagus, chronic pain, and viral hepatitis C. These results highlight the cross-disorder predictive capability of the transdiagnostic physical illness factor in the PMBB.”\nAnd p. 15:\n“Phewas of Physical Illness Factor in PMBB\nFirst, a polygenic score weight file was generated from the summary statistics by using the “auto” setting of LDpred2 and a linkage disequilibrium reference matrix for HapMap3+ variants56. Principal component-normalized scores were then calculated for participants of the Penn Medicine BioBank (PMBB) using pgsc_calc57. The PMBB is a biobank that links to the electronic health records of consenting participants who receive care within the Penn Medicine health system58. Finally, a phenome-wide association study (PheWAS) was performed in individuals genetically similar to a European reference population in the Penn Medicine BioBank (PMBB) using age at enrollment, sex, and the first 10 principal components of genetic ancestry as covariates, testing for associations with disease phecodes.”\nGiven its predictive validity, we believe that the public availability of these summary statistics will be a valuable addition to the field and aid in future cross-disorder work.\n9. Since Genomic SEM builds on the structure of pairwise genetic correlations between traits, it would be natural to display these correlations, for example in a supplementary table or a heat map. Furthermore, the results that account for Fig. 3 are not provided anywhere in the manuscript and are thus difficult to interpret.\nWe thank the reviewer for highlighting how including this information can aid in the interpretability and downstream applications of the present findings. We now include a table of the results and a figure for pairwise genetic correlations in Supplementary Fig. 9.\n10. The authors name their depression phenotype MDD, but the sample is not based on clinically ascertained cases. Instead, this phenotype should be called major depression or depression like the original GWAS (PMID: 30718901). Please also note that there are publicly available and larger MDD GWAS datasets that the authors could use instead, at least if they want to use the largest data available.\nWe have updated our language to reflect phenotypic ascertainment in MD GWAS on ps. 4, 7, 8, & 9.\n11. It would be good if the specific psychiatric factors were specified in the introduction to help readers unfamiliar with previous Genomic SEM work.\nWe agree that this will add helpful context to readers unfamiliar with previous Genomic SEM work. We now write on p. 1:\n“Here, we apply Genomic SEM and introduce and validate a new extension, Genomic Exploratory SEM (E-SEM), to explore the relationships across psychiatric (compulsive, psychotic/thought, neurodevelopmental, internalizing, and substance use) and newly identified latent genomic factors of physical illness.”\n12. There are a few instances of double spaces in the text. e.g. page 11, and misspellings, e.g. genomic ssem on page 6\nWe appreciate the reviewer pointing out these errors and have now extensively copy-edited the manuscript.\n13. Fig. 6. It would be good to display the categories on y-axis instead of in a box\nWe agree this would increase the interpretability of Fig. 7 (previously Fig. 6) and now include the psychiatric systems on the y-axis.\nReviewer #2 (Remarks to the Author):\nDr. Lawrence et al., explored the shared genetic liability across psychiatric and physical illness. They evaluated the genetic overlap between latent genomic factors from 76 physical outcomes, encompassing 1.9 million cases across eight domains. They found that internalizing, neurodevelopmental, and substance use factors share genetic risk with factors across all physical illness domains.\nThis study is of interest but the authors must address the issues outlined below.\n1. Which diagnostic criteria for the psychiatric disorders has been used in the present study? According to the ICD or DSM diagnostic criteria, common psychiatric disorders including anxiety disorder, major depressive disorder, obsessive compulsive disorder, schizophrenia, bipolar disorder, alcohol and drug use disorders, etc.\nWe thank the reviewer for pointing out the lack of clarity in the ascertainment protocol across the utilized psychiatric GWAS data. We now write on p. 10:\n“Encompassing compulsive, psychotic/thought, neurodevelopmental, internalizing, and substance use disorders, the psychiatric factors represent groupings of genome-wide genetic liability across 13 major psychiatric disorders, ascertained with internationally recognized diagnostic criteria and participant self-report of symptom presentation. Supplementary Table 7 lists the disorders defining each psychiatric factor, with specific ascertainment methods presented in the corresponding paper of each discovery GWAS.”\n2. Why did the authors focus on the compulsive, psychotic/thought, neurodevelopmental, internalizing, and substance-use disorders, which may not be the typical types of the psychiatric disorders? How did the authors assess the above-mentioned psychiatric characteristics?\nWe agree that modeling the shared signal across psychiatric disorders may miss disease-specific risk pathways. We now include bivariate analyses that display the pairwise genetic correlations between all psychiatric and physical traits. In addition, we now write on ps. 4-5:\n“Providing broad characterization for the specific psychiatric disorders driving heterogeneous factor relationships, we conducted bivariate analyses between psychiatric and physical outcomes (Supplementary Fig. 9). Two findings of note here are that: (i) compared to all other psychiatric disorders in psychiatric domains with multiple Q significant relationships, PTSD, ADHD, and Factor MD had the highest average r with physical illnesses and (ii) perhaps most surprising was that g the average r between these disorders and physical illnesses was higher than the average pairwise g r between physical disorders and other physical disorders (Fig. 3). This pattern of results collectively indicates that the genetic signal for PTSD, ADHD, and MD indexes an expansive set of risk pathways that transcend psychiatric and physical illness diagnostic boundaries that are likely driving the observed instances of heterogeneous factor relationships (as indexed by Q ).”\nFactor\nWe present all bivariate genetic correlations in Supplementary Fig. 9. We focus on the shared signal captured by the inter-factor correlations across medical or psychiatric factors (e.g., compulsive) as the bivariate associations indicates widespread risk sharing across clusters of disorders that are more parsimoniously understood and modeled in an SEM framework. As this is one of the motivating frameworks for the paper, we sought to significantly revise the Introduction to make this explicit early on. We now write on p. 2:\n“Pervasive comorbidity indicates that there are likely transdiagnostic risk pathways shared beyond single psychiatric-physical disorder pairings. For example, literature investigating patterns of physical comorbidity among patients with depression has found an increased risk for coronary artery disease14. However, it may be that the risk pathways shared between coronary artery disease and depression are shared broadly across other internalizing disorders (e.g., anxiety) that phenotypically1 and genetically15 overlap. However, a hindrance to multivariate applications of modeling broader systems of comorbidity patterns (e.g., structural equation modeling) has historically been the practical limitations of obtaining large samples that encompass a variety of psychiatric and physical outcomes. The current study leverages cutting-edge multivariate genomic methods applied to a wide range of physical and psychiatric datasets to circumvent these pragmatic barriers and evaluate whether genetic overlap across pairs of psychiatric and physical disorders cuts across diagnostic boundaries.\nGenomic Structural Equation Modeling (SEM16) extends bivariate frameworks for estimating genome-wide genetic correlations, such as LD-score regression (LDSC17), into a multivariate approach designed to model genome-wide genetic correlations across multiple traits.\nAs genome-wide association study (GWAS) summary statistics need not come from the same participant sample, Genomic SEM offers the unique opportunity to examine genetic overlap across a broad range of rare and potentially mutually exclusive outcomes. Here, we apply Genomic SEM and introduce and validate a new extension, Genomic Exploratory SEM (E-SEM), to explore the relationships across psychiatric (compulsive, psychotic/thought, neurodevelopmental, internalizing, and substance use) and newly identified latent genomic factors of physical illness.\nGenomic E-SEM is the first software developed for genomic exploratory factor analyses that enables examining the multivariate genomic architecture of complex traits for which there are no a priori model specifications or when theoretically derived model specifications do not provide a reasonable fit to the observed data. The combination of these exploratory and confirmatory genomic modeling approaches allows us to include a broad range of psychiatric and physical health phenotypes ascertained from different participant samples in the same statistical model, conduct exploratory and confirmatory factor analyses to identify genomic factors indexing shared risk across groupings of traits, and delineate transdiagnostic risk pathways from more circumscribed shared risk pathways that can reflect overlap unique to a pair of disorders. These results provide a critical context for interpreting the consistent observations of psychiatric-physical comorbidity in the extant literature.”\n3. The authors used the Genomic SEM and Genomic E-SEM software for analyses. Genomic E-SEM is the first software developed for genomic exploratory factor analyses that enables examining the multivariate genomic architecture of complex traits. How about the statistic power for the genomic E-SEM when the compulsive, psychotic/thought, neurodevelopmental, internalizing, and substance-use disorders, are not independent diagnostic units?\nWe agree that the field standard of treating disorders defining psychiatric factors as independent diagnostic units may not reflect the ground-truth organization of symptoms, with theoretical diagnostic models potentially reflecting historical approaches to disease conceptualization. We hope that GWAS of disease symptoms becomes publicly available for assumptions embedded in utilized diagnostic systems to be formally tested, allowing for more parsimonious models of symptom clustering. We now write on p. 8:\n“Further, as GWAS data for psychiatric and physical symptoms become publicly available, future work should evaluate the validity of conceptualizing disorders as independent diagnostic units.”\nRegarding the statistical power of E-SEM, we validate this approach in a series of simulations, demonstrating its capability to draw valid inferences from differentially powered GWAS data. We write on p. 12:\n“Simulations. We validated key parameters of Genomic E-SEM using simulated LDSC objects for 10% heritable phenotypes with a univariate LDSC intercept of 1.05 to model slight levels of population stratification. These simulations validated the sandwich-corrected SEs and determined the power to discern different population structures in an EFA conducted with Genomic E-SEM. For each model, we simulated participant samples of 50,000 or 500,000. We validated sandwich-corrected SEs on simulated population sizes of 50,000 and the power to discern different population structures on both sets of simulated population sizes. The first population-generating structure was a two-factor model, where each factor had three indicators with a standardized loading of 0.71. We ran three iterations of this model using inter-factor correlations of 0.3, 0.5, and 0.7, with 100 simulations per iteration (e.g., 300 simulations). The second was a single-factor model using standardized factor loadings across all indicators of either 0.32, 0.5, or 0.71 with 100 simulations per iteration (e.g., 300 simulations). The power to detect the ground-truth structures across the population-generating models is presented in Supplementary Fig. 9, with the validation of sandwich-corrected SEs presented in Supplementary Fig. 10. Across simulations, we find valid sandwich-corrected SEs, with ratios of sandwich-corrected SEs to empirical SEs (calculated as the standard deviation of the parameter estimates across the 100 simulations) close to 1. In addition, we find adequate power to discern population structures, with power scaling in the expected direction. That is, we observed greater power to detect the ground-truth factor structure as factor loadings increased in the case of a single-factor generating model and as inter-factor correlations increased in the case of a two-factor generating model. At a simulated population size of 50,000, power was ≥ 78% when the standardized inter-factor correlation was 0.5 in a two-factor generating model and when the standardized factor loading was 0.71 (~50% of the variance explained by the factor) in a common-factor generating model. When the simulated population size reached 500,000, power was ≥ 99% across all generating structures.”\n4. The authors conducted a confirmatory factor analysis (CFA) on the 76 curated physical illness traits using Genomic SEM. Have the authors tried to classified 76 curated physical illness traits as several system diseases of the digestive system, cardiovascular system, and respiratory system, etc., before that they explored the shared genetic liability across 13 psychiatric and 76 physical illness?\nWe thank the reviewer for pointing out that organ-based systems potentially create artificial diagnostic categories that don’t reflect genome-wide levels of risk-sharing that can exist in disorders across organ systems. Regarding this concern, we conduct a systems-agnostic analysis, modeling all 73 physical illnesses in a common factor that captures disease risk across all eight physical domains. In addition, we now include a PheWAS of the common-factor physical-illness GWAS, which provides substantial predictive validity across several pertinent disease domains. Further, the common factor exhibited substantial genome-wide relationships with psychiatric factors that exhibit similar broad risk-sharing patterns observed in systems-level analyses.\nWe now write on p. 5:\n“We applied Genomic SEM to conduct a multivariate GWAS on the physical illness factor, which estimates the effects of single-nucleotide polymorphisms (SNPs) on the shared liability across the 21 physical disorders defining the factor. We removed genetic loci that were genome-wide significant (p < 5 ×10−8) for Q a heterogeneity metric that is conceptually similar to the SNP, Q metric described above, that functions to remove SNPs that are unlikely to operate via the Factor factor (e.g., if the SNP is only associated with one of the 21 physical disorders). After removing these Q loci, we identified 27 genome-wide significant risk loci (Supplementary Table 4 and SNP Fig. 5).\nWe then conducted a phenome-wide association study (PheWAS) in the Penn Medicine BioBank (PMBB), testing for associations between multivariate GWAS results for the transdiagnostic physical illness factor and a wide range of disease phecodes (see Fig. 6). Providing external validity for our risk loci, we found substantial effects across multiple physical illness systems (e.g., chronic airway obstruction, type 2 diabetes, obesity, hypertension, and ischemic heart disease). Additional significant associations emerged for sleep-related breathing disorders (sleep apnea, shortness of breath), emphysema, GERD/diseases of the esophagus, chronic pain, and viral hepatitis C. These results highlight the cross-disorder predictive capability of the transdiagnostic physical illness factor in the PMBB.\nWe then modeled the common physical factor alongside the five psychiatric factors utilized in the primary analysis, using the same Bonferroni corrected threshold (p < 4.5 ×10−4) to identify significant inter-factor correlations. Consistent with the large, pervasive genetic correlations observed in the systems-level analysis (see Fig. 7; Supplementary Table 5). We identified substantial inter-factor correlations with substance use (r = 0.63, SE = 0.04), internalizing (r = g g 0.61, SE = 0.03), and neurodevelopmental factors (r = 0.68, SE = 0.05). In addition, we identified g significant Q results with the psychotic/thought disorders factor.”\nFactor\nAnd on ps. 14-15:\n“The physical illness domains were chosen based on prior epidemiological work on psychiatric and physical comorbidity6. These domains largely reflect physical illness grouped by the body system where the cardinal symptom or disease manifests, with the implicit assumption that these disorders will, in turn, have a shared etiology. Given prior associations with psychiatric illness observed in the phenotypic space using these illness domains, we retained these broad system-based categories for the primary analysis to produce a set of results that are comparable to extant phenotypic literature. In addition, we also conducted a follow-up analysis to explore the factor structure across all curated physical illness traits (e.g., largest publicly available GWAS; h2 Z- SNP statistic > 4) and subsequent relationships with the five psychiatric factors utilized in the primary analysis.\nOur data-driven follow-up analysis involved multiple stages. First, we employed the acceleration factor, optimal coordinates, and Kaiser rule to determine the number of factors. Since each index suggested a different number of factors, we chose the most parsimonious factor structure—a common factor solution indicated by the acceleration factor. Next, we conducted a confirmatory factor analysis (CFA) on the 73 curated physical illness traits using Genomic SEM and applied a standardized loading cutoff of 0.6 to decide which traits to include in further modeling. We then fit two subsequent CFAs: (i) to assess the model fit for our data-driven common physical illness factor, and (ii) to evaluate the relationship between the common physical illness factor and the five psychiatric factors used in the primary analysis.\nWe applied a standard set of QC filters to the individual-level GWAS using the sumstats function in GenomicSEM, which involved filtering SNPs with a minor allele frequency > 1% present in the 1000 Genomes Phase 3 reference panel, aligning to the same reference allele, and ensuring that univariate GWAS results were all partially standardized with respect to the variance in the outcome (the physical illness) but not the predictor (the SNP). The userGWAS function in Genomic SEM16 was then used to perform a multivariate GWAS on the physical illness factor (expected sample size50 (𝑁̂) = 106,118.7) to assess the effects of 6,867,052 SNPs on the shared signal across the 21 physical disorders defining the factor. We utilized FUMA v1.5.251 to identify genomic risk loci with a significance threshold of p < 5 ×10−8. A genetic locus was defined as genome-wide significant SNPs not in linkage disequilibrium (LD) with other significant SNPs at an R² > 0.1 based on LD structure from the 1000 Genomes Phase 3 European ancestry reference panel. In addition, genomic loci within 250 kilobases (kb) of one another were merged into a single locus as these proximal SNPs may not have been identified as in LD due to a mismatch between the physical illness GWAS samples and reference panel. As a quality control step, we also estimated heterogeneity metrics for each SNP (Q 52), which evaluate if SNP effects do not SNP conform with the identified factor structure (e.g., if a SNP is only associated with one of the physical illness that defined the common factor). We pruned SNPs with a Q p-value < 5 × SNP 10−8 and any SNPs within a one-megabase window upstream or downstream of a Q -significant SNP SNP.\nWe then applied MAGMA v1.0853 to conduct gene property analyses, identifying tissue-specific gene expression, and gene set analyses, identifying previously established sets of genes with shared functional properties or associations with external traits. Gene expression data were obtained from GTEx v854 for 54 body tissues, and gene set data were obtained from MSigDB v755. No tissue-specific gene expression or gene sets passed our multiple-testing corrected significance threshold.\nFirst, a polygenic score weight file was generated from the summary statistics by using the “auto” setting of LDpred2 and a linkage disequilibrium reference matrix for HapMap3+ variants56. Principal component-normalized scores were then calculated for participants of the Penn Medicine BioBank (PMBB) using pgsc_calc57. The PMBB is a biobank that links to the electronic health records of consenting participants who receive care within the Penn Medicine health system58. Finally, a phenome-wide association study (PheWAS) was performed in individuals genetically similar to a European reference population in the Penn Medicine BioBank (PMBB) using age at enrollment, sex, and the first 10 principal components of genetic ancestry as covariates, testing for associations with disease phecodes.”\n5. When the authors explored the shared genetic liability across 13 psychiatric and 76 physical illness, how about the multiple testing corrections?\nWe thank the reviewer for pointing out the importance of multiple testing corrections. We apply a Bonferroni correction that reflects the number of inter-factor relationships at systems-level analyses (.05/110). To make this more clear to readers, we write on p. 4:\n“Across all eight physical illness domains, we identified 46 significant inter-factor correlations with latent psychiatric factors with a Bonferroni corrected threshold of p < 4.5 ×10−4 (0.05/110 inter-factor correlations) and not Q significant at the same threshold.”\nFactor\nAnd on p. 5:\n“We then modeled the common physical factor alongside the five psychiatric factors utilized in the primary analysis, using the same Bonferroni corrected threshold (p < 4.5 ×10−4) to identify significant inter-factor correlations.”"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have significantly improved their manuscript and addressed most of my comments. I only have a few more suggestions:\nI agree that the phewas using the Penn Medicine BioBank provides external validation for the transdiagnostic physical illness factor. However, as shown in Fig. 6, there were no significant associations with any mental phenotype, unlike the numerous significant associations observed with physical traits across multiple domains. How come? This seems somewhat at odds with the reported shared genetic risk between psychiatric and physical illness factors.\nI appreciate the more nuanced discussion in the updated manuscript, but I still believe that the claim in the abstract and discussion that the study calls into question the bifurcation of psychiatric and physical conditions should be toned down. Although the study clearly demonstrates shared transdiagnostic genetic signal between mental and physical disorders, disease classification also reflects other dimensions like clinical presentation, affected organ systems, etc. As discussed by the authors, the observed correlations could e.g. be driven by personality traits that affect risk of multiple disorders across different domains, like risky behavior. This doesn’t necessarily imply that mental and physical disorders are equivalent in kind.\nMinor:\n- I couldn’t find information about the sample size for the Penn Biobank\n- Consider replacing the «.» with «,» on bottom of page 6?: «Consistent with the large, pervasive genetic correlations observed in the systems-level analysis (see Fig. 7; Supplementary Table 5). We…»\n- Fig. 3 replace MDD with MD", "reviewer_2": "The authors have faithfully responded the reviewer's comments as well as they can. Therefore, this manuscript should be acceptable for the publication in Nature Communications.\nTwo little suggestions:\n1. Should the author consider merging 7 figures into 5 figures?\n2.The authors found shared genetic liability across psychiatric and physical illness. In the limitation part, the authors should also talk about the validation of biological mechanism."}, "author_response": "We thank the reviewer for pointing out the lack of significant PheWAS associations for mental phenotypes with instruments derived from the common physical illness risk factor. These null associations are likely a reflection of power issues, which is supported by two pieces of evidence. First, we find that while the common medical factor of 21 diseases (implied sample size (𝑁̂) = 106,118.7), incorporated in the last round of revisions, is not predictive of psychiatric illness in the PheWAS, the prior, better-powered common factor of 34 indicators (𝑁̂ =129,650.2) yields significant associations with psychiatric phenotypes in the same PMBB cohort. Second, the ratio of cases to controls was much smaller for many of the psychiatric disorders, thereby yielding a smaller effective sample size in the predicted psychiatric outcomes in PMBB and consequently lower power to detect associations relative to the other medical phenotypes. We have updated the manuscript to incorporate results from the 34-indicator model for the PMBB PheWAS and clarify these power-based interpretations, now writing on ps. 5-6 of the manuscript:\n“We then conducted a phenome-wide association study (PheWAS) in the Penn Medicine BioBank (PMBB; N = 29,883), testing for associations between multivariate GWAS results for the transdiagnostic physical illness factor (expected sample size [𝑁̂] = 106,118.7) and a wide range of disease phecodes (see Fig. 4C). Providing external validity for our risk loci, we found substantial effects across multiple physical illness systems (e.g., chronic airway obstruction, type 2 diabetes, obesity, hypertension, and ischemic heart disease). Additional significant associations emerged for sleep-related breathing disorders (sleep apnea, shortness of breath), emphysema, GERD/diseases of the esophagus, chronic pain, and viral hepatitis C. PheWAS results revealed no significant associations with psychiatric outcomes. This result may be due to more limited power to predict outcomes, such as psychiatric disorders, not used to define the 21-disorder common physical factor in Genomic SEM. Further, the ratio of cases to controls was much smaller for many of the psychiatric disorders in the PMBB, resulting in a much smaller effective sample size, relative to many of the physical outcomes. We provide a summary of PheWAS results and the number of cases for all associations with a p < .05 in Supplementary Table 5.\nTo investigate this power-based interpretation of null psychiatric associations, we went on to examine better powered 34-indicator common medical factor (𝑁̂ =129,650.2), expanded to include a combination of symptoms and disease proxies (Supplementary Fig. 10a). We find that this updated medical factor significantly predicted psychiatric illnesses (see Supplementary Fig. 10b). However, the lack of associations for the 21-disorder model should still be kept in mind when interpreting genetic correlations reported below.”\nI appreciate the more nuanced discussion in the updated manuscript, but I still believe that the claim in the abstract and discussion that the study calls into question the bifurcation of psychiatric and physical conditions should be toned down. Although the study clearly demonstrates shared transdiagnostic genetic signal between mental and physical disorders, disease classification also reflects other dimensions like clinical presentation, affected organ systems, etc. As discussed by the authors, the observed correlations could e.g. be driven by personality traits that affect risk of multiple disorders across different domains, like risky behavior. This doesn’t necessarily imply that mental and physical disorders are equivalent in kind.\nWe thank the reviewer for pointing out additional refinements to our arguments, calling for additional nosological frameworks based on statistically inferred systems of risk-sharing. We refine our argument in the abstract, discussion, and conclusion, and now write on p. 1:\n“Our results reveal pervasive risk sharing between specific groups of psychiatric and physical conditions and suggest the need for additional nosological frameworks.”\np. 8:\n“Although this organ- and symptom-based classification remains essential for biological classification and mechanistic insight, it does not exhaust the array of useful nosologies. Our findings provide further context to this historical dichotomy, highlighting an avenue for disease organization based on systems of shared risk from statistically inferred factors.”\n“This framework need not imply that psychiatric and physical disorders are identical in clinical expression or management, but sheds light on common etiological factors that can inform prevention strategies targeting liability rather than its diverse downstream manifestations.”\nand p. 9:\n“This study enriches our knowledge of the multivariate systems of relationships underlying complex multimorbidity, with implications for interpreting emergent genetic signals and considering additional nosological frameworks that capture the systems of interconnected risk between two theoretically dichotomized sets of disorders.”\nMinor:\n- I couldn’t find information about the sample size for the Penn Biobank\nWe now provide the total N (29,883) on p. 5 and a supplementary table with PheWAS results and case/control information for all outcomes with p < .05. Further information on the Penn Biobank can be found at https://pmbb.med.upenn.edu.\n- Consider replacing the «.» with «,» on bottom of page 6?: «Consistent with the large, pervasive genetic correlations observed in the systems-level analysis (see Fig. 7; Supplementary Table 5). We…»\nWe have now replaced the «.» with «,» on bottom of page 6 and thank the reviewer for finding this misplacement.\n- Fig. 3 replace MDD with MD\nWe have now updated Fig. 3 to display MD rather than MDD.\nReviewer #2 (Remarks to the Author):\nThe authors have faithfully responded the reviewer's comments as well as they can. Therefore, this manuscript should be acceptable for the publication in Nature Communications.\nTwo little suggestions:\n1. Should the author consider merging 7 figures into 5 figures?\nWe think this is a great suggestion and have now merged Figs. 4-6 into a single Fig. 4 for a total of five main text figures.\n2.The authors found shared genetic liability across psychiatric and physical illness. In the limitation part, the authors should also talk about the validation of biological mechanism.\nWe agree that this is an important area of future research. We now write in the Limitations and Future Directions section on p. 9:\n“Future work should integrate cross-ancestry fine-mapping with functional genomics and multi-omic approaches to prioritize causal mechanisms, followed by causal mediation analyses to adjudicate vertical vs. horizontal pleiotropy.”"}]

Decision letter