Overview

Our initial aim was to estimate the effect of F ROH on 45 quantitative traits and to assess whether any of these effects differed significantly from zero. Previous work7,11 has shown that inbreeding coefficients are low in most human populations, and that very large samples are required to reliably estimate the genetic effects of inbreeding13. To maximise sample size, a collaborative consortium (ROHgen6) was established, and research groups administering cohorts with SNP chip genotyping were invited to participate. To ensure that all participants performed uniform and repeatable analyses, a semi-automated software pipeline was developed and executed locally by each research group. This software pipeline required cohorts to provide only quality-controlled genotypes (in plink binary format) and standardised phenotypes (in plain-text) and used standard software (R, PLINK12,32, KING33) to perform the analyses described below. Results from each cohort were returned to the central ROHgen analysts for meta-analysis.

During the initial meta-analysis, genotypes were released for >500,000 samples from the richly phenotyped UK Biobank (UKB)10. It was therefore decided to add a further 34 quantitative phenotypes and 21 binary traits to the ROHgen analysis. Many of these additional traits were unique to UKB, although 7 were also available in a subset of ROHgen cohorts willing to run additional analyses. In total, the effect of F ROH was tested on 100 traits and therefore experiment-wise significance was defined as 5 × 10−4 (=0.05/100).

Cohort recruitment

In total, 119 independent, genetic epidemiological study cohorts were contributed to ROHgen. Of these, 118 were studies of adults and contributed multiple phenotypes, while 1 was a study of children and contributed only birth weight. To minimise any potential confounding or bias caused by within-study heterogeneity, studies were split into single-ethnicity sub-cohorts wherever applicable. Each sub-cohort was required to use only one genotyping array and be of uniform ancestry and case-status. For example, if a study contained multiple distinct ethnicities, sub-cohorts of each ancestry were created and analysed separately. At minimum, ancestry was defined on a sub-continental scale (i.e. European, African, East Asian, South Asian, West Asian, Japanese, and Hispanic were always analysed separately) but more precise separation was used when deemed necessary, for example, in cohorts with large representation of Ashkenazi Jews. In case-control studies of disease, separate sub-cohorts were created for cases and controls and phenotypes associated with disease status were not analysed in the case cohort: for example, fasting plasma glucose was not analysed in Type 2 diabetes case cohorts. Occasionally, cohorts had been genotyped on different SNP genotyping microarrays and these were also separated into sub-cohorts. There was one exception (deCODE) to the single microarray rule, where the intersection between all arrays used exceeded 150,000 SNPs. In this cohort the genotype data from all arrays was merged since the correspondence between F ROH for the individual arrays and F ROH the intersection dataset was found to be very strong (\(\beta _{{\mathrm{merged}},{\mathrm{hap}}}\) = 0.98, r2 = 0.98; \(\beta _{{\mathrm{merged}},{\mathrm{omni}}}\) = 0.97, r2=0.97). Dividing studies using these criteria yielded 234 sub-cohorts. Details of phenotypes contributed by each cohort are available in Supplementary Data 4.

Ethical approval

Data from 119 independent genetic epidemiology studies were included. All subjects gave written informed consent for broad-ranging health and genetic research and all studies were approved by the relevant research ethics committees or boards. PubMed references are given for each study in Supplementary Data 2.

Genotyping

All samples were genotyped on high-density (minimum 250,000 markers), genome-wide SNP microarrays supplied by Illumina or Affymetrix. Genotyping arrays with highly variable genomic coverage (such as Exome chip, Metabochip, or Immunochip) were judged unsuitable for the ROH calling algorithm and were not permitted. Imputed genotypes were also not permitted; only called genotypes in PLINK binary format were accepted. Each study applied their own GWAS quality controls before additional checks were made in the common analysis pipeline: SNPs with >3% missingness or MAF <5% were removed, as were individuals with >3% missing data. Only autosomal genotypes were used for the analyses reported here. Additional, cohort-specific, genotyping information is available in Supplementary Data 2.

Phenotyping

