Mapping rare protein-coding variants on multi-organ imaging traits

Domain: biology
Published: 2025-12-23
Editorial outcome: minor_revision

Abstract
Human organ structure and function are important endophenotypes for clinical outcomes. Genome-wide association studies (GWAS) have identified numerous common variants associated with phenotypes derived from magnetic resonance imaging (MRI) of the brain and body. However, the role of rare protein-coding variations affecting organ size and function is largely unknown. Here we present an exome-wide association study that evaluates 596 multi-organ MRI traits across over 50,000 individuals from the UK Biobank. We identified 107 variant-level associations and 224 gene-based burden associations (67 unique gene-trait pairs) across all MRI modalities, including PTEN with total brain volume, TTN with regional peak circumferential strain in the heart left ventricle, and TNFRSF13B with spleen volume. The singleton burden model and AlphaMissense annotations contributed 8 unique gene-trait pairs including the association between an approved drug target gene of KCNA5 and brain functional activity. The identified rare coding signals elucidate some shared genetic effects across organs, prioritize previously identified GWAS loci, and are enriched for drug targets. Overall, we demonstrate how rare variants enhance our understanding of genetic effects on human organ morphology and function and their connections to complex diseases.

Introduction
Magnetic resonance imaging (MRI)-derived traits enable us to study the structure, function, and abnormalities of human organs in vivo. Many of these traits serve as established endophenotypes implicated in complex diseases and related traits. Therefore, it is of great interest to uncover genetic effects using imaging data to better understand the biology of human organs in health and disease. Recent genome-wide association studies (GWAS) have successfully identified common variants associated with multi-organ imaging traits, including brain structural, diffusion, and functional MRI1–8, as well as cardiovascular magnetic resonance imaging (CMR)9–11 and abdominal MRI12–14. However, a limitation of common variant signals identified by GWAS is that they often reside in non-coding regions and exhibit small effect sizes, which complicates the direct derivation of biological insights or the identification of causal genes15–19. By focusing on rare variants in protein-coding regions of the genome, whole exome sequencing (WES) studies aim to directly identify genes of interest.

Although existing exome-wide association studies16,17,20,21 (ExWAS) have identified associations for some imaging traits, our knowledge of the rare variant genetic architectures of human organs and their roles in various diseases remains substantially limited. Specifically, most previous ExWAS have focused on a single organ and/or a small set of imaging traits, lacking a multi-organ perspective that simultaneously explores the genetic effects of the human brain and body. For example, Park et al.20 analyzed CT imaging-derived hepatic fat, Haas et al.14 studied liver fat based on abdominal MRI using machine learning, Jurgen et al.21 studied several CMR traits, and Backman et al.16 and Karczewski et al.17 included brain MRI traits in the UK Biobank (UKB) WES studies of a wider range of phenotypes. In addition, these studies might be underpowered because of the relatively limited sample sizes, which are critical for rare variant association studies. Therefore, analyzing more and refined imaging traits, spanning multiple organs, in larger sample sizes will yield deeper insights into the genetic effects of rare variants across the whole body and shared genetic effects across organs10.

Here, we conducted ExWAS for 596 MRI traits derived from the brain, heart, liver, kidney, and lung (Supplementary Data S1) of over 50,000 participants from the UKB study. We used an internal discovery-replication design to make the best use of available data resources to identify rare variants and genes associated with human organ structure and function (Fig. 1). Specifically, all the primary findings and the corresponding interpretations are based on the joint sample for higher statistical power and more stable summary statistics. The joint sample was then split into the internal discovery sample and validation sample to partially verify the robustness of the findings within the UKB study (Fig. 1, Methods, and Supplementary Note). We evaluated various functional annotation approaches for missense variants in gene-level set-based association testing, including conventional methods such as SIFT22, PolyPhen2 HDIV23, PolyPhen2 HVAR23, LRT24, and MutationTaster25, as well as the deep learning-based method AlphaMissense26. We compared our results with previous GWAS on the same set of MRI traits in order to provide shared common and rare variant evidence for a gene’s involvement in a trait. We also used burden heritability regression (BHR)27 tests to characterize the genetic architecture for ultra-rare coding variants (minor allele frequency [MAF] < 1 × 10−4) across functional classes. In summary, our study identified exome-wide associations for multi-organ structure and function, providing a valuable source of evidence that could be useful in drug discovery and clinical therapeutics. These findings may also enhance our understanding of the complex interrelations between the human organs, health, and disease.

Results
Overview of variant-level associations with multi-organ MRI traits
Based on European individuals from UKB phases 1 to 5 MRI data (released up to late 2023, average n = 40,038; Fig. 1), we identified 107 rare (MAF < 0.01) variant-level associations (Supplementary Data S2) using a conservative P-value threshold of 2.8 × 10−10 (Bonferroni adjusted for all variant-trait tests with minor allele count [MAC] > 5 as 0.05/178,280,016, Methods) (Fig. 2A). Only predicted loss-of-function variants (pLoF) and missense variants were included in our study. There were 75 gene-trait pairs between 24 unique genes and 62 MRI traits, including 3 abdominal MRI traits, 7 CMR traits of the heart and the aorta, 10 brain structural MRI traits (regional brain volumes), as well as 42 brain diffusion MRI traits (diffusion tensor imaging [DTI] parameters) (Fig. 2B). Additionally, 138 associations (92 unique genes and 107 MRI traits) showed suggestive evidence, surviving a more liberal P-value threshold of 1 × 10−8 (Supplementary Data S3, Methods). These associations spanned across all categories of MRI traits (Supplementary Data S2). As expected, larger effect sizes were linked with lower MAF, which revealed the process of negative selection17,28,29 (Fig. 2C).

Due to the lack of independent data sources and the inherent rarity of the variants, identifying rare variant associations typically requires large sample sizes, and findings are difficult to fully replicate. To assess the reliability of our findings derived from the joint sample and maximize the use of available data, we adopted an internal replication procedure30 (Fig. 1 and Methods). Specifically, we split our joint sample to the internal discovery dataset and validation dataset for internal replication. The discovery dataset only included individuals from phases 1 to 3 of MRI release (average n = 30,739), while the validation dataset is independent of it by design and contained individuals from phases 4 and 5 of the MRI release (Method). For the purpose of replication analysis, summary statistics from the discovery sample will be separately compared with results from the validation sample and the joint sample (Method and Supplementary Note). When we restricted our study to the discovery sample, which is the subset of the joint sample, 59 associations passed the Bonferroni threshold of P < 3.5 × 10−10 (Bonferroni adjusted for all variant-level tests with MAC > 5 in this subset, 0.05/142,935,657, Methods) (Supplementary Data S4) and were used for internal replication. Among these 59 associations, 48 had mutations with MAC > 5 suitable for association effect = 0.028 L, effect = 0.42 s.d. units, 95% CI = [0.34, 0.49], P = 2.0 × 10−28), while the latter was the top hit among all variant-level associations. TNFRSF13B is a well-known risk gene of common variable immunodeficiency and mutations in TNFRSF13B typically occur in patients who developed splenomegaly31–33. In addition, the missense variant rs72553883 was associated with a series of blood-related traits including platelet, myeloid white cell, and lymphoid white cell indices in a previous GWAS34 using genotyping array data while similar associations between rare variants in TNFRSF13B and multiple blood biomarkers were also observed in Karczewski et al.17 using exome data (Fig. 2D). Consistent with these results on spleen abnormalities derived from immunodeficiency and the central role of the spleen in blood filtration and blood cell turnover, our finding points out a direct association between spleen volume and missense mutations in TNFRSF13B. Another missense variant in the manganese transporter SLC30A10 (rs188273166, MAF = 0.09%, χ2 = 94.79, 1 df, effect = 63.72 ms35,36 [ms], effect = 1.18 s.d. units, 95% CI = [0.94, 1.42], P = 2.1 × 10−22) was associated with higher corrected T1 liver iron. SLC30A10 is known to be involved in maintaining the manganese level37. This association aligned with recent GWAS that reported the role of rs188273166 in hypermagnesemia symptoms and a SLC30A10-targeted study found its positive effect on corrected T1 liver iron38,39. Interestingly, we also identified a missense variant in DDX51 associated with the fat-free muscle volume of right posterior thigh (rs200735214, MAF = 0.01%, χ2 = 41.06, 1 df, effect = −0.87 L, effect = −1.09 s.d. units, 95% CI = [−1.42, −0.76], P = 1.5 × 10−10), which constituted the only signal that passed our stringent P value threshold for human muscle/fat composition traits. Notably, at a more permissive P-value threshold at 1 × 10−8 as our suggestive evidence, we observed that rs200735214 in DDX51 was also associated with the fat-free muscle volume for total thigh (MAF = 0.01%, χ2 = 36.02, 1 df, effect = −2.26 L, effect = −0.90 s.d. units, 95% CI = [−1.19, −0.60], P = 1.95 × 10−9) among other signals for human muscle/fat composition traits such as rs534069616 in PLIN4 and posterior thigh muscle fat infiltration (MAF = 0.2%, χ2 = 36.94, 1 df, effect = −0.95%, effect = −0.40 s.d. units, 95% CI = [−0.53, −0.27], P = 1.2 × 10−9 for the left; MAF = 0.2%, χ2 = 34.52, 1 df, effect = −0.91%, effect = −0.38 s.d. units, 95% CI = [−0.51, −0.26], P = 4.2 × 10−9 for the right). DDX51, a member of the DEAD-box helicase family, plays a crucial role in ribosomal RNA processing and is essential for ribosome biogenesis40,41. Although no direct link between DDX51 and muscle function has been reported previously, the strong association identified in our exome data highlights its potential involvement in muscle composition traits. Meanwhile, PLIN4 encodes a protein known as Perilipin4, which belong to the perilipin family of proteins that coat lipid droplets and involve in regulating lipid metabolism42. Perilipins are essential for the proper storage and release of lipids within cells. In particular, PLIN4 has been shown to be highly expressed in skeletal muscle and was found at periphery of skeletal muscle fibers42, which is consistent with our findings.

Variant-level tests identified associations for brain and heart MRI traits
We highlight some rare variants associations with regional brain volumes, brain DTI parameters, and CMR traits of the heart and aorta. A missense variant in ADRA1A was associated with the volume of the left ventral diencephalon (rs771722367, MAF = 0.01%, χ2 = 42.67, 1 df, effect = 683.56 mm3, effect = 1.28 s.d. units, 95% CI = [0.90, 1.67], P = 6.5 × 10−11), though the P value for this variant and right ventral diencephalon was P = 3.6 × 10−8 (χ2 = 30.37, 1 df), which did not pass our liberal threshold at 1 × 10−8. ADRA1A was previously reported to be associated with schizophrenia43,44. The ventral diencephalon area was also associated with a missense variant in PKD1 (rs1181041827, MAF = 0.01%, χ2 = 45.71, 1 df, effect = 923.88 mm3, effect = 1.73 s.d. units, 95% CI = [1.21, 2.24], P = 5.2 × 10−11 for right ventral diencephalon; MAF = 0.01%, χ2 = 43.11, 1 df, effect = 945.24 mm3, effect = 1.77 s.d. units, 95% CI = [1.26, 2.29], P = 1.4 × 10−11 for left ventral diencephalon). For DTI parameters, two missense variants (rs2652098 and rs143368552) in Versican (VCAN) contributed the largest number of associations (63 in total, effect range = [−0.48, 0.42] s.d. units, P < 2.8 × 10−10). Figure 2E visualizes the associations between these two missense variants and 12 mean diffusivity traits. Common variants in VCAN have been extensively associated with white matter traits in previous GWAS2,5,8,45–47. Notably, Elliot et al.5 identified an intronic variant rs67827860 with strong associations with DTI traits. Thus, we conducted conditional analysis for all the 63 associations between these two variants in VCAN and DTI traits by conditioning on this nearby common variant as a covariate (Methods). We found that 85.7% of these associations still survived our stringent threshold at 2.8 × 10−10 while all the associations passed our suggestive threshold at 1 × 10−8, which support the identified missense variants in VCAN as independent signals in addition to the previously identified common variant (Supplementary Data S5). VCAN plays a pivotal role in various neural processes, which may influence the pathophysiology of neurological disorders such as multiple sclerosis48–50. Another missense variant rs201680145 in NOTCH3 was also associated with multiple DTI parameters (effect range = [−1.50, 1.70] s.d. units, P < 2.5 × 10−10). Previous studies have revealed the role of NOTCH3 in white matter hyperintensities and several neurodegenerative diseases51–53.

We found associations between a missense variant rs189569984 in RBM20 and three CMR traits, including left ventricular end-systolic volume (LVESV), right ventricular end-systolic volume (RVESV) (MAF = 0.87%, χ2 = 64.60, 1 df, effect = −3.77 mL, effect = −0.20 s.d. units, 95% CI = [−0.25, −0.15], P = 9.1 × 10−16 for LVESV; MAF = 0.87%, χ2 = 44.88, 1 df, effect = −3.37 mL, effect = −0.16 s.d. units, 95% CI = [−0.21, −0.11], P = 2.1 × 10−11 for RVESV), and left ventricular ejection fraction (LVEF) (MAF = 0.87%, χ2 = 56.37, 1 df, effect = 1.49 %, effect = 0.24 s.d. units, 95% CI = [0.18, 0.30], P = 6.0 × 10−14). Mutations in RBM20 were previously known to be related to cardiovascular diseases including heart failure and dilated cardiomyopathy54–56, and were associated with LVEF in a previous GWAS on CMR traits9. Our rare variant analysis prioritized RBM20 and provided additional evidence for its impact on heart structure and function. Furthermore, we found that ascending aorta maximum/minimum areas were associated with a missense variant in ANO1 (rs201870990, MAF = 0.61%, χ2 = 42.35, 1 df, effect = 41.95 mm2, effect = 0.22 s.d. units, 95% CI = [0.15, 0.29], P = 7.6 × 10−11 for maximum; MAF = 0.61%, χ2 = 45.61, 1 df, effect = 42.05 mm2, effect = 0.23 s.d. units, 95% CI = [0.16, 0.29], P = 1.4 × 10−11 for minimum). ANO1 was among the loci identified by a recent GWAS57 on ascending aorta diameter (see Fig. 1 in their study57) but was not pointed out or discussed explicitly. Indeed, ANO1 was also reported to be effective in preventing cardiac fibrosis and maybe a potential target for therapy58,59. However, future research might elucidate its role in aortic development and/or geometric remodeling. In summary, our analysis of rare variants directly prioritized and implied a small set of genes related to human organs. Due to the rarity of these variants and their lack of linkage disequilibrium (LD) with common variants, they could not be effectively studied in previous GWAS that focused on common variants.

