key: cord-0863388-nn7egvkm
authors: Schmitz, Lauren L.; Goodwin, Julia; Miao, Jiacheng; Lu, Qiongshi; Conley, Dalton
title: The impact of late-career job loss and genetic risk on body mass index: Evidence from variance polygenic scores
date: 2021-04-07
journal: Sci Rep
DOI: 10.1038/s41598-021-86716-y
sha: db5ed04191f7ebde5425f701a0e54e018d39c8a2
doc_id: 863388
cord_uid: nn7egvkm

Unemployment shocks from the COVID-19 pandemic have reignited concerns over the long-term effects of job loss on population health. Past research has highlighted the corrosive effects of unemployment on health and health behaviors. This study examines whether the effects of job loss on changes in body mass index (BMI) are moderated by genetic predisposition using data from the U.S. Health and Retirement Study (HRS). To improve detection of gene-by-environment (G × E) interplay, we interacted layoffs from business closures—a plausibly exogenous environmental exposure—with whole-genome polygenic scores (PGSs) that capture genetic contributions to both the population mean (mPGS) and variance (vPGS) of BMI. Results show evidence of genetic moderation using a vPGS (as opposed to an mPGS) and indicate genome-wide summary measures of phenotypic plasticity may further our understanding of how environmental stimuli modify the distribution of complex traits in a population.

within-person variability is hindered by a lack of large datasets with genotype data and longitudinal phenotypic data on participants, researchers have developed methods that can detect population-level variance effects that are not driven by mean effects, referred to herein as vGWAS 16, [19] [20] [21] [22] [23] . In this study, we apply summary statistics from both mGWAS and vGWAS to construct whole genome PGSs for BMI that capture mean effects (mPGS) and variance effects (vPGS) 24 . Evaluating both measures in a G × E framework is necessary because environmental shifts may moderate individuals' propensity for higher or lower BMI, and/or their propensity towards changes in BMI or BMI plasticity 21, 24 .

Our data on older workers aged 50-70 come from the U.S. Health and Retirement Study (HRS). The HRS is a nationally representative, longitudinal study with genotype data and over twenty years of sociodemographic data on respondents, including individual-level exposures to involuntary job losses from business closures. We focus specifically on business closures because they are typically the byproduct of external, firm level decisions to restructure or relocate businesses and are therefore considered more exogenous than layoffs or firings, which may be correlated with unobserved health or worker characteristics that could bias G × E estimates [25] [26] [27] . The majority of G × E interaction studies use endogenous measures of the environment that cannot address the non-random distribution of genes across environments. This is important because G × E interactions can, in that case, be proxying a different, unmeasured E that is interacting with G, G × G interactions (i.e., epistasis), or even E × E (if the measured genes proxy other environments). Specifically, if the cause of job loss is endogenous, such a measure could be intertwined with a host of unobserved genetic or environmental influences that are associated with health and changes in BMI [28] [29] [30] .

To address this, our empirical strategy interacts business closures with, respectively, an mPGS and vPGS in a regression-adjusted semiparametric difference-in-differences (DiD) propensity score matching framework that compares the BMI of those before and after an involuntary job loss with a control group that was not laid off. Combining propensity score matching with DiD estimation makes the model more robust to selection on observables and unobservables with time invariant effects (e.g., ability or worker preferences) 31 . This is necessary because although business closures are plausibly more exogenous than layoffs or firings, it is still possible that workers with unhealthy behaviors or poor health, for example, could select into more vulnerable or volatile industries 32 . To date, we are aware of only one other study that has leveraged a vPGS and a quasi-natural experiment (education reform in the UK) to detect G × E interaction effects on BMI and educational attainment 24 . Results from this study found evidence of mPGS and vPGS interaction effects, indicating that both forms of moderation need to be tested in G × E interaction studies.

In the context of older workers in the U.S., we focus on changes in BMI for two reasons. First, BMI is an inexpensive, non-invasive proxy measure of adiposity that is available for all HRS waves and is predictive of metabolic syndrome and other more difficult to measure anthropomorphic measures like abdominal adiposity that increase risk for cardiovascular disease and type 2 diabetes 33 . In older adults, unintentional weight loss or frailty can also be harmful and indicative of decreased resistance to stressors, resulting in greater vulnerability to disease and disability [34] [35] [36] [37] [38] . Past studies have found slight increases in BMI from unemployment 31, 39 , particularly in middle-aged workers, but there is no consensus, and some studies that look at the causal impact of job loss on BMI or related health outcomes have been unable to locate an average treatment effect [40] [41] [42] . Using the HRS, Salm finds no causal effect of business closures on various measures of physical and mental health, whereas Gallo et al. find involuntary job loss is associated with increased depressive symptoms and risk of stroke but not myocardial infarction 6, 40, 43 . A few studies have explored the possiblity that the health impacts of job loss vary across the distribution of health status 39, 42 . For example, using finite mixture models, Deb et al. found increases in drinking and BMI were concentrated among workers who were already pursuing unhealthy behaviors pre-job loss, indicating the effects of job loss may be especially problematic for high-risk individuals. However, because genetic data have only recently become available in large population studies, past research was not able to explore the possibility that the impact of job loss on intra-individual fluctuations in BMI may vary across the spectrum of genetic risk.

