Rare genetic variants explain missing heritability in smoking

Common genetic variants explain less variation in complex phenotypes than inferred from family-based studies, and there is a debate on the source of this ‘missing heritability’. We investigated the contribution of rare genetic variants to tobacco use with whole-genome sequences from up to 26,257 unrelated individuals of European ancestries and 11,743 individuals of African ancestries. Across four smoking traits, single-nucleotide-polymorphism-based heritability (hSNP2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h^2_{\mathrm{SNP}}$$\end{document}) was estimated from 0.13 to 0.28 (s.e., 0.10–0.13) in European ancestries, with 35–74% of it attributable to rare variants with minor allele frequencies between 0.01% and 1%. These heritability estimates are 1.5–4 times higher than past estimates based on common variants alone and accounted for 60% to 100% of our pedigree-based estimates of narrow-sense heritability (hped2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h^2_{\mathrm{ped}}$$\end{document}, 0.18–0.34). In the African ancestry samples, hSNP2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h^2_{\mathrm{SNP}}$$\end{document} was estimated from 0.03 to 0.33 (s.e., 0.09–0.14) across the four smoking traits. These results suggest that rare variants are important contributors to the heritability of smoking. The team of authors led by Seon-Kyeong Jang use whole-genome sequencing data and show that rare genetic variants explain much of the ‘missing heritability’ in smoking behaviours. These results help address a long-standing mystery in behavioural genetics.

use are strong indicators of addiction to nicotine. For example, the number of cigarettes smoked per day is genetically highly correlated (r = 0.95) with dependence on nicotine 13 and other commonly used substances 14,15 . Tobacco use is also considered a part of externalizing dimensions of personality and psychopathology along with other substance use and hyperactive behaviours 16 . The heritability of smoking behaviours has been estimated at approximately 50% (s.e., 5%) 17 in twin studies, comparable to many other complex behavioural traits.
At the same time, estimates of tobacco use heritability from GWAS of single nucleotide polymorphisms (SNPs) have routinely found much lower SNP-based heritability (h 2 SNP ) estimates 18,19 . Such analyses to date have been based on common variants with MAF > 1% from GWAS of imputed microarrays. In a recent GWAS of tobacco use in up to 1.1 million individuals, Liu et al. 18 reported h 2 SNP estimates ranging between 5% and 11%, with smoking initiation and age of smoking initiation showing the highest and the lowest estimates, respectively. Even more recently, Evans et al. 20 reported h 2 SNP estimates of 5-18% for smoking traits in UK Biobank imputed genotypes of up to 323,068 individuals. In this latter study, the estimated contribution of rare variants to the heritability was minimal, probably due to poor imputation of rare variants with MAF < 1%. Similar to results for other complex traits 17,21 , some but far from all of the twin-based heritability of smoking behaviours can be attributed to common variants obtained through the imputation of microarray genotypes.
There is an extensive literature regarding the possible contributors to missing heritability, including inflated family-based heritability estimates 22 , increased phenotypic heterogeneity in GWAS compared with family-based studies 23 , epistasis 24 and genetic variants not in linkage disequilibrium (LD) with common variants, including structural 21 and rare variants 25 . Rare genetic variants are one compelling explanation for the missing heritability of fitness-related traits, as one expects negative selection to force deleterious alleles to low frequencies 26,27 . Common variants, in contrast, are expected to explain most trait variance under a neutral model where most mutations have little selective effect 27,28 . The contribution of rare variants can inform competing explanations for the missing heritability as well as competing population genetics models, but this requires large samples of extremely well-imputed microarrays or whole-genome sequencing (WGS).
To date, SNP-based heritability estimates of complex behaviours have been based on a few million common variants. With imputed variants, the quality of imputation depends on the reference panel used 29,30 , and even the best imputation strategies perform poorly for variants with MAF < 1% in population-based samples 31,32 of unrelated individuals. With the advent of relatively affordable deep WGS, it is now possible to directly genotype variants of lower frequencies in larger samples. While genetic association studies may be underpowered to detect an association between a given single rare variant and a complex trait 33 , recent extensions of mixed-effects models allow one to estimate a random effect representing the aggregate contribution of rare variants to phenotypic variance over and above common variation 31,33,34 . To date, a small number of recent WGS studies have reported evidence that rare variants account for a substantial part of the heritability of anthropometric, transcriptomic and medical phenotypes [35][36][37][38] (but see also refs. 39,40 for counterexamples). Notably, rare variants, especially those in regions of low LD, capture a large part of the missing heritability for height and body mass index, albeit with large standard errors 41 .
Here we used deep WGS (mean depth 30×) from the Trans-Omics for Precision Medicine (TOPMed) programme to estimate the heritability explained by variants as rare as MAF of 1 in 10,000 for individuals of European ancestry (up to 26,257 individuals) and 1 in 1,000 for African admixed individuals (up to 11,743 individuals). We performed extensive sensitivity analyses to test the influence of rare-variant population structure along with various analytic parameters on our SNP-based heritability estimates for tobacco use. The study protocol was approved by the Institutional Review Board at the University of Minnesota and the TOPMed consortium. Informed consent was obtained originally by the participating studies.

