Genetic architecture of four smoking behaviors using partitioned h²_SNP

Introduction

Cigarette smoking is a leading cause of premature death worldwide(1), and many smokers struggle to quit, despite interest and numerous attempts(2). Smoking prevalence has decreased in recent decades due to public health efforts(3); however, rates of alternative forms of nicotine use (e.g., vaping) have grown rapidly during this time(4), demonstrating a pressing need to characterize the underlying biology of nicotine use and smoking to reduce subsequent premature death.

A key aspect of that underlying biology is the genetic architecture(5) of smoking behaviors, including the relative contribution of rare vs. common variants, functional annotation of associated loci, and characterization of the neutral and selective forces shaping that architecture. Numerous(6-10) twin, adoption, and family studies have demonstrated that up to 50% of the variance in nicotine dependence and individual smoking behaviors, such as quantity, is attributable to genetic influences. Recent genome-wide association studies (GWAS) have improved our understanding of this genetic basis by identifying over two hundred conditionally independent loci associated with these traits to date(11-14). This genetic signal is enriched in loci that influence the epigenome and within specific brain regions such as the hippocampus, providing a more nuanced interpretation of specific class(es) of variants, candidate brain regions, and potential causal mechanisms that influence smoking(11). Together, this body of work strongly indicates a highly polygenic architecture to smoking behaviors. Nonetheless, significantly associated loci collectively explain only a small proportion of the family-based genetic variance, leaving many additional loci undiscovered and the majority of the genetic variance unexplained.

While additional common variants of very small effect are likely to be identified as sample sizes grow, some of the unexplained variability undoubtedly arises from uncommon and rare variants (MAF<0.01), though their relative contribution is uncertain. The most recent large GWASs(11, 15) of smoking behaviors and nicotine dependence, using LD score regression (LDSC), estimate the SNP-based heritability (h²_SNP) due to common variants as 0.05-0.09 across traits(16). A related exome sequencing study(17) estimated that rare coding variants explained approximately 1-2% of the phenotypic variance. However, given that the majority of identified associations are intergenic(11), exome-based studies are unlikely to identify most rare variants influencing these behaviors. Thus, the rare-variant contribution to smoking behaviors may yet be substantial when assessed with methods that can account for the aggregated influence of common and rare variation.

Additionally, the contribution of non-additive genetic variance to these smoking behaviors is poorly understood. Twin-based studies have typically evaluated ACE models(9), which estimate additive genetic (A), common environment (C) and unique environment (E) variances using twin correlations, implicitly assuming zero dominance genetic variance. Extended twin kinship models can estimate dominance genetic variance and shared environmental effects simultaneously, and the only such model to evaluate smoking initiation found no evidence of dominance genetic variance(18). Although several h²_SNP estimates of smoking behaviors have been to published, to our knowledge only one estimate of SNP-based dominance genetic variance (δ²_SNP) has been reported, which found δ²_SNP of smoker status indistinguishable from zero(19). Furthermore, the allele frequency spectrum and contribution of functional annotations related to LD, allele frequency, recombination and related genomic features for smoking behaviors has not been fully explored. One study applied partitioned h²_SNP approaches to evaluate tissue-specific effects, with results indicating that genes expressed in the cerebellum are enriched in their contribution to nicotine dependence(15). The only published work has examined a single trait, smoking status, finding contributions of low-LD and -MAF variants consistent with negative or purifying selection(20, 21). Whether these same patterns exist for other smoking behaviors, such as quantity of use or cessation, is unknown.

A comprehensive evaluation of the frequency spectrum, the influence of dominance genetic variance, and the contributions of functional annotations is needed to provide a more complete picture of the genetic architecture underlying complex smoking behaviors. Here, we use recently developed methods(19, 21-26) to evaluate these heritable contributions and characterize the genetic architecture of four smoking behaviors: smoking initiation (whether an individual has ever been a regular smoker), age of initiation of regular smoking, cigarettes per day (evaluated with different data encodings), and smoking cessation.

Methods

Results

Using GREML-LDMS-I with unrelated individual, we estimated smoking initiation h²_SNP(SE)=0.176(0.007), smoking cessation h²_SNP=0.119(0.018), cigarettes per day h²_SNP=0.098(0.011), and age of initiation h²_SNP=0.055(0.011) (Figure 1, Table S1). MAF- and LD-partitioned heritability estimates differed across traits. Common variants (MAF>0.05) contributed substantially to all traits, particularly common variants with relatively low LD (Figure 1). Uncommon variants (MAF 0.01-0.05) with low LD, but not high LD, contributed to all traits. Alternatively, uncommon (MAF<0.01) variants contributed significantly only to smoking initiation and age of initiation, and rare (MAF<0.001) variants did not contribute significantly to any trait.

• Download figure
• Open in new tab

Figure 1.

h²_SNP estimates (+/- standard error) across four smoking behaviors, partitioned using GREML-LDMS-I. Note that twin-based estimates are roughly 50% across these smoking traits.