Second, we focus on changes in BMI because BMI is currently the most well studied phenotypic trait in vGWAS. Previous meta-analyses of mGWAS have identified more than 100 genome-wide significant loci associated with BMI [44] [45] [46] [47] [48] . The largest cluster of highly significant genetic variants is located in the FTO (fat mass and obesity associated) gene region on chromosome 16. Studies suggest FTO polymorphisms increase obesity risk through subtle changes in food intake and preference and affect pathways in the central nervous system that regulate appetite 45, 49 . In particular, the SNP rs1421085 underlies the association between the FTO locus and obesity via activation of IRX3 and IRX5, which play a role in the differentiation of adipocyte subtypes 50 . Recent vGWAS have found evidence for loci with variance effects on BMI located in genes responsible for adipocyte differentiation (PPARG ) and genes implicated in the pathology of obesity, diabetes, atherosclerosis, and cancer (FTO, PPARG , CCNL1, TCF7L2, ZNF668, GIPR) 16 . Most BMI-associated loci have their largest impact early in life or during adolescence 51 , although a few loci, which have also been associated with type 2 diabetes or coronary artery disease, exhibit stronger effects in older adults 52 . Past studies have also found genetic moderation of social aspects of the environment that affect BMI, including lifetime socioeconomic status (SES), social norms, birth cohort, and institutional policies [53] [54] [55] .

In this study, we used a quasi-natural experiment to investigate whether the effect of job loss on BMI-a stressful and often debilitating lifetime event-is moderated by genetic predisposition. To incorporate genetic effects associated with the population variance in BMI, we used summary statistics from recently developed GWAS methods that can separate SNP mean and variance effects 16 to construct vPGS for BMI 24 . Prior to constructing the vPGS in the HRS sample, we validated its performance in the UK Biobank (see Fig. 1 for an overview of our analytic workflow). Similar to past studies that have used business closures to examine the effects of job loss on BMI, we classified workers who have been laid off due to a business closure as being "treated" and compared these individuals to a "control" group who reported working for the same employer the entire time they were in the HRS sample. Our sample is limited to individuals of European descent because comparable mGWAS and www.nature.com/scientificreports/ vGWAS in other ancestral populations are currently unavailable. Overall, we hypothesized that a vPGS would be better able to capture downstream differences in the genotype-phenotype relationship between treated and control groups, and that highly plastic individuals, or those with a higher vPGS, would adapt more quickly to changes in their environment, and therefore maintain a more stable weight in the face of job loss.

Matching quality and summary statistics. Table 1 displays the means and matching statistics for covariates by treatment group, control group, and matched control group. Detailed variable descriptions can be found in Supplementary Table S1, and additional descriptive statistics and matching statistics are reported  in Supplementary Tables S2 and S3 , respectively. After matching, covariates should be balanced with little to no significant differences remaining. We included both the standardized bias and two-sample t-tests for equality of the means to check for significant differences between covariates for both groups (Methods) 56 . Before matching, individuals affected by a business closure have lower socioeconomic standing, were less likely to have health insurance, and reported worse mental health and health behaviors than continuously employed individuals. Labor statistics show they were more likely to work part time, for smaller firms, in the agriculture/fishing/farming, construction/mining, manufacturing, or trade industries, were more likely to be blue collar, and have lower job tenure than workers in the control group. Significant differences in the mean or variance of BMI between treatment and control groups are not apparent before or after matching.

After matching, covariates are more balanced overall, and the standardized biases for the majority of variables are at or below 5% and/or the t-test p value > 0.05, which indicates that mean differences between the treatment and control group are small and the balancing procedure was effective 57 . Notable exceptions include mean differences in education, industry, household income, access to health insurance, and smoking behavior. To minimize any remaining differences between groups and increase the precision of our treatment effect estimates, we controlled for all covariates in our empirical model 58 . Importantly, we do not see any significant differences in the mPGS or vPGS between treatment and control groups before or after matching, indicating the absence of gene-environment correlation (rGE), or genetic selection into the treatment group.

Construction and predictive performance of vPGS. Because we are incorporating mPGS and vPGS into our G × E interaction model, it is important to use a vPGS that captures variance effects that are distinct from mean effects 24 . To decorrelate the mean and variance effects, Young et al. proposed a dispersion effects test that can identify differences in the variance of the GWAS sample as a whole that are not driven by mean effects at the SNP level 16 . We used dispersion weights from Young et al. to construct vPGS in the HRS (Methods).

The predictive performance of mPGS for BMI in the HRS and other European ancestry population-based samples has been well studied [53] [54] [55] 59, 60 , and previous GWAS have reported predicted R 2 values that range between 5 and 10% 47 . Consistent with these studies, the mPGS explains 7.2% of the variance in BMI in the HRS European ancestry sample (Methods). To evaluate the predictive performance of the vPGS, we fit a Double Generalized www.nature.com/scientificreports/ Linear Model (DGLM) that allowed us to assess the association between the vPGS and the between-individual variance in BMI in a UK Biobank (UKB) test sample that is independent of our HRS testing sample (Methods). Table 2 reports associations from the DGLM with and without mPGS adjustment (Models 1 and 2, respectively). The dispersion vPGS is significantly associated with the population variance in BMI in the UKB test sample (p = 1.01E−04), and this association holds after controlling for the mPGS (p = 1.44E−04). Since DGLM is a loglinear variance model, the effect size estimate of 0.019 means that a one unit increase in the standardized vPGS results in an approximately 1% increase in the standard deviation of BMI on the inverse normal transformed scale (i.e., exp(0.019) − 1 ≈ 1%).