Gene-based burden tests identified complementary signals
Gene-based burden tests enable us to study the collective effects of rare variants within a gene, thus boosting power. In burden tests, the variants of interest defined by annotation and MAF (Methods) are collapsed into a burden mask. Then, the burden score is calculated as the maximum genotype dosage across the single variant sites in the mask, and this score is tested for association (Methods). However, the involved burden models pose challenges to the adjustment of multiple testing since they have unknown correlated structure. Thus, we prioritized a P-value threshold at 1 × 10−9 based on two empirical null distributions60,61 (Methods). A total number of 17,474 genes are available for burden tests in the joint sample and discovery sample. Using the same joint sample as in variant-level tests, we identified 224 significant associations in gene-based tests (P < 1 × 10−9, Methods; Supplementary Data S6 for all the significant results). As different burden models or MAF cutoffs may implicate the same gene-trait associations (Methods), we further summarize the nonredundant results of 67 unique gene-trait pairs60 (involving 26 genes and 57 MRI traits) in Table 1 (a more comprehensive version is provided in Supplementary Data S7). For suggestive evidence, we additionally put all the associations that passed a more relaxed P-value threshold at 1 × 10−8 in Supplementary Data S8. Manhattan plot for all the gene-level associations were presented in Fig. 3A, across all MRI phenotype categories. Three gene-trait pairs overlapped across variant-level and gene-level tests (two associations between ANO1 and max/min aorta area, and an association between TNFRSF13B and spleen volume), indicating the robustness of these gene-trait pairs as both a single rare variant in this region and the joint effects of variants across the region are strong enough to drive an association. Consistent with observations in previous ExWAS17,61, we identified 21 additional genes, highlighting the power of gene-based burden tests in rare variant association studies. In particular, we found 3 genes and 18 gene-trait associations for brain functional MRI (fMRI) traits, which did not have any signals in variant-level tests. To evaluate the reportability of our findings, we wanted to follow a similar procedure in variant-level tests to perform an internal replication. However, many rare mutations, particularly those with a MAF of less than 0.01%, accounted for a large proportion of our discoveries but were not observed in the smaller independent dataset from phases 4 to 5 (average n = 8,989). Therefore, we only made use of the internal discovery sample and examined whether the associations in the combined phases 1 to 5 sample had smaller P-values than those in the phases 1 to 3 data, similar to internal-replication analysis for variant-level tests (Methods and Supplementary Note). When we restricted the discovery sample to individuals in phases 1 to 3, there were 96 significant associations for 25 unique gene-trait pairs (P < 1 × 10−9, Supplementary Data S9). Within these 96 associations, 79.2% (76/96, 19 gene-trait pairs) had smaller P-values in the combined phases 1 to 5 sample, all of which had concordant directions of effects (Supplementary Data S9). The detailed summary statistics and results of internal replication analysis for gene-based associations have been provided in Supplementary Data S9. We discuss some of the 67 gene-trait pairs in the following two sections.

Genes associated with abdominal, brain, and heart MRI traits in burden tests
As shown in Table 1, we identified 5 genes for abdominal MRI traits, 7 for CMR traits, 7 for regional brain volumes, 8 for brain DTI parameters, and 3 for brain fMRI traits. Figure 3B highlights the genes with large effect sizes. Below, we highlight several interesting findings for each trait category.

We found that 4 out of 5 signals of abdominal MRI traits from the combined phases 1 to 5 sample also present in the phases 1 to 3 data (Supplementary Data S4). The highest hit among all the gene-level signals was the association between spleen volume and TNFRSF13B when aggregating the effects of pLoF and Alpha damaging (Methods) missense variants (MAF = 1.4%, χ2 = 246.24, 1 df, effect = 0.027 L, effect = 0.40 s.d. units, 95% CI = [0.35, 0.45], P = 1.7 × 10−55) (Fig. 2D). In addition to TNFRSF13B, we observed that SH2B3 was also associated with spleen volume (MAF = 0.4%, χ2 = 63.40, 1 df, effect = 0.025 L, effect = 0.38 s.d. units, 95% CI = [0.29, 0.48], P = 1.7 × 10−15), where spleen abnormalities (e.g. splenomegaly) and associated blood traits have been previously reported62–65. For muscle measurements, we found an association between fat-free muscle volume of the posterior thigh and TANC1 (MAF = 0.01%, χ2 = 37.63, 1 df, effect = −1.06 L, effect = −1.32 s.d. units, 95% CI = [−1.74, −0.90], P = 8.5 × 10−10). TANC1 had genetic overlaps with the identified loci in previous GWAS on heel bone mineral density66 and its role in muscle development and rhabdomyosarcoma has been discussed in previous studies67,68.

As shown in Fig. 4A, TTN was associated with 8 CMR traits, including LVESV (MAF = 0.39%, χ2 = 109.60, 1 df, effect = 7.34 mL, effect = 0.40 s.d. units, 95% CI = [0.32, 0.47], P = 1.2 × 10−25), LVEF (MAF = 0.39%, χ2 = 140.48, 1 df, effect = −3.48%, effect = −0.56 s.d. units, 95% CI = [−0.66, −0.47], P = 2.1 × 10−32), as well as global and regional peak circumferential strain measurements (effect range = [0.32, 0.60] s.d. units, P range = [1.2 × 10−21, 2.4 × 10−11]). Figure 4C visualizes the original effect sizes (in percentage) of TTN on different regional peak circumferential strain measurements and a global peak circumferential strain measurement. We observed that they generally produced similar and consistent results. These results make sense as TTN is a well-known gene associated with cardiac structure9 and cardiovascular diseases such as heart failure69, dilated cardiomyopathy70,71, atrial fibrillation72, supraventricular tachycardia, and mitral valve disease21. In addition to the previously known exome-wide association between TTN and LVESV21, our findings suggest a broader influence of this gene on cardiac structure and function, using a larger sample size and an expanded set of CMR traits.

We would like to highlight our exome-wide associations with CMR traits of aorta. Notably, three genes (ANO1, COL21A1, GEM) were associated with both maximum and minimum area of ascending aorta, while PLCE1 was associated with maximum and minimum area of descending aorta (Fig. 4B,C). In particular, both COL21A1 and GEM lower the size of the ascending aorta (effect range = [−61.16, −54.95] mm2, P range = [1.7 × 10−13, 6.4 × 10−13] for COL21A1 and effect range = [−87.68, −87.06] mm2, P range = [1.3 × 10−12, 1.06 × 10−11] for GEM, respectively). Rare variants in COL21A1 and GEM were reported to be linked to higher pulse blood pressure in a previous rare variant association study of blood pressure73. A smaller aorta leads to a wider pulse pressure because it increases aortic characteristic impedance, which is highly dependent on aortic diameter. Our results on these aorta area-lowering genes further provide the evidence for the underlying mechanisms of previous observations. Moreover, GEM is involved in cardiac activities74 and is related to heart failure74,75. Another study76 pointed out that such effects of GEM could potentially be used as a gene therapy target in treating ventricular arrhythmias and heart failure. Our discovery indicated that GEM may also impact the structural integrity of the ascending aorta, thus improving our understanding of the role of GEM in cardiovascular health. In contrast, missense variants in ANO1 have positive effect. The role of ANO1 has been discussed in our variant-level analysis and were captured these associations again by aggregating the missense variants in ANO1 using burden test (effect range = [32.43, 32.79] mm2, P range = [1.7 × 10−13, 6.4 × 10−13]), suggesting the similar effects and aligned directions across a group of missense variants within ANO1. In addition, TAGLN was also associated with descending aorta minimum area (MAF = 0.05%, χ2 = 38.71, 1 df, effect = 62.28 mm2, effect = 0.67 s.d. units, 95% CI = [0.46, 0.88], P = 4.91 × 10−11) while its association with descending aorta maximum area was observed at a more liberal P value threshold at 1 × 10−8 (Supplementary Data S8). At the bottom of Fig. 4C, we found that the effects (in mm2) of different genes on the aorta area could largely vary in terms of both directions and the sizes. GOLM1 was found to be associated with descending aorta distensibility (MAF = 0.01%, χ2 = 44.23, 1 df, effect = 2.85 mmHg, effect = 2.5 s.d. units, 95% CI = [1.77, 3.24], P = 2.93 × 10−11). A recent study demonstrated that GOLM1 may present a potential therapeutic target for treat sepsis-induced cardiac dysfunction in animal models77. GOLM1 was also related to Alzheimer’s disease and the implied cognitive deficits78. However, its role in aortic wall remodeling requires further study. Using GTEx data resources79, we found that four (ANO1, COL21A1, GEM, TAGLN) out of five genes (ANO1, COL21A1, GEM, TAGLN, PLCE1) associated with ascending/descending aorta areas had moderate to high expression and overexpressed in aorta while PLCE1 had high expression in artery tissues in general, which further supports their roles in affecting the structure of aorta.

For brain MRI traits, we highlight in Fig. 3C, D the strongest association between PTEN and total brain volume with a large effect size when aggregating the pLoF and missense variants (MAF = 0.01%, effect = 386,703.40 mm3, effect = 2.77 s.d. units, 95% CI = [2.22, 3.32], P = 5.8 × 10−23). The distribution of total brain volume in cm3 for PTEN mutation carriers versus non-carriers is provided in Fig. 3C. Notably, as illustrated in Fig. 3D, different burden masks would yield associations with various strengths, and we observed that the signal became stronger if we aggregate ultra-rare (MAF < 1 × 10−4) pLoF variants and damaging missense variants in PTEN, which may not be observed and testable in the previous ExWAS16,17 with a sample size smaller than ours due to the low MAF. In a recent ExWAS80, PTEN was associated with autism spectrum disorder. Indeed, PTEN is a well-known gene that has been consistently linked to autism spectrum disorder81–84, brain development85–87, and brain clinical phenotypes including brain overgrowth88–90, with an enlarged brain volume representing an associated phenotype in these conditions. To our knowledge, no previous GWAS or ExWAS has linked PTEN to any brain MRI measurements, thus our analysis provided direct evidence connecting rare variants in PTEN with total brain volume. Another gene OMA1 was observed to be associated with 7 regional brain volume traits including superior frontal gyrus, lateral orbitofrontal cortex, and precentral gyrus for both sides of the brain and right superior parietal lobule (effect range = [−524.83, −171.02] mm3 units, P range = [1.77 × 10−14, 6.35 × 10−10]). Notably, PTEN and OMA1 regulate the PTEN-induced kinase 1 (PINK1) and thus may be useful in preventing epileptogenesis91 while a protective role of OMA1 in neurodegeneration92–94 was previously reported. The remarkably strong association PTEN demonstrate that even within a general healthy cohort (as opposed to disease-focused or case-control studies), we can identify biological links between genes associated with neurodevelopmental disorders using MRI data.

For DTI parameters, in contrast to variant-level analysis where missense variants in VCAN accounted for the most associations, our gene-based burden tests revealed a more diverse range of effect on white matter, with 8 distinct genes accounting for 10 gene-trait associations. For brain fMRI traits, PLCE1, GLUL, and KCNA5 were identified to be associated with 18 phenotypes, where PLCE1 contributed to 16 of them. We found that many gene-trait pairs would not have been identified if singleton burden models or the recent deep learning-based AlphaMissense26 burden models were not included. We will provide more detailed discussions about these interesting associations in the following sections.

Overall, gene-based burden tests examined the aggregated effects of rare variations within genes and advanced our understanding of the underlying genetic dispositions across the abdomen, heart, and brain organs. Our findings prioritized a set of rare genes previously unknown to be associated with human organs. We further discussed examples of genetic overlaps between these MRI signals and other health-related phenotypes, which may help to understand complex diseases by providing exome-wide insights into the genetic mechanisms underlying multiple human organs including abdominal, cardiac, and neurological pathologies.

Pleiotropic effects of PLCE1 and COL21A1 on brain and heart MRI traits
In our gene-based burden tests, PLCE1 was associated with the largest number of phenotypes including 16 brain fMRI traits (including 5 traits from parcellation-based approach95 and 11 traits from whole brain spatial independent component analysis [ICA]96–98) and 2 CMR traits (descending aorta maximum/minimum areas) (Fig. 4B). In addition, COL21A1 was associated with both brain cerebrospinal fluid (CSF) volume and ascending aorta maximum/minimum areas (Fig. 4B). We found that a recent GWAS99 discovered the associations of descending and ascending aorta areas separately with PLCE1 and COL21A1 while PLCE1 was also found in previous GWAS3,5,100 for brain fMRI. PLCE1 encodes a phospholipase (PLCε) that catalyzes the hydrolysis of phosphatidylinositol-4,5-bisphosphate to generate inositol 1,4,5-triphosphate and diacylglycerol. They are common second messengers regulating multiple cellular processes, including cell activation, growth, differentiation, and gene expression101. Interestingly, prior mechanistic study identified a causal role for PLCε on the development of thoracic dilation102 and dissection103. Similarly, COL21A1 encodes the alpha-1 chain of collagen 21, a known extracellular matrix component in the arterial wall, secreted by vascular smooth muscle cells. While the previous GWAS99 suggested distinct genetic bases between the ascending and descending aorta and linked aortic traits to brain small vessel disease, it did not establish direct associations of these genes with other brain MRI traits. Our results provide exome-wide evidence highlighting the associations of PLCE1 and COL21A1 with both brain and heart MRI traits. Additionally, a previous study suggested that the observed GWAS effects of PLCE1 on fMRI traits may be blood-derived100. Our exome-wide evidence supporting PLCE1’s role in both aortic and fMRI traits aligns with this hypothesis. Therefore, it might be helpful to integrate cardiac features and cardiac-related genetic information when analyzing fMRI data, given that fMRI phenotypes are blood oxygen level-dependent and are influenced by blood flow.

