Main

GWAS have identified thousands of genetic variants associated with complex phenotypes1,2, typically using samples of non-closely related individuals3. GWAS associations can be interpreted as estimates of direct individual genetic effects, that is, the effect of inheriting a genetic variant (or a correlated variant) on a phenotype4,5,6. However, there is growing evidence that GWAS associations for some phenotypes estimated from samples of unrelated individuals also capture effects of demography7,8 (assortative mating9,10,11 and population stratification12) and indirect genetic effects of relatives13,14,15,16,17,18,19 (Fig. 1). For example, Lee et al.14 found that within-sibship GWAS estimates for educational attainment variants were around 40% lower than estimates from unrelated individuals, indicating the presence of demographic and indirect genetic effects. These nondirect sources of genetic associations are themselves of interest for estimating parental effects13,18, understanding human mate choice9,10,11 and genomic prediction14,19. However, they can also impact downstream analyses using GWAS summary data such as biological annotation, heritability estimation20,21,22, genetic correlations23, Mendelian randomization (MR)7,24,25 and polygenic adaptation tests26,27,28,29.

Fig. 1: Demographic and indirect genetic effects.
figure 1

Population stratification: population stratification is defined as the distortion of associations between a genotype and a phenotype when ancestry A influences both genotype G (via differences in allele frequencies) and the phenotype X. Principal components and linear mixed model methods control for ancestry but they may not completely control for fine-scale population structure. Assortative mating: assortative mating is a phenomenon where individuals select a partner based on phenotypic (dis)similarities. For example, tall individuals may prefer a tall partner. Assortative mating can induce correlations between causes of an assorted phenotype in subsequent generations. If a phenotype X is influenced by two independent genetic variants G1 and G2 then assortment on X (represented by effects of X on mate choice M) will induce positive correlations between G1 in parent 1 and G2 in parent 2 and vice versa. Parental transmission will then induce correlations between otherwise independent G1 and G2 in offspring. These correlations can distort genetic association estimates. Indirect genetic effects: indirect genetic effects are effects of relative genotypes (via relative phenotypes and the shared environment) on the index individual’s phenotype. These indirect effects influence population GWAS estimates because relative genotypes are also associated with genotypes of the index individual. Indirect genetic effects of parents on offspring are of most interest because they are likely to be the largest. However, indirect genetic effects of siblings or more distal relatives are also possible.

Within-family genetic association estimates, such as those obtained from samples of siblings, can provide less biased estimates of direct genetic effects because they are unlikely to be affected by demographic and indirect genetic effects of parents7,17,30,31,32,33,34. GWAS using siblings (within-sibship GWAS) (Fig. 2) have been previously limited by available data, but are now feasible by combining well-established family studies with recent large biobanks that incidentally or by design contain thousands of sibships35,36,37,38,39.

Fig. 2: Population GWAS estimate the association between raw genotypes G and phenotypes X.
figure 2

As outlined in Fig. 1, estimates from population GWAS may not fully control for demography (population stratification and assortative mating) and may also capture indirect genetic effects of relatives. For simplicity we use N to represent all sources of associations between G and X that do not relate to direct effects of G. Circles indicate unmeasured variables and squares indicate measured variables. If parental genotypes are known, G can be separated into nonrandom (determined by parental genotypes) and random (relating to segregation at meiosis) components. Within-sibship GWAS include the mean genotype across a sibship (GF) (a proxy for the mean of the paternal and maternal genotypes GP, M) as a covariate to capture associations between G and X relating to parents. The within-sibship estimate is defined as the effect of the random component: that is, the association between family-mean-centered genotype GC (that is, G − GF) and X. Demography and indirect genetic effects of parents (N) will be captured by GF. The association between GC and X will not be influenced by these sources of association but could be affected by indirect effects of the siblings themselves, which are not controlled for.

Here, we report findings from a within-sibship GWAS of 25 phenotypes using data from 178,076 siblings from 19 studies, the largest GWAS conducted within sibships to date (Fig. 3). Our results are broadly consistent with previous studies comparing population and within-sibship genetic effect estimates in smaller sample sizes13,14,19,40. We found that within-sibship meta-analysis GWAS estimates are smaller than population estimates for seven phenotypes (height, educational attainment, age at first birth, number of children, cognitive ability, depressive symptoms and smoking). We show that these differences in GWAS estimates, which are likely to partially reflect demographic and indirect genetic effects, can affect downstream analyses such as estimates of heritability, genetic correlations and MR. However, we find that genetic associations with most clinical phenotypes, such as lipids, are less strongly affected. We found strong evidence of polygenic adaption on taller human height using within-sibship data. Our study illustrates the importance of collecting genome-wide data from families to understand the effects of inherited genetic variation on phenotypes that are affected by assortative mating, population stratification and indirect genetic effects.

Fig. 3: A flowchart of analyses undertaken in this project.
figure 3

We started by performing quality control and running GWAS models in 19 individual cohorts. We then meta-analyzed GWAS data from 18 of these cohorts with European-ancestry individuals. We then used the European meta-analysis data for downstream analyses including LDSC, MR and polygenic adaptation testing. We performed analyses in the China Kadoorie Biobank separately. QC, quality control.

Results

Within-sibship and population-based GWAS comparison

For GWAS analyses we used data from 178,076 individuals (with one or more genotyped siblings) from 77,832 sibships in 19 studies. Sample sizes for individual phenotypes ranged from 13,375 to 163,748 (median: 82,760, mean: 79,794). More information on sample sizes from individual cohorts and for each phenotype is contained in Supplementary Table 1. We used within-sibship models which use deviations of the individual’s genotype from the mean genotype within the sibship (that is, all siblings in the family present in the study). For example, in a sibling pair where one sibling has two risk alleles and the other sibling has one risk allele, the mean sibship genotype is 1.5 risk alleles and the individual’s deviations are +0.5 and −0.5, respectively. The within-sibship model includes the mean sibship genotype as a covariate to capture the between-family contribution of the SNP14. For comparison, we also applied a standard population GWAS model; a covariate-adjusted linear regression of the outcome on raw genotype, which does not account for the mean sibship genotype. Standard errors were clustered by sibship. Age, sex and principal components were included as covariates in both models. All GWAS analyses were performed in individual cohort studies separately using R (v.3.5.1) and meta-analyses were conducted across these using summary data. Amongst the phenotypes analyzed, the largest available sample sizes in a meta-analysis of European cohorts were for height (N = 149,174), body mass index (BMI) (N = 140,883), educational attainment (N = 128,777), ever smoking (N = 124,791) and systolic blood pressure (SBP) (N = 109,588) (Supplementary Table 2). We also report stratified results from non-European samples including 13,856 individuals from the China Kadoorie Biobank. Sample sizes here refer to the number of individuals across all sibships.