We used a propensity score matched DiD model to evaluate the effect of job loss from a business closure (Methods). Table 3 shows separate propensity score adjusted DiD results from specifications with and without the mPGS and vPGS interactions. To control for multiple hypothesis testing, we adjusted p values using the Benjamini-Hochberg false discovery rate (FDR), which controls for the proportion of falsely rejected null hypotheses among all those rejected 61 . Results with FDR < 5% are indicated with an asterisks. The results in Column 1 are similar in magnitude and direction to the HRS results reported by Deb et al., which find a positive, but insignificant main effect from business closures in the full HRS sample on changes in BMI in all ancestry groups 39 . Among European ancestry respondents between the ages of 50 and 70 who reported being in the labor force, the main effect is still positive (Column 2). In the genotyped European ancestry sample, the main effect of business closures is insignificant but negative (Column 3). Columns 4-6 show that the inclusion of both the mPGS and vPGS is necessary to uncover a genetic main effect and an interaction effect: the mPGS captures a significant main effect of genotype on BMI (p = 0.001), while the vPGS captures a significant G × E effect (p = 0.001). Graphically, this can be seen in Fig. 2 , which used the estimated parameters from the DiD regression models in Table 3 to predict BMI at different values of the mPGS and vPGS for treated workers in the wave following a business closure and for control workers that were matched to treated individuals in the same HRS wave (Methods). The regression coefficients that were used to estimate predicted BMI varied depending on treatment status and the value of the mPGS and vPGS (Methods). In the mPGS figure, which plots predicted BMI based on coefficients from the mPGS interaction model in Column 5, we see differences in predicted BMI by mPGS, but no significant differences between groups. Conversely, in the vPGS figure, which plots predicted BMI based on coefficients from the vPGS interaction model in Column 6, there are no differences in the predicted level of BMI by vPGS, but as indicated by the cross-over shape of the interaction, there is suggestive evidence of a G × E interaction, or evidence of environmental moderation by vPGS for the treatment group relative to the control group in the post-treatment wave. Significant differences between treatment and control groups can only be seen in the lower half of the vPGS distribution; workers below the mean in the treatment group appear more likely to lose weight as a result of a business closure relative to control workers with similar vPGS scores.

We conducted an event time study (ETS) to assess the validity of our findings and to show the evolution in BMI by vPGS for the treatment and control groups up to four years post job loss (Methods). The assumption underlying the DiD research design is that in the absence of an involuntary job loss, BMI would have evolved similarly for the treatment and control groups (i.e., the "parallel trends" assumption). Figure 3 plots the coefficient estimates from the ETS model, which can be interpreted as the difference in BMI between treatment and control groups (Supplementary Table S4 ). The first panel of Fig. 3 indicates the presence of parallel trends in BMI prior to a business closure for the full sample-i.e., the difference between treatment and control groups is close to zero and not statistically significant. We then estimated separate ETS regressions for respondents in the top and bottom 50% of the vPGS distribution to compare trajectories in BMI for treatment and control groups by vPGS. Similar to the results in Fig. 2 , we found suggestive evidence that individuals in the bottom 50% of the vPGS distribution have a lower BMI on average compared to the control group (p = 0.043). These effects did not persist in the next HRS interview wave at t + 2, or up to four years post job loss. 

Gene-environment interplay is a fundamental biological process that influences the diversity of outcomes we observe in human populations 62 . However, because genetic differences are tightly interwoven with environmental differences, it is challenging to identify genomic and environmental factors underlying phenotypic plasticity. The search for interaction effects is further complicated by the fact that the majority of GWAS methods are unable to separate genetic effects on phenotypic variability from effects on the mean or level of trait values 16 . As a result, most PGS × E interaction studies cannot detect interaction effects that are driven by loci that affect plasticity 24 .

Research suggests that SNPs associated with the variance of BMI (vQTLs) are highly enriched for G × E interactions, and that these vQTLs play an important role in cellular response to the external environment 16, 22 .

Results from this study indicate that an unexpected job loss did not alter the mPGS-BMI relationship; individuals with higher mPGS had higher levels of BMI regardless of whether they were in the treatment or control group. However, we do see suggestive evidence of genetic moderation by vPGS. Genetic moderation is particularly pronounced in the lower half of the vPGS distribution; less-plastic individuals in the bottom 50% appear to adjust more slowly to environmental changes, resulting in minor weight loss compared to similarly matched individuals in the control group. Results from an ETS analysis show that changes in BMI were detectable in the wave immediately preceding job loss but did not persist in subsequent HRS waves. In the context of the job loss Table 3 . Effect of job loss from a business closure on BMI with and without mPGS and vPGS interactions for workers aged 50-70 in the Health and Retirement Study (HRS). *FDR corrected p value < 0.05. Abbreviations: SE, standard error; CI, confidence interval. Robust standard errors in parentheses. Each column in the table shows separate propensity score matched models for workers aged 50-70. All specifications adjust for BMI in the previous wave, or BMI (t − 2), and for the conditioning variables used in the propensity score matching that are reported in Table 1 and defined in detail in Supplementary Table S1. Columns 1-3 do not include mPGS or vPGS in the matching procedure. Column 1 includes all workers, regardless of ancestral background or the availability of genotype data and includes additional controls for race and Hispanic ethnicity in the matching procedure. Column 2 reports results for all European ancestry workers, regardless of the availability of genotype data. Column 3 reports results for the European ancestry sample with genotype data prior to matching on the mPGS and vPGS. Columns 4-5 include the mPGS in the matching procedure, and Column 6 includes the mPGS and the vPGS in the matching procedure. Individuals in the control group can have multiple observations. In the analytic sample, unique N(control) = 3564; unique N(treated) = 375. 63 , genetic predisposition may be another avenue through which health inequalities emerge and deepen within a population. As a result, a better understanding of the extent to which biological forces act as a mechanism between worker displacement and workers' health may improve our estimation of the short-and long-term effects of job loss. It is important to note that we deployed, by necessity, a vPGS that was constructed from weights that were trained in a discovery sample to predict variation between individuals net of mean effects. However, we used this measure in analysis that examined within-subject variation. This is an important distinction that may inform the interpretation of our results. Namely, our theory is that highly plastic individuals are better able to adapt to changing environmental contexts and thus maintain a more stable weight in the face of job loss over time. We classified individuals as plastic or non-plastic based on their score from a genetic model that predicts whether an individual scores higher on a cross-individual model of dispersion (independent of mean effects). In a sense, this means that someone with a higher vPGS has more "noise" in their prediction than another individual's score does-that is, s/he is less well-predicted from a levels (mean values) regression than someone with a lower score. Someone who is low on the plasticity score has a BMI that is better predicted by their levels effects than someone who is high on the score. This, in turn, we think is indicative of someone whose phenotype is more affected by non-additive genetic effects (i.e., epistasis) as well as by environmental effects. That is, imagine two groups: One with a low vPGS and one with a high vPGS. It is the lower-vPGS group that may have a narrow range around a population mean BMI of 25 (say, SD = 2 units), while those with a high vPGS may display the same mean BMI in their group (25) but have a wider dispersion (SD = 4 units). Thus, high vPGS individuals do a better job of buffering differences in environments they encounter and, as a result, their phenotypes vary more widely. In one sense of "plasticity" they are higher, as their phenotype varies more.