Results
Heritability estimates. We estimated the heritability of four smoking phenotypes in samples of European ancestries: (1) age of smoking initiation (AgeSmk, N = 14,709), which measures the age at which an individual started regularly smoking; (2) cigarettes smoked per day (CigDay, N = 15,384), an index of heaviness of smoking and the average number of cigarettes smoked per day as a current or former smoker, grouped into five bins, with higher numbers indicating greater use; (3) smoking cessation (SmkCes, N = 17,827), a binary variable indicating whether a person is a current smoker; and (4) smoking initiation (SmkInit, N = 26,257), a binary variable indicating people who ever smoked regularly (that is, over 100 cigarettes) in their life. We stratified SNPs and short indels by their MAF and LD so that we had a total of six bins (MAF 5-50% and high LD, MAF 5-50% and low LD, MAF 1-5% and high LD, MAF 1-5% and low LD, MAF 0.1-1%, and MAF 0.01-0.1%). Variants were assigned to the high-LD group when their LD scores were higher than the median of the variants in the same MAF bin (Table 1). After we adjusted for the fixed effects of ten common variant principal components (PCs), ten rare variant PCs, age, sex and a random effect of study, the total SNP-based heritability (ĥ 2 SNP ) from variants of all MAFs was 0.226 (s.e., 0.116) for AgeSmk, 0.134 (0.095) for CigDay, 0.283 (0.127) for SmkCes and 0.225 (0.096) for SmkInit (Fig. 1). Heritability estimates from the four common-variant bins (that is, MAF 1-50%)-which include two bins for variants with high LD scores and two bins for variants with low LD scores-were summed to compute the total heritability attributable to common variants (ĥ 2 common ). Rare-variant heritability (ĥ 2 rare ) was computed likewise from the two MAF 0.01-1% bins, which were not grouped by LD score because most rare variants have low LD (Methods). Across the four smoking phenotypes, rare variants accounted for 35-74% of the  AgeSmk  14,709  3,092,517  3,092,534  1,342,734  1,342,736  5,392,813  28,280,118   CigDay  15,384  3,092,240  3,092,269  1,341,881  1,341,882  5,435,505  20,415,037   SmkCes  17,827  3,092,593  3,092,594  1,340,764  1,340,771  5,413,019  23,483,166   SmkInit  26,257  3,092,454  3,092,475  1,341,068  1,341,071  5,395,579  21,108,704 total heritability. The majority of the heritability attributable to rare variants was in the rarest frequency bin (MAF 0.01-0.1%), except for CigDay, where most ĥ 2 rare was attributed to variants with MAF 0.1-1% (Supplementary Table 1). When we further partitioned rare variants according to their functional impact, most of the heritability was localized to lower-functional-impact (non-protein-altering) variants (Supplementary Table 2).