Previous studies have found that association estimates of height and educational attainment genetic variants are smaller in within-family models13,14,40. We aimed to investigate whether similar shrinkage in association estimates is observed for other phenotypes by comparing within-sibship and population genetic association estimates for 25 phenotypes that were widely available in family-based studies. We observed the largest within-sibship shrinkage (% decrease in association estimates from population to within-sibship models) for genetic variants associated with number of children (67%; 95% confidence interval (95% CI) 4%, 130%), age at first birth (52%; 30%, 75%), depressive symptoms (50%; 18%, 82%) and educational attainment (47%; 41%, 52%). We also found evidence of shrinkage for cognitive ability (22%; 6%, 37%), ever smoking (19%; 9%, 30%) and height (10%; 8%, 12%). In contrast, within-sibship association estimates for C-reactive protein (CRP) were larger than population estimates (−9%; −15%, −2%). We found limited evidence of within-sibship differences for the remaining 17 phenotypes, including BMI and SBP (Fig. 4 and Supplementary Table 3).

Fig. 4: Estimates of shrinkage between population and within-sibship models with corresponding 95% CIs.
figure 4

Shrinkage is defined as the % decrease in association between the relevant weighted score and phenotype when comparing the population estimate with the within-sibship estimate. Shrinkage was computed as the ratio of two weighted score association estimates with standard errors derived using leave-one-out jackknifing. The number of individuals contributing to each phenotype ranged from n = 149,174 for height to n = 13,375 for age at menopause. Further information on the sample sizes of each phenotype is contained in Supplementary Table 2. SG, weighted score at genome-wide significance (P < 5 × 10−8); SL, weighted score at more liberal threshold (P < 1 × 10−5); Education, educational attainment; EverSmk, ever smoking; WHR, waist-to-hip ratio; Alcohol, weekly alcohol consumption; Menarche, age at menarche; AFB, age at first birth; Children, number of biological children; Menopause, age at menopause; Cognition, cognitive ability; Depressive, depressive symptoms; PA, physical activity; CPD, cigarettes per day; LDL, low-density lipoprotein-cholesterol; HDL, HDL-cholesterol; TG, triglycerides; eGFR, estimated glomerular filtration rate; FEV1, forced expiratory volume; FEV1FVC, ratio of FEV1/forced vital capacity; HbA1c, hemoglobin A1C.

We investigated possible heterogeneity in shrinkage for height and educational attainment genetic variants across variants and between cohorts. Using the meta-analysis results, we did not observe strong evidence of heterogeneity in shrinkage across variants that were strongly associated with height and educational attainment. This suggests that shrinkage may be largely uniform across these signals for these phenotypes. We also found limited evidence of cohort heterogeneity in shrinkages for height (heterogeneity P = 0.89) and educational attainment (P = 0.40) across the European-ancestry cohorts (Extended Data Figs. 1 and 2). In contrast, there was limited evidence for shrinkage on height in the China Kadoorie Biobank (shrinkage −3%; 95% CI −13%, 7%; heterogeneity with European meta-analysis P = 0.006) but some evidence of shrinkage on ever smoking (shrinkage = 134%; 10%, 258%) (Extended Data Fig. 3).

Within-sibship SNP heritability estimates

Linkage disequilibrium (LD) score regression (LDSC) can use GWAS data to estimate SNP heritability, the proportion of phenotypic variation explained by common SNPs20,23. We used simulations to investigate the applicability of LDSC when using within-sibship GWAS data, finding evidence that LDSC can estimate SNP heritability using both population and within-sibship model GWAS data if effective sample sizes (based on standard errors) are used to account for differences in power between the models (Methods).

To evaluate the impact of controlling for demographic and indirect genetic effects, we compared LDSC SNP heritability estimates based on population and within-sibship effect estimates for 25 phenotypes. Theoretically, within-sibship shrinkage in GWAS estimates will also lead to attenuations in within-sibship SNP heritability estimates (Methods). The within-sibship SNP heritability point estimate for educational attainment attenuated by 76% from the population estimate (population h2: 0.13; within-sibship h2: 0.04; difference P = 5.3 × 10−26), with attenuations also observed for cognitive ability (population h2: 0.24; within-sibship h2: 0.14; attenuation 44%; difference P = 0.011), ever smoking (population h2: 0.10; within-sibship h2: 0.07; attenuation 25%; difference P = 0.029) and height (population h2: 0.41; within-sibship h2: 0.34; attenuation 17%; difference P = 1.6 × 10−3). The observed attenuations were consistent with theoretical expectation (Supplementary Table 4), suggesting that the lower within-sibship SNP heritability estimates are explained by genetic association estimate shrinkage. Across the 21 additional phenotypes, population and within-sibship SNP heritability estimates were relatively consistent (Fig. 5 and Supplementary Table 5). SNP heritability estimates using SumHer21 with the LDAK-Thin model (expected heritability contribution of each SNP is dependent on allele frequencies and local LD) provided consistent evidence for within-sibship attenuations in SNP heritability for height, educational attainment and cognitive ability (Supplementary Table 6 and Extended Data Fig. 4).

Fig. 5: LDSC SNP h2 estimates for 25 phenotypes using population and within-sibship meta-analysis data with corresponding 95% CIs.
figure 5

The number of individuals contributing to each phenotype ranged from n = 149,174 for height to n = 13,375 for age at menopause. BMI, body mass index; Education, educational attainment; EverSmk, ever smoking; SBP, systolic blood pressure; WHR, waist-hip ratio; Alcohol, weekly alcohol consumption; Menarche, age at menarche; AFB, age at first birth; Children, number of biological children; Menopause, age at menopause; Cognition, cognitive ability; Depressive, depressive symptoms; PA, physical activity; CPD, cigarettes per day; LDL, LDL cholesterol; HDL, HDL cholesterol; TG, triglycerides; CRP, C-reactive protein; eGFR, estimated glomerular filtration rate; FEV1, forced expiratory volume; FEV1FVC, ratio of FEV1/forced vital capacity; HbA1c, Haemoglobin A1C. Further information on the sample sizes of each phenotype is contained in Supplementary Table 2.

Within-sibship r g with educational attainment

We used LDSC23 to estimate cross-phenotype genome-wide genetic correlations (rg) between educational attainment and 20 phenotypes with sufficient heritability (population/within-sibship h2 point estimate > 0) and statistical power. To determine the effects of demographic and indirect genetic effects on rg, we compared estimates of rg using population and within-sibship estimates.

There was strong evidence using population estimates that educational attainment is negatively correlated with BMI (rg = −0.32; −0.37, −0.26), ever smoking (rg = −0.41; −0.49, −0.34) and CRP (rg = −0.46; −0.67, −0.25). However, these correlations attenuated towards zero when using within-sibship estimates: BMI (rg = −0.05; −0.22, 0.12), ever smoking (rg = −0.14; −0.42, 0.14) and CRP (rg = −0.06; −0.43, 0.30), with some evidence at nominal significance for differences between population and within-sibship rg estimates (BMI difference P = 5.3 × 10−4, ever smoking difference P = 0.040, CRP difference P = 0.039). These attenuations indicate that genetic correlations between educational attainment and these phenotypes from population estimates may be inflated by demographic and indirect genetic effects (Fig. 6 and Supplementary Table 7).

Fig. 6: LDSC rg estimates between educational attainment and 20 phenotypes using population and within-sibship meta-analysis data with corresponding 95% CIs.
figure 6