Turning to our within-person analysis in the HRS, low vPGS-scoring individuals may be more "stable" in their weight from year to year, given the smaller "error" term from the levels regression. Indeed, when we simply compare the standard deviation of BMI within individuals in our HRS sample across waves by a quartile split in the vPGS, irrespective of treatment status, we find that individuals in the lowest quartile of the vPGS score display a (non-significantly) lower within-person standard deviation in BMI than individuals in the highest quartile (Supplementary Table S5 ). Thus, in the absence of a specific, measured environmental shock, an individual with a low vPGS score derived from the cross-person discovery analysis has more phenotypic stability in the within-person analysis (i.e., less plasticity). This may be one reason why we are able to detect significant differences between treated and control individuals up to two years post job loss in the lower half of the vPGS distribution-i.e., in the presence of an unexpected job loss, it takes less plastic individuals longer to recover to their pre-job loss weight. On the other hand, it could be that less plastic individuals display a more stable weight trajectory overall-regardless of environmental forces, rendering our ability to see a "specific" environmental effect more clearly. Finally, it is also possible that higher vPGS individuals are less affected by environmental exposures overall because their BMI is more under polygenic control than those with lower vPGS.

Either way, with only 375 individuals in our treatment group we are likely underpowered to detect precise G × E effects between treatment and control groups, particularly in the upper part of the vPGS distribution. Specifically, because high vPGS individuals display a higher within-person standard deviation in BMI, it may be harder to detect differences between treatment and control groups in smaller samples because their weight oscillates more between waves, independent of any particular treatment effect, as mentioned above. Conversely, it's possible that the true shape of the interaction does not display a crossover effect at higher levels of the vPGSperhaps because of the larger within-or between-person standard deviation in BMI. Due to a lack of detailed job loss data in other population studies that also collect genetic data on participants, we were unable to pursue replication of our quasi-experimental approach in other samples. Thus, we caution that our results are suggestive, and we cannot draw any definitive conclusions about the short-and long-term dynamics of changes in BMI as they relate to job loss from this study.

There are several limitations of the HRS data, all of which may bias our estimates downwards or reduce the precision of our estimates. First, we only observe BMI in the HRS every two years, which makes it difficult to assess stress-related changes in BMI that are more proximal to the timing of the event, or in the months immediately pre-and post-job loss. It is entirely possible that high vPGS individuals gained or lost more weight than low vPGS individuals in the months following a job loss but they bounced back quicker to their pre-job loss weight, which would make it more difficult to detect differences between high vPGS treatment and control individuals in the subsequent HRS wave. Second, to obtain the largest sample of treated individuals, we were limited to using self-reports of BMI, which may induce measurement error in our estimates. In 2006, the HRS did start collecting in person, objective measures of BMI; however, these measures are only collected at every other wave, or every four years, and are not available for all participants. Third, because the HRS is a sample of older individuals who were genotyped in 2006, 2008, or 2010, our results may be subject to mortality selection. To reduce the potential of mortality selection, we limited our analyses to individuals born after 1930 64 .

In addition, there is significant complexity surrounding obesity and aging such that differences in BMI may not indicate an actual change in body fat. Higher BMI at midlife is a risk factor for age-related disease and early mortality; however, at older ages it might be somewhat protective of mortality because age-related diseases and aging itself are wasting conditions that stimulate weight loss. Therefore, while incrementally higher BMI in midlife is more likely a measure of risk for disease, later in life it may actually signal the absence of disease. In addition, individuals generally lose muscle mass with increasing chronological age, meaning older individuals could maintain a constant BMI while simultaneously losing lean body mass and gaining a greater portion of adiposity 65 . Thus, any increases in BMI from a job loss may be offset by these other countervailing trends that Scientific Reports | (2021) 11:7647 | https://doi.org/10.1038/s41598-021-86716-y www.nature.com/scientificreports/ are inherent to the aging process, which may also explain the null findings we report for more plastic individuals in the top half of the vPGS distribution. Furthermore, the relatively nominal findings we report may reflect a greater culmination of environmental and lifestyle factors on adiposity in older adults that overwhelm any genetic effects. The genomic influence on BMI has been shown to both weaken over the life course and increase in magnitude since the current obesity epidemic began in the mid-1980s 55, 59, 60, 66 . Finally, a significant limitation of this study is we were limited to conducting analyses in individuals of European decent. We focus on individuals of European decent because comparable GWAS in other ancestral populations are currently unavailable. Estimates from a European ancestry GWAS are not necessarily accurate or valid in other ancestral populations, and PGSs constructed from European ancestry GWAS summary statistics will not have the same predictive power for individuals from other ancestral backgrounds 67, 68 . Restricting our analysis to one ancestral group is also important because SNPs within regions of interest may tag different causal variants if the underlying linkage disequilibrium (LD) structure varies across ancestral groups 68, 69 . Thus, we caution that PGS constructed from European ancestry GWAS cannot be generalized to other ancestral populations. On the environmental side, limiting our analysis to white, non-Hispanic HRS respondents restricts the scope of potential job loss effects that we can observe. Race powerfully shapes structural-and institutionally-derived differences in occupational sorting and occupational opportunities across the life course 70, 71 . For example, white HRS respondents were more likely to work in higher status jobs with better working conditions than their Black counterparts 72 , and following a job loss, they were more likely to be reemployed or have additional economic resources to buffer stressful declines in income [73] [74] [75] , all of which may further bias effects from this study downwards.