The connection between the heart and brain has increasingly garnered attention10,104. Moreover, the role of hemodynamic role of the aorta involved the cushioning of pressure and flow pulpability caused by the intermittent ventricular ejection, and an impaired cushioning function for changes in aortic stiffness and/or diameter has been identified as a key mechanism for target organ damage, including the brain and the heart105. The associations we identified, involving rare coding mutations in the same genes (PLCE1 and COL21A1) linked to heart, aortic and brain structure and function, suggest that similar biological mechanisms may underlie the shared pathways of these organs that influence both cardiac and neurological health.

Associations uniquely identified by singletons and AlphaMissense
Singleton variants are those observed only once in the study cohort. AlphaMissense26 is a recent deep learning-based method that offers a novel approach to annotating damaging missense variants. The Singleton burden model captured the rarest category in our association tests and may be of special interest16,106, and the incorporation of AlphaMissense offered additional genetic findings in our study. Therefore, we would like to highlight the contributions of these two groups of burden models in this section. In our gene-based burden tests, 10 gene-trait pairs would not be detected without the singleton burden models or including AlphaMissense as part of our annotation resources. Among these signals, 8 of them did not even have any counterparts that passed the relaxed P < 1 × 10−8 threshold (Supplementary Data S10).

Singleton damaging missense variants in glutamate-ammonia ligase GLUL, which encodes the glutamine synthetase protein, were associated with functional connectivity of the somatomotor and secondary visual networks (MAF = 0.01%, χ2 = 37.54, 1 df, effect = 2.48 s.d. units, 95% CI = [1.68, 3.27], P = 8.9 × 10−10). GLUL knockout mice were reported to be neuronally affected in multiple regions including somatosensory and visual cortices107. Our singleton analysis provides consistent evidence in human genetics. In addition, a recent study suggested that visual-somatosensory integration may be a biomarker for preclinical Alzheimer’s disease108. Together, rare mutations in GLUL may play a potential role in human neurodegenerative diseases109–112. Another gene, KCNA5 related to potassium voltage-gated channel, was associated with functional activity in the subcortical-cerebellum network when we combined the rare pLoF variants and Alpha damaging missense variants (MAF = 0.1%, χ2 = 44.21, 1 df, effect = 0.59 s.d. units, 95% CI = [0.41, 0.76], P = 2.9 × 10−11). Notably, voltage-gated potassium channels are essential for neurons and cardiac activities and previous GWAS have discovered the common variants mapped to KCNA5 associated with cortical surface7 and thickness113. In addition, a previous study on the rat cerebellum model revealed strong KCNA5 immunoreactivities in the cerebellar nuclei114. Consistent with these findings, we provided additional evidence for KCNA5’s role in brain function. Moreover, some studies also found close relationships between KCNA5 and heart diseases including atrial fibrillation115 and pulmonary arterial hypertension116. Indeed, KCNA5 is an approved drug target for the treatment of cardiac arrhythmias. Inspired by this additional context, we investigated whether KCNA5 was associated with any CMR traits. Interestingly, we found that when aggregating the singleton pLoF and damaging missense variants, KCNA5 was associated with a regional peak circumferential strain measurement (MAF = 0.02%, χ2 = 21.69, 1 df, effect = 5.11%, effect = 1.08 s.d. units, 95% CI = [0.63, 1.54], P = 3.2 × 10−6; “Ecc_AHA_4”, Supplementary Data S1) at a P-value threshold of 1.74 × 10−5 (0.05/2870, Bonferroni adjusted for 82 CMR phenotypes across all the 7 variant function classes and 5 MAF classes). These discussions further suggested that multi-organ imaging genetic studies would bring about insights to the complex interplay across the brain-heart system. In summary, our results underscore the efficacy of our gene-based burden test approach and highlight the importance of using multiple MAF cutoffs and innovative annotation tools.

Comparison with previous ExWAS for brain MRI and CMR traits
In this section, we discuss the links between our gene-based rare variant signals and existing ExWAS results on MRI traits. To our knowledge, there were no rare variant associations reported for brain fMRI or abdominal MRI traits. The studies most similar to ours are those by Backman et al.16, which included regional brain volumes and DTI parameters among a broad range of phenotypes, and Jurgen et al.21, which focused on CMR traits, while Pirruccello et al.117 declared no findings of rare variant associations for aorta traits. Both studies used data from the UKB cohort but with smaller sample sizes compared to our current study. Importantly, we analyzed different brain MRI and DTI parameters, extracted from raw images using our own pipeline, compared to those investigated by Backman et al.16. Additionally, we incorporated a broader set of CMR traits, including strain and thickness metrics, expanding upon the traits analyzed by Jurgen et al.21.

For regional brain volumes and DTI parameters (corresponding to STR and dMRI traits in Backman et al.16), 5 genes (AMPD3, HTRA1, MYCBP2, RBL1, SCUBE2) whose associations passed our stringent P-value threshold (P < 1 × 10−9) were also significantly associated with their brain MRI traits. For example, associations between AMPD3 and mean MD of splenium of corpus callosum passed the stringent P value threshold both in our study and in Backman et al.16. SCUBE2 was significantly associated with the regional volume of the left cerebellum exterior in our analysis, while its association with the volume of the cerebellum cortex in the left hemisphere was reported in their result as suggestive evidence. MYCBP2 was associated with the external capsule in both our analysis and Backman et al.16. Additionally, several associations between PLEKHG3 and multiple dMRI traits were identified in Backman et al.16 and two gene-trait pairs within these associations passed our relaxed threshold (P < 1 × 10−8). For CMR traits, we replicated the only exome-wide significant association between TTN and LVEF in Jurgen et al.21, while their suggestive association between TTN and LVESV also passed our stringent threshold. In addition, we found 7 associations between TTN and 7 peak circumferential strain metrics. We summarized the overlapping associations in these previous ExWAS studies in Supplementary Data S11.

Concordant evidence with GWAS signals
Associations with rare coding variants could prioritize genes among the numerous loci identified in GWAS for polygenic complex traits by providing concordant evidence16,21. In this section, we leveraged the previous GWAS summary statistics1–3,10,118 of all the same 596 imaging traits to compare the identified signals between common variants and rare coding variants (Methods). As expected, we observed convergent evidence for all categories of MRI traits. Briefly, more than half of the 174 rare coding signals (53.4% = 93/174, that is, 64/107 variant-trait associations and 29/67 gene-trait pairs) were within the 1 Mb range of the independent GWAS signals, which is consistent with the observation in a previous large-scale UKB phenotypes screening16. Based on our results, we found that variant-level associations with regional brain volumes and CMR traits were not within the neighborhood of GWAS signals, while all categories of MRI traits had at least one rare coding association within the 1 Mb of GWAS signal for signals from gene-level tests.

We zoomed into these shared signals between GWAS and ExWAS to provide more detailed insights. The association between SH2B3 and spleen volume identified in our gene-level burden test was further supported by an independent GWAS signal rs2239194, which is an expression quantitative trait locus (eQTL) in spleen tissue79. For DTI parameters, two missense variants (rs2652098 and rs143368552) in VCAN contributed the most to the variant-level associations and all these associations were within 1 Mb of the GWAS signals. However, similar to findings from a previous GWAS study47, we did not identify any brain tissue-related eQTLs. Associations between AMPD3 and DTI parameters of splenium of corpus callosum were also accompanied by signals from GWAS. These common variants are eQTLs of multiple brain tissues including cerebellum and spinal cord where brain white matter also presents. For fMRI traits, all the associations between PLCE1 and 18 MRI traits (16 brain fMRI traits and 2 CMR traits) were consistent with previous GWAS results, though there were no eQTLs in brain cortex tissues. On the other hand, we found that the only signal for brain volumes that fell into the GWAS loci was the association between COL21A1 and CSF volume. The corresponding common variants rs3857615 and rs9475654 were eQTLs in cerebellum and cerebellum hemisphere.

In addition, we identified gene-trait associations that were not reported in our previous GWAS of the same MRI traits but appeared in other studies (e.g., RBM20 and LVEF in Pirruccello et al.9). We also observed that some rare coding associations in our study were linked to GWAS loci of related traits within the same phenotype cluster. For example, though both ExWAS and GWAS revealed the associations between TTN and multiple regional peak circumferential strain measurements, the phenotype for the signals did not exactly match. In summary, rare variant associations with MRI traits not only identified novel signals but also prioritized genes among thousands of GWAS loci. These associations provide new insights and context for understanding the genetic basis of structure and function in the human brain and body. We provided these exome-wide signals that showed convergent evidence with previous GWAS signals in Supplementary Data S12 and S13.

Multi-organ imaging genetics and drug targets
We used the Therapeutic Target Database (TTD)119 to query potential drug target genes (Methods). Among the 26 unique genes identified by gene-based burden tests, 6 genes were approved drug targets (KCNA5 and ANO1) or in clinical trials (TTK, GEM, LPAR3, and TNFRSF13B). The identified genes in gene-based burden tests were significantly enriched for potential drug targets (6 out of 26, compared with 1735 out of 17,448 genes, OR = 2.7, P = 0.0391, Supplementary Data S14). Moreover, when including a broader set of potential drug genes that were documented in TTD to be reported in the literature, the enrichment became stronger (13 out of 26, compared with 3468 out of 17,448 genes, OR = 4.0, P = 5.737 × 10−4, Supplementary Data S14). Consistent with the procedure in a previous study61 and observations in other ExWAS16,61,120, our results suggest that rare coding variants identified in imaging genetics may facilitate the discovery of promising drug target genes121,122. Given the complicated nature of human organ interplay and the corresponding biological pathways involved, we would emphasize that such enrichment should be interpreted carefully as they might not directly point to the corresponding diseases and related traits for existing drugs. However, these associations with human organs could provide novel drug targets if they are supported by further evidence, and in turn these MRI traits themselves might be useful to validate the new drugs if the underlying mechanisms are well established123. Moreover, As endophenotypes for complex diseases, MRI traits can offer additional information on known drug targets, such as potential off-targets and/or side effects121,124, and also aid in drug repurposing (e.g., ANO1 and cardiac disease as discussed in the previous section58,59).

Burden heritability and genetic correlation
To investigate the rare coding genetic architecture for multi-organ MRI traits, we applied burden heritability regression (BHR)27 to estimate the heritability for pLoF and damaging missense variants (i.e., “int1 missense variants”, Methods), while the heritability for synonymous variants is also calculated, and it serves as negative control. Designed for rare variants with MAF < 1 × 10−3, BHR stratifies variants based on functional annotation and MAF to allow for different variant classes to have different mean effect sizes. Given our relatively small sample size (average n = 40,038), we restricted our analysis to the ultra-rare MAF bin (MAF < 1 × 10−4) which contained most of the coding variants in our analysis (Supplementary Data S15). As expected, ultra-rare pLoF variants consistently demonstrated higher heritability than damaging missense variants across all MRI categories (Fig. 5A). Supplementary Data S16 and S17 show the complete list of burden heritability for ultra-rare pLoF variants and damaging missense variants. In contrast, there is no evidence that synonymous variants have significant heritability, showing that the estimates from our results for rare variants are well calibrated (Supplementary Data S18). Furthermore, the highest average heritability of ultra-rare pLoF variants were observed in the DTI parameters category. For ultra-rare pLoF variants, the trait had the highest heritability was the caudal anterior cingulate (h² = 0.012, SE = 0.0035) among regional brain volumes, the mean axial diffusivity of the body of the corpus callosum (h² = 0.020, SE = 0.0048) among DTI parameters, functional connectivity between the second visual and auditory networks (h² = 0.020, SE = 0.0047) among brain fMRI traits, regional peak circumferential strain (h² = 0.016, SE = 0.0053) among CMR traits, and pancreas iron (h² = 0.019, SE = 0.0053) among abdominal MRI traits. We note that these burden h² estimates are relatively large, which is similar to the high common variant heritability for many imaging traits. This in general, support an endophenotype model where genes have more direct connections with imaging phenotypes than downstream health outcomes.

We further explored the burden genetic correlation pattern for ultra-rare pLoF variants between MRI traits (MAF < 1 × 10−4) and 186 complex traits and diseases (MAF < 1 × 10−5) downloaded from Genebass17 (Supplementary Data S19). We focused on 67 heritable MRI traits that passed the Bonferroni adjusted threshold at P < 8.4 × 10−5 (= 0.05/596). Consequently, we estimated 12,462 (67 × 186) pairs of burden genetic correlations. We found 3 pairs that passed the MRI-multiple testing threshold at P < 7.5 × 10−4 (0.05/67) and 591 pairs passed a nominal P-value threshold at P < 0.05. Although no pairs survived the most stringent P-value threshold adjusted for both the imaging traits and Genebass traits at P < 4.0 × 10−6 (0.05/12,462), we observed that there were interesting pairs in the top-ranking subset for continuous traits from Genebass (Supplementary Data S20). The results for binary traits can be found in Supplementary Data S21. Figure 5B shows the selected top hits of burden genetic correlation pairs, including insular volume and reaction time125, regional radial strain and systolic blood pressure, left atrium maximum volume and triglycerides to high-density lipoprotein cholesterol ratio, and abdominal fat ratio and birth weight. These burden genetic correlations between MRI traits and health-related traits may reveal a shared genetic basis in ultra-rare pLoF variants. While the current MRI sample size was slightly underpowered, future larger samples may enable a more robust quantification and confirmation of these genetic links.