The number of individuals contributing to the educational attainment GWAS was n = 128,777 with sample sizes for outcomes ranging from n = 149,174 for height to n = 27,638 for cognitive ability. BMI, body mass index; Education, educational attainment; EverSmk, ever smoking; SBP, systolic blood pressure; WHR, waist-hip ratio; Alcohol, weekly alcohol consumption; Menarche, age at menarche; AFB, age at first birth; Children, number of biological children; Menopause, age at menopause; Cognition, cognitive ability; Depressive, depressive symptoms; CPD, cigarettes per day; LDL, LDL cholesterol; HDL, HDL cholesterol; TG, triglycerides; CRP, C-reactive protein; eGFR, estimated glomerular filtration rate; FEV1, forced expiratory volume; HbA1c, Haemoglobin A1C. Further information on the sample sizes of each phenotype is contained in Supplementary Table 2.

Within-sibship MR (WS-MR): effects of height and BMI

MR uses genetic variants as instrumental variables to assess the causal effect of exposure phenotypes on outcomes24,41. MR was originally conceptualized in the context of parent–offspring trios where offspring inherit a random allele from each parent24. However, with limited family data, most MR studies have used data from unrelated individuals. WS-MR is largely robust against demographic and indirect genetic effects that could distort estimates from nonfamily designs7,25. Here, we used population MR and WS-MR to estimate the effects of height and BMI on 23 phenotypes. These provide a useful comparison as we find evidence of shrinkage in GWAS estimates for height but little evidence of shrinkage for BMI, and both height and BMI have large sample sizes.

WS-MR estimates for height and BMI on the 23 outcome phenotypes were largely consistent with population MR estimates for height based on the slope of a regression of the WS-MR and population MR estimates (−3%; 95% CI −16%, 10%) and BMI (−5%; 95% CI −14%, 4%). However, in agreement with the genetic correlation analyses, we observed differences between population MR and WS-MR estimates of height and BMI on educational attainment. Population MR estimates provided strong evidence that taller height and lower BMI increase educational attainment (0.06 s.d. increase in education per s.d. taller height; 95% CI 0.04, 0.07; 0.19 s.d. decrease in education per s.d. higher BMI; 0.16, 0.22). In contrast, WS-MR estimates for these relationships were greatly attenuated (height: 0.02 s.d. increase; −0.01, 0.04; difference P = 1.2 × 10−3; BMI: 0.05 s.d. decrease; 0.01, 0.09; difference P = 2.8 × 10−7). We also observed similar attenuation from population MR and WS-MR estimates for BMI on age at first birth (difference P = 2.3 × 10−3) and cognitive ability (difference P = 0.020); phenotypes highly correlated with education. These differences illustrate instances where population-based MR estimates might be distorted by demographic and indirect genetic effects or other factors (Table 1).

Table 1 WS-MR: effects of height and BMI on 23 phenotypes

Polygenic adaptation

Polygenic adaptation is a process via which phenotypic changes in a population over time are induced by small shifts in allele frequencies across thousands of variants. One method of testing for polygenic adaptation is to compare Singleton Density Scores (SDS), measures of natural selection over the previous 2,000 years (ref. 28), with GWAS P values. However, this approach is sensitive to population stratification as illustrated by recent work using UK Biobank data which showed that population stratification in GWAS data likely confounded previous estimates of polygenic adaptation on height26,27. Within-sibship GWAS data are particularly useful in this context as they are robust against population stratification26,27,29. Here, we recalculated Spearman’s rank correlation (r) between tSDS (SDS aligned with the phenotype-increasing allele) and our population/within-sibship GWAS P values for 25 phenotypes, with standard errors estimated using jackknifing over blocks of genetic variants.

We found strong evidence for polygenic adaptation on taller height in the European meta-analysis GWAS using both population (r = 0.022; 95% CI 0.014, 0.031) and within-sibship GWAS estimates (r = 0.012; 0.003, 0.020) (Extended Data Figs. 5 and 6). These results were supported by several sensitivity analyses: (1) evidence of enrichment for positive tSDS (mean = 0.18, s.e. = 0.06, P = 0.003) amongst 310 putative height loci from the within-sibship meta-analysis results (Extended Data Fig. 7); (2) positive LDSC rg between height and tSDS in the meta-analysis results (Supplementary Table 8); and (3) evidence for polygenic adaptation on taller height when meta-analyzing correlation estimates from eight individual studies (for example, SDS using only UK Biobank GWAS summary data) for population (r = 0.013; 0.010, 0.015) and within-sibship (r = 0.004; 0.002, 0.007) estimates (Fig. 7). There was also some putative within-sibship evidence for polygenic adaptation on increased number of children (P = 0.024) and lower high-density lipoprotein (HDL)-cholesterol (P = 0.024) (Extended Data Fig. 5).

Fig. 7: Spearman rank correlation estimates and corresponding 95% CIs between tSDS (SDS aligned with height-increasing alleles) and absolute height Z scores.
figure 7

Positive correlations indicate evidence of historical positive selection on height-increasing alleles. The pooled estimate is a meta-analysis of the correlation estimates from the individual studies shown above while the European meta-analysis estimate is the correlation estimate using the meta-analysis GWAS data. The number of individuals in the meta-analysis estimate was n = 149,174 with the sample sizes for the displayed individual studies ranging from n = 40,068 in UK Biobank to 4,708 in the Netherlands Twin Register. Further information on available height data in each phenotype is contained in Supplementary Table 1. QIMR, Queensland Institute of Medical Research.

Discussion

Here, we report results from the largest within-sibship GWAS to date which included 25 phenotypes and combined data from 178,076 siblings. Consistent with previous studies13,14,19,40, we found that GWAS results and downstream analyses of behavioral phenotypes (for example, educational attainment, smoking behavior) as well as some anthropometric phenotypes (for example, height, BMI) are affected by demographic and indirect genetic effects. However, we found that most analyses involving more molecular phenotypes, such as lipids, were not strongly affected. This suggests that the best strategy for gene discovery and polygenic prediction for these phenotypes remains to maximize sample sizes using unrelated individuals. For phenotypes sensitive to demographic and indirect genetic effects, such as educational attainment, family-based estimates are likely to provide less biased estimates of direct genetic effects.

A key aim of GWAS is to estimate direct genetic effects on phenotypes, but other sources of genetic associations can be extremely informative. For example, knowledge of indirect genetic effects can be used to elucidate maternal effects15,42 or the extent to which health outcomes are mediated by family environments13,18. Future family-based GWAS could also provide further estimates of indirect genetic effects6,18,43.

We found little evidence of heterogeneity in shrinkage estimates at genetic variants strongly associated (P < 1 × 10−5) with height and educational attainment, although power was limited by available samples. The limited detectable heterogeneity could indicate that the observed shrinkage is largely driven by assortative mating or indirect genetic effects. Both of these tend to influence associations proportional to the direct effect, whereas population stratification is likely to have larger effects on ancestrally informative markers. Notably, twin studies have indicated effects of the common environment on many of the phenotypes for which we observed shrinkage, such as educational attainment44, cognitive ability45 and smoking46, potentially consistent with indirect genetic effects of parents. In contrast, twin studies do not find strong evidence for common environmental effects on height, where shrinkage is more likely to be a consequence of assortative mating10,46,47.