These limitations are counterbalanced by several strengths of our study. The use of a large, nationally representative cohort of individuals from the same ancestry group is an advantage in that it both increases our power to detect effects while also minimizing the presence of ascertainment bias and other selection issues. Having access to detailed, longitudinal job loss data in the HRS also allowed us to exploit a quasi-experimental research design that limited the treatment group to individuals who lost their job due to a business closure while also creating a control group that is matched on a rich suite of pre-job loss characteristics. Current G × E interaction studies that utilize population data are often unable to separate gene-environment correlation (rGE) from G × E effects, which limits our understanding of social-environmental effects on health 27 . Finally, to our knowledge, this is one of the first studies to integrate genetic measures that can separately capture phenotypic mean and variance effects into PGS × E interaction analysis.

erty of biological systems that impacts how species adapt to environmental changes [76] [77] [78] . Incorporating vPGS measures into G × E interaction research may further our understanding of how and to what extent environmental stimuli modify the distribution of anthropomorphic traits in a population. In particular, sizable unemployment shocks from the Great Recession and the COVID-19 pandemic have highlighted the importance of understanding the short-and long-run health consequences of business cycles. Future studies that are able to observe the biology underlying these types of large, social-environmental effects on physiological changes that precede disease promises to inform new opportunities for effective intervention 79 .

Standardized bias estimates. The standardized bias compares the distance between the marginal distributions, or the difference in sample means between the treated X T and matched control X C subsamples as a percentage of the square root of the average of the sample variances in both groups for a covariate X 56 :

Health and Retirement Study (HRS) data. The HRS is a nationally representative, longitudinal panel study of individuals over the age of 50 and their spouses that is sponsored by the National Institute on Aging (NIA U01AG009740) and conducted by the University of Michigan 80, 81 . Launched in 1992, the HRS introduces a new cohort of participants every six years and interviews around 20,000 participants every two years. To maximize sample size, we compiled data from 13 waves (1992-2016). Information on job loss and smoking behavior was obtained from the 1992-2016 Public Use Core Files; demographic and socioeconomic data came from the RAND Data File (version P).

Genotype data on ~ 15,000 participants was collected from a random subset of the ~ 26,000 total participants that were selected to participate in enhanced face-to-face interviews and saliva specimen collection for DNA in 2006, 2008, and 2010. We restricted our sample to individuals of European ancestry who were between the ages of 50-70 who reported working part-time or full-time in the previous wave and who were not self-employed. The final sample consists of 3939 workers with 11,934 observations.

Research (CIDR) in 2011, 2012, and 2015 (RC2 AG0336495, RC4 AG039029). Full quality control details can be found in the Quality Control Report 82 . Genotype data on over 15,000 participants was obtained using the llumina HumanOmni2.5 BeadChips (HumanOmni2.5-4v1, HumanOmni2.5-8v1), which measures ~ 2.4 million SNPs. Genotyping quality control was performed by the Genetics Coordinating Center at the University of Washington, Seattle, WA. Individuals with missing call rates > 2%, SNPs with call rates < 98%, HWE p value < 0.0001, www.nature.com/scientificreports/ chromosomal anomalies, and first-degree relatives in the HRS were removed. Imputation to 1000G Phase I v3 (released March 2012) was performed using SHAPEIT2 followed by IMPUTE2. The worldwide reference panel of all 1092 samples from the Phase I integrated variant set was used. These imputation analyses were performed and documented by the Genetics Coordinating Center at the University of Washington, Seattle, WA. All positions and names are aligned to build GRCh37/hg19. Principal component (PC) analysis was performed on a selected set of independent SNPs to identify population group outliers and to provide sample eigenvectors as covariates in the statistical model to adjust for possible population stratification and were provided by the HRS. The European ancestry sample included all respondents that had PC loadings within ± one standard deviations for eigenvectors one and two in the PC analysis of all unrelated study subjects and who self-identified as White on survey data. A second set of principal components was then calculated within the European ancestry sample to further account for any population stratification within the sample. The genotype sample has been defined by the HRS and is available on dbGaP 83 .

Mean polygenic score (mPGS) construction and performance. We calculated a linear mPGS for the HRS sample based on a GWAS of 457,824 European ancestry individuals in the UK Biobank 48 . Imputed HRS genotype data were accessed through dbGap (phs000428). The mPGS BMI score was constructed in PRSice 84 by taking a weighted sum across the number of SNPs (n) of the number of reference alleles x (zero, one, or two) at that SNP multiplied by the effect size for that SNP (β):

GWAS summary statistics were pruned for linkage disequilibrium (LD) using the clumping procedure in PLINK (R 2 = 0.1, range = 1000 kb) 85, 86 . Since these GWAS summary statistics were pre-clumped, no LD-clumping or p value threshold was implemented in PRSice. After LD clumping was applied, 90,326 SNPs were used to construct the BMI mPGS. The mPGS was standardized to have a mean of zero and a standard deviation of one for all analyses.

To verify the performance of the mPGS in the HRS European ancestry sample (N = 10,550), we leveraged the longitudinal nature of the BMI data in the HRS. We first fitted a multilevel linear growth curve model on BMI and age:

where Y it and Age it denote the BMI and age of respondent i at time point t , respectively ( i = 1, . . . , n and t = 1, . . . , T i ). We included linear and quadratic terms for age to reflect the non-linear age-dependent trajectory of BMI. Next, we used linear regression to evaluate predictive performance of mPGS on individual intercepts (i.e., β 0i ) estimated in the level one model described above. We also adjusted for the effect of sex and the first 10 genetic principal components. The mPGS has a predictive R 2 of 7.2%.