In total, results are reported for 79 quantitative traits and 21 binary traits. These traits were chosen to represent different domains of health and reproductive success, with consideration given to presumed data availability. Many of these traits have been the subject of existing genome-wide association meta-analyses (GWAMA), and phenotype modelling, such as inclusion of relevant covariates, was copied from the relevant consortia (GIANT for anthropometry, EGG for birth weight, ICBP for blood pressures, MAGIC for glycaemic traits, CHARGE-Cognitive, -Inflammation & -Haemostasis working groups for cognitive function, CRP, fibrinogen, CHARGE-CKDgen for eGFR, CHARGE-ReproGen for ages at menarche and menopause, Blood Cell & HaemGen for haematology, GUGC for urate, RRgen, PRIMA, QRS & QT-IGC for electrocardiography, GLGC for classical lipids, CREAM for spherical equivalent refraction, Spirometa for lung function traits, and SSGAC for educational attainment and number of children ever born). Further information about individual phenotype modelling is available in Supplementary Note 1 and Supplementary Data 8.

ROH calling

Runs of homozygosity (ROH) of >1.5 Mb in length were identified using published methods6,11. In summary, SNPs with minor allele frequencies below 5% were removed, before continuous ROH SNPs were identified using PLINK with the following parameters: homozyg-window-snp 50; homozyg-snp 50; homozyg-kb 1500; homozyg-gap 1000; homozyg-density 50; homozyg-window-missing 5; homozyg-window-het 1. No linkage disequilibrium pruning was performed. These parameters have been previously shown to call ROH that correspond to autozygous segments in which all SNPs (including those not present on the chip) are homozygous-by-descent, not chance arrangements of independent homozygous SNPs, and inbreeding coefficient estimates calculated by this method (F ROH ) correlate well with pedigree-based estimates (F PED )11. Moreover, they have also been shown to be robust to array choice6.

Calculating estimators of F

For each sample, two estimates of the inbreeding coefficient (F) were calculated, F ROH and F SNP . We also calculated three additional measures of homozygosity: F ROH<5Mb , F ROH>5Mb and F SNP_outsideROH .

F ROH is the fraction of each genome in ROH >1.5 Mb. For example, in a sample for which PLINK had identified n ROH of length l i (in Mb), i ϵ {1..n}, then F ROH was then calculated as

$$\begin{array}{*{20}{c}} {F_{{\mathrm{ROH}}} = \frac{{\mathop {\sum }

olimits_{i = 1}^n l_i}}{{3Gb}}} \end{array},$$ (1)

where F ROH<5Mb and F ROH>5Mb are the genomic fractions in ROH of length >5 Mb, and in ROH of length <5 Mb (but >1.5 Mb), respectively, and the length of the autosomal genome is estimated at 3 gigabases (Gb). It follows from this definition that

$$\begin{array}{*{20}{c}} {F_{{\mathrm{ROH}}} = F_{{\mathrm{ROH}} > 5{\mathrm{Mb}}} + F_{{\mathrm{ROH}} < 5{\mathrm{Mb}}}} \end{array}.$$ (2)

Single-point inbreeding coefficients can also be estimated from individual SNP homozygosity without any reference to a genetic map. For comparison with F ROH , a method of moments estimate of inbreeding coefficient was calculated34, referred to here as F SNP , and implemented in PLINK by the command–het.

$$\begin{array}{*{20}{c}} {F_{{\mathrm{SNP}}} = \frac{{O\left( {{\mathrm{{HOM}}}} \right) -E\left( {{\mathrm{{HOM}}}} \right)}}{{N - E\left( {{\mathrm{{HOM}}}} \right)}}} \end{array},$$ (3)

where O(HOM) is the observed number of homozygous SNPs, E(HOM) is the expected number of homozygous SNPs, i.e. \(\mathop {\sum}

olimits_{i = 1}^N {\left( {1 - 2p_iq_i} \right)}\), and N is the total number of non-missing genotyped SNPs.

F ROH and F SNP are strongly correlated, especially in cohorts with significant inbreeding, since both are estimates of F. To clarify the conditional effects of F ROH and F SNP , an additional measure of homozygosity,F SNPoutsideROH , was calculated to describe the SNP homozygosity observed outside ROH.

$$\begin{array}{*{20}{c}} {F_{{\mathrm{SNP}}_{{\mathrm{outsideROH}}}} = \frac{{O\prime \left( {\mathrm{{HOM}}} \right) - E\prime \left( {\mathrm{{HOM}}} \right)}}{{N\prime - E\prime \left( {\mathrm{{HOM}}} \right)}}} \end{array},$$ (4)