Notably, we estimated significantly different (non-overlapping 95% CI) total and binned h²_SNP for different CPD encodings. Total h²_SNP ranged from 0.092(0.011) for the raw CPD count to 0.289(0.038) when CPD was dichotomized into heavy(CPD>20)/light(CPD<=10) smokers (Figures 2 & S2-S3, Tables S1-3). All dichotomized CPD total h²_SNP estimates (except using the median) were >0.2. We found differences in partitioned estimates, where common variants (MAF>0.05) contributed to substantially higher h²_SNP of heavy/light CPD than the other CPD encodings. Rarer (MAF 0.001-0.01) variant contribution was also higher, though the smaller sample size of the dichotomized data led to larger standard errors.

• Download figure
• Open in new tab

Figure 2.

h²_SNP estimates (+/- standard error) of CPD using different phenotype encodings, partitioned using GREML-LDMS-I. Heavy vs. light is dichotomized with Light: CPD<=10 and Heavy: CPD>20; estimated h²_SNP shown on the liability scale using a prevalence of 0.42.

We estimated the contribution of dominance variance. For all traits, the 95% CI of δ²_SNP estimates overlapped zero (Table S4).

The relatedness threshold strongly influenced estimated h²_SNP when using PCGC, but not when using GREML (Tables S1-3, Figs. S2-S6). Specifically, the PCGC estimates were considerably higher than GREML estimates when applying a relatedness<0.05 cutoff with smoking initiation and smoking cessation, but dropped and had overlapping 95% CIs at lower relatedness thresholds. The higher estimates when using PCGC with relatedness<0.05 were driven by a much greater contribution of rare variant h²_SNP (MAF<0.0001; Figure S5-S6).

We applied partitioned LDSC to assess contribution of functional annotations and the role of LD and selective constraint in smoking behaviors. Across smoking behaviors, we found that SNPs that were highly conserved, that had lower MAF-adjusted LD or lower MAF quantiles (MAF>0.001 in Liu et al.(11)), and that were in areas of high CpG content and low recombination rate contributed significantly to heritable genetic variation (Figures 3 & S7, Table S5).

• Download figure
• Open in new tab

Figure 3.

Partitioned LDSC regression coefficient p-values for all annotations with at least one significant coefficient across all traits. See Supplemental Figure S5 & Table S3 for all annotations.

Discussion

We estimated h²_SNP and δ²_SNP across four key smoking behaviors, and partitioned variance among rare vs. common variants and functional annotations. Our h²_SNP estimates are more than double the previously reported(11) LDSC-based and single-component GREML-based estimates for smoking initiation (0.18 vs. 0.08 and 0.12) and smoking cessation (0.12 vs. 0.05 and 0.06), but are nearly identical for binned CPD (0.1). Our estimate of age of smoking initiation h²_SNP=0.05 is nearly identical to the LDSC-based estimate, but is much lower than the previous single-component GREML estimate of 0.11. The difference in age of initiation h²_SNP may be due to including all variants in a single GRM when the causal variants are relatively common(22). Partitioned estimates of common, well-tagged variants are similar to the LDSC-based estimates(11) across all four traits, consistent with expectations, as LDSC estimates variance due to common, well-tagged variants(16, 22). The higher h²_SNP estimates for smoking initiation and cessation results from larger contributions of low-LD and low-frequency variants (MAF<0.01), suggesting that for these traits, a non-trivial portion of the genetic variance is due to rarer variants and those that are poorly tagged by surrounding SNPs. This contribution is likely underestimated in the current study, because even with HRC-imputed data, these sites are typically poorly imputed, which leads to a downward bias in h²_SNP estimates(22, 23).

Alternative CPD encodings led to different estimates, wherein total h²_SNP for dichotomized heavy/light smoker status was over twice that of other encodings. This may be explained by one or more possible phenomenon that occur after restricting the analyses to phenotypic extremes, i.e., removing the center of the distribution. First, the extremes of the CPD distribution may be capturing a phenotype more closely approximating physical dependence on nicotine. Tolerance and withdrawal may index severity of nicotine dependence(36), a construct for which we do not have formal diagnoses, but which is highly heritable. In this case, though lacking other important aspects of the clinical presentation such as craving or loss of control, the dichotomized heavy/light phenotype is comparing individuals who may find overnight abstinence less aversive and start smoking later in the day, and endorse lower levels of nicotine dependence (light) to those who meet criteria for severe nicotine dependence (heavy), whereas the standard continuous CPD encoding includes intermediate levels of smoking heaviness that may or may not correlate with clinical presentations of nicotine dependence. Our GREML-based estimate of common, well-tagged h²_SNP (∼0.09) is approximately the same as one recently reported LDSC-based estimate of nicotine dependence(15), consistent with this hypothesis.

Alternatively, the dichotomized phenotype may reflect lower environmental variance and result in higher h²_SNP, if, for example, environmental effects such as reduced access to cigarettes or regular use of nicotine replacement therapy lead to intermediate values of CPD. Such differences in variance cannot be tested when either trait is dichotomous because the liability underlying the dichotomous trait must be assumed to have unit variance. Ongoing work will seek to distinguish between these two possibilities, and determine whether variants that contribute to heavy/light CPD and other smoking behaviors examined here also contribute to nicotine dependence liability or severity.