The weak evidence for within-sibship shrinkage in the association between BMI genetic variants and BMI is in contrast to the strong evidence from MR analyses (and genetic correlation analyses) that the association between BMI genetic variants and educational attainment does attenuate. These results indicate cross-trait shrinkage in association estimates for BMI genetic variants even in the absence of same-trait shrinkage.

Within-sibship GWAS data can be useful for validating results from larger samples of unrelated individuals. Here, we showed that population MR and WS-MR estimates of the effects of height and BMI were generally consistent for 23 outcome phenotypes. However, we observed differences between within-sibship and population MR estimates of height (on educational attainment) and BMI (on educational attainment, cognitive ability and age at first birth). This suggests the MR assumptions do not hold for these relationships in samples of unrelated individuals. In subsequent studies, WS-MR could be used as a sensitivity analysis when including phenotypes likely to be affected by demographic or indirect genetic effects7,25.

We used non-European data from the China Kadoorie Biobank to evaluate whether demographic and indirect genetic effects influence GWAS analyses conducted in the Chinese population. In this sample, we found minimal evidence of shrinkage for height genetic variants but—consistent with the European meta-analysis—suggestive evidence of shrinkage for variants associated with smoking initiation. The absence of shrinkage for height suggests that demographic effects such as assortative mating may differ between populations. Larger within-family studies in non-European populations could be used to evaluate population differences in demographic and indirect effects.

We also used the within-sibship GWAS data to evaluate evidence for recent selection. A previous study reporting polygenic adaptation on height in the UK population was found to be biased by population stratification in the Genetic Investigation of ANthropometric Traits (GIANT) consortium26,27,28. Previous evidence for adaptation on height using siblings in UK Biobank was suggestive of some adaptation, but statistically inconclusive26. Here, using within-sibship GWAS estimates from a larger (~4-fold) sample of siblings, we found strong evidence of polygenic adaptation on increased height and some evidence of adaptation on number of children and HDL-cholesterol. We anticipate that future studies on human evolution will benefit from using large within-family datasets such as our resource.

Within-family GWAS are limited by both available family data and statistical inefficiency (homozygosity within families). To help address this issue, future population-based biobanks could recruit the partners, siblings and offspring of study participants. In contrast, conventional population GWAS designs sampling unrelated individuals are likely to be the optimal approach to maximize statistical power for discovery GWAS for genetic associations. Indeed, we found that many genotype–phenotype associations from population GWAS models were also observed in within-sibship GWASs, albeit sometimes attenuated towards zero. A notable limitation of within-sibship models is that they do not control for indirect genetic effects of siblings, that is, effects of sibling genotypes on the shared environment. Sibling effects have been estimated to be modest compared with parental effects6,48 but could have impacted our GWAS estimates. Another limitation is that while assortative mating is unlikely to affect within-sibship GWAS estimates, it can bias within-sibship estimates of heritability downwards49 and so may have affected our LDSC SNP heritability and genetic correlation estimates. However, the within-sibship shrinkage in GWAS estimates and LDSC heritability estimates were largely consistent, suggesting any such bias is unlikely to have impacted our conclusions. Our findings are also limited to adult phenotypes. Within-family GWAS (for example, using parent–offspring trios) could use data from children to evaluate if childhood phenotypes are more strongly affected by indirect genetic effects.

Methods

Study participants

Nineteen cohorts contributed data to the overall study (Supplementary Table 1). These cohorts were selected on the basis of having at least 500 genotyped siblings (an individual with 1 or more siblings in the study sample) with at least 1 of the 25 phenotypes that were analyzed in the study. Phenotypes were selected based on available data and to include a range of different phenotypes. Detailed information on genotype data, quality control and imputation processes are provided in the Cohort Descriptions in the Supplementary Materials. Individual cohorts defined each phenotype based on suggested definitions from an analysis plan (see the Phenotype Definitions in the Supplementary Materials).

GWAS analyses