Sensitivity analyses.
We used a variety of sensitivity analyses to evaluate the influence of residual population structure and cryptic relatedness. First, we evaluated the influence of increasing the number of PCs from 20 to 40 and 100 (half of which derived from common and rare variants) on ĥ 2 SNP to test whether greater numbers of PCs account for potential residual population structure.
The results for AgeSmk, SmkCes and SmkInit were unaffected (  Table 3). We observed model convergence issues when correcting CigDay for a large number of PCs, resulting in negative ĥ 2 SNP (for example, −0.12) attributable to the rarest MAF bin, which we deemed well outside the plausible range.
Second, we tested whether ĥ 2 SNP , especially that of rare variants, was driven by individuals sharing long identity-by-descent (IBD) segments. Using FastSMC 42 , we detected about 7.5 million IBD segments longer than 2 cM and shared among at least two individuals (mean length, 2.50 cM; s.d., 0.50 cM). After we removed up to 178 individuals whose shared IBD segments stretched longer than 2.5% of the total genome length, ĥ 2 SNP did not show meaningful differences, with changes in ĥ 2 SNP ranging from −0.025 to 0.031 (Fig. 2  and Supplementary Table 3). Third, to assess the influence of distant The error bars represent standard errors. 'rare' is the sum of the mAF 0.1-1% and mAF 0.01-0.1% bins. 'Common' is the sum of the other mAF bins. 'Total' is the sum of 'rare' and 'Common'. All estimates were adjusted for demographic variables and 20 PCs (half of them from rare variants) as fixed effects along with a random effect of cohort, except for CigDay, which was adjusted for five common PCs to allow model convergence.
Fifth, we obtained the approximate null distribution of ĥ 2 SNP under the assumption that population structure explains the observed SNP-based heritability (that is, the genotypes of ancestrally close individuals are exchangeable). This was done by estimating ĥ 2 SNP in samples randomly permuted within individuals close in ancestry as determined by a weighted sum of the top ten PCs with weights proportional to their corresponding eigenvalues 44 . We computed ĥ 2 null from 100 permuted trials and found limited evidence of inflation, as mean heritability estimates from the permuted trials were near zero across all bins of the four smoking phenotypes (Supplementary Table 4).
Finally, we estimated ĥ 2 SNP in samples with decreasing control of population stratification by ascertaining samples with less family structure on heritability estimates, we adjusted for the top 20 PCs from an IBD-based genetic relatedness in addition to the 20 PCs from rare and common variants and found that ĥ 2 SNP was not affected, with the largest change in ĥ 2 SNP being 0.013 ( Fig. 2 and Supplementary Table 3).
Fourth, previous studies indicated that rare variants are more likely to show geographical specificity 43 . While such geographic specificity was at least partially captured by our inclusion of study as a random effect, we went further by including a random effect of recruitment site (which is nested within a given study) in the model. After adjustment of site, ĥ 2 rare of SmkCes was about 50% of the original estimate, while ĥ 2 SNP of the other smoking traits showed little change ( Fig. 2 Table 3). Using related individuals in TOPMed across cohorts (N up to 21,546), narrow-sense heritabilities for AgeSmk, CigDay, SmkCes and SmkInit were estimated at 0.343 (0.075), 0.175 (0.066), 0.179 (0.067) and 0.288 (0.057), respectively. To check whether ĥ 2 ped was downwardly estimated due to a relatively large proportion of unrelated pairs, we also estimated pedigree-based heritability in the Framingham Heart Study (FHS) (N up to 3,024), which is composed of families as part of the study design. Similar estimates (AgeSmk, 0.225; CigDay, 0.244; SmkCes, 0.195; SmkInit, 0.299) were obtained in the FHS as in the larger analysis of all available families, albeit with larger standard errors. The estimates from the pedigree analysis are presented in Fig. 3 and Supplementary Table 5.

and Supplementary
Heritability estimates in African ancestries. We explored ĥ 2 SNP attributable to variants with MAF down to 0.1-1% for SmkInit (N = 11,743) and common variants for the other smoking phenotypes (N = 6,796-7,549) in samples of African ancestries. After we adjusted for 100 PCs, admixture proportions of five ancestries, local ancestry and the same demographic covariates as in European ancestry, common-variant-based SNP heritability was estimated as follows: 0.098 (0.091) for AgeSmk, 0.028 (0.113) for CigDay, 0.075 (0.120) for SmkCes and 0.085 (0.10) for SmkInit. Variants from MAF 0.1-1% additionally accounted for about 24% of the phenotypic variance of SmkInit, leading to total SNP-based heritability (combining both common and low-frequency variants) of SmkInit of 0.329 (0.144). Heritability estimates for individual bins are reported in Table 2 and Supplementary Table 6.