We found no evidence of dominance genetic variance for any phenotype, though we note that the power to detect δ²_SNP is lower relative to h²_SNP(19) and therefore sample size may be a limiting factor to detect low, but non-zero δ²_SNP. Our findings are consistent with those of Zhu et al.(19), who reported low δ²_SNP across 79 traits and δ²_SNP∼0 for one smoking phenotype, smoking status. For the four smoking phenotypes in the current study, dominance genetic variance likely contributes little or not at all to the phenotypic variance. We note that dominance interactions of alleles within individual loci may still be contributing to these traits, but as this contributes to additive genetic variance (i.e., h²_SNP), its contribution to δ²_SNP can be limited, particularly for low-frequency variants(37). Alternatively, interactions between, rather than within, loci may lead to epistatic genetic variation underlying smoking behaviors, and such effects could not be tested using the current approach.

We identified several functional annotations related to LD, MAF, and sequence conservation that significantly contribute to h²_SNP (Figure 3, Table S5). In addition, GREML-LDMS-I h²_SNP analyses identified higher contribution of poorly tagged variants relative to well-tagged variants within the same MAF range across all four traits, and also identified nominally significant (95% CI>0) contribution of rare variants (MAF<0.01) for smoking initiation, raw CPD count, and age of initiation. Across the four traits analyzed, rare variants accounted for between 10 and 20% of total h²_SNP (Table S1). This suggests a role of low-frequency SNPs in low LD with surrounding regions, consistent with purifying and background selection acting to remove mutations with deleterious effects. Given that tobacco use in high concentrations, such as found in cigarettes, is evolutionarily novel for humans, it is unlikely that negative selection acted directly on these smoking behaviors, but rather mutations that today influence nicotine related behaviors may have pleiotropic effects on other traits that were subject to negative selection across evolutionary time(20).

Our h²_SNP estimates are still considerably lower than twin-based estimates, which range from 50%-80% for dependence, smoking initiation and quantity of use(6-10), suggesting that additional still-missing heritability remains. This is unlikely to be explained by common causal variants, which are well-tagged in current imputation reference panels and from which we expect little downward bias in h²_SNP estimates(22). Further work will be required to fully characterize non-additive genetic variance, such as epistasis or gene-environment interaction. Rare variants are a likely source of the still-missing heritability. The SE of the rarest MAF partitions were substantially larger than the common variant partition SE, indicating that increased sample size will improve the precision of estimates of rare-variant contribution. Overall, estimates are still generally low compared with those attributable to common variants, and even with large reference panels such as the HRC, rare variants are expected to be poorly imputed, resulting in downwardly-biased h²_SNP(22, 23). Further work through deep sequencing of large samples(38) or using those deeply-sequenced individuals as an improved imputation reference panel is needed to obtain less-biased estimates of rare-variant h²_SNP. For example, height and BMI h²_SNP estimates using whole genome sequencing have approached twin-based heritability estimates; rare variants account for a substantial proportion of the heritability(39).

Beyond the limitation of rare variant imputation, our study highlights several key issues in h²_SNP estimation. First, while we used the largest relatively homogenous sample available, even larger samples will be needed for more precise estimation of rare variant contribution, as demonstrated by the much smaller SE of h²_SNP estimates of traits with larger sample sizes. Second, estimates are sensitive to the estimation method, i.e., H-E Regression-based vs. GREML, which may be due to how environmental confounding differentially influences estimates across methods. GREML-based estimates were relatively stable across relatedness thresholds (Table S1). However, PCGC-based estimates were quite sensitive to relatedness thresholds, being much higher than GREML-based estimates at a .05 threshold and declining with lower thresholds. Although a full assessment of performance of estimators is beyond the scope of this study, it will be important to assess the potential for environmental confounding. As with the possibility of rare variant-environment confounding in GWAS(40), environmental confounding is particularly relevant to estimates of rare variant h²_SNP because very rare variants are more likely to be shared by individuals sharing recent common ancestors and who may therefore be more likely to share environmental influences. Models that incorporate environmental sharing of families, partners, and close relatives or geography (e.g.,(41, 42)) are a possible avenue to address confounding. To this effect, we note that a full extended twin family design found a lower and possibly sex-dependent estimate of common additive genetic variance, as well as strong environmental influences(18).

In conclusion, though our h²_SNP estimates of the four different smoking behaviors were generally modest, they are higher than previously published estimates for smoking initiation and cessation and indicate that additional genetic variance may be explained by low- and rare-frequency variants, which may be due to the impact of purifying selection on genes involved in these highly polygenic traits. Quantity of use, as measured by CPD, may also be modestly heritable, but as it depends on the encoding of the variable, additional characterization of the phenotype and its relationship with nicotine dependence is required. All estimates will be improved by the use of complete whole genome sequencing of large numbers of individuals(38), including the contribution of rare variants to smoking behaviors.