We calculated BMI vPGS for HRS participants of European ancestry. SNP weights in the vPGS were based on dispersion effects estimated in the UKB using the heteroskedastic linear mixed model (HLMM) approach 16 . Pre-pruned HLMM summary statistics were obtained from Young et al. 16 . We did not perform additional LD-clumping or p value thresholding to filter variants. A total of 242,870 SNPs remained in the vPGS model after overlapping the HLMM summary statistics and HRS genotype data. The vPGS was constructed in PRSice 84 and standardized to have a mean of zero and a standard deviation of one for all analyses.

European ancestry UKB samples identified from genetic PCs (data field 22,006). To avoid overfitting, HRS samples were not used for model validation. Quality control procedures for the UKB genetic data have been described elsewhere 87 . We excluded participants recommended by UKB (data field 22,010), those with conflicting genetically-inferred (data field 22,001) and self-reported sex (data field 31), and those who withdrew from the study. We randomly apportioned UKB participants (N = 406,873) into training (N = 325,498) and testing sets (N = 81,375), with an 80-20 split. We applied the HLMM approach to estimate the dispersion effect of each SNP on BMI using samples in the training set, controlling for sex, age, age 2 , age 3 , age × sex, age 2 × sex, age 3 × sex, genotyping array, and the first 40 genetic PCs. Following Young et al., we analyzed related and unrelated samples in the training set separately and performed fixed-effect meta-analysis to combine the results 16 . Related samples were inferred from genetic kinship (third-degree relatives or higher; data field 22,021). Random effects were included to account for genetic relatedness in the analysis of related samples. We then pruned SNPs following Young et al. and used dispersion effect estimates to generate vPGS for samples in the testing set.

We then fitted a Double Generalized Linear Model (DGLM) to associate the vPGS with the between-individual BMI variance in testing samples 88 . The DGLM takes the form of where BMI i denotes the inverse normal-transformed BMI of individual i , G i is the vPGS of individual i , X i is the vector of covariates including sex, age, age 2 , age 3 , age × sex, age 2 × sex, age 3 × sex, genotyping array, and the first 40 genetic principal components. Here, α 1 quantifies the effect of vPGS on the variability of BMI and is the www.nature.com/scientificreports/ parameter of interest in this analysis. The vPGS was standardized to have a mean of zero and a variance of one for all analyses. We fitted DGLM using the dglm package 89 in R.

To assess the performance of vPGS after adjusting for the effect of mPGS, we performed a standard, meaneffect GWAS of BMI on the training set and used the effect estimates to generate mPGS for the testing samples. GWAS summary statistics were pruned for LD using the clumping procedure in PRSice (R2 = 0.1, range = 250 kb) when calculating the mPGS 84 . We then fitted the same DGLM model as above with the mPGS added to the vector of covariates.

Treatment and control groups. For each observation, we used information from two waves-before and after treatment. Before treatment (t − 2), all respondents were working for pay either full-or part-time. At the following HRS interview two years later (t), respondents in the treatment group report they were no longer working for their previous-wave employer. These respondents were asked why they left their employer. Possible answers included 'business closed' , 'laid off/let go' , 'poor health/disabled' , 'quit' 'family care' , 'better job' , 'retired' , 'family moved' , 'strike' , 'divorce/separation' , 'transportation/distance to work' , and 'early retirement incentive/ offer' . Respondents could report up to three reasons. Our definition of exogenous job loss includes observations that reported being laid off due to a business closure. We excluded workers who also stated that they quit or left for health reasons but included workers who stated as a second reason being laid off, retiring, family care, better job, or ownership change because these circumstances could have occurred concurrently with a business closure (e.g., a worker may have retired because the business closed) 40 . Of those who experienced a job loss, 347 gave business closure as the sole reason for leaving their job. Of the remaining 28 respondents, 13 cited 'laid off/let go' , 8 cited 'retired' , 3 cited 'family care' , 3 cited 'better job' , and 1 cited 'ownership change' as a second reason.

For the control group, we used individuals who reported working for the same employer the entire time they were in the sample-i.e., we did not include individuals in the control group if they ever quit their job or were laid off for any reason. Treated individuals are only in the analytic sample for two waves, or pre-and post-job loss. For control individuals, we used matching with replacement to increase the average quality of matching, which reduces bias 57 . As a result, control individuals can be in the analytic sample for multiple HRS waves. The control group consisted of 11,559 observations on 3564 workers.

Difference-in-differences (DiD) approach. We used DiD estimation combined with nonparametric kernel matching to estimate the average treatment effect on the treated (ATT) by genotype, or the change in BMI by genotype brought about by the job loss of those who actually lost their job 31, 90 . This approach compares individuals who have been laid off due to a business closure with a group of similar individuals who are still working for their same employer. To construct a control group with a similar distribution of covariates as the treatment group, the kernel-based matching estimator uses a distance-weighted average of all propensity scores in the control group to construct a counterfactual outcome for each individual in the treatment group. These weights were applied to the DiD regression model to obtain a balanced sample of treated and untreated individuals. The coefficients from the DiD regression were then used to estimate the ATT by mPGS and vPGS.

A traditional DiD setting assumes that after conditioning on a vector of observables X , the BMI of individuals in the treatment group would have evolved similarly over time to the BMI of individuals in the control group if they had never been laid off:

where BMI it − BMI it−2 refers to the change in BMI before and after the treatment, BC denotes the treatment group indicator (i.e., whether an individual lost their job due to a business closure), and i ′ denotes an individual in the control group with the same characteristics as individual i in the treatment group. While conditioning on genotype and a rich set of covariates minimizes the possibility of violating this assumption, other systematic differences between the treated and control groups may remain even after conditioning on observables.