Discussion
Using up to 26,257 whole-genome sequences, we found that about 13-28% of the phenotypic variance of four smoking phenotypes in individuals of European ancestry could be attributed to genetic variants with minimum MAFs as low as 1 in 10,000. This is about 1.5-4.5 times larger than the previous SNP-based heritability estimates based on common variants alone.
This increase was largely driven by the inclusion of rare variants with MAF 0.1-0.01%, which accounted for 35-74% of the estimated SNP-based heritability of smoking phenotypes in samples of European ancestries. Smoking cessation showed the highest ĥ 2 rare (0.209; s.e., 0.123), followed by age of smoking initiation (0.153; s.e., 0.114), smoking initiation (0.079; s.e., 0.091) and cigarettes per day (0.054; s.e., 0.092). Different smoking phenotypes would in principle vary in their genetic architectures if they differentially relate to fitness or fitness-related traits. However, the standard errors associated with these estimates prevent us from drawing strong conclusions about their relative magnitudes.
The contributions of rare variants observed here seem consistent with the action of negative selection. Under negative selection, harmful mutations are maintained at low frequency in the population, whereas under the neutral model, common variants are expected to explain the majority (>90%) of the heritability 28 . The same selective pressures may apply to the heritability patterns found for common variants in this study. We found that most of the heritability from common variants was attributed to low-LD variants in the European ancestry samples. Consistent with this, Gazal et al. 45 reported that common variants with low LD in low-recombination-rate regions had larger per-SNP heritability than those with high LD. One possible explanation for this finding is that low-LD common variants are more likely to have arisen recently than high-LD variants. The low-LD variants have therefore had less time to be pushed to low frequency or fixation by selection pressure, comprising a major source of the heritability attributed to common variants 45 .
Our suggestion that smoking traits have been under negative selection might seem surprising. While nicotine itself has been pervasive in the environment for millions of years, tobacco in its stringent PC-based ancestry filtering thresholds (Methods). The largest change in ĥ 2 SNP was observed for SmkCes, where the estimates were up to 0.10 higher in samples with less stringent PC-based filtering than in the primary analysis, possibly reflecting sampling variations and residual population structure (Supplementary Fig. 1 and Supplementary Table 3).
In summary, in 36 comparisons (4 phenotypes × 9 conditions), ĥ 2 SNP showed relatively large changes (ĥ 2 SNP ± −0.10) only for SmkCes when we adjusted for geographical site or relaxed ancestry filtering. Our primary results presented in Fig. 1 are from the most stringent PC-filtering condition.
Pedigree-based heritability estimates. We removed closely related individuals when estimating the heritabilities described above. However, we were able to leverage thousands of related family members in TOPMed to independently estimate the total narrow-sense heritability of our four smoking phenotypes in pedigrees (ĥ 2 ped ). These heritability estimates were derived using a single GRM that included all individuals with at least one close relative in TOPMed. The relatedness of closely related individuals was estimated from common variants, which, when estimated for close family members, serves as a proxy for rare and common variant sharing. The relatedness of classically unrelated individuals (relatedness estimate (π) < 0.05) was set to be zero (same as K IBS>t in Zaitlen et al. 22 ). This procedure provides an upper bound on the narrow-sense heritability estimate, to which we can compare the estimates of ĥ 2 SNP described in the previous section and quantify any of the remaining missing heritability. cigarette form represents an evolutionarily novel environment for humans, and the selection against such traits is likely to be weak due to its recency. Nevertheless, genetic variants influencing smoking can affect multiple different traits more directly related to evolutionary fitness. For example, smoking is highly genetically correlated with age at onset of reproductive behaviours (especially age of smoking initiation 46 ) and with many social and health outcomes such as education and metabolic diseases. Selection pressure on these correlated traits can therefore influence the genetic architecture of smoking phenotypes 46 . More studies are needed to understand the mechanisms that maintain the contribution of rare genetic variants to the heritability of smoking phenotypes in the population.
We found that almost all the contribution from rare variants to ĥ 2 SNP was attributed to genetic variants classified as having relatively low expected protein-altering consequences. This may suggest an importance of rare regulatory genetic variation in the etiology of smoking behaviours. Despite limited contributions to heritability from protein-altering variants, it is still possible that they explain higher heritability individually 41 given the smaller number of such variants (about 0.2 million protein-altering versus 20 million non-protein-altering variants in MAF 0.01-0.1%) and the presumably higher selective pressure against them. Unfortunately, high standard errors and overall small heritabilities from the protein-altering variants bin prevented us from drawing strong conclusions from the present data.
Our SNP-based heritability estimates were largely robust to our sensitivity checks. Heritability changed little after we increased the number of PCs and tested finer-scale population structure via permutation. The estimates also largely remained unaffected by removing individuals that shared several long IBD segments or adjusting for the top 20 PCs from the IBD relatedness matrix. Controlling for site of recruitment led to about a 50% reduction in the rare-variant heritability of smoking cessation, while the other three smoking phenotypes were little affected. One explanation is that the genetic risk of smoking cessation is stratified geographically 47,48 , and rare variants are sensitive to capturing such ecological clustering. However, it is not clear why the same pattern is not observed for other closely related measures of smoking. Alternatively, for practical reasons, studies recruited samples from sites with systematic differences in current smoking rates or disease status (for example, clinical samples versus community controls), leading to the confounding of site and smoking status. We therefore consider ĥ 2 SNP of smoking cessation obtained after adjusting for geographical site as a conservative lower bound.
There is a long, ongoing debate on the source of missing heritability for virtually all complex traits. We found that current ĥ 2 SNP accounted for 60% to 100% of the pedigree-based, narrow-sense heritability estimates (ĥ 2 ped ) across the four phenotypes, closing the missing heritability gap for smoking phenotypes. Note that our ĥ 2 ped values are estimated from pedigrees from multiple cohorts, who also may span different generations (for example, parents-offspring) and share markedly different smoking environments. We consider our ĥ 2 ped as the most relevant benchmark by which to compare ĥ 2 SNP from our genome-based restricted maximum likelihood (GREML) results and past GWAS studies, which also have varying sources of sample heterogeneity, including gene-environment interaction 49 and varying heritability 50 . Twin studies, in contrast, usually consist of large, homogeneous cohorts with well-defined pedigrees sharing most of the demography (that is, monozygotic and dizygotic twins). For individuals of African ancestries, we found that common SNPs and indels (MAF > 1%) accounted for 3-15% of the phenotypic variance of smoking phenotypes, similar in magnitude to existing studies 18,20 . We observed an additional contribution of rare variants (MAF 0.1-1%) for smoking initiation, accounting for about 25% of smoking initiation liability. Note that estimating heritability in samples with complex population structure is in its infancy, and how recent admixture influences causal allele frequency spectra is not well understood. More research on this topic will be welcome.
Our findings should be interpreted in light of several limitations. First, although we used the largest WGS sample reported to date for heritability analysis, even larger datasets will be important. Larger sample sizes would allow greater precision in estimation and a more comprehensive assessment of the genetic architecture of these complex traits by finer partitioning by MAF and functional annotations. Second, even with the use of deep sequences, we did not fully reach the trait heritability as estimated using available pedigrees or twin heritability reported in the literature. There remain many explanations, including ultra-rare variants (MAF ≤ 0.01%), other types of genetic variations (for example, copy number variations) that are not in LD with the variants used here and still other sources of heritable variation. Cigarettes per day showed the lowest SNP-based heritability and often encountered model convergence issues, especially with a large number of covariates. Since it is a binned variable, its heritability may be more robustly estimated on a liability scale 51 .
Third, although we carefully assessed the potential influence of population structure in extensive sensitivity analyses, existing methods may not be entirely capturing the population stratification of rare variants 52,53 . More studies are needed to understand the extent and the best ways to handle it, although we are also sensitive to the possibility that such corrections may be counterproductive, removing meaningful heritable variance. Fourth, our pedigree-based heritability estimates may be inflated by shared environment, as we were unable to model genetic similarity and environmental similarity separately. The pedigree estimates therefore may best be interpreted as an upper bound on the SNP-based narrow-sense heritability. However, we expect the magnitude of inflation due to shared environment to be smaller than in typical twin studies given that our pedigree sample includes multiple classes of family relations who share much less of their environments than twins would.
Fifth, there is ample documentation of spousal similarity in smoking behaviours [54][55][56] , and one of the potential sources of this similarity, assortative mating, can artefactually inflate marker-based heritability estimates 57 . While recent studies have not found evidence for primary assortative mating on smoking traits in the genome 58 , more research is needed on this topic.
Sixth, genetic association signals of human behavioural traits can arise from multiple sources, not only from variants' downstream biological influences on a phenotype of interest but also from its genetically correlated traits (for example, education and income) along with indirect effects of parents' and siblings' genetics [59][60][61] . Our SNP-based heritability estimates include influences from multiple sources and could overlap with polygenic signals of socio-economic factors given their substantial genetic correlations with smoking (|r| = 0.26-0.55) 18 . With bigger sample sizes and studies of close families, future research can localize the origins of SNP-based heritability and examine the genetic architecture of shared and distinct genetic components of smoking and its correlated phenotypes (such as education and income) 62 .
Seventh, smoking phenotypes were measured by one or two questions and were limited to those commonly collected in biomedical studies such as those in TOPMed. This allows accumulations of large sample sizes across multiple independently collected samples, but each measure is noisy and may reflect idiosyncrasies of the participating studies. We attempted to account for the latter through including a random effect of study.
In conclusion, the present study expanded our understanding of the genetic architecture of tobacco use by showing that much of the missing heritability of smoking phenotypes can be explained by rare genetic variants. The use of WGS allows less biased estimation of SNP-based heritability by including virtually all potential causal SNPs down to a MAF of 0.01% in individuals of European ancestries and 0.1% in those of African ancestries, rather than relying on the degree to which imputed variants tag causal variants. The current study informs the genetic etiology of nicotine addiction and provides a benchmark for the future study of other complex behavioural traits.