where

$$\begin{array}{*{20}{c}} {O\prime \left( {{\mathrm{{HOM}}}} \right) = O\left( {{\mathrm{{HOM}}}} \right) - N_{{\mathrm{SNP}}\_{\mathrm{ROH}}}} \end{array},$$ (5)

$$\begin{array}{*{20}{c}} {E\prime \left( {\mathrm{{HOM}}} \right) = \left( {\frac{{N - N_{{\mathrm{ROH}}}}}{N}} \right) \ast E\left( {\mathrm{{HOM}}} \right)} \end{array},$$ (6)

$$N\prime = N - N_{{\mathrm{ROH}}}$$ (7)

And N SNP_ROH is the number of homozygous SNPs found in ROH. Note that:

$${{F_{\mathrm{{SNP}}}}}_{{\mathrm{outsideROH}}} \approx F_{{\mathrm{SNP}}} - F_{{\mathrm{ROH}}}$$ (8)

A further single point estimator of the inbreeding coefficient, described by Yang et al.20 as \(\widehat F^{{\mathrm{{III}}}}\), is implemented in PLINK by the parameter –ibc (Fhat3) and was also calculated for all samples.

$$F_{{\mathrm{GRM}}} = \widehat F^{\mathrm{{III}}} = \frac{1}{N}\mathop {\sum}

olimits_{i = 1}^N {\frac{{\left( {x_i^2 - \left( {1 + 2p_i} \right)x_i + 2p_i^2} \right)}}{{2p_i\left( {1 - p_i} \right)}}},$$ (9)

where N is the number of SNPs, pi is the reference allele frequency of the ith SNP in the sample population and x i is the number of copies of the reference allele.

Effect size estimates for quantitative traits

In each cohort of n samples, for each of the quantitative traits measured in that cohort, trait values were modelled by

$$\begin{array}{*{20}{c}} {y = \beta _{F_{{\mathrm{ROH}}}} \ast {\mathbf{F}}_{{\mathbf{ROH}}} + {\mathbf{{Xb}}} +{\boldsymbol{\varepsilon}}} \end{array},$$ (10)

where y is a vector (n × 1) of measured trait values, \(\beta_{{F}_{\rm{ROH}}}\) is the unknown scalar effect of F ROH on the trait, F ROH is a known vector (n × 1) of individual F ROH , b is a vector (m × 1) of unknown fixed covariate effects (including a mean, μ), X in a known design matrix (n × m) for the fixed effects, and ε is an unknown vector (n × 1) of residuals.

The m fixed covariates included in each model were chosen with reference to the leading GWAMA consortium for that trait and are detailed in Supplementary Data 8. For all traits, these covariates included: age (and/or year of birth), sex, and at least the first 10 principal components of the genomic relatedness matrix (GRM). Where necessary, additional adjustments were made for study site, medications, and other relevant covariates (Supplementary Data 3).

For reasons of computational efficiency, it was decided to solve Eq. (10) in two steps. In the first step, the trait (y) was regressed on all fixed covariates to obtain the maximum likelihood solution of the model:

$$\begin{array}{*{20}{c}} {\mathbf{y}} = {\mathbf{{Xb}}} + \boldsymbol{\varepsilon}\prime \end{array}.$$ (11)

All subsequent analyses were performed using the vector of trait residuals ε′, which may be considered as the trait values corrected for all known covariates.

In cohorts with a high degree of relatedness, mixed-modelling was used to correct for family structure, although, because ROH are not narrow-sense heritable, this was considered less essential than in Genome-Wide Association Studies. Equation (11) becomes

$$\mathbf{y} = {\mathbf{{Xb}}} + {\mathbf{u}} + \boldsymbol{\varepsilon}\prime,$$ (12)

where u is an unknown vector (n × 1) of polygenic effects with multivariate normal distribution of mean 0 and covariance matrix σ g 2A, where A is the genomic relationship matrix (GRM). In these related cohorts, a GRM was calculated using PLINK v1.9 and Grammar+ residuals of Eq. (12) were estimated using GenABEL35. These Grammar+ residuals (ε′) were used in subsequent analyses.