To minimize potential confounding from unobservable characteristics, we used the weights from propensity score matching (W) to reduce unmeasured differences between the treatment and control groups that could bias estimates:

Covariates used to estimate the propensity score, or the probability of treatment, were also included in the DiD regression model. Thus, coefficients from the regression-adjusted semiparametric DiD matching estimator are considered "doubly robust" because the estimator is consistent if the regression model or the propensity score model is correctly specified 91, 92 . As a result, the DiD matching estimator accounts for selection on observable and unobservable variables with time invariant effects, or the model allows for systematic differences between treatment and control groups even after conditioning on observables 93 .

Difference-in-differences (DiD) empirical strategy. Our empirical strategy can be broken down into three parts. First, we estimated propensity scores using a probit regression that regresses business closures on the mPGS and vPGS, as well as a rich set of covariates that are both standard in the job loss literature and satisfy the conditional independence assumption-i.e. they influence job loss and/or changes in BMI 90, 94 . In addition, we only conditioned on observables that were unaffected by job loss (or the anticipation of it), or variables that were either fixed over time or measured in t − 2 57 . A complete list of covariates can be found in Supplementary Tables S1 and S2. To avoid losing observations with missing information on a covariate, we set missing values equal to zero and included an additional dichotomous variable that is equal to one if the observation is missing. www.nature.com/scientificreports/ As a result, matching is not only on observed values but also on the missing data pattern 31, 95 . Throughout, we restricted our analysis to the region of common support, or the subset of individuals in the control group that were comparable to individuals in the treatment group 94 . Specifically, we dropped treatment observations whose propensity score was greater than the maximum or less than the minimum propensity score of the controls. We used the estimates from the probit regression to compute the weights for the control group with kernel matching, a nonparametric matching estimator that uses the weighted averages of all observations on common support to construct the counterfactual outcome 90, 93 . Specifically, the weight given to a non-treated individual j was in proportion to the closeness of their observables to treated individual i:

where P is the propensity score for individual i or j in the treated or control group, respectively, K[·] is the kernel function, and b is the bandwidth parameter. We used the program psmatch2 96 in Stata 14 to compute w i, j with the Epanechnikov kernel function and a bandwidth of 0.06 90 . In addition, when computing the weights, we performed exact matching on survey year and sex in t − 2. This ensured 1) individuals who were laid off were matched with controls from the same time period, and 2) treated individuals were grouped with same-sex nontreated individuals.

In the final step, we incorporated the weights from propensity score matching into the DiD regression model:

where BC is an indicator for job loss due to a business closing in the years between HRS survey waves, or between t − 2 and t for individual i , X is a vector of observable time invariant and variant covariates measured at t − 2, including the first 10 principal components of the genetic data. We also include BMI t−2 to control for baseline BMI, or to estimate deviations in BMI between t − 2 and t. All regressions were estimated with robust standard errors.

Estimating average treatment effects by genotype. Estimated parameters from the DiD regression model were used to estimate the conditional mean or predicted BMI for treated and untreated individuals at various values of the mPGS and vPGS (Fig. 2) . For example, the BMI for treated (BC = 1) and untreated (BC = 0) individuals with hypothetical mPGS and vPGS values at 0 and 1, respectively, would be estimated as follows:

From here, the ATT can be estimated by taking the difference in E[BMI t |W(X)] between treated and nontreated individuals with the same mPGS and vPGS values: Event time study analysis. We estimated an event time study (ETS) model for individuals in the top and bottom 50% of the vPGS distribution using the following specification:

This model is similar to the DiD model outlined above except the business closure term is replaced by a series of event terms that are the product of indicators for each HRS survey year y relative to the survey year the respondent reported a job loss t * i , I t − t * i = y , and their treatment status (BC i ) . The omitted category is the survey year prior to treatment y = −2 . We also present ETS results for the full sample that includes controls for the vPGS. Each estimate of y gives the difference in BMI for treated individuals compared to non-treated individuals relative to the excluded year. If outcomes were evolving similarly for treated and untreated individuals prior to a business closure, the coefficient estimates for y < 0 should be close to zero and not statistically significant.

Health and Retirement Study (HRS) phenotypic data is publicly available on the HRS website: https:// hrs. isr. umich. edu/ data-produ cts. HRS genotype data is available through the NCBI Database of Genotypes and Phenotypes (dbGaP): https:// www. ncbi. nlm. nih. gov/ gap/. UK Biobank phenotype and genotype data is publicly available through their Access Management System (AMS) after applying for access: https:// www. ukbio bank. ac. uk/ enable-your-resea rch/ apply-for-access. 

Unemployment rates by age, sex, and marital status, seasonally adjusted. Labor Force Statistics from the Current Population Survey

A first in nearly 50 years, older workers face higher unemployment than midcareer workers

Can unemployed older workers find work