Sample.
We considered individuals of European ancestries in TOPMed (freeze, 8; mean depth, >30) 34 measured for at least one of four smoking phenotypes for inclusion. Our data access application was approved by relevant cohorts in TOPMed. We determined European ancestries in two steps. First, we identified an initial ancestry-inclusive set by projecting TOPMed genotypes (N = 137,977) onto genetic PC axes from the 1000 Genomes project 8 (1000 G) and then used a k-nearest-neighbour method to assign the ancestry of TOPMed individuals with 1000 G as a reference set. More specifically, we used online augmentationdecomposition-transformation to calculate the PC scores of TOPMed individuals, which implements Procrustes transformation with an augmented dataset (that is, combining TOPMed and 1000 G reference genomes; https://github.com/ daviddaiweizhang/fraposa) 63 . Then, for a given TOPMed sample, we chose the top 20 reference individuals in 1000 G who were the closest in terms of the Euclidean distance of 20 PC scores and assigned European ancestry when at least 87.5% of the reference individuals had European ancestry (Supplementary Note). This resulted in 38,915 individuals initially classified as having European ancestry who also had at least one smoking phenotype. Second, after visually inspecting PCs 1-4 of the selected individuals, we suspected residual population heterogeneity, especially in PC 4 ( Supplementary Fig. 2). We then further restricted the samples to those whose summed Euclidean distance of PCs 1-4 fell within the 1 IQR of the European sample (N = 38,915) identified in the first step. We additionally created samples using 1.5, 2 and 3 IQR and reserved them for sensitivity analysis (Supplementary Figs. 1 and 2). After IQR filtering, we retained only unrelated individuals. Relatedness was estimated with HapMap3 variants (Hardy-Weinberg equilibrium P > 10 −6 , MAF > 0.01) using GCTAv1.92 to obtain a list of nominally unrelated individuals with pairwise relatedness (π < 0.025). We additionally removed  individuals who showed very high rare-variant sharing (π > 0.25 in the MAF 0.01-0.1% bin) who seemed to have inherited small IBD segments from distant non-European ancestors. This resulted in the final sample size per phenotype shown in Table 1 (N ranging from 14,709 to 26,257).
Phenotypes. TOPMed is a consortium of independent studies, where DNA samples were sequenced and called in a unified way. Smoking phenotypes had previously been collected independently in each of the constituent TOPMed studies. Four smoking phenotypes, each representing self-report questions assessing different stages of tobacco use, were available across most TOPMed studies. We used the same definition and coding scheme as Liu et al. 18 . AgeSmk measures the age at which an individual started regularly smoking. CigDay is the average number of cigarettes smoked per day as a current or former smoker and is grouped into five bins, with higher numbers indicating greater use. For both AgeSmk and CigDay, lifelong non-smokers are excluded (set to missing). SmkCes and SmkInit are binary variables indicating former versus current smokers and never versus ever regular smokers, with case defined as current and ever smoker, respectively. These four variables were correlated but not redundant, each measuring distinct aspects of smoking (Supplementary Note and Supplementary Table 7). Descriptive statistics for each phenotype across cohorts can be found in Supplementary Table 8.
Genotypes, LD scores, GRMs and GREML-LDMS-I. GREML estimates heritability by comparing phenotypic similarity to observed genetic similarity among distantly related individuals using a linear mixed model 64 . GREML yields biased estimates when causal variants are unevenly distributed as a function of LD and MAF 31 . To mitigate this bias, the GREML-LDMS-I method partitions SNPs into different LD × MAF bins. We initially considered ~710 million genotypes that passed strict quality filters described elsewhere 34 . We additionally removed 95,750 variants with Hardy-Weinberg equilibrium P values less than 10 −6 in the European ancestry sample (N = 38,915). Then, we calculated allele frequency separately for each phenotype using plink1.9 in a final sample that went through PC, IQR and relatedness filtering. We stratified variants by MAF and then further stratified by LD scores using the median LD scores within the two most common MAF bins. This  Table 9). This process resulted in approximately 35 million SNPs and indels (Table 1).
For each phenotype, we performed GREML-LDMS-I with the GRMs for the above-mentioned six bins and a cohort indicator matrix as random effects. The cohort matrix was an N × N matrix indicating whether a given pair of individuals belonged to the same study (1, otherwise 0). The variance explained by cohort was included in the total phenotypic variance, which could yield somewhat smaller heritability estimates than when the heritability is conditioned on cohort. We regressed out age, squared age, sex and their two-way interaction terms 65-67 from the quantitative phenotypes (that is, AgeSmk and CigDay) and used residuals. The residuals of CigDay were further inverse-rank normalized to facilitate model convergence by smoothing its distribution. GREML-LDMS-I results from other types of transformations for AgeSmk and CigDay are presented in Supplementary Table 3. We calculated 10 PCs from common variants (MAF ≥ 1%) and another 10 PCs from rare variants (0.01% < MAF < 1%) after LD pruning (plink -indep-pairwise 100 5 0.2 for common and 2,000 10 0.02 for rare variants). These 20 PCs and the sequencing centre were entered as fixed effects in the linear mixed model. The same demographic variables and 20 PCs were entered as fixed effects for the binary variables. We used GCTAv1.92 for both the construction of GRMs and the GREML-LDMS-I analysis. We replaced the original GRM diagonals with 1+ inbreeding coefficients (Fhat3 computed from the GCTA software) because the original formula tended to yield far higher variance of diagonals for rare variants than expected in theory (Supplementary Table 10) 68 . We allowed the estimates to be negative to obtain unbiased estimates. Heritabilities for binary phenotypes were analysed under a liability threshold model 69 . Population prevalence was set at 0.15 and 0.42 for SmkCes and SmkInit, respectively, on the basis of smoking prevalence in the UK Biobank dataset, to allow for ready comparison with this publicly available and widely used dataset. For all traits, total heritability was calculated by adding the heritability estimates of the six bins with standard errors approximated by the delta method 70 .
Partitioning rare-variant heritability. To further interrogate sources of rare-variant heritability, we divided rare-variant bins into protein-altering versus non-protein-altering variant bins 71 . Functional impact of variants was assessed by snpEff 4.3 with 'HIGH' and 'MODERATE' impact categorized as protein-altering and 'LOW' and 'MODIFIER' categorized as non-protein-altering 72 . 'HIGH' includes variants expected to have a disruptive impact on the protein, such as protein truncation and loss of function. 'MODERATE' includes variants that are expected to influence the effectiveness of the protein, such as missense and splice region variants. 'MODIFIER' and 'LOW' include variants that influence non-coding genes or are located in non-coding regions and are usually considered harmless to protein behaviours. The details of this classification can be found at https://pcingola.github.io/SnpEff/se_inputoutput/. Variance components were then estimated for a total of eight bins.