To estimate \(\beta_{{F}_{\rm{ROH}}}\) for each trait, trait residuals were regressed on F ROH to obtain the maximum likelihood (ML) solution of the model

$${\boldsymbol{\varepsilon }}\prime = \mu + \beta _{F_{\mathrm{ROH}}} \ast {\mathbf{F}}_{{\mathbf{ROH}}} + {\boldsymbol{\varepsilon}}.$$ (13a)

The sex-specific estimates of \(\beta_{{F}_{\rm{ROH}}}\) (Supplementary Data 12) were obtained from Eq. (13) applied to the relevant sex.

For all traits, a corresponding estimates of \(\beta_{{F}_{\rm{SNP}}}\) and \(\beta_{\rm{F}_{\rm{GRM}}}\) were obtained from the models

$${\boldsymbol{\varepsilon }}\prime = \mu + \beta _{F_{\mathrm{SNP}}} \ast {\mathbf{F}}_{{\mathbf{SNP}}} + {\boldsymbol{\varepsilon}},$$ (13b)

$${\boldsymbol{\varepsilon }}\prime = \mu + \beta _{F_{\mathrm{GRM}}} \ast {\mathbf{F}}_{{\mathbf{GRM}}} + {\boldsymbol{\varepsilon}}$$ (14)

and the effects of different ROH lengths and of SNP homozygosity (Fig. 4b) were obtained from the model

$$\begin{array}{l}{\boldsymbol{\varepsilon }}\prime = \mu + \left( {\beta _1 \ast {\mathbf{F}}_{{\mathbf{SNP}}_{\mathrm{outsideROH}}}} \right) + \left( {\beta _2 \ast {\mathbf{F}}_{{\mathbf{ROH}} < 5{\mathbf{Mb}}}} \right)\\ {\hskip -62pt}+ \left( {\beta _3 \ast {\mathbf{F}}_{{\mathbf{ROH}} > 5{\mathbf{Mb}}}} \right) + {\boldsymbol{\varepsilon }} \end{array}.$$ (15)

Effect size estimates for binary traits

Binary traits were analysed by two methods. The primary estimates of \(\beta_{{F}_{\rm{ROH}}}\) (Fig. 3b and Supplementary Data 10) were obtained from full logistic models:

$$\begin{array}{*{20}{c}} {g\left( {E\left[ {\mathbf{y}} \right]} \right) ={\mathbf{{Xb}}}} \end{array},$$ (16)

where g() is the link function (logit), and where F ROH and all applicable covariates (Supplementary Datas 3, 8) were fitted simultaneously. Mixed modelling for family structure was not attempted in the logistic models since an accepted method was not apparent.

For all subsequent results, y was scaled by \(1/\sigma _y^2\) and analysed by linear models, as for quantitative traits, including mixed-modelling where appropriate for family studies. This method of estimating binary traits with simple linear models gives asymptotically unbiased estimates of \(\beta_{{F}_{\rm{ROH}}}\) and se(\(\beta_{{F}_{\rm{ROH}}}\)) on the ln(Odds-Ratio) scale36. For all significant binary traits, a comparison of \(\widehat \beta _{F_{\mathrm{ROH}}}\) from the full model with \(\widehat \beta _{F_{\mathrm{ROH}}}\) from the linear model approximation is presented in Supplementary Fig. 8.

To give \(\widehat \beta _{F_{\mathrm{ROH}}}\) a more tangible interpretation, effect estimates are frequently quoted in the text as β 0.0625 , i.e. the estimated effect in the offspring of first cousins, where 6.25% of the genome is expected to be autozygous.

Religiosity and educational attainment as additional covariates

To assess the importance of potential social confounders, proxy measures of socio-economic status and religiosity were separately included in Eq. (13) as additional covariates. The modified effect estimates (\(\widehat {\beta \prime }_{F_{\mathrm{ROH}}}\)) were tested for significance (Supplementary Data 20) and compared to the uncorrected estimates (\(\beta_{{F}_{\rm{ROH}}}\)) (Supplementary Fig. 10a, b).

Since Educational Attainment (EA) was measured in many cohorts, this was chosen as the most suitable proxy for socio-economic status. However, since F ROH is known to affect EA directly6 any change in \(\beta_{{F}_{\rm{ROH}}}\) when conditioning on EA cannot be assumed to be entirely due to environmental confounding.