GWAS analyses were performed uniformly across individual studies using automated scripts and a preregistered analysis plan (https://github.com/LaurenceHowe/SiblingGWAS). Scripts checked strand alignment, imputation scores and allele frequencies for the genetic data as well as missingness for covariates and phenotypes. Scripts also summarized covariates and phenotypes and set phenotypes to missing for sibships if only one individual in the sibship had nonmissing phenotype data. To harmonize variants for meta-analysis, genetic variants were renamed in a format including information on chromosome, base pair and polymorphism type (SNP or INDEL: insertion or deletion). The automated pipeline restricted analyses to common genetic variants (minor allele frequency (MAF) > 0.01) and removed poorly imputed variants (INFO: information score < 0.3). Analyses were restricted to include individuals in a sibship, that is, a group of two or more full siblings in the study. Monozygotic twins were included if they had an additional sibling in the study.

GWAS analyses involved fitting both population and within-sibship models to the same samples. The population model is synonymous with a conventional principal component adjusted model, and was fit using linear regression in R (v.3.5.1). The within-sibship model is an extension of the population model including the mean sibship genotype (the mean genotype of siblings in each sibship) as a covariate to account for family structure, with each individual’s genotype centered around the mean sibship genotype7,14. Age, sex and up to 20 principal components (10 principal components were included in smaller studies at the discretion of study co-authors) were included as covariates in both models. The pipeline used imputed ‘best guess’ genotype calls rather than dosage data.

For individual j in sibship i with ni > 2 siblings:

Population model:

$${{{\mathrm{Phenotype}}}}_{{{{ij}}}}\sim {{{\mathrm{G}}}}_{{{{ij}}}} + {{{\mathrm{Sex}}}}_{{{{ij}}}} + {{{\mathrm{Age}}}}_{{{{ij}}}} + {{{\mathrm{PC}}}}1_{{{{ij}}}} + {{{\mathrm{PC}}}}20_{{{{ij}}}}$$

Within-sibship model:

$${{{\mathrm{Phenotype}}}}_{{{{ij}}}}\sim {{{\mathrm{G}}}}_{{{{ij}}}}^{{{\mathrm{C}}}} + {{{\mathrm{G}}}}_{{{i}}}^{{{\mathrm{F}}}} + {{{\mathrm{Sex}}}}_{{{{ij}}}} + {{{\mathrm{Age}}}}_{{{{ij}}}} + {{{\mathrm{PC}}}}1_{{{{ij}}}} + {{{\mathrm{PC}}}}20_{{{{ij}}}}$$

where\({{{\mathrm{G}}}}_{{{i}}}^{{{\mathrm{F}}}} = \frac{{\mathop {\sum }\nolimits_1^{{{n}}} {{{\mathrm{G}}}}_{{{{ij}}}}}}{{{{n}}}}\,{{{\mathrm{and}}}}\,{{{\mathrm{G}}}}_{{{{ij}}}}^{{{\mathrm{C}}}} = {{{\mathrm{G}}}}_{{{{ij}}}} - {{{\mathrm{G}}}}_{{{i}}}^{{{\mathrm{F}}}}\)

Gij, genotype of sibling j in sibship i; \({{{\mathrm{G}}}}_{{{i}}}^{{{\mathrm{F}}}}\), mean family genotype for sibship i over n siblings; \({{{\mathrm{G}}}}_{{{{ij}}}}^{{{\mathrm{C}}}}\), genotype of sibling j in sibship i centered around \({{{\mathrm{G}}}}_{{{i}}}^{{{\mathrm{F}}}}\); PC, principal component.

Standard errors from both estimators were clustered over families at the sibling level to account for nonrandom clustering of siblings within families. Note that this clustering accounts for sibling relationships but does not account for further relatedness present in each sample. For example, a sibling pair could be related to another sibling pair (that is, two pairs of siblings who are first-cousins). We performed simulations, described below, confirming that such relatedness can lead to underestimating standard errors in the population model and has no effect on the standard errors of the within-sibship model.

GWAS models were performed in individual studies, harmonized and then meta-analyzed for each phenotype using a fixed-effects model in METAL50 with population and within-sibship data meta-analyzed separately. We performed meta-analyses using only samples of European ancestry. We used data from 13,856 individuals from the China Kadoorie Biobank separately in downstream analyses. Information on sample sizes for individual phenotypes is contained in Supplementary Table 2. Information on further quality control performed before meta-analysis is detailed in the Supplementary Methods.

Meta-analysis

Phenotypes were harmonized between studies using phenotypic summary data on means and standard deviations. GWAS of study-specific phenotypes that did not conform to analysis plan definitions (for example, binary instead of continuous) were excluded from meta-analyses. GWAS presented in different continuous units (for example, not standardized) were transformed before meta-analysis by dividing association estimates and standard errors by the standard deviation of the phenotype as measured in the cohort. Meta-analyses for 25 phenotypes were performed using a fixed-effects model in METAL50.

Within-sibship and population-based GWAS comparison

Overview

We hypothesized that the within-sibship estimates would differ compared with population-based estimates due to the exclusion of effects from demographic and familial pathways. In general, these effects have been shown to inflate (rather than shrink) population-based estimates, so we estimated within-sibship shrinkage (the % difference from population to within-sibship estimates). To estimate this shrinkage, we required estimates of the associations with a phenotype from each within-sibship and population-based analysis that was not affected by winner’s curse. Hence, we adopted a strategy where we used an independent reference dataset to select the variants associated with a phenotype. Using the meta-analysis results to obtain association estimates for these variants, we generated summary-based weighted scores of those association estimates in the within-sibship and population-based analyses and estimated the ratio of those scores. We used the UK Biobank dataset excluding sibling data as the independent reference dataset.

GWAS in independent reference discovery dataset

We performed GWAS in an independent sample of UK Biobank (excluding siblings) for each phenotype using a linear mixed model as implemented in BOLT-LMM51. We started with a sample of 463,006 individuals of ‘European’ ancestry derived using in-house k-means cluster analysis performed using the first four principal components provided by UK Biobank with standard exclusions also removed52. To remove sample overlap, we then excluded the sibling sample (N = 40,276), resulting in a final sample of 422,730 individuals. To model population structure in the sample, we used 143,006 directly genotyped SNPs, obtained after filtering on MAF > 0.01; genotyping rate > 0.015; Hardy–Weinberg equilibrium P < 0.0001; and LD pruning to an r2 threshold of 0.1 using PLINK v.2.0 (ref. 53). Age and sex were included in the model as covariates.

All 25 phenotypes (conforming to our phenotype definition) were available in UK Biobank data except for a continuous measure of depressive symptoms. For depressive symptoms, we performed a GWAS of binary depression which was excluded from the meta-analysis (see definition in Supplementary Materials). Using the BOLT-LMM UK Biobank GWAS summary data, we performed strict LD clumping in PLINK v.2.0 (ref. 53) (r2 < 0.001, physical distance threshold = 10,000 kb) using the 1000 Genomes Phase 3 EUR reference panel54 to generate independent variants associated with each phenotype at genome-wide significance (P < 5 × 10−8) and at a more liberal threshold (P < 1 × 10−5).

Summary-based weighted scores

For a particular phenotype the sets of independent variants obtained from the independent UK Biobank GWAS were used to generate a summary-based weighted score using an inverse variance weighting (IVW) approach55,56:

$${{{S}}} = \frac{{\mathop {\sum }\nolimits_{{{k}}}^{{{M}}} \frac{{{{{w}}}_{{{k}}}\beta _{{{k}}}}}{{\sigma _{{{k}}}^2}}}}{{\mathop {\sum }\nolimits_{{{k}}}^{{{M}}} \frac{{{{{w}}}_{{{k}}}^2}}{{\sigma _{{{k}}}^2}}}}$$

with standard error

$$\sigma _{{{S}}} = \sqrt {\frac{1}{{\mathop {\sum }\nolimits_{{{k}}}^{{{M}}} \frac{{{{{w}}}_{{{k}}}^2}}{{\sigma _{{{k}}}^2}}}}}$$

Here, the score S represents the weighted average of the association estimates of the M variants on a phenotype, where β and σ represent the beta coefficients and standard errors from the within-sibship (W) or population-based (P) meta-analysis results. The discovery association estimates from the UK Biobank GWAS were used as weights (w). The set of M variants were determined using either the genome-wide significance (G) or the more liberal threshold (L). Hence, depending on which model is used to determine the association estimates and which set of SNPs are used, four scores can be calculated for each phenotype—SP,G, SP,L, SW,G and SW,L.

These sets of scores were obtained for each of the 25 phenotypes with weights for binary depression used as a substitute for depressive symptoms because a suitable measure was unavailable in UK Biobank. The scores were strongly associated with the set of phenotypes in the meta-analysis data based on determining P values from their Z scores. The SW,L scores were nominally associated at P < 0.05 for 24 of 25 (exception: number of children) of the phenotypes, with the SP,L scores associated with all 25 phenotypes at this threshold (Supplementary Table 9).

Estimating shrinkage from population to within-sibship estimates

We used the within-sibship and population-based scores to calculate the average shrinkage (δ, that is, proportion decrease) of genetic variant–phenotype associations

$$\delta = 1 - \frac{{{{{S}}}_{{{{\mathrm{W}}}},}}}{{{{{\mathrm{S}}}}_{{{{\mathrm{P}}}},}}}$$

The standard errors of δ could be estimated using the delta method as below using the standard errors of the scores and the covariance between the scores Cov(Sw, SP,):

$$\sigma _{\delta}\sim \left( {\frac{{{{{S}}}_{{{{\mathrm{W}}}},}}}{{{{{S}}}_{{{{\mathrm{P}}}},}}}} \right)\sqrt {\left( {\frac{{\sigma _{{{{S}}}_{{{{\mathrm{W}}}},}}^2}}{{{{{S}}}_{{{{\mathrm{W}}}},}^2}} + \frac{{\sigma _{{{{S}}}_{{{{\mathrm{P}}}},}}^2}}{{{{{S}}}_{{{{\mathrm{P}}}},}^2}}} \right) - \frac{{2{{{\mathrm{Cov}}}}\left( {{{{S}}}_{{{{\mathrm{w}}}},},{{{S}}}_{{{{\mathrm{P}}}},}} \right)}}{{{{{S}}}_{{{{\mathrm{W}}}},}{{{S}}}_{{{{\mathrm{P}}}},}}}}$$

However, we do not have an estimate of this covariance term because the two GWAS were fit in separate regression models. We therefore used the jackknife to estimate \(\sigma _{\delta }\). For a score of M variants, we removed each variant in turn and repeated IVW and shrinkage analyses as above, extracting the shrinkage point estimate in each of the M iterations. We then calculated \(\sigma _{\delta }\) as follows:

$$\sigma _{\delta } = \sqrt {\frac{{{{{M}}} - 1}}{{{{M}}}}\mathop {\sum}\limits_{{{k}}}^{{{M}}} {(\sigma _{\delta _{,{{{k}}}}} - \mu )^2} }$$

where

$$\mu = \frac{{\mathop {\sum }\nolimits_{{{k}}}^{{{M}}} \sigma _{\delta _{,{{{k}}}}}}}{{{{M}}}}$$

As a sensitivity analysis, we investigated the effects of positive covariance between the population and within-sibship models on the shrinkage standard errors using individual-level participant data from UK Biobank. Analyzing shrinkage on height, we used seemingly unrelated regression to estimate the covariance term between the population and within-sibship estimators. We found that standard errors for shrinkage estimates decreased by around 15% when the covariance was modeled (Supplementary Table 10). Seemingly unrelated regression standard errors were consistent with the jackknife approach standard errors.

As the primary analysis, we reported shrinkage results using the liberal threshold (P < 1 × 10−5), with results using the genome-wide threshold (P < 5 × 10−8) reported as a sensitivity analysis. In the main text, we report the shrinkage estimates that reach nominal significance (P < 0.05). We presented shrinkage estimates in terms of % (multiplying by 100).

As a sensitivity analysis, we also presented study-level shrinkage estimates for height and educational attainment and tested for heterogeneity. These phenotypes were chosen because of previous evidence for shrinkage on these phenotypes and available data.

Heterogeneity of shrinkage across variants within a phenotype

We used results of the within-sibship and population-based meta-analyses to estimate whether shrinkage estimates were consistent across all variants within a phenotype, using an estimate of heterogeneity. As above, we only evaluated heterogeneity for height and educational attainment because of previous evidence and available data. For each variant we estimated the Wald ratio of the shrinkage estimate

$${{{s}}}_{{{k}}} = \frac{{\beta _{{{{\mathrm{P}}}},{{{k}}}}}}{{\beta _{{{{\mathrm{W}}}},{{{k}}}}}}$$

The heterogeneity estimate was obtained as

$${{{Q}}} = \mathop {\sum }\limits_{{{k}}}^{{{M}}} {{{w}}}_{{{k}}}^2\left( {{{{s}}}_{{{k}}}-{{{S}}}} \right)^2$$

where

$${{{w}}}_{{{k}}} = \sqrt {\frac{{{{{S}}}^2}}{{\sigma _{{{{\mathrm{W}}}},{{{k}}}}^2 + {{{S}}}^2\sigma _{{{S}}}^2}}}$$

Applying LDSC to within-sibship data

LDSC is a widely used method that can be applied to GWAS summary data to estimate heritability and genetic correlation20,23. The LDSC ratio, a function of the LDSC intercept unrelated to statistical power, is a measure of the proportion of association signal that is due to confounding. In this work, we apply LDSC to estimate SNP heritability and genetic correlation using the population and within-sibship GWAS data, so we investigated the LDSC intercept/ratio estimates from these data. Further detail is contained in the Supplementary Methods.

LDSC confounding estimates varied across the 25 phenotypes in the within-sibship model. Confounding estimates were modest for height (10%; 95% CI 6%, 14%) and BMI (9%; 2%, 16%), while the estimate for educational attainment was imprecise (35%; 12%, 57%). Across all phenotypes in the within-sibship data, the median confounding estimate was 21% (Q1–Q3: 10%, 28%), but stronger conclusions are limited by imprecise estimates (Supplementary Table 11 and Extended Data Fig. 8). The LDSC confounding estimates were higher using the population GWAS data (median 42%: Q1–Q3, 35%, 56%) than both the within-sibship model and previous studies (Supplementary Table 12). For example, the population model LDSC ratio estimates were higher for height (23%; 21%, 26%), BMI (22%; 19%, 25%) and educational attainment (41%; 37%, 45%).

The observed nonzero confounding in the within-sibship model was unexpected because of the intuition that the within-sibship GWAS models are unlikely to be confounded. The LDSC ratios in the population GWAS were also higher than previous studies. We followed up these findings by evaluating the effects of LD score mismatch and cryptic relatedness on the LDSC ratios.

Evaluation of LD score mismatch

A large proportion of samples in the meta-analysis were from UK-based studies such as UK Biobank and Generation Scotland, for which the LD scores, generated using 1000 Genomes project (phase 3) European samples (CEU, TSC, FIN, GBR), have been shown to fit reasonably well20. However, a large number of samples were from Scandinavian populations (HUNT study, FinnTwin), where LD mismatch leading to elevated LDSC intercept/ratios has been previously discussed20. We investigated this possibility using empirical and simulated data.

We investigated variation in LDSC ratios across populations by comparing ratios for height across well-powered individual studies (N > 5,000): UK Biobank, HUNT, the China Kadoorie Biobank (using default East Asian LD scores), Generation Scotland, DiscovEHR, Queensland Institute of Medical Research (QIMR) study and FinnTwin. We found some evidence of heterogeneity between studies: ratio estimates were higher in Scandinavian studies compared with UK-based studies (Extended Data Fig. 9). We also calculated within-sibship ratio estimates for BMI, SBP and educational attainment using UK Biobank summary data. UK Biobank estimates were largely consistent with zero confounding although confidence intervals were wide (Supplementary Table 13).

We also performed simulations to evaluate potential mismatch between the Norwegian HUNT study and the default LD scores, which were generated using 1000 Genomes data, finding evidence of LD score mismatch between the 1000 Genomes LD scores and HUNT. The simulation setup and results are detailed in the Supplementary Methods.

The combined findings from the empirical and simulated analyses suggest that LD score mismatch with the 1000 Genomes LD scores in the Norwegian HUNT study and other studies likely contributed to inflated LDSC ratios in both population and within-sibship GWAS models.

Cryptic relatedness

One source of inflation in GWAS associations is cryptic relatedness: nonindependence between close relatives in the study sample results which leads to inflated precision. In sibling GWAS models we clustered standard errors over sibships, but this clustering does not account for nonindependence between related sibships, for example, uncle/mother and two offspring. Inflated signal relating to cryptic relatedness may result in confounded signal, which is detected by the LD score intercept/ratio. In conventional population GWAS, either close relatives are removed or a mixed model is used to account for close relatives. We performed empirical and simulated analyses detailed in the Supplementary Methods to investigate the effect of cryptic relatedness on the population and within-sibship models.

The results suggest that the standard errors in the within-sibship model are not underestimated because of cryptic relatedness relating to common environmental effects shared between relatives. Thus, cryptic relatedness likely inflated LDSC ratios in the population models but not in the within-sibling data.

Within-sibship SNP heritability estimates

LDSC was used to generate SNP heritability estimates for 25 phenotypes using the LDSC harmonized (see above) meta-analysis summary data. The summary data were harmonized using the LDSC munge_sumstats.py function, and we used the precomputed European LD scores from 1000 Genomes Phase 3.

LDSC requires a sample size parameter N to estimate SNP heritability. For this parameter, we used the effective sample size for each meta-analysis phenotype, equivalent to the number of independent observations. This was estimated as follows using GWAS standard errors, minor allele frequencies and the phenotype standard deviations (after adjusting for covariates).

$${{{\mathrm{Effective}}}}\,{{{N}}} = \frac{1}{{{{{\mathrm{s.e.}}}}^2}}\frac{{{{{\mathrm{s.d.}}}}\_{{{\mathrm{Resid}}}}^2}}{{2 \times {{{\mathrm{MAF}}}} \times (1 - {{{\mathrm{MAF}}}})}}$$

s.e., GWAS model standard error; MAF, MAF of the variant; s.d._Resid, standard deviation of the regression residual.

Effective sample size was estimated for each individual study GWAS and each model (for example, UK Biobank population GWAS of height). To reduce noise from low-frequency variants, we restricted to variants with MAF between 0.1 and 0.4 (from 1000 Genomes EUR). At the meta-analysis stage, the effective sample size for each variant was calculated as the sum of sample sizes of studies in which the variant was present. Simulations evaluating the use of effective sample sizes are detailed in the Supplementary Methods.

In empirical analyses, we decided to focus on the differences between the population model (\({{{h}}}_{{{{\mathrm{Pop}}}}}^2\)) and within-sibship model \(({{{h}}}_{{{{\mathrm{WS}}}}}^2)\) SNP heritability estimates. If we assume that biases affect the estimates equally then the difference between the two estimates will be unbiased. We estimated the difference between the heritability estimates (\({{{h}}}_{{{{\mathrm{Diff}}}}}^2\)) using a difference-of-two-means test57 as below.

$${{{h}}}_{{{{\mathrm{Diff}}}}}^2 = {{{h}}}_{{{{\mathrm{Pop}}}}}^2 - {{{h}}}_{{{{\mathrm{WS}}}}}^2$$
$${{{\mathrm{s.e.}}}}\left( {{{{h}}}_{{{{\mathrm{Diff}}}}}^2} \right)\sim \sqrt {{{{\mathrm{s.e.}}}}({{{h}}}_{{{{\mathrm{Pop}}}}}^2)^2 + {{{\mathrm{s.e.}}}}({{{h}}}_{{{{\mathrm{WS}}}}}^2)^2 - 2{{{\mathrm{Cov}}}}({{{h}}}_{{{{\mathrm{Pop}}}}}^2,{{{h}}}_{{{{\mathrm{WS}}}}}^2)}$$

To estimate \({{{\mathrm{Cov}}}}({{{h}}}_{{{{\mathrm{Pop}}}}}^2,{{{h}}}_{{{{\mathrm{WS}}}}}^2)\), we computed the cross-GWAS LDSC intercept between the population and within-sibship GWAS data (for the same phenotype) which is an estimate of \({{{\mathrm{Cor}}}}({{{h}}}_{{{{\mathrm{Pop}}}}}^2,{{{h}}}_{{{{\mathrm{WS}}}}}^2)\). The estimates of this term were ~0.40 across phenotypes. We then calculated the covariance term as follows:

$${{{\mathrm{Cov}}}}({{{h}}}_{{{{\mathrm{Pop}}}}}^2,{{{h}}}_{{{{\mathrm{WS}}}}}^2) = {{{\mathrm{Cor}}}}({{{h}}}_{{{{\mathrm{Pop}}}}}^2,{{{h}}}_{{{{\mathrm{WS}}}}}^2) \times {{{\mathrm{s.e.}}}}({{{h}}}_{{{{\mathrm{Pop}}}}}^2) \times {{{\mathrm{s.e.}}}}({{{h}}}_{{{{\mathrm{WS}}}}}^2)$$

We used the difference Z score (that is, \(\frac{{{{{h}}}_{{{{\mathrm{Diff}}}}}^2}}{{{{{\mathrm{s.e.}}}}\left( {{{{h}}}_{{{{\mathrm{Diff}}}}}^2} \right)}}\)) to generate a P value for the difference between \({{{h}}}_{{{{\mathrm{Pop}}}}}^2\) and \({{{h}}}_{{{{\mathrm{WS}}}}}^2\). In the text, we report differences reaching nominal significance (difference P < 0.05).

We calculated the expected effect of shrinkage on LDSC SNP heritability estimates. LDSC heritability estimates (h2) are derived from the formulation below20:

$$\chi ^2\sim \frac{{{{{Nh}}}^2{{{l}}}_{{{j}}}}}{{{{M}}}} + {{{Na}}} + 1$$

where χ2 is the square of the GWAS Z score, N is the sample size, M is number of variants such that \(\frac{{{{{h}}}^2}}{{{{M}}}}\) is the average heritability for each variant, lj is the LD score of variant j and a is the effect of confounding biases.

Uniform shrinkage across the genome would lead to GWAS Z scores being multiplied by a factor (1 − k), where k is the shrinkage coefficient, and χ2 statistics being multiplied by (1 − k)2. As above, we have used effective sample size to account for differences in N between the population and within-sibship models. Therefore, assuming all other coefficients remain consistent, the expectation of \({{{h}}}_{{{{\mathrm{WS}}}}}^2\) can be written as a function of k and \({{{h}}}_{{{{\mathrm{Pop}}}}}^2\).

$${{{h}}}_{{{{\mathrm{Pop}}}}}^2 = {{{y}}}$$
$${{{h}}}_{{{{\mathrm{WS}}}}}^2 = (1 - {{{k}}})^2{{{y}}}$$

To evaluate the sensitivity of our results to assumptions of heritability models, we also estimated SNP heritability using SumHer21, which allows the use of different heritability models with regard to how local LD and allele frequencies affect the heritability contributions of individual SNPs. In SumHer analyses, we followed the same procedure as above for LDSC using effective sample sizes and estimating SNP heritability for all 25 phenotypes. We used the LDAK-Thin model with the precomputed tagging file over the BLD-LDAK model because of the limited power of our datasets (the BLD-LDAK model includes additional parameters so generates less precise estimates).

Within-sibship r g with educational attainment

We used LDSC to estimate rg between educational attainment and other phenotypes using both population and within-sibship data. LDSC requires nonzero heritability to generate meaningful rg estimates, so we restricted analyses to the 22 phenotypes with SNP heritability point estimates greater than zero in both population and within-sibship models (that is, omitted physical activity and ratio of forced expiratory volume (FEV1)/forced vital capacity (FEV1FVC)). We estimated only pairwise genetic correlations between educational attainment and all other phenotypes because of previous evidence that educational attainment is influenced by demographic and indirect genetic effects and, given the limited statistical power, to reduce the multiple testing burden. Estimates failed to converge for genetic correlation analyses involving age at first birth and age at menopause, so these phenotypes were not analyzed here. We estimated the difference between the population (rg,Pop) and within-sibship (rg,WS) estimates (rg,Diff) using a difference-of-two-means test57.

$${{{r}}}_{{{{\mathrm{g}}}},{{{\mathrm{Diff}}}}} = {{{r}}}_{{{{\mathrm{g}}}},{{{\mathrm{Pop}}}}} - {{{r}}}_{{{{\mathrm{g}}}},{{{\mathrm{WS}}}}}$$

We used the jackknife to estimate the standard error of the difference, \({{{\mathrm{s.e.}}}}({{{r}}}_{{{{\mathrm{g}}}},{{{\mathrm{Diff}}}}})\). After restricting to ~1.2 million Hapmap 3 variants present in the 1000 Genomes LD scores, we ordered variants by chromosome and base pair and separated variants into 100 blocks. We removed each block in turn and computed \({{{r}}}_{{{{\mathrm{g}}}},{{{\mathrm{Diff}}}}}\) using LDSC 100 times. We then calculated \({{{\mathrm{s.e.}}}}({{{r}}}_{{{{\mathrm{g}}}},{{{\mathrm{Diff}}}}})\) across the 100 iterations as follows:

$${{{\mathrm{s.e.}}}}({{{r}}}_{{{{\mathrm{g}}}},{{{\mathrm{Diff}}}}}) = \sqrt {\frac{{99}}{{100}}\mathop {\sum}\nolimits_1^{100} {({{{r}}}_{{{{\mathrm{g}}}},{{{\mathrm{Diff}}}}\,{{{k}}}} - \mu )^2} }$$

where

$$\mu = \frac{{\mathop {\sum }\nolimits_1^{100} {{{r}}}_{{{{\mathrm{g}}}},{{{\mathrm{Diff}}}},{{{k}}}}}}{{100}}$$

\({{{r}}}_{{{{\mathrm{g}}}},{{{\mathrm{Diff}}}},{{{\mathrm{k}}}}}\) is the rg estimate in the kth iteration and μ is the mean rg estimate across all 100 iterations.

We used the difference Z score (that is, \(\frac{{{{{r}}}_{{{{\mathrm{g}}}},{{{\mathrm{Diff}}}}}}}{{{{{\mathrm{s.e.}}}}\left( {{{{r}}}_{{{{\mathrm{g}}}},{{{\mathrm{Diff}}}}}} \right)}}\)) to generate a P value for heterogeneity between \({{{r}}}_{{{{\mathrm{g}}}},{{{\mathrm{Pop}}}}}\) and \({{{r}}}_{{{{\mathrm{g}}}},{{{\mathrm{WS}}}}}\). In the text, we report differences reaching nominal significance (heterogeneity P < 0.05).

WS-MR: effects of height and BMI

We performed MR analyses using the within-sibship meta-analysis GWAS data to estimate the effect of two exposures (height and BMI) on 23 outcome phenotypes. For the exposure instruments, we used 803 and 418 independent genetic variants for height and BMI, respectively. These variants were identified by LD clumping in PLINK (r2 < 0.001, physical distance threshold = 10,000 kb, P < 5 × 10−8) using the 1000 Genomes Phase 3 EUR reference panel54. We then performed an MR-IVW analysis using the within-sibship meta-analysis data to estimate the effect of the exposure on the outcome as

$$\beta _{{{{\mathrm{MR}}}}} = {\sum} {\frac{{\beta _{{{{\mathrm{Exp}}}}} \ast \beta _{{{{\mathrm{Out}}}}}}}{{(\sigma _{{{{\mathrm{Out}}}}})^2}}} /{\sum} {\frac{{(\beta _{{{{\mathrm{Exp}}}}})^2}}{{(\sigma _{{{{\mathrm{Out}}}}})^2}}}$$

where βExp is the association estimate from exposure GWAS, βOut is the association estimate from outcome GWAS and σOut is the standard error from outcome GWAS.

We also performed MR analyses using the population meta-analysis GWAS data for comparison. We estimated differences between population MR and WS-MR estimates using the difference-of-two-means test57:

$$\beta _{{{{\mathrm{MR}}}},{{{\mathrm{Diff}}}}} = \beta _{{{{\mathrm{MR}}}},{{{\mathrm{Pop}}}}} - \beta _{{{{\mathrm{MR}}}},{{{\mathrm{WS}}}}}$$

We used the jackknife to estimate the standard error of the difference, s.e.(βMR,Diff). With n genetic instruments, we removed each variant from the analysis in turn and then computed βMR,Diff, storing the estimate from the n iterations. We then calculated s.e.(βMR,Diff) as follows:

$${{{\mathrm{s.e.}}}}(\beta _{{{{\mathrm{MR}}}},{{{\mathrm{Diff}}}}}) = \sqrt {\frac{{{{{n}}} - 1}}{{{{n}}}}\mathop {\sum }\limits_1^{{{n}}} (\beta _{{{{\mathrm{MR}}}},{{{\mathrm{Diff}}}},{{{k}}}} - \mu )^2}$$

where

$$\mu = \frac{{\mathop {\sum }\nolimits_1^{{{n}}} \beta _{{{{\mathrm{MR}}}},{{{\mathrm{Diff}}}},{{{k}}}}}}{{{{n}}}}$$

n is the number of genetic variants used as instruments, βMR, Diff, k is the βMR, Diff estimate in the kth iteration and μ is the mean βMR, Diff estimate across all n iterations.

We used the difference Z score (that is, \(\frac{{\beta _{{{{\mathrm{MR}}}},{{{\mathrm{Diff}}}}}}}{{{{{\mathrm{s.e.}}}}\left( {\beta _{{{{\mathrm{MR}}}},{{{\mathrm{Diff}}}}}} \right)}}\)) to generate a P value for heterogeneity between βMR, Pop and βMR,WS. In the text, we report differences reaching nominal significance (heterogeneity P < 0.05).

Polygenic adaptation

Polygenic adaptation was estimated using similar methods to a previous publication28. Precomputed SDS were downloaded for UK10K data from https://web.stanford.edu/group/pritchardlab/. Genomic regions under strong recent selection (MHC chr6: 25,892,529–33,436,144; lactase chr2: 134,608,646–138,608,646) were removed and SDS were normalized within each 1% allele frequency bin.

SDS were merged with GWAS meta-analysis data for 25 phenotypes. Variants with low effective sample sizes (<50% of maximum) were removed for each phenotype. SDS were transformed to tSDS such that the reference allele was the phenotype-increasing allele.

Spearman’s rank test was used to estimate the correlation between tSDS and the absolute value of GWAS Z scores from the population and within-sibship models. Standard errors were estimated using the jackknife. The genome was ordered by chromosome and base pair and divided into 100 blocks. Correlations were estimated 100 times with each kth block removed in turn. The standard error of the correlation estimate, s.e.(Cor), was calculated as follows:

$${{{\mathrm{s.e.}}}}\left( {{{{\mathrm{Cor}}}}} \right) = \sqrt {\frac{{99}}{{100}}\mathop {\sum }\limits_1^{100} ({{{\mathrm{Cor}}}}_{{{k}}} - \mu )^2}$$

where

$$\mu = \frac{{\mathop {\sum }\nolimits_1^{100} {{{\mathrm{Cor}}}}_{{{k}}}}}{{100}}$$

Cork is Spearman’s rank correlation estimate in the kth iteration and μ is the mean correlation estimate across the 100 iterations.

Given previous concerns26,27, we performed several sensitivity analyses for the height analysis detailed in the Supplementary Methods.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.