Sensitivity analyses.
To evaluate the degree of sensitivity of our results to residual population structure (especially that of rare variants) and various analytic decisions, we conducted extensive sensitivity analyses. First, we performed GREML-LDMS-I with 40 PCs (20 common and 20 rare) and 100 PCs (50 common and 50 rare) as fixed effects. Second, we estimated heritability after removing individuals whose total length of shared IBD segments covered more than 2.5% of the genome (N = 121-177 depending on the phenotype) to test whether heritability attributed to rare variants is primarily driven by a subset of samples sharing recent ancestors and presumably higher environmental similarity. Third, we also adjusted for 20 PCs from the IBD-based GRM, in addition to the 10 PCs from common variants and 10 PCs from rare variants (40 total PCs). IBD segments were estimated using 632,008 common variants (LD pruned by plink -indep-pairwise 1,000 50 0.6) with FastSMC 68 with the default settings, except we set the minimum IBD length to ≥2 cM. We used decoding files precomputed by the authors of FastSMC and only retained IBD segments with quality score >0.1. The IBD-based GRM was created by summing the length of the shared IBD segments in a given pair of individuals and dividing it by the length of the diploid human genome (3,608 × 2 cM). We then derived PCs from the IBD-based GRM. Fourth, we added 53 recruitment sites nested within the cohorts as random effects. This additional control is meant to better account for rare variants that may cluster geographically. Fifth, we performed GREML-LDMS-I with permuted samples within genetically close neighbours 36,44,73 . For this, we created an N × N distance matrix populated by scaled Euclidian distance of PCs 1-10 calculated from LD-pruned common variants for every pair of individuals. We then randomly exchanged genotypes of a given individual with one of their 100 nearest neighbours using sampling algorithms from LocPerm, which was developed to control for population structure in the rare-variant association test 44 . We estimated heritability in 100 replicates of permuted genotypes with the same demographic covariates and cohort used in the main analysis. The mean (ĥ 2 null ) and s.d. of heritability estimates from 100 replicates were calculated for each bin and were tested against zero using a one-sided Z-test, as heritability is bounded above zero (Supplementary Note). Last, we estimated heritability in samples with increasingly less strict PC-based ancestry filtering mentioned earlier (1.5, 2 and 3 IQR). The greater the IQR threshold, the more ancestral variation can be present in the resulting sample, the larger the sample and the greater the chance of observing effects of population stratification.
Pedigree-based heritability. We created a GRM with all available samples, including related ones but excluding pairs related greater than 0.80 to exclude identical twins and duplicates. To aid in model identification, we included cohorts that had at least ten first-degree relatives defined as relatedness greater than 0.375. A list of cohorts included in this analysis is presented in Supplementary Table  11. We created a GRM using common variants (MAF > 5%) and set GRM entries of pairs related less than 0.05 as zero and retained the rest of the entries as they were. This GRM was fitted together with a cohort matrix and the same set of fixed-effect covariates used in the primary analysis. To test whether the resulting ĥ2 ped is underestimated due to the relatively low level of relatedness structure in the sample, we repeated the analysis without a random effect of cohort, using only the FHS, which has high proportions of related individuals.