Methods
Ethics and data resources
This study made use of the data from UK Biobank (UKB) involving ~500,000 participants aged from 40 to 69 when recruited between 2006 and 2010 (https://www.ukbiobank.ac.uk/). The UKB study received the ethics approval from the North West Multi-centre Research Ethics Committee (reference number: 11/NW/0382) with informed consent obtained by all the participants. The present study was under UKB application number 76139.
Whole-exome sequencing (WES) data were obtained from the UK Biobank and the detailed sequencing protocols can be found elsewhere18. By following a unified OQFE protocol121, raw sequencing reads were aligned to the GRCh38 reference genome, and small variants were called on a per-sample basis using DeepVariant. Per-sample gVCFs were then joint-genotyped using GLnexus to produce multi-sample project-level VCFs (pVCF) suitable for association analysis. We made use of the OQFE WES files in BGEN format provided by UK Biobank after these centrally pre-processing steps.

Multi-organ imaging phenotypes
The information of imaging phenotypes can be found in previous studies1–4,10. Briefly, we studied 596 imaging traits encompassing brain MRI traits1–4, heart CMR traits10, and abdominal MRI traits. Three major modalities of brain MRI were included. First, we included 101 regional brain volume traits1 derived from structural MRI. Second, we used 110 tract-averaged parameters2 from DTI capturing the microstructure of brain white matter. Third, for fMRI traits, we involved 76 node amplitude traits and six global functional connectivity traits based on ICA3 as part of resting fMRI traits; we also included 180 (90+90) parcellation-based4,95 resting fMRI traits and task fMRI traits. Regarding CMR traits, we made use of 82 traits10 including regional and local measurements from cardiac chambers and the aorta. For abdominal MRI traits, we contained 41 imaging traits based on MRI data of the liver, kidneys, lungs, pancreas, spleen, and body muscle/fat composition. For each phenotype, values greater than five times the median absolute deviation from the median value were removed as in our previous work1–4,10. A full list of these 596 traits can be found in Supplementary Data S1.

Internal replication and joint analysis design
Given the limited sample size (average n = 40,038 across imaging traits with non-missing data) and lack of independent data sources, we prioritized our signals through joint analysis30 and adopted an internal validation procedure to evaluate the robustness of our discoveries. Specifically, we first restricted our sample to include only participants from phases 1 to 3 of MRI release (average n = 30,739) and then excluded any individuals related to phases 1 to 3 participants in the rest of the sample from phases 4 to 5 of MRI release (average n = 8,989, relatives of the phases 1 to 3 were removed). The joint sample included all the individuals from phases 1 to 5 and was used to report significant results. To investigate the robustness of variant-level associations, we tested the associations on phases 1 to 3 individuals and phases 4 to 5 individuals separately. Then, focusing on the significant level and effect sign, we investigated (i) whether signals identified from the phases 1 to 3 sample could be replicated in the independent phases 4 to 5 sample and (ii) whether the signals identified from phases 1 to 3 sample had stronger evidence (i.e., smaller P values and concordant effect directions) in the joint sample. For gene-based burden tests, we noted that many rare mutations, especially those variants with MAF < 1 × 10−4, accounted for a large proportion of our discoveries but were not observed in the smaller independent dataset from phases 4 to 5 (average n = 8,989; even within a burden model with MAF cutoff at 0.01, the lack of those rare variants in the sample from phase 4 to phase 5 would make the burden associations incomparable). Thus, we only investigated whether the identified associations from phases 1 to 3 dataset had stronger evidence in the joint sample. Consequently, based on this internal replication procedure, we were able to expect the level to which our discoveries from the joint analysis could be replicated.

Quality control for UK Biobank exome data
We used Plink v. 2.0132 to conduct the quality control steps, restricting the sample to individuals with imaging data (n = 54,365). For exome sequencing data, we included all the variants that had a MAF below 0.01. Variant-level quality control excluded all the variants that had missing rate larger than 10% or had Hardy–Weinberg equilibrium P < 1 × 10−15. Sample-level quality control excluded any sample with missing rate over 10%. Consequently, no individual was excluded from the quality control process, while 8,127,841 variants remained for the downstream analyzes before annotation.

Functional annotation for protein-coding variants
We used Variant Effect Predictor133 (VEP v. 108) for variant annotation. Each variant was mapped to the most severe consequence across the canonical transcripts. We defined the loss-of-function variants using the Loss-of-Function Transcript Effect Estimator134 (LOFTEE) plugin, which further collapsed stop-gained, essential splice, and frameshift variants into high-confidence predicted loss-of-function variants (hcpLoF or pLoF hereafter) or low-confidence predicted loss-of-function variants (lcpLoF). Missense variants were prioritized by dbNSFP135,136 (v. 4.5a) plugin using five prediction algorithms16: SIFT22, PolyPhen2 HDIV23, PolyPhen2 HVAR23, LRT24, and MutationTaster25. Missense variants were defined as “int5 damaging missense variants” if predicted damaging or possible damaging by the intersection of all the five algorithms and “int1 damaging missense variants” if predicted damaging or possible damaging by any one of the five algorithms. Parallelly, we also used AlphaMissense26 to assign the missense variants to be “Alpha damaging” if predicted as “pathogenic”. Synonymous variants served as empirical null control to support the study-wise P-value threshold for burden tests and the method would be detailed below in the section for determining study-wise significance level. Predicted loss-of-function variants (including both hcpLoF and lcpLoF) and missense variants were included in our downstream association tests. They were referred to as non-synonymous variants or coding variants of interest (n = 2,143,707).

Exome-wide association tests
The covariate adjustment for brain MRI, DTI parameters, brain fMRI, and heart CMR traits was consistent with our previous GWAS1–3,10,118. In short, we adjusted a set of basic covariates for all the imaging traits including age (at imaging), age-squared, sex, age-sex interaction, age-squared-sex interaction, imaging site, and the top 40 genetic PCs (for the phases 4 to 5 sample replication analysis, we only adjusted for top 10 genetic PCs). For brain structural MRI traits, we additionally adjusted for total brain volume (for traits other than itself). For brain fMRI traits, we additionally adjusted for the effects of volumetric scaling, head motion, head motion-squared, brain position, and brain position-squared. For heart CMR traits, we additionally adjusted for the effects of standing height and weight. For abdominal MRI traits137, we additionally adjusted for the effects of standing height and body mass index.
We used REGENIE (v. 3.1.3)138 to conduct the exome-wide association tests. Common variants from genotyping array data (MAF > 0.01, genotype missing rate < 10%, Hardy–Weinberg equilibrium P > 1 × 10−15, LD pruning with r2 < 0.9) were included in the step 1 of REGENIE to capture genome-wide polygenic effects. Then, the predictors obtained from step 1 were used in step 2 for both variant-level tests and gene-based burden tests. For all the variant-level tests and gene-based burden tests, we followed the default setting of REGENIE that only variants/burden masks with MAC > 5 were tested. For variant-level tests, we tested all the coding variants of interest as defined above. For gene-based burden tests, we included three general types of burden masks: (i) hcpLoF variants only, (ii) damaging missense variants only, and (iii) the combinations hcpLoF, lcpLoF, and damaging missense variants. For damaging missense variants, we further defined three categories, that is “int5 damaging”, “int1 damaging”, and “Alpha damaging” missense variants as described above, which resulted in seven finer variant sets for each gene: (i) hcpLoF variants only, (ii) “Alpha damaging” missense variants only, (iii) “int5 damaging” missense variants only, (iv) “int1 damaging” missense variants only, (v) hcpLoF variants and “Alpha damaging” missense variants, (vi) hcpLoF variants and “int5 damaging” missense variants, and (vii) hcpLoF variants, lcpLoF variants, and “int1 missense” variants. Additionally, for each gene set, we considered four levels of variants based on the alternative allele frequency (REGENIE used alternative allele frequency for separation, but we would keep the notation MAF as they were mostly concordant in our analysis): singletons only, MAF ≤ 1 × 10−4, MAF ≤ 1 × 10−3, and MAF ≤ 1 × 10−2. Then for each given gene set, the score is calculated as the maximum of the genotype dosage across all the variants in this set (the “max rule” as the default implementation in REGENIE), and the score is tested for association. We noted that some sets may not be testable due to lack of qualified variants (for example, considering a gene that only had variants with MAF > 1 × 10−3) while some sets would produce repetitive results (for example, considering a gene that only had variants with MAF ≤ 1 × 10−4). To account for these issues, we adopted an empirical null based P-value threshold for multiple testing adjustment as described below.

Study-wise significance level for association tests
For variant-level associations, we applied the Bonferroni correction at 0.05 level, which yields a conservative P-value cutoff at 2.8 × 10−10 (= 0.05/178,280,016, that is 0.05 adjusted for all the variant-level tests for coding variants of interest, see Supplementary Data S24). The variant-level tests include all the coding variants of interest (n = 2,143,707) across all 596 phenotypes, but only those with sufficient MAC (MAC > 5) were tested and the MAC was calculated for each phenotype respectively by REGENIE, so a large number of tests of variants with MAC < 5 were ignored (Supplementary Note). The summary count and more details are available in Supplementary Data S24. Briefly, 2,143,707 coding variants of interest were included for variant-level association tests, and an average of 299,128 tests were conducted for each phenotype in joint sample and similarly an average of 239,825 tests per phenotype were conducted for the discovery sample (Supplementary Data S24 and Supplementary Note). We note that our primary Bonferroni-type P-value threshold is quite stringent as the phenotypes are not entirely independent, while we also provided the results for the variant-level association that passed the rare variant genome-wide significance threshold at 1 × 10−8 as suggestive evidence139,140. The P-value cutoff for phases 1 to 3 sample was 3.5 × 10−10 (= 0.05/142,935,657, Supplementary Data S24).
For gene-based burden tests, Bonferroni correction for all the tests conducted may be too strict and inappropriate here since distinct burden models and MAF cutoffs were highly correlated and may produce repetitive results. Moreover, the dependence structure made it unclear whether the Benjamini and Hochberg procedure could provide valid false discovery rate control. Alternatively, we used resampling-based methods to derive the empirical distributions, which suggested setting 1 × 10−9 as our conservative P-value threshold for burden tests. Specifically, we investigated the empirical null distributions of P-values from both permutation test60,61 and synonymous burden test61,141 separately. We conducted the permutation for the imaging traits once for every phenotype while preserving the genetic structure and the burden models60,61. At the tail of 243,691,511 P-values derived from the permutation-based null test, we only observed 5 results that had P-values smaller than 1 × 10−9 (Supplementary Data S25). Therefore, based on this permutation-derived threshold at 1 × 10−9, the expected false discoveries would be 5 out of 224 significant associations (3 out of 67 significant gene-trait pairs) across all burden models and imaging phenotypes.
Furthermore, we found that under this permutation-based P-value cutoff, only synonymous variants in NOSTRIN and HIGD1B were observed to be associated with the mean diffusivity of the cingulum (cingulate gyrus) tract in brain white matter and the volume of left rostral middle frontal region, respectively, at the tail of the distribution of 40,545,911 P-values from synonymous burden tests (Supplementary Data S26). These two genes were previously reported to be associated with white matter142,143 and other brain-related traits144,145, although the biological significance of the roles of synonymous variants remains unclear. Nevertheless, the very few significant associations further supported the validity of our choice of this P-value cutoff through a complementary perspective as synonymous variants generally would not contribute to the gene-trait associations61,141. We observed similar patterns of the tail distributions of P-values from permutation test and synonymous burden test from the phases 1 to 3 sample (Supplementary Data S27 and S28). Thus, the P-value cutoff for burden tests of this sample (used in our internal replication) was set to be 1 × 10−9. We also provided a more liberal P-value cutoff at 1 × 10−8 as suggestive evidence (approximately empirical false discovery proportion 0.1, Supplementary Data S25 and S26) for gene-based burden test.

Cross-reference with previous GWAS
For the comparison with GWAS signals, we leveraged the summary statistics derived from previous studies1–4,10 on the same study cohort. The P-value cutoff for GWAS signals was set to be the Bonferroni adjustment for all the traits (P < 5 × 10−8/trait numbers for each category) as described previously. Moreover, we focused on reported independent variants146 for GWAS signals and examined whether our exome-wide signals fell within the 1 Mb range of such independent GWAS signals to prioritize the effector genes among thousands of identified GWAS loci. The genetic build for previous GWAS was GRCH37, so we did liftover147 for our exome-wide results back from GRCH38 to GRCH37 and made comparisons based on the same genetic coordinates.

Imaging genetics and drug detection
We downloaded Therapeutic Target Database (TTD, last updated by January 10th, 2024)119 and evaluated all the entries documented in TTD that were claimed to be approved drug targets, in clinical trials, or supported by literature. We performed enrichment tests to see whether the identified genes in our gene-based burden tests were enriched in TTD drug target genes. We conducted Fisher’s exact tests for two classes of drug target genes in TTD. Specifically, we first included all drug targets that are either approved drug targets or in clinical trials (documented as “Successful” or “Clinical trial” in TTD), and then additionally included drug targets that were supported by literature (documented as “Literature-reported” in TTD).

Ultra-rare burden heritability and genetic correlation with health-related traits
We used the recently proposed BHR27 method to quantify the burden heritability and the burden genetic correlation for rare variants. BHR required one to divide the rare variants into multiple subgroups based on MAF cutoff and annotation outcome. Following the practical guide and previous examples, we investigated the genetic architecture of ultra-rare variants (MAF < 1 × 10−4) in our sample across 596 imaging traits. We used “univariate” mode and ran BHR to estimate the burden heritability for ultra-rare pLoF variants and “int1 damaging” missense variants. Synonymous variants were used as negative control. The corresponding P-values were manually calculated using one-sided test. Moreover, we downloaded 186 health-related traits or diseases from Genebass17 to estimate the burden genetic correlation with imaging traits in terms of the ultra-rare pLoF variants. Among the 596 imaging traits examined, those that showed significant heritability after multiple testing adjustment (P = 8.4 × 10−5 = 0.05/596, that is Bonferroni adjusted for 596 imaging traits) were included in burden correlation analysis. Then, we ran “bivariate” mode of BHR to calculate the burden genetic correlation for the selected imaging traits and downloaded phenotypes.

Rank-based inverse normal transformation and conditional analysis
To verify the robustness of the identified associations, rank-based inverse normal transformation (RINT) was applied by calling “--apply-rint” in both step 1 and step 2 of REGENIE for all identified associations in our study, while we followed the same testing procedure. The summary statistics produced after applying RINT were compared with the original ones. Conditional analysis was performed by specifying the nearby common variants as covariates to assess the independence of the identified signals. The conditional analysis is restricted to step 2 in REGENIE, while the genome-wide predictors in step 1 were reused.

Statistics and reproducibility
No statistical method was used to predetermine sample size. We made full use of the sample with both imaging phenotypes and exome sequencing data available. No data were excluded from the analyzes after quality control step as described previously. Specific details on study design and association analyzes have been provided in the previous sections. The experiments were not randomized, and the investigators were not blinded to allocation during analysis, because this is an observational association analysis of existing cohort data.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript by Yijun Fan et al. presents an exome-wide association study assessing rare coding variants associated with 596 multi-organ MRI traits in over 50,000 individuals from the UK Biobank. The study uses both variant-level and gene-based burden associations to explore how these variants contribute to the structure and function of human organs. While the analysis is extensive and well-documented, I have concerns regarding the reliability of the novel loci identified.\nMajor comments:\n1. The authors have chosen joint analysis and internal replication to validate their findings. Given that meta-analysis is a robust method for integrating results across multiple GWAS to identify novel loci, could the authors clarify why this approach was not considered?\n2. In Page 5, Line 18, the manuscript states that 107 rare variant-level associations were identified using a conservative P-value threshold. However, the derivation of these 107 associations from the reported 59 that passed the Bonferroni threshold is unclear. Furthermore, it is mentioned that 48 of these 59 associations had mutations with a Minor Allele Count (MAC) > 5. It seems there might be a discrepancy in the numbers reported or a potential misunderstanding in how these figures relate. Could the authors clarify how the total of 107 was reached and confirm whether all these associations also underwent internal replication and joint analysis?\n3. For a clearer understanding of the impact and reliability of the reported associations, please provide detailed statistical information (such as Beta, P-values) for the 107 rare variant-level associations across different phases (1-3, 4-5) and in the joint analysis.\n4. Similar to the variant-level associations, the details surrounding the 244 gene-based associations are not entirely clear. Including a detailed workflow, possibly in the supplementary materials, might help clarify the process by which these associations were identified and analyzed.\n5. On Page 6, the comparison between effect directions and P-values is noted, but there is no mention of comparing effect sizes. Did the authors evaluate the correlation of effect sizes across studies? Including this analysis would provide deeper insights into the consistency of the associations identified.\n6. Considering the manuscript addresses both variant-level and gene-based associations, it would be beneficial to understand if there are overlapping signals between these two levels of analysis. Such overlaps could signify particularly robust associations.\nMinor comments:\n1. The manuscript mentions that 8,127,841 variants remained after quality control for downstream analyses (Page 38, Line 5). However, there seems to be a discrepancy, as Page 5, Line 19 refers to a Bonferroni adjustment based on 178,280,016 variant-trait tests. Can the authors clarify these numbers? Specifically, how do the 8,127,841 variants relate to the much larger number used for the Bonferroni correction?\n2. Please provide clearly variants number and/or gene number in each stage.\n3. There are numerous NA entries in the BETAorg (95% CI) and Unit columns of Table 1. Could the authors explain the reason for these missing values?", "reviewer_2": "In this work Fan et al. investigate the link between rare genetic variation measured using exome sequence and MRI-derived phenotypes for 3 ‘organs’ (brain, heart, abdomen). They filtered down from millions of rare variants by restricting the analysis to damaging variants as predicted by various tools (classic ones as well as AlphaMissense). The analysis is restricted to the UK Biobank. Since there were no external validation datasets (as is common in genetics) the authors split the UK Biobank into a discovery and a validation cohort. Imaging releases 1 to 3 were used for discovery and releases 4 to 5 were used for validation. However, the nature of rare variants means that some identified SNPs were not present in the validation dataset (with about 10k people in the validation you can expect two carriers of the variant at MAF 0.0001). In such cases the authors considered a confirmatory analysis if the direction was maintained and the p-value improved.\nIn addition to the variant-level analysis the authors conducted a series of gene-burden tests using various groupings of the rare variants based on the algorithms that rated them as damaging. These tests uncovered a few additional gene-trait associations that were not picked up by single variant tests.\nThe authors then discuss some of the gene-trait associations grouped by organ and by analysis technique. Focusing on three different organs renders these discussions quite long and probably only interesting to researchers who are working on the specific organ. While others might simply skip these parts. The synergy between organs appears to be limited. There is a link between cardiac associations and fMRI-traits, which is sensible and has been commented on before.\nThere are a few novel findings, some of them support gene-trait associations that have been uncovered with GWAS in the same cohort. Overall, this is a large dataset and it is somewhat surprising that not more associations were revealed by this analysis.\nOverall, this manuscript is well written and easy to follow. The authors conducted a lot of analyses (e.g., it is nice to evaluate the effects of different SNP annotations and the drug target database enrichment). In my opinion the statistics are appropriate although there remain some questions about the detailed methodology.\nMajor\n1. It is unclear how the exome sequencing data was pre-processed. Presumably the authors did not conduct the analysis from scratch and relied on the VCF files provided by UK Biobank. However, it would still be great to provide a summary of the processing details from raw sequences to called variants.\n2. It is not very clear how exactly the authors arrive at the multiple testing threshold. Even the Supplementary table 21 does not quite reveal the underlying principle sufficiently. It does not even become clear whether the threshold is adjusted for all the different imaging traits. For instance for the volume trait there are 101 phenotypes and 31096939\nVariants tested. Does that mean that there are 307890.5 variants tested per phenotype? If that is the case, then the authors should also emphasize that their adjustment is quite stringent because the phenotypes are not independent. Please, can the authors clarify this aspect of the analysis.\n3. All the explored traits are of a quantitative nature. Especially with rare variants where there may only be a hand full of carriers, outliers in these traits can drive the associations and lead to inflated type 1 error (see https://pmc.ncbi.nlm.nih.gov/articles/PMC4970685/). Have the authors considered any mitigations for these cases? If not, it would be valuable to verify the identified associations after applying inverse normal transformation (INT) of the quantitative data.\n4. For the validation analysis. I don’t think that the argument that the direction was maintained in the larger cohort and that the p-values improve is a sufficiently strong one (e.g., P6 Line 11). As outlined in the previous comment, outliers may be driving the association and they may still be the driving force in the larger cohort. To make this a valid argument, the authors should show that the p-value improved beyond what would be expected from the increased sample size alone.\n5. The authors evaluate the variant hits in terms of existing GWAS hits. In particular for DTI traits they uncover a rare variant in the VCAN gene. In an earlier GWAS (with comparable sample size) an intronic SNP in the gene was strongly associated with DTI traits (rs67827860 (P = 4.06 × 10−37; Elliot et al. ). What is the interpretation here given the common variant shows a stronger association? Wouldn’t it be more likely that the rare variants are only in LD with the variation causing driving the association?\n6. When reading the manuscript it does not become quite clear which gene based burden test the authors are using. There are multiple different ones implemented in regenie. Please, can the authors provide detailed information (also in the results section on page 9 around line 20, not just the methods section).\n7. I recommend adding to the discussion also the authors’ perspective on the value of rare-variant analyses. Compared to GWAS with the same sample size, much fewer loci-trait associations have been identified. Does this mean that perhaps one should focus cost and effort to other types of variation?\nMinor\n1. P4 Line 28: “whole body and shared genetic regulation across organs.” I think the reasoning is a bit inconsistent in that the rare variants (which are limited to damaging, protein altering variations) are viewed as regulators for processes across organs. Often the focus on rare variation is to illuminate a more direct “mutation”-effect association rather than uncovering a regulatory process. I recommend to reword this.\n2. When reporting SNP effect sizes and association strengths it would be useful to have the MAF or MAC in the main text\n3. P8 Line 11-13: I do not think the reasoning/causality here is one to support the association (i.e., the gene is associated with Schizophrenia (common) and that Schizophrenia is associated with enlarged ventricles etc) to support the link between the gene and brain volumes… Suggest rephrasing this.\n4. P38 Line 21 “Synonymous variants served as empirical null control to support the study-wise P-value threshold for burden tests”. The statement is a bit vague and unclear. The method behind this is only explained further down in the methods section. I recommend to add a forward point to the relevant section to clarify that this will be explained below.\n5. Regarding the Hardy-Weinberg filtering: I recommend to omit this step. Especially with rare variants checking the HW may not be useful. Previous work that examined the HW QC step suggested to retain variants and only examine the HW after associations have been found: https://www.medrxiv.org/content/10.1101/2024.02.07.24301951v2\n6. It appears as if genetics that influences blood flow affects the fMRI phenotypes. This appears to be a bias in the analysis and interpretation of fMRI signals. Perhaps the authors can comment on the need to incorporate cardiac features/genetics into the analysis of fMRI data?"}, "author_response": "Response to the Reviewer 1:\nThank you so much for your careful review and constructive suggestions! Based on your\nsuggestions, we have made substantial improvements and updated our manuscript\naccordingly, with all changes highlighted in yellow. Here, we provide our responses to each\nof your points. For your convenience, we first state your comments in italic and then\nprovide our responses.\nReviewer #1 (Remarks to the Author):\nThe manuscript by Yijun Fan et al. presents an exome-wide association study assessing\nrare coding variants associated with 596 multi-organ MRI traits in over 50,000 individuals\nfrom the UK Biobank. The study uses both variant-level and gene-based burden\nassociations to explore how these variants contribute to the structure and function of\nhuman organs. While the analysis is extensive and well-documented, I have concerns\nregarding the reliability of the novel loci identified.\nResponse: Thank you very much for reviewing our paper and your constructive comments!\nWe appreciate that you have highlighted the strengths of our study, including the extensive\nanalysis across muti-organ MRI traits, both variant-level and gene-based analyses, and the\nwell-documented steps. Following your comments, we conducted additional analyses and\nprovided clarifications to improve the reliability of the newly identified loci. We have tried\nour best to improve this work by closely following your specific questions and suggestions\nbelow.\nMajor comments:\n1. The authors have chosen joint analysis and internal replication to validate their findings.\nGiven that meta-analysis is a robust method for integrating results across multiple GWAS\nto identify novel loci, could the authors clarify why this approach was not considered?\nResponse: Thank you so much for your insightful question! We agree that meta-analysis is a\nrobust and convenient method for integrating results across multiple studies with only\nsummary statistics is available, especially for common variant analysis in GWAS. In this UK\nBiobank-based WES study, we had full access to individual-level data from 54,365 European\nparticipants, enabling us to directly perform a joint analysis on the pooled sample to\nestimate summary statistics. Therefore, we first reported the association analysis based on\nthe joint sample, which maximizes statistical power by leveraging all available data. To make\nthe best use of this data resource for evaluating reproducibility, we further split the sample\ninto independent internal discovery and validation subsets for an internal replication\nanalysis. This line of reasoning led us to favor a joint analysis over a meta-analysis in our\nstudy design.\nsome variants (for example, with frequencies as low as MAF < 1 × 10⁻⁴) poses unique\nchallenges, making meta-analysis particularly challenging. For example, in gene-based\nburden tests, rare variants are aggregated within a gene based on functional annotation and\nMAF thresholds. However, given a specific gene region and a predefined burden mask, the\nsets of variants identified in the discovery and validation samples may differ, making the\nresults not fully comparable and complicating gene-level meta-analysis of rare variants.\nBased on these considerations, using a joint analysis for the primary results, alongside an\ninternal replication strategy via sample splitting, was motivated by several factors. First,\ngiven the low MAF of rare mutations, joint analysis enhances the stability of summary\nstatistics and may increase statistical power. Second, this design allows us to make full use\nof the available data to evaluate the robustness of associations identified through both\nvariant-level and gene-level tests in a conventional manner for validation analysis.\nFurthermore, in the absence of non-UK Biobank WES datasets of comparable size suitable\nfor replicating associations with complex imaging phenotypes, our approach also provides\nbroader insights into the replicability of such phenotypes, for example, by quantifying the\nproportion of signals that can be confirmed given the available sample size.\nAs noted by Reviewer 2, the lack of external validation datasets is a common limitation in\ncurrent WES analyses. When additional external data sources become available, it will be\nhighly valuable to conduct similar rare variant association studies for imaging traits, as done\nin our UK Biobank study. Meta-analysis using summary statistics from multiple studies can\nthen be performed to integrate and better validate the findings. Following your comments,\nwe have clarified the considerations that we adopted this internal discovery-replication\ndesign within the UK Biobank study and mentioned the benefits of meta-analysis in the\nfuture studies to integrate results from multiple studies. In the revised manuscript, we have\nincluded these details in main text, discussion, and Supplementary Note:\n“Specifically, all the primary findings and the corresponding interpretations are based on the\njoint sample for higher statistical power and more stable summary statistics. Then the joint\nsample was split into the internal discovery sample and validation sample to verify the\nrobustness of the findings within the UKB study (Fig. 1, Methods, and Supplementary\nNote).” (Page 5, Lines 2-6)\n“Consequently, this led us to use joint analysis with a stringent P-value cutoff, and an\ninternal validation procedure to assess the robustness of our findings” (Page 25, Lines 20-22)\n“When future data resources become available, it will be valuable to replicate our findings\nand extend the scope of WES studies across multiple cohorts. Meta-analysis may be a robust\nand efficient approach for identifying and replicating associations using only summary\nstatistics.” (Page 25, Lines 28-31)\n2. In Page 5, Line 18, the manuscript states that 107 rare variant-level associations were\nidentified using a conservative P-value threshold. However, the derivation of these 107\nassociations from the reported 59 that passed the Bonferroni threshold is unclear.\nFurthermore, it is mentioned that 48 of these 59 associations had mutations with a Minor\nAllele Count (MAC) > 5. It seems there might be a discrepancy in the numbers reported or a\npotential misunderstanding in how these figures relate. Could the authors clarify how the\ntotal of 107 was reached and confirm whether all these associations also underwent\ninternal replication and joint analysis?\nResponse: Thank you so much for pointing this out! We appreciate your careful review on\nthe details of our association results. The “107 rare variant-level associations” were\nidentified based on the joint sample (phase 1 to phase 5) with a stringent P-value cutoff at\n2.8 × 10-10. The “59 associations” were identified using the internal discovery sample (phase\n1 to phase 3, a subset of the joint sample) specifically for the internal replication analysis.\nThat is, all 107 associations were identified through the joint analysis, leveraging the full\nsample size to maximize discovery power. However, for the internal replication analysis, we\nsplit the data, and only 59 associations derived from a subset of the sample were evaluated.\nSpecifically, we first conducted the joint analysis to report the 107 significant associations\nand then split the joint sample into the internal discovery sample and validation sample as\nmentioned in our manuscript and the response to your Comment 1. The P-value cutoff for\nthe discovery sample was achieved by the same procedure as the joint sample, resulting in\n59 significant variant-level associations. For these 59 associations, we found that 48 of them\nhave MAC > 5 in the validation sample (phase 4 to phase 5) and were suitable for\nassociation testing. During this revision, we have provided additional information, such as\nthe MAC for all identified rare variants in the validation sample, which may be helpful for\nclarification.\nFollowing your suggestions, we have carefully reviewed the results and added additional\ncolumns to Table S4 for the variant-level tests in the validation sample, including MAC\nvalues and indicators for the variants excluded from the validation tests as well as for\nwhether the variants passed the internal replication steps. In the revised manuscript, we\nhave also clarified procedures for the variant-level tests and explicitly stated the number of\ndiscoveries and their relationships. These clarifications are included in the main text,\nSupplementary Note, and the updated Table S4 to avoid any potential misunderstandings:\n“To assess the reliability of our findings derived from the joint sample and make full use of\navailable data, we adopted an internal replication procedure30. (Fig. 1, Methods, and\nSupplementary Note). Specifically, we split the joint sample to the internal discovery dataset\nand validation dataset after joint analysis for internal replication. The discovery dataset\nincluded individuals from phases 1 to 3 of MRI release (average n = 30,739) and is a subset of\nthe joint sample, while the validation dataset is independent of it by design and contained\nindividuals from phases 4 to 5 of the MRI release (Method, average n = 8,989). For the\npurpose of replication analysis, summary statistics from the discovery sample will be\nseparately compared with results from the validation sample and the joint sample (Method\nand Supplementary Note). When we restricted our study to the discovery sample, 59\nassociations passed the Bonferroni threshold of P < 3.5 × 10-10 ... and were used for internal\nreplication.” (Page 6, Lines 6-18)\n“...48 had mutations with MAC > 5 suitable for association tests (Table S4) in the\nindependent validation sample…” (Page 6, Lines 18-20)\n3\n“Consequently, this led us to use joint analysis with a stringent P-value cutoff, and an\ninternal validation procedure to assess the robustness of our findings” (Page 25, Lines 20-22)\n3. For a clearer understanding of the impact and reliability of the reported associations,\nplease provide detailed statistical information (such as Beta, P-values) for the 107 rare\nvariant-level associations across different phases (1-3, 4-5) and in the joint analysis.\nResponse: Thank you so much for your suggestions! Following your comment (including\nprevious ones), we have provided a clearer summary of our study design and the details of\nassociation tests conducted across different phases.\nSpecifically, as explained in our main text and the responses to your previous comments, the\n107 associations identified in the joint sample are our primary findings and all the summary\nstatistics can be found in Table S2. They include both the associations restricted to the\ndiscovery sample and associations unique to the joint sample.\nWhen restricted to the discovery sample, 59 associations were identified for replication\nanalysis. The summary statistics (including Beta and P-values) for these 59 associations from\nthe internal discovery sample, those available in validation sample, and those available in\nthe joint sample have already been provided in Table S4. Following your suggestions, we\nadded several additional columns including MAC, number of individuals, indicator of the\nresults for clearer understanding of the impact and reliability of the associations and\nstatistical information. We also explicitly pointed out these details in the main text:\n“More details of the results from the internal replication analysis and summary statistics of\nvariant-level associations across different stages can be found in Table S4.” (Page 7, Lines 4-\n6)\n4. Similar to the variant-level associations, the details surrounding the 244 gene-based\nassociations are not entirely clear. Including a detailed workflow, possibly in the\nsupplementary materials, might help clarify the process by which these associations were\nidentified and analyzed.\nResponse: Thank you so much for your suggestions! Similar to variant-level tests, these 224\ngene-based associations were based on the joint sample (phase 1 to phase 5) and they are\nour primary findings. For internal replication, we only conducted the confirmation analysis\nbetween the discovery sample (phase 1 to phase 3) and the joint sample as gene-based\ntesting based on validation sample is not comparable with the discovery sample or joint\nsample due to the limited sample size mentioned in the main text and in the response to\nyour Comment 1. Following your comments, we have provided similar details and statistical\ninformation in Table S9. We also clarified the process in the main text and in the\nSupplementary Note:\n“Therefore, we only made use of the internal discovery sample and examined..., similar to\ninternal-replication analysis for variant-level tests (Methods and Supplementary Note)”\n(Page 11, Lines 18-21)\n4\n“The detailed summary statistics and results of internal replication analysis for gene-based\nassociations have been provided in Table S9.” (Page 11, Lines 25-27)\n5. On Page 6, the comparison between effect directions and P-values is noted, but there is\nno mention of comparing effect sizes. Did the authors evaluate the correlation of effect\nsizes across studies? Including this analysis would provide deeper insights into the\nconsistency of the associations identified.\nResponse: Thank you so much for your insightful suggestion! Following your suggestion, we\nhave added the correlation metric in the validation stage for variant-level associations to\nverify the consistency of the associations identified. We agree that the correlation of effect\nsizes between discovery sample and validation sample is interpretable as they are from\nindependent samples.\nSpecifically, for all the variant-level associations identified in the phase 1 to phase 3\ndiscovery sample under the primary p-value cutoff at 3.5 × 10-10 as well as a more\npermissive cutoff at 1 × 10-8, we retrieved the effect sizes of 48 associations available in the\nvalidation sample for the correlation test (with MAC > 5). Then, we conducted Pearson\ncorrelation test for those associations with summary statistics available in both stages. The\ncorrelation showed that the associations across the two stages are consistent. Specifically,\nat the primary threshold 3.5 × 10-10, there are 48 (out of 59) associations with both stages of\nsummary statistics available, and we have observed great correlation of effect sizes (r =\ng\n0.935, P = 2.8 × 10-22). At a more lenient cutoff 1 × 10-8, there were 83 (out of 154)\nassociations available for both stages and showed moderate correlation (r = 0.653, P = 2.2 ×\ng\n10-11). Overall, the associations identified showed consistency across the discovery and\nvalidation sample, while it also verified that our primary threshold gave a credible set of\nassociations. In this revision, we have added the results of correlation test in the main text:\n“We conducted the Pearson correlation test for the effect sizes of the identified associations\nin the discovery sample and the validation sample. Among these 48 associations with\nsummary statistics available in both stages, we observed a substantial correlation (r = 0.935,\ng\nP = 2.8 × 10-22), which shows a high level of agreement between the associations identified in\nthe discovery sample and validation sample. When we relaxed the P-value cutoff to 1 × 10-8,\nthere was still a moderate to high correlation (r = 0.653, P = 2.2 × 10-11).” (Page 6, Lines 23-\ng\n29)\n6. Considering the manuscript addresses both variant-level and gene-based associations, it\nwould be beneficial to understand if there are overlapping signals between these two\nlevels of analysis. Such overlaps could signify particularly robust associations.\nResponse: Thank you so much for your insightful comments! We agree that evaluating\noverlapping signals between variant-level and gene-based tests is interesting, as these tests\nexamine associations from different perspectives/assumptions and can yield\ncomplementary insights. Therefore, the overlapping signals indicate that both a single rare\nvariant in this region and the joint effects of multiple variants across the region are both\nstrong enough to drive a significant association. In our study, we found that ANO1 with\n5\nmax/min aorta area and TNFRSF13B with spleen volume are the associations shared in both\nvariant-level and gene-based signals. We discussed these associations respectively in our\noriginal manuscript. In the revised manuscript, we highlighted these associations and\nemphasized the fact that they are overlapping signals in the main text following your\nsuggestion:\n“Three gene-trait pairs overlapped across variant-level and gene-level tests (two\nassociations between ANO1 and max/min aorta area, and an association between\nTNFRSF13B and spleen volume), indicating that both a single rare variant in this region and\nthe joint effects of variants across the region are strong enough to drive an association.”\n(Page 11, Lines 5-10)\nMeanwhile, we would like to note that the WES literature generally expects variant-level\nand gene-level tests to capture distinct signals (such as Wang et al., Nature PMID: 34375979;\nBackman et al., Nature PMID: 34662886; Karczewski et al., Cell Genom. PMID: 36778668).\nFor example, gene-based burden test is more powerful for identifying signals driven by pLoF\nvariants as these variants are rare individually but have consistent effect directions\n(Karczewski et al., Cell Genom. PMID: 36778668), while the pattern for missense variants in\nvariant-level tests and gene-level tests is less obvious (Backman et al., Nature PMID:\n34662886; Karczewski et al., Cell Genom. PMID: 36778668). These observations are\nconsistent with our findings. In the revised manuscript, we have also discussed this point in\nthe discussion:\n“Analyses at the variant level and gene-based burden level may yield complementary\ninsights into gene associations. For example, gene-based burden testing is more powerful for\nidentifying signals driven by pLoF variants as these variants are rare individually but have\nconsistent effect directions17, while the pattern for missense variants signals in variant-level\ntests and gene-level tests is less obvious16,17,61. In our study, we found that pLoF signals were\nexclusively detected by gene-based tests, whereas three gene–trait pairs were observed\nacross both variant- and gene-level analyses. This pattern aligns with previous findings and\nunderscores the importance of incorporating both approaches. Moreover, although gene-\nbased burden tests can enhance signal detection16,17, annotating the deleterious variants is\nstill an active research area.” (Page 23, Lines 6-15)\nMinor comments:\n1. The manuscript mentions that 8,127,841 variants remained after quality control for\ndownstream analyses (Page 38, Line 5). However, there seems to be a discrepancy, as\nPage 5, Line 19 refers to a Bonferroni adjustment based on 178,280,016 variant-trait tests.\nCan the authors clarify these numbers? Specifically, how do the 8,127,841 variants relate\nto the much larger number used for the Bonferroni correction?\nResponse: Thank you so much for your careful review! The 8,127,841 was the number of\nvariants that were filtered after the basic quality control steps before functional annotation\n(Page 41, Line 8). After annotation, only coding variants of interest were included in the\nassociation tests (specifically predicted loss-of-function variants and missense variants,\n6\ndefined in the annotation section in Methods). This step would further filter out a subset of\n2,143,707 variants of interest for association tests (Page 41, Line 28).\nThe key reason we use a much larger number for the Bonferroni correction is that we\nadditionally adjust for all imaging phenotypes. Specifically, the Bonferroni correction was\nadjusted for the number of all the variant-level tests conducted based on these coding\nvariants of interest across all the 596 phenotypes. This would give us a quite stringent\nthreshold to determine the significant variants. Additionally, REGENIE computed the MAC\nfor each phenotype, and following its suggested and default settings, only variants or gene\nmasks with MAC > 5 were used for association testing, as described in our original\nmanuscript. Thus, a great number of rare variants would be omitted due to insufficient MAC\nand the number of tests omitted was reported in the log files for each phenotype category.\nFor example, in the variant-level tests of 101 brain volume phenotypes, 690,562,178 tests in\ntotal were ignored due to low MAC as documented in the log files, which means 6,837,249\nvariants were ignored on average for each brain volume phenotype. Meanwhile, there were\n130,349,763 results (including all the variants after quality control) observed in the variant-\nlevel tests for these 101 brain volume phenotypes and we focused on 31,096,939 results\namong them for coding variants of interest based on the annotation results. This means\n1,290,592 variants (including both coding variants of interest and other variants) are present\nfor each brain volume phenotype and among them 307890.5 variants were coding variants\nof interest on average. Finally, the number of all categories of variants present and the\nnumber of variants ignored on average for each brain volume phenotype would add up to\nthe total number of variants after quality control, that is 8,127,841 (=1,290,592+6,837,249).\nBut we would like to emphasize that the number 31,096,939 would be used for multiple\ntesting adjustment for the part of 101 brain volume phenotypes (corresponding to 307890.5\non average for each of the brain volume phenotypes) as they represent the number of\ncoding variants.\nBased on these facts, we manually retrieved the number of coding variants of interest\nincluded in the variant-level tests for each of the 596 phenotypes by customized code for\nprecise multiple testing adjustment. Finally, the number of tests for each organ category of\nphenotypes (such as 31,096,939 for brain volume phenotypes) was summed up and this\nwould yield 178,280,016 variant-trait tests in total as documented in Table S24. We also\nincluded the total number of variant-level tests for each phenotype category in Table S24,\nwhich is helpful to track the number of variants.\nFollowing your question, we have carefully checked our procedure and put this example\ninto Supplementary Note for illustration. We have also clarified it in the revised manuscript\nand Table S24:\n“The variant-level tests include all the coding variants of interest (n = 2,143,707) across all\n596 phenotypes, but only those with sufficient MAC (MAC > 5) were tested and the MAC was\ncalculated for each phenotype respectively by REGENIE, so a large number of tests of\nvariants with MAC < 5 were ignored (Supplementary Note). The summary count and more\ndetails are available in Table S24.” (Page 43, Lines 11-16)\n7\n2. Please provide clearly variants number and/or gene number in each stage.\nResponse: Thank you so much for your suggestion! Following your suggestion, we added the\nnumber of variants and genes for association tests in the joint sample. The validation\nanalysis was only performed for the identified associations from the discovery sample for\nvariant-level tests and was not conducted for gene-based burden tests. We also provided\nthe average number of variants tested in the internal discovery sample and joint sample, as\nthis would give a clearer intuition of the number of tests affected by the sample size. In the\nrevised manuscript, we have made such changes:\n“The variant-level tests include all the coding variants of interest (n = 2,143,707) across all\n596 phenotypes, but only those with sufficient MAC (MAC > 5) were tested and the MAC was\ncalculated for each phenotype respectively by REGENIE, so a large number of tests of\nvariants with MAC < 5 were ignored (Supplementary Note). The summary count and more\ndetails are available in Table S24. Briefly, 2,143,707 coding variants of interest were included\nfor variant-level association tests, and an average of 299,128 tests were conducted for each\nphenotype in joint sample and similarly an average of 239,825 tests per phenotype were\nconducted for the discovery sample (Table S24 and Supplementary Note).” (Page 43, Lines\n11-20)\n“A total number of 17,474 genes are available for burden tests in the joint sample and\ndiscovery sample. Using the same joint sample as in variant-level tests...” (Page 10, Lines 28-\n31)\n3. There are numerous NA entries in the BETAorg (95% CI) and Unit columns of Table 1.\nCould the authors explain the reason for these missing values?\nResponse: Thank you very much for pointing this out! These NA entries in the BETAorg (95%\nCI) and Unit columns of Table 1 correspond to DTI traits and fMRI traits. For DTI and fMRI\nphenotypes, the physical unit is typically missing or not interpretable. Instead, these kinds of\nphenotypes may stand for summary signals extracted from the measurements to indicate\nsome specific types of information. By contrary, structural traits are more interpretable in\nterms of physical unit (for example, volume), so the original units would offer more intuitive\ninformation regarding the effect sizes for these traits. For clarification, we have mentioned\nit in the caption of Table 1 in the revised manuscript:\n“Additional original units were provided for structural traits including regional brain volumes,\nabdominal MRI traits, and CMR traits for intuitive interpretations of the effect sizes.” (The\ncaption of Table 1, Page 50, Lines 17-19)\nThank you so much for your careful review and constructive suggestions! Based on your\nsuggestions, we have made substantial improvements and updated our manuscript\naccordingly, with all changes highlighted in yellow. Here, we provide our responses to each\nof your points. For your convenience, we first state your comments in italic and then\nprovide our responses.\nReviewer #2 (Remarks to the Author):\nIn this work Fan et al. investigate the link between rare genetic variation measured using\nexome sequence and MRI-derived phenotypes for 3 ‘organs’ (brain, heart, abdomen). They\nfiltered down from millions of rare variants by restricting the analysis to damaging\nvariants as predicted by various tools (classic ones as well as AlphaMissense). The analysis\nis restricted to the UK Biobank. Since there were no external validation datasets (as is\ncommon in genetics) the authors split the UK Biobank into a discovery and a validation\ncohort. Imaging releases 1 to 3 were used for discovery and releases 4 to 5 were used for\nvalidation. However, the nature of rare variants means that some identified SNPs were\nnot present in the validation dataset (with about 10k people in the validation you can\nexpect two carriers of the variant at MAF 0.0001). In such cases the authors considered a\nconfirmatory analysis if the direction was maintained and the p-value improved.\nIn addition to the variant-level analysis the authors conducted a series of gene-burden\ntests using various groupings of the rare variants based on the algorithms that rated them\nas damaging. These tests uncovered a few additional gene-trait associations that were not\npicked up by single variant tests.\nThe authors then discuss some of the gene-trait associations grouped by organ and by\nanalysis technique. Focusing on three different organs renders these discussions quite long\nand probably only interesting to researchers who are working on the specific organ. While\nothers might simply skip these parts. The synergy between organs appears to be limited.\nThere is a link between cardiac associations and fMRI-traits, which is sensible and has\nbeen commented on before.\nThere are a few novel findings, some of them support gene-trait associations that have\nbeen uncovered with GWAS in the same cohort. Overall, this is a large dataset and it is\nsomewhat surprising that not more associations were revealed by this analysis.\nOverall, this manuscript is well written and easy to follow. The authors conducted a lot of\nanalyses (e.g., it is nice to evaluate the effects of different SNP annotations and the drug\ntarget database enrichment). In my opinion the statistics are appropriate although there\nremain some questions about the detailed methodology.\nResponse: Thank you very much for reviewing our paper and for your insightful comments!\nWe completely agree that including extensive analyses such as incorporating multiple\nannotation strategies, conducting both variant-level and gene-level tests, and leveraging\n9\nother data sources such as drug target database could be interesting. We also appreciate\nthat you have highlighted our study design inspired by the limitations of rare variant\nassociation studies. For example, the nature of rare mutation requires relatively large\nsample to investigate, and the lack of external dataset poses a challenge on validation\nanalysis. We have tried our best to improve this work by carefully following your\nconstructive comments and questions below.\nMajor\n1. It is unclear how the exome sequencing data was pre-processed. Presumably the\nauthors did not conduct the analysis from scratch and relied on the VCF files provided by\nUK Biobank. However, it would still be great to provide a summary of the processing\ndetails from raw sequences to called variants.\nResponse: Thank you so much for your suggestion! We have provided the relevant\ninformation and citations (Hout et al. Nature, PMID: 33087929, Szustakowski et al. Nat\nGenet. PMID: 34183854) about the pre-processing steps of exome sequencing data from UK\nBiobank in the revised manuscript:\n“Whole-exome sequencing (WES) data were obtained from the UK Biobank and the detailed\nsequencing protocols can be found elsewhere18. By following a unified OQFE protocol121, raw\nsequencing reads were aligned to the GRCh38 reference genome, and small variants were\ncalled on a per-sample basis using DeepVariant. Per-sample gVCFs were then joint-\ngenotyped using GLnexus to produce multi-sample project-level VCFs (pVCF) suitable for\nassociation analysis. We made use of the OQFE WES files in BGEN format provided by UK\nBiobank after these centrally pre-processing steps.” (Page 39, Lines 16-22)\n2. It is not very clear how exactly the authors arrive at the multiple testing threshold. Even\nthe Supplementary table 21 does not quite reveal the underlying principle sufficiently. It\ndoes not even become clear whether the threshold is adjusted for all the different imaging\ntraits. For instance for the volume trait there are 101 phenotypes and 31096939\nVariants tested. Does that mean that there are 307890.5 variants tested per phenotype? If\nthat is the case, then the authors should also emphasize that their adjustment is quite\nstringent because the phenotypes are not independent. Please, can the authors clarify this\naspect of the analysis.\nResponse: Thank you so much for your careful review! The Bonferroni correction was\nadjusted for all the variant-level tests conducted based on these coding variants of interest\nand across all the 596 phenotypes. Notably, REGENIE suggested only including variants/gene\nmasks with MAC > 5 for association tests and set it in default for stable results. The MAC\nwas calculated for each phenotype respectively by REGENIE. Therefore, a large number of\nrare variants were omitted due to low MAC and the exact number of tests omitted was\nreported in the log files for each phenotype category. In addition, for implementation, we\nran REGENIE for all the variants in the original bfile but only included coding variants of\ninterest for analysis.\ntotal were ignored due to low MAC as documented in the log file, which means 6,837,249\nvariants were ignored on average for each brain volume phenotype. Meanwhile, there were\n130,349,763 results in total observed in the variant-level tests for these 101 brain volume\nphenotypes and we focused on 31,096,939 results among them for coding variants of\ninterest. That is, 1,290,592 variants (including both coding variants of interest and other\nvariants) were present for each brain volume phenotype and among them 307890.5\nvariants were coding variants of interest on average. Our multiple testing adjustment would\nuse the number 31,096,939 here for these 101 brain volume phenotypes. The number of\ncoding variants of interest included in the variant-level tests for each of the 596 phenotypes\nwas retrieved by customized code for precise adjustment. Finally, the number of tests for all\nphenotypes (for example, 31,096,939 for 101 brain volume phenotypes) was summed up\nand this would yield 178,280,016 variant-trait tests in total as documented in Table S24. We\nalso included the total number of variant-level tests for each phenotype category Table S24,\nwhich is helpful to track the number of variants.\nWe completely agree that this threshold is quite stringent, and we would like to use this\nstringent threshold to ensure the reliability of our primary findings, given that we lack large-\nscale independent dataset suitable for external replication analysis. For example, we used\nthe term ‘conservative’ to describe this P-value threshold in the manuscript (Page5, Line 24).\nFollowing your suggestion, we have clarified these aspects of analysis in the revised\nmanuscript:\n“The variant-level tests include all the coding variants of interest (n = 2,143,707) across all\n596 phenotypes, but only those with sufficient MAC (MAC > 5) were tested and the MAC was\ncalculated for each phenotype respectively by REGENIE, so a large number of tests of\nvariants with MAC < 5 were ignored (Supplementary Note). The summary count and more\ndetails are available in Table S24.” (Page 43, Lines 11-16)\n“Consequently, this led us to use joint analysis with a stringent P-value cutoff, and an\ninternal validation procedure to assess the robustness of our findings” (Page 25, Lines 20-22)\n“We note that our primary Bonferroni-type P-value threshold is quite stringent as the\nphenotypes are not entirely independent, while we also...” (Page 43, Lines 20-21)\n3. All the explored traits are of a quantitative nature. Especially with rare variants where\nthere may only be a hand full of carriers, outliers in these traits can drive the associations\nand lead to inflated type 1 error (see\nhttps://pmc.ncbi.nlm.nih.gov/articles/PMC4970685/). Have the authors considered any\nmitigations for these cases? If not, it would be valuable to verify the identified\nassociations after applying inverse normal transformation (INT) of the quantitative data.\nResponse: Thank you so much for your insightful feedback and suggestions!\nIn terms of outliers, all the phenotype data present in our study were pre-processed to have\nthose entries with values greater than five times the median absolute deviation from the\nmedian value removed. This step was part of our pipeline for generating phenotype data as\n11\nin our previous work. In the revised manuscript, we have included this information in the\nMethod and cited our previous work accordingly following your question:\n“For each phenotype, values greater than five times the median absolute deviation from the\nmedian value were removed as in our previous work1-4,10.” (Page 40, Lines 4-5)\nFollowing your suggestions, we further verified our identified rare-variant associations by\napplying rank-based inverse normal transformation (RINT). Specifically, we applied RINT by\ncalling “--apply-rint” in both step 1 and step 2 in REGENIE as suggested in their guidance\n(https://rgcgithub.github.io/regenie/options/). We obtained all the summary statistics after\napplying RINT for the previously identified associations that passed our primary threshold in\nboth variant-level tests and gene-level tests following the same pipeline. The new results\nlargely supported our original findings. Specifically, all associations except one (LPAR3 with\nmean MD of uncinate fasciculus) still showed significance and passed the nominal threshold\nadjusted for the number of signals undergoing RINT analysis at 4.7 × 10-4 (= 0.05/107) for\nvariant-level tests and 2.2 × 10-4 (= 0.05/224) for gene-based burden tests, respectively.\nUsing our original stringent threshold, 75.7% (81/107) variant-level associations still\nsurvived our conservative threshold at 2.8 × 10-10, and the proportion became 85.0%\n(91/107) under the suggestive threshold at 1 × 10-9. The proportion of gene-based\nassociations passing the primary threshold at 1 × 10-9 is 87.0% (195/224), and became 94.2%\n(211/224) at a more lenient cutoff 1 × 10-8.\nIn the revised manuscript, we added two additional supplementary tables for the summary\nstatistics of RINT (Tables S22-S23) and discussed these results in the Discussion section and\nprovided the details in the Methods and Supplementary Note:\n“There is also ongoing research regarding the best practice of performing rare variant\nassociation tests. For example, Auer et al.127 pointed out that rare variant associations for\nquantitative traits might be susceptible to outliers and non-normality and ignoring such\nissues could result in false discoveries. Although outliers were removed in our pre-processing\nsteps for phenotype data (Methods), there is still potential concern about non-normality. We\nverified all the identified associations by applying rank-based inverse normal transformation\n(RINT) (Methods). The results largely supported our original findings (Tables S22-S23 and\nSupplementary Note). Briefly, all associations except one (LPAR3 with the mean MD of\nuncinate fasciculus) passed the nominal threshold adjusted for the number of signals\nundergoing RINT analysis (4.7 × 10-4 = 0.05/107 for variant-level tests and 2.2 × 10-4 =\n0.05/224 for gene-based burden tests, respectively). The majority of these associations also\nremained significant under our primary stringent threshold.” (Page 23, Lines 25-32; Page 24,\nLines 1-4)\n“To verify the robustness of the identified associations, rank-based inverse normal\ntransformation (RINT) was applied by calling ‘--apply-rint’ in both step 1 and step 2 of\nREGENIE for all identified associations in our study, while we followed the same testing\nprocedure. The summary statistics produced after applying RINT were compared with the\noriginal ones.” (Page 46, Lines 2-6)\n12\n4. For the validation analysis. I don’t think that the argument that the direction was\nmaintained in the larger cohort and that the p-values improve is a sufficiently strong one\n(e.g., P6 Line 11). As outlined in the previous comment, outliers may be driving the\nassociation and they may still be the driving force in the larger cohort. To make this a valid\nargument, the authors should show that the p-value improved beyond what would be\nexpected from the increased sample size alone.\nResponse: Thank you so much for your comments! As mentioned in the response to your\nComment 3, we have considered the potential issues caused by outliers in our preprocessing\nstep for phenotype data and results from rank-based inverse normal transformation (RINT)\nsupported our findings. In internal replication, we compared the P-values from the discovery\nsample and joint sample. For testing a null hypothesis of no effect, the summary statistics\nfrom the joint sample are sufficient. For the purpose of confirming the associations by\ninternal replication based on sample splitting, we intended to investigate whether signals\nderived from the discovery sample were attenuated in the enlarged sample (which might\nindicate some evidence for false discovery). However, it is not direct to quantify the\nconfidence in this setting. If we assume that the effect sizes are constant, we could expect\nthat the Z-scores observed in the discovery sample (or the corresponding P-values) would\nimprove based on the enlarged sample at a scale of square root of the sample size.\nConsistent with your comments, if the summary statistics (i.e. Z-scores and P-values)\nobserved in the joint sample have similar behaviors with the expected values, then such\nassociations would be reasonable. However, it is not direct to quantify the randomness in\nthis step, which led us to conduct a heuristic comparison between the analytical values and\nobserved values.\nFollowing your suggestions, we calculated the expected Z-scores (and the corresponding P-\nvalues) by multiplying the original Z-scores from the discovery sample by the square root of\nthe ratio of sample size of the joint sample and the discovery sample for each association.\nWe compared them with the observed Z-scores in the joint sample. Overall, they showed\nagreement and similar behaviors were observed. Specifically, among the 50 associations\nhaving smaller P-values in the joint sample, 34 of them have larger Z-scores than the\nexpected ones (or smaller P-values than the expected ones), while 40 of them have the\ncorresponding P-values at the same scale of the expected ones in terms of logarithm.\nNotably, the potentially different effect sizes across stages, the instable estimates from rare\nvariants in the smaller sample, and the sampling randomness in the realization could still\ncomplicate the interpretations here. In the revised manuscript, we have included them into\nTable S4, and discussed them in the Supplementary Note:\n“A comparison between the observed summary statistics in the internal analysis and the\nexpected improvements of these values in the enlarged joint sample is discussed in\nSupplementary Note.” (Page 7, Lines 2-4)\n5. The authors evaluate the variant hits in terms of existing GWAS hits. In particular for\nDTI traits they uncover a rare variant in the VCAN gene. In an earlier GWAS (with\ncomparable sample size) an intronic SNP in the gene was strongly associated with DTI\ntraits (rs67827860 (P = 4.06 × 10−37; Elliot et al. ). What is the interpretation here given\n13\nthe common variant shows a stronger association? Wouldn’t it be more likely that the rare\nvariants are only in LD with the variation causing driving the association?\nResponse: Thank you so much for raising this point! As noted, the common variant showed\na stronger association with DTI traits (has a smaller P-value) in a previous GWAS with\ncomparable sample size than two rare variants (rs2652098 and rs143368552) identified in\nour study. Typically, common variants, which often rely on LD tagging, would have greater\nstatistical power as more carriers of such mutations are observed in the cohort. This fact can\npartially explain why a smaller P-value was observed in the previous study when this\ncommon variant alone was included in the test. Conversely, rare variants are expected to\ncapture more ‘direct’ effect and have independent signals that are potentially causal\ncompared with common variants in GWAS (for example, this point has been mentioned in\nDeBoever et al. Nat Commu., PMID: 29691392, Wang et al. Nature, PMID: 34375979, and\nhas been observed in Backman et al. Nature, PMID: 34662886). However, the power might\nbe limited due the nature of rare variations as there are less carriers.\nTo clarify the relationships between these signals, we checked the LD and conducted\nconditional analysis for DTI traits. We found both of these two rare variants are in low LD\nwith rs67827860 (R2 = 0.008 for both), calculated in the European ancestry from LDLink\n(https://ldlink.nih.gov/?tab=home). This suggests that the rare variant is not simply tagging\nthe common variant.\nFurthermore, we conducted conditional analysis by conditioning on this nearby common\nvariant in REGENIE (the conditional analysis is restricted to step 2 only; the genome-wide\npredictors in step 1 were reused). Among 63 variant-level associations between these two\nvariants in VCAN and DTI traits (Table S2), 85.7% of these associations (54 out of 63)\nsurvived the stringent threshold at 2.8 × 10-10, while all of these 63 associations remain\nsignificant at a relaxed threshold 1 × 10-8 and have similar effect sizes.\nTaken together, these results support the interpretation that two missense variants\nidentified in our study likely represent additional, independent associations around VCAN\nlocus in addition to the previously reported GWAS hit (or the causal variation in LD with\nrs67827860). Following your suggestion, we have added the details of conditional analysis in\nthe method, discussed the results in the main text, and reported the summary statistics of\nconditional analysis in an additional supplementary table Table S5:\n“Notably, Elliot et al. identified an intronic variant rs67827860 with strong associations with\nDTI traits. Thus, we conducted conditional analysis for all the 63 associations between these\ntwo variants in VCAN and DTI traits by conditioning on this nearby common variant as a\ncovariate (Methods). We found that 85.7% of these associations still survived our stringent\nthreshold at 2.8 × 10-10 while all the associations passed our suggestive threshold at 1 × 10-8,\nwhich supports the newly identified missense variants in VCAN as independent signals in\naddition to the previously identified common variant (Table S5).” (Page 9, Lines 14-21)\n“Conditional analysis was performed by specifying the nearby common variants as\ncovariates to assess the independence of the identified signals. The conditional analysis is\n14\nrestricted to step 2 in REGENIE, while the genome-wide predictors in step 1 were reused.”\n(Page 46, Lines 6-9)\n6. When reading the manuscript it does not become quite clear which gene based burden\ntest the authors are using. There are multiple different ones implemented in regenie.\nPlease, can the authors provide detailed information (also in the results section on page 9\naround line 20, not just the methods section).\nResponse: Thank you so much for your careful review! For all the gene-based tests in this\nstudy, we used the burden test following the ‘max rule’ (default implementation in REGENIE,\nsee https://rgcgithub.github.io/regenie/options/#gene-based-testing). Specifically, the ‘max\nrule’ scores each person based on the most damaging genotype they carry in terms of\ndosage across the single variant sites of interest defined by annotation as a burden\nmask/score, and then this burden score is tested for association similar to the manner of the\nsingle variant test (also see https://rgcgithub.github.io/regenie/options/#gene-based-\ntesting). In the revised manuscript, these details are now provided in the main text as well\nas Methods:\n“In burden tests, the variants of interest defined by annotation and MAF (Methods) are\ncollapsed into a burden mask. Then, the burden score is calculated as the maximum\ngenotype dosage across the single variant sites in the mask, and this score is tested for\nassociation (Methods).” (Page 10, Lines 22-25)\n“Then for each given gene set, the score is calculated as the maximum of the genotype\ndosage across all the variants in this set (the ‘max rule’ as the default implementation in\nREGENIE), and the score is tested for association.” (Page 42, Lines 31-32; Page 43, Lines 1-2)\n7. I recommend adding to the discussion also the authors’ perspective on the value of rare-\nvariant analyses. Compared to GWAS with the same sample size, much fewer loci-trait\nassociations have been identified. Does this mean that perhaps one should focus cost and\neffort to other types of variation?\nResponse: Thank you so much for your suggestions! The fewer associations identified in our\nstudy compared with GWAS partially result from the limited sample size and the nature of\nrare mutations, though rare variants have stronger effect sizes. We also recognize that the\nsignals from rare variants and from common variants would provide complementary\ninformation, which together contributes to a more complete human genetic landscape.\nFollowing your recommendations, we have discussed the value of rare variant associations\nin the Discussion section in the revised manuscript:\n“Compared with GWAS of similar sample size (at the scale of 50,000) to our study, the\nidentified associations are much fewer partially due to the limited statistical power induced\nby the low frequency of the protein-coding variants17,18,62. However, rare coding variants\noffer more direct and stronger effects with interpretable biological implications, which would\ncontribute to complementary genetic information with common variants, prioritize the\neffector genes, and elucidate the underlying mechanisms between gene-trait pairs.” (Page\n26, Lines 5-11)\n1. P4 Line 28: “whole body and shared genetic regulation across organs.” I think the\nreasoning is a bit inconsistent in that the rare variants (which are limited to damaging,\nprotein altering variations) are viewed as regulators for processes across organs. Often\nthe focus on rare variation is to illuminate a more direct “mutation”-effect association\nrather than uncovering a regulatory process. I recommend to reword this.\nResponse: Thank you so much for your thoughtful suggestions to improve the wording here!\nWe totally agree with this point as our study focuses on rare variants annotated as\npredicted loss-of-function or damaging missense with interpretable biological consequences\nat protein-level. In this revision, we have rephrased this sentence:\n“whole body and shared genetic effects across organs.” (Page 4, Line 28)\nWe have also changed the wording in other places when suitable. The following are a few\nexamples:\n“elucidate some shared genetic effects across organs.” (Page 3, Line 14)\n“...provided additional evidence for its impact on heart structure and function.” (Page 10,\nLine 5)\n“Our novel discovery indicated that GEM may also impact the structural integrity of the\nascending aorta...” (Page 13, Lines 16-17)\n“The consequences of these rare coding variants to human organs were largely\nunderstudied.” (Page 22, Lines 25-26)\n2. When reporting SNP effect sizes and association strengths it would be useful to have\nthe MAF or MAC in the main text\nResponse: Thank you so much for your helpful comments! Following your suggestion, we\nwent through our paper and added the MAF within the parentheses when reporting the\neffect sizes (together with the information of P value and 95% CI) in the main text. The\nfollowing are a few examples:\n“...we identified associations between spleen volume and two missense variants in\nTNFRSF13B (rs72553883, MAF = 0.55%, effect = 0.025 L, effect = 0.37 s.d. units, 95% CI =\norg\n[0.29, 0.45], P = 5.6 × 10-19, and rs34557412, MAF = 0.66%, effect = 0.028 L, effect = 0.42\norg\ns.d. units, 95% CI = [0.34, 0.49], P = 2.0 × 10-28), ...” (Page 7, Lines 20-23)\n“...the strongest association between PTEN and total brain volume with extraordinary effect\nsize when aggregating the pLoF and missense variants (MAF = 0.01%, effect = 386703.40\norg\nmm3, effect = 2.77 s.d. units, 95% CI = [2.22, 3.32], P = 5.8 × 10-23).” (Page 14, Lines 7-10)\n16\n“KCNA5 was associated with a regional peak circumferential strain measurement (MAF =\n0.02%, effect = 5.11%, effect = 1.08 s.d. units, 95% CI = [0.63, 1.54], P = 3.2 × 10-6;\norg\n“Ecc_AHA_4”, Table S1) ...” (Page 17, Lines 27-29)\n3. P8 Line 11-13: I do not think the reasoning/causality here is one to support the\nassociation (i.e., the gene is associated with Schizophrenia (common) and that\nSchizophrenia is associated with enlarged ventricles etc) to support the link between the\ngene and brain volumes… Suggest rephrasing this.\nResponse: Thank you so much for pointing this out! We intended to point out fact that the\nidentified gene associated with both the brain volume and the mental disorder. In the\nrevised manuscript, we have simplified the wording and rephrased accordingly to avoid\npotentially overclaiming the causal mechanism:\n“ADRA1A was previously reported to be associated with schizophrenia43,44.” (Page 9, Lines 4-\n5)\n4. P38 Line 21 “Synonymous variants served as empirical null control to support the study-\nwise P-value threshold for burden tests”. The statement is a bit vague and unclear. The\nmethod behind this is only explained further down in the methods section. I recommend to\nadd a forward point to the relevant section to clarify that this will be explained below.\nResponse: Thank you so much for your thoughtful comments and suggestion on the writing!\nFollowing your suggestion, we have added a forward point to the Method section in the\nmain text:\n“Synonymous variants served as empirical null control to support the study-wise P-value\nthreshold for burden tests and the method would be detailed below in the section for\ndetermining study-wise significance level.” (Page 41 Lines 23-25)\n5. Regarding the Hardy-Weinberg filtering: I recommend to omit this step. Especially with\nrare variants checking the HW may not be useful. Previous work that examined the HW QC\nstep suggested to retain variants and only examine the HW after associations have been\nfound:\nhttps://www.medrxiv.org/content/10.1101/2024.02.07.24301951v2\nResponse: Thank you so much for the suggestion and bringing our attention into this paper\n(Greer et al.). Our sample size is around 50,000 and our HW filter was set at 1 × 10-15, which\nwould not remove too many variants as discussed in Greer et al. Specifically, in our previous\nquality control (QC), only 8,272 variants were filtered out due to the HW QC step but passed\nother QC steps. We carefully investigated the impact of these variants. Among these\nvariants, 2,202 variants were coding variants of interest based on our annotation method.\nThen, we conducted the association tests for these new variants and reused the results from\nstep 1 of REGENIE. For variant-level associations, none of them were significant under the\nsuggestive P-value cutoff at 1 × 10-8. For gene-based burden test, 270 unique genes were\naffected according to our definition for burden masks. When incorporating these variants\n17\ninto the gene-based burden test, there is no new signals identified using the suggestive\nthreshold. In summary, the impact of the HW step in our QC has minimal impact on our\nfindings.\nFollowing your suggestions, we mentioned this aspect of HW filtering in the Discussion\nsection and Supplementary Note:\n“In addition, Greer et al.128 suggested retaining all the variants before conducting\nassociation tests and examining the Hardy-Weinberg equilibrium afterward to avoid the\nover-removal of variants included in the study. We examined the variants filtered out by\nHardy-Weinberg equilibrium step in our quality control step and no new discoveries were\nidentified, suggesting that this step had minimal impact on our study (Supplementary\nNote).” (Page 24, Lines 4-10)\n6. It appears as if genetics that influences blood flow affects the fMRI phenotypes. This\nappears to be a bias in the analysis and interpretation of fMRI signals. Perhaps the\nauthors can comment on the need to incorporate cardiac features/genetics into the\nanalysis of fMRI data?\nResponse: Thank you so much for your comment! In the revised manuscript, we have\ncommented on this aspect of the analysis of fMRI data in the main text:\n“Therefore, it might be helpful to integrate cardiac features and cardiac-related genetic\ninformation when analyzing fMRI data, given that fMRI phenotypes are blood oxygen level-\ndependent and are influenced by blood flow.” (Page 16, Lines 11-13)\n18"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors appropriately addressed all comments.", "reviewer_2": "I would like to thank the authors for addressing all my comments. I have no further comments."}, "author_response": ""}]

Decision letter