The analysis of religiosity was only carried out in UKB, where a rough proxy was available. Although no direct questions about religious beliefs were included, participants were asked about their leisure activities. In response to the question Which of the following do you attend once a week or more often? (You can select more than one), 15.6% of UKB participants selected Religious Group from one of the seven options offered. In the models described, religiosity was coded as 1 for those who selected Religious Group and 0 for those who did not. Although this is likely to be an imperfect measure of actual religious belief it is currently the best available in a large dataset.

Assortative mating

Humans are known to mate assortatively for a number of traits including height and cognition37, and so we sought to investigate if this could influence our results, for example, by the trait extremes being more genetically similar and thus the offspring more homozygous. We see no evidence for an effect of assortative mating on autozygosity, however. Firstly, a polygenic risk score for height (see Supplementary Note 1), which explains 18.7% of the phenotypic variance in height, was not associated with F ROH (p = 0.77; Supplementary Fig. 5). Secondly, linear relationships between traits and autozygosity extend out to very high F ROH individuals (Supplementary Figs. 9a–f). Samples in the highest F ROH group are offspring of genetically similar parents, very likely first or second degree relatives and, for example, the height of these samples is on average 3.4 cm [95% CI 2.5–4.3] shorter than the population mean. Assortative mating would suggest this height deficit has been inherited from genetically shorter parents, but this would require an implausibly strong relationship between short stature and a propensity to marry a very close relative. Thirdly, the sex-specific effects we observe could only be explained by assortative mating if the additive heritability of these traits also differed by gender.

Average trait values in groups of similar F ROH

In each cohort individuals were allocated to one of ten groups of similar F ROH . The bounds of these groups were the same for all cohorts, specifically {0, 0.002, 0.0041, 0.0067, 0.0108, 0.0186, 0.0333, 0.06, 0.10, 0.18 and 1.0}. Within each group the mean trait residual (ε′) and mean F ROH were calculated, along with their associated standard errors. Within each cohort the expectation of ε′ is zero at the mean F ROH , however as mean F ROH varies between cohorts (Fig. 2, Supplementary Data 5) it was necessary to express ε′ relative to a common F ROH before meta-analysis. Hence, for this analysis only, the trait residuals (ε′) were expressed relative to the F ROH = 0 intercept, i.e. by subtracting μ from Eq. (13).

Effect of F ROH within adoptees

We compared \(\beta_{{{F}_{\rm{ROH}\_{\rm{ADOPTEE}}}}}\) to cross-cohort \(\beta_{{F}_{\rm{ROH}}}\), not that from UKB alone, as we consider the latter to be a noisy estimate of the former; estimates in UKB are consistent with those from meta-analysis.

Effect of F ROH within full-sibling families

In a subset of cohorts, with substantial numbers of related individuals, further analyses were performed to investigate the effect of F ROH within full-sibling families. In each of these cohorts, all second-degree, or closer, relatives were identified using KING (parameters:–related–degree 2). Full-siblings were then selected as relative pairs with genomic kinship >0.175 and IBS0 >0.001. This definition includes monozygotic twins, who were intentionally considered as part of full-sibling families. Although monozygotic twins are expected to have identical F ROH , they may not have identical trait values, and including additional trait measurements decreases the sampling error of the within-family variance estimate, hence increasing statistical power. Dizygotic twins were also included.

For each individual (j) with identified siblings, the values of F ROH and trait residual (ε′) were calculated relative to their family mean (and called F j ROH_wSibs and ε j wSibs, respectively), i.e. for individual j with n full-siblings S k where k ϵ {1..n}

$$\begin{array}{*{20}{c}} {F_j^{{\mathrm{ROHwSibs}}} = F_j^{{\mathrm{ROH}}} - \frac{1}{{\left( {n + 1} \right)}}\mathop {\sum}

olimits_{i{\it{\epsilon }}\left\{ {j,S_k} \right\}} {F_i^{{\mathrm{ROH}}}} } \end{array},$$ (17)