SNP-based heritability in African ancestries.
We initially had 19,788 individuals who had one of smoking phenotypes and were classified as having African ancestry by the ancestry assignment procedure described earlier. We selected unrelated individuals (π < 0.025) by applying the pcairPartition function in the GENESIS package on kinship coefficients from pcrelate 74 , where familial relatedness is estimated with MAF (>1%) and LD-pruned (|LD| < 0.32) variants after accounting for population structure. We additionally excluded 58 individuals who showed more than two-way admixture as determined by global ancestry proportions of either South-Central Asia, East Asia, Native America or the Middle East greater than 10%. This resulted in the sample sizes in Table 2. The resulting samples were predominantly of African ancestry with recent European ancestry admixture (mean African ancestry proportion, 83%; s.d., 8.9%; mean European ancestry proportion, 15%; s.d., 8.4%). The ancestry proportions of each individual were estimated from local ancestry inference by RFMix 75 , which estimates the ancestry of admixed individuals at each genomic segment of two homologous chromosomes, using the Human Genome Diversity Panel as a reference panel 76 (Supplementary Fig. 3). Given the smaller sample size of African ancestry samples, we restricted the variants to be analysed to MAF 0.1-1% for SmkInit (N = 11,744) and common variants for the other three phenotypes (N = 6,796-7,549). We applied the same analytic procedure as in the European ancestry samples: (1) we selected genetic variants with minor allele count > 2 and Hardy-Weinberg equilibrium P > 10 −6 , (2) we partitioned the variants into MAF (0.05, 0.5] and high LD, MAF (0.05, 0.5] and low LD, MAF (0.01, 0.05] and high LD, MAF (0.01, 0.05] and low LD, and additionally MAF (0.001, 0.01] for SmkInit (Table 2), and (3) we calculated the GRM for each bin and performed GREML-LDMS-I with the same set of demographic and cohort covariates as in the European ancestry analyses. We also included 50 common-and 50 rare-variant-based PCs and global ancestry proportions of five continents (sub-Saharan Africa, South-Central Asia, East Asia, Native America and the Middle East) as fixed effects and local ancestry kinship (Supplementary Note) as a random effect to account for complex population structure in the admixed sample 42,43 . Pedigree heritability was not estimated for individuals of African ancestries, as SNP-based relatedness is no longer proportional to IBD in admixed populations 77 .
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
Phenotypes are available through an authorized access portal in dbgap (https://dbgap.ncbi.nlm.nih.gov/) or direct request to TOPMed principal investigators. Accession numbers and email addresses of the principal investigators are presented in the Supplementary Note. Genetic data are available through the dbgap TOPMed exchange area.

Code availability
All software used is publicly available and can be found at the references cited.