Older Workers, Retirement, and the Great Recession 1-7 (Russell Sage Foundation

What makes retirees happy

The persistence of depressive symptoms in older workers who experience involuntary job loss results from the health and retirement survey

Retirement and subjective well-being

Intermittent lack of health insurance coverage and use of preventive services

Rough passage: affordable health coverage for near-elderly Americans

Earnings losses of displaced workers

The effects of job displacement on job quality: findings from the Wisconsin longitudinal study

Job loss in the United States

Unemployment and change of tobacco habits: a study of young people from 16 to 21 years of age

Stress-related eating and drinking behavior and body mass index and predictors of this behavior

Current challenges and new opportunities for gene-environment interaction studies of complex diseases

Identifying loci affecting trait variability and detecting interactions in genome-wide association studies

Interactions between polygenic scores and environments: methodological and conceptual challenges

Power and predictive accuracy of polygenic risk scores

Recent developments in statistical methods for detecting genetic loci affecting phenotypic variability

FTO genotype is associated with phenotypic variability of body mass index

A sibling method for identifying vQTLs

Leveraging phenotypic variability to identify genetic interactions in human phenotypes

Genotype-by-environment interactions inferred from genetic effects on phenotypic variability in the UK

Polygenic scores for plasticity: a new tool for studying gene-environment interplay

The promise and challenges of incorporating genetic data into longitudinal social science surveys and research. Biodemogr

The challenge of causal inference in gene-environment interaction research: leveraging research designs from the social sciences

Modeling gene-environment interactions with quasi-natural experiments

Toward a better estimation of the effect of job loss on health

Job loss from poor health, smoking and obesity: a national prospective survey in France

Losing life and livelihood: a systematic review and meta-analysis of unemployment and all-cause mortality

Does job loss make you smoke and gain weight?

Job loss, firm-level heterogeneity and mortality: evidence from administrative data

Comparison of the associations of body mass index and measures of central adiposity and fat mass with coronary heart disease, diabetes, and all-cause mortality: a study using data from 4 UK cohorts

Frailty in older adults: evidence for a phenotype

Frailty and failure to thrive

Toward an understanding of frailty

Frailty: an emerging public health priority

Frailty consensus: a call to action

The effect of job loss on overweight and drinking

Does job loss cause ill health

The effect of job loss on health: evidence from biomarkers

Quantile treatment effects of job loss on health

Involuntary job loss as a risk factor for subsequent myocardial infarction and stroke: findings from the health and retirement survey

Association analyses of 249,796 individuals reveal 18 new loci associated with body mass index

Genetic studies of body mass index yield new insights for obesity biology

New genetic loci link adipose and insulin biology to body fat distribution

Meta-analysis of genome-wide association studies for height and body mass index in ∼700000 individuals of European ancestry

Mixed-model association for biobank-scale datasets

The bigger picture of FTO-the first GWAS-identified obesity gene

FTO obesity variant circuitry and adipocyte browning in humans

Polygenic prediction of weight and obesity trajectories from birth to adulthood

The influence of age and sex on genetic associations with adult body size and shape: a large-scale genome-wide interaction study

Gene-environment interactions related to body mass: school policies and social context as environmental moderators

Lifetime socioeconomic status, historical context, and genetic inheritance in shaping body mass in middle and late adulthood

The genome-wide influence on human BMI depends on physical activity, life course, and historical period

The bias due to incomplete matching

Some practical guidance for the implementation of propensity score matching

Large sample properties of matching estimators for average treatment effects

Changing polygenic penetrance on phenotypes in the 20th century among adults in the US population

Association of a genetic risk score with body mass index across different birth cohorts

Controlling the false discovery rate: a practical and powerful approach to multiple testing

Phenotypic plasticity and the origins of diversity

The far-reaching impact of job loss and unemployment

Mortality selection in a genetic sample and implications for association studies

Association of a genetic risk score with BMI along the life-cycle: evidence from several US cohorts

Generalization and dilution of association results from european GWAS in populations of non-European ancestry: the PAGE study

Human demographic history impacts genetic risk prediction across diverse populations

Genome-wide association studies in diverse populations

The structure of disadvantage: individual and occupational determinants of the black-white wage gap

The Extent of occupational segregation in the United States: differences by race, ethnicity, and gender

Structural racism in the workplace: does perception matter for health inequalities?

Jobs lost, jobs regained: an analysis of black/white differences in job displacement in the 1980s

After the rainy day: how private resources shape personal trajectories following job loss and amplify racial inequality

Family adaptations to income and job loss in the U.S

Biological robustness

Phenotypic plasticity's impacts on diversification and speciation

Biological robustness: paradigms, mechanisms, and systems principles

Social determinants of health and survival in humans and other animals

An overview of the health and retirement study

Cohort profile: the health and retirement study (HRS)

Quality control report for genotypic data

Health and retirement study (HRS) dbGaP study accession: phs000428.v2.p2

PRSice: polygenic risk score software

Second-generation PLINK: rising to the challenge of larger and richer datasets

The UK Biobank resource with deep phenotyping and genomic data

Genetic heterogeneity of residual variance -estimation of variance components using double hierarchical generalized linear models

Matching as an econometric evaluation estimator: evidence from evaluating a job training programme

Toward a curse of dimensionality appropriate (CODA) asymptotic theory for semi-parametric models

Nonparametric estimation of average treatment effects under exogeneity: a review

Does matching overcome LaLonde's critique of nonexperimental estimators

Causal effects in nonexperimental studies: reevaluating the evaluation of training programs

Matching methods for causal inference: a review and a look forward

psmatch2: stata module to perform full Mahalanobis and propensity score matching, common support graphing, and covariate imbalance testing

This research was funded by the following awards from the National Institute on Aging (NIA): K99 AG056599, R00 AG056599, P30 AG012846, T32 AG000221 (Schmitz), and P30 AG017266 (Lu and Schmitz). The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. The authors would like to thank Jason Fletcher, Rebecca Johnson, and Ramina Sotoudeh for helpful comments and feedback on earlier versions of this draft.

L.L.S. and D.C. conceived the project and designed the study. L.L.S. and J.G. performed the data analysis in the HRS J.M. constructed the mPGS and vPGS in the HRS and UKB and carried out the validation analysis in the UKB. Q.L. oversaw all analysis in the UKB. L.L.S. drafted the manuscript and all authors contributed to critically reviewing and revising the draft. All authors read and approved the final manuscript before submission.

The authors declare no competing interests.

Supplementary Information The online version contains supplementary material available at https:// doi. org/ 10. 1038/ s41598-021-86716-y.Correspondence and requests for materials should be addressed to L.L.S.Reprints and permissions information is available at www.nature.com/reprints.Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.