$$\begin{array}{*{20}{c}} {\varepsilon _j^{{\mathrm{wSibs}}} = \varepsilon _j^\prime - \frac{1}{{\left( {n + 1} \right)}}\mathop {\sum}

olimits_{i{\it{\epsilon }}\left\{ {j,S_k} \right\}} {\varepsilon \prime _i} } \end{array}.$$ (18)

The effect of F ROH within-full-siblings (\(\beta_{{F}_{\rm{ROH}}\_{\rm{wSibes}}}\)) was estimated by linear regression of εwSibs on FROH_wSibs.

Importantly, the variation of F ROH within full-siblings is entirely caused by differences in Mendelian segregation, and is therefore completely independent of all possible confounders. Hence, the effect estimates obtained by this method are estimates of the genetic effects of F ROH , unbiased by any possible confounder. Since confounding by social factors is a major concern in this field, methods that can definitively exclude this possibility are of critical importance.

Between-cohort meta-analysis

As is typical in genome-wide association meta-analyses (GWAMA), genetic effects were estimated within single-ethnicity sub-cohorts, and meta-analysis of the within-cohort effect sizes was used to combine results38. This established method eliminates any potential confounding caused by between-cohort associations between F ROH and traits.

Each cohort returned estimates and standard errors of: \(\beta_{{F}_{\rm{ROH}}}\), \(\beta_{{F}_{\rm{SNP}}},\beta_{{F}_{\rm{ROH} > {\rm{Mb}}}},\beta_{{F}_{\rm{ROH} < {\rm{Mb}}}},\beta_{{F}\_{\rm{outsideROH}}}, \beta_{{F}_{\rm{ROH}\_{\rm{wSibs}}}}\), as well as trait means (\(\overline {\varepsilon \prime }\)) and standard errors within each of 10 F ROH bins. The between-cohort mean of each of these 16 estimates was then determined by fixed-effect, inverse-variance meta-analysis using the R package metafor39. Results shown in Figs. 3–5 are meta-analysed averages of the within-cohort effects.

The meta-analysis was also run for various subsets of cohorts, stratified by ancestry as defined in Supplementary Data 18. Meta-analysis estimates from these groupings are shown in Supplementary Fig. 1.

Median and 95% CI of a ratio

In the analyses of adoptees (Fig. 5c), siblings (Fig. 5d) and potential confounders (Supplementary Figs. 10a, b) we wished to compare the effect estimates (\(\beta_{{F}_{\rm{ROH}}}\)) from two different methods across a wide range of traits. The units of \(\beta_{{F}_{\rm{ROH}}}\) differ by trait so, to allow comparison across all traits, the unitless ratio of effect size estimates was calculated (for example \(\beta_{{F}_{\rm{ROH}\_{\rm{wSibs}}}}\): \(\beta_{{F}_{\rm{ROH}}}\)). Figure 5c, d and Supplementary Figs. 10a, b show the medians and 95% CI of these ratios. These were determined empirically by bootstrap since, although formulae exist for the mean and standard error of a ratio40, the assumption of normality is violated when \(\beta_{{F}_{\rm{ROH}}}\)/se(\(\beta_{{F}_{\rm{ROH}}}\)) is not large.

Genetic correlations in UK Biobank

Genetic correlations were calculated using LD-Score Regression41, implemented in LDSC v1.0.0 (https://github.com/bulik/ldsc). Summary statistics were parsed using default parameters in the LDSC ‘munge_sumstats.py’ script, extracting only variants present in the HapMap 3 reference panel.

Accuracy of F ROH measures of inbreeding effects

A recent paper suggested that ROH may overestimate inbreeding effects by as much as 162%42; however, this could only be the case if F ROH underestimates excess homozygosity at the causal loci by at least 162%. We do not believe this to be the case since the maximum F ROH measured in many cohorts is around 0.25 (the expectation in the offspring off first-degree relatives), and the effect size estimates from these samples are consistent with the overall estimates (Fig. 5c, d and Supplementary Fig. 9a–f). We note that Yengo et al. applied the ROH calling parameters used here to imputed data. These parameters have been validated for called genotype data6 but not, to our knowledge, for the higher SNP density and error rate of imputed data (see also Supplementary Note 4). The simple method for detecting ROH used here was well suited to our study, since it could be easily implemented on over one million samples, and most of the variation in F ROH is caused by easily-identified long ROH.43,44,45

Reporting summary

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