The impact of late-career job loss and genetic risk on body mass index: Evidence from variance polygenic scores.

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.

Introduction

Recent unemployment shocks from the COVID-19 pandemic have left millions of older workers unemployed. In the U.S. alone, the seasonally adjusted unemployment rate for adults aged 55 and over jumped from 2.6 in February 2020 to as high as 13.6 in April 2020[1]. Recent evidence indicates unemployment rates for workers 55 and older have exceeded those of mid-career workers since the pandemic began—the first time in nearly 50 years that older workers have faced higher unemployment than mid-career workers[2]. For older workers in particular, the scarring effects of unexpected job loss could be severe. Job loss at older ages has been associated with longer periods of unemployment than any other age group[3,4], higher rates of depression and anxiety[5-7], and a sharp increase in the need for medical care due to heightened stress levels and gaps in health insurance coverage[8,9]. Further, when reemployed, older workers suffer significant wage penalties and lower levels of employer-offered pension and health insurance[10-12]. All these factors could trigger chronic stress and adverse changes in health and health behaviors[13,14].

This study expands on past work by examining the degree to which underlying genetic predisposition moderates changes in body mass index (BMI) after a job loss. Genotype-by-environment (G × E) interaction is a fundamental component of population variance for complex traits like BMI, but there has been limited success in identifying G × E effects in human populations due to several challenges, including the complexity of measuring environmental exposures, a need for statistical methods that can screen for genetic effects on phenotypic variability, and inadequate power to detect small G × E effects at loci across the genome[15-17]. To overcome the third challenge of lower power to detect individual effects in the context of multiple hypothesis testing, researchers have used whole-genome polygenic scores (PGSs) constructed from well-powered genome-wide association studies (GWAS) that summarize the genomic contribution to a trait or disease across common variants in the genome[18]. That is, PGSs aggregate thousands of genome-wide genetic influences on a phenotype into a single index using results from GWAS that estimate the association between genetic variants and the conditional mean of a phenotype, which we refer to herein as mGWAS. However, PGSs constructed from mGWAS may not capture the impact of loci that contribute to within-individual variance in an outcome that are more responsive to environmental stimuli (i.e., variance quantitative loci or vQTLs). Since estimating genetic contributions to 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-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-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-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-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-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-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-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 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.

Figure 1

Analytic workflow in the UK Biobank and the Health and Retirement Study (HRS).

Results

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].

Table 1

Before treatment means of treated workers who lost their job due to a business closure, control workers, and matched controls.

	Match status	Means		% bias	t-test p value	V(T)/V(C)
		Treated	Control
BMI	U	27.35	27.76	− 7.9	0.137	0.96
	M		27.53	− 3.5	0.631	1.01
BMI mPGS	U	0.03	0.00	3.4	0.510	1.06
	M		0.05	− 1.2	0.873	1.07
BMI vPGS	U	0.01	0.00	0.9	0.863	1.00
	M		0.04	− 3.5	0.640	0.94
Femalea	U	0.54	0.56	− 3.9	0.453
	M		0.54	− 0.2	0.974
Age (years)	U	57.43	57.15	6.3	0.212	1.15
	M		57.48	− 1.2	0.869	1.15
No degreea	U	0.11	0.06	18.2	< 0.0001
	M		0.11	− 0.2	0.981
High school degreea,b	U	0.67	0.52	31.3	< 0.0001
	M		0.61	12.3	0.088
College degreea,b	U	0.22	0.42	− 44.3	< 0.0001
	M		0.28	− 13.0	0.061
Household income (log)b	U	11.08	11.34	− 26.3	< 0.0001	1.82
	M		11.17	− 9.1	0.243	1.30
Household wealth ($100 k)	U	3.10	3.64	− 5.2	0.453	0.12
	M		3.33	− 2.2	0.763	0.11
Works part timea	U	0.19	0.11	21.5	< 0.0001
	M		0.18	3.5	0.655
White Collara	U	0.64	0.71	− 15.7	0.002
	M		0.65	− 1.8	0.808
Blue Collara,b	U	0.24	0.18	15.1	0.002
	M		0.22	5.9	0.434
Servicea	U	0.11	0.09	8.3	0.095
	M		0.12	− 4.0	0.613
Job tenure (years)b	U	11.36	15.25	− 35.6	< 0.0001	1.03
	M		12.54	− 10.7	0.135	1.11
Health excellent/very gooda	U	0.62	0.64	− 4.1	0.437
	M		0.61	0.8	0.914
Health insurancea,b	U	0.65	0.82	− 38.2	< 0.0001
	M		0.69	− 10.0	0.203
Exercise vigorously 3 + times/weeka	U	0.35	0.36	− 1.6	0.759	0.99
	M		0.36	− 1.3	0.859	0.99
Ever smoke cigarettesa	U	0.62	0.53	17.6	< 0.0001
	M		0.60	5.0	0.487
Cigarettes per dayb	U	5.30	2.75	27.4	< 0.0001	1.94
	M		4.58	7.8	0.328	1.27
Drinks per week	U	2.24	2.35	− 2.0	0.713	0.78
	M		2.10	2.7	0.709	0.78
CES-D	U	0.99	0.75	15.1	0.002	1.33
	M		0.95	2.5	0.747	1.00

aIndicates a variable is binary. bIndicates a variable is not balanced after matching (p value < 0.05 or standardized bias > 5%). Abbreviations: U, unmatched; M, matched; % bias: percent standardized bias; CES-D, Center for Epidemiological Studies-Depression 8 item scale. The t-test p value refers to the p value from the difference in means between the treated and control groups before and after matching. The V(T)/V(C) column presents the ratio of a continuous variable’s variance for the treatment group over the variance for the control group. Additional covariates in the matching procedure: marital status, firm size, industry, survey year, regional Census division, additional categories for variables with missing values, and the first 10 principal components of the European ancestry genetic data. We used kernel-based propensity score matching with a bandwidth of 0.06. Unmatched control observations = 11,629; Unmatched treated observations = 399; Matched control observations = 11,559; Matched treated observations = 375.

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-55,59,60], and previous GWAS have reported predicted R2 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 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 log-linear 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., ).

Table 2

Variance polygenic score (vPGS) validation in the UK Biobank.

	Model 1	Model 2
vPGS	Beta (SE)	Beta (SE)
	95% CI	95% CI
	p value	p value
	0.019 (0.005)	0.019 (0.005)
	[0.009, 0.029]	[0.009, 0.029]
	1.01E−04	1.44E−04

Abbreviations: SE, standard error; CI, confidence interval. Results show the association between the vPGS and the between-individual variance in BMI in a 20% holdout sample of the UK Biobank. Model 2 adjusts for the effect of the mean PGS (mPGS). Unique N = 81,375.

Gene-by-environment (G × E) interaction results

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).

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).

	(1)	(2)	(3)	(4)	(5)	(6)
	Full sample	European ancestry sample	Genotyped European ancestry sample
			No PGS	With mPGS	With mPGS interaction	With mPGS & vPGS interaction
	Beta (SE)	Beta (SE)	Beta (SE)	Beta (SE)	Beta (SE)	Beta (SE)
	95% CI	95% CI	95% CI	95% CI	95% CI	95% CI
	p value	p value	p value	p value	p value	p value
Business closure (BC)	0.137 (0.088)	0.051 (0.103)	− 0.099 (0.105)	− 0.091 (0.104)	− 0.091 (0.103)	− 0.095 (0.102)
	[− 0.036, 0.310]	[− 0.151, 0.252]	[− 0.305, 0.106]	[− 0.294, 0.113]	[− 0.294, 0.112]	[− 0.296, 0.106]
	0.121	0.623	0.344	0.383	0.379	0.352
mPGS				0.175* (0.051)	0.170* (0.048)	0.163* (0.047)
				[0.075, 0.275]	[0.076, 0.263]	[0.071, 0.255]
				0.001	0.001	0.001
Business closure × mPGS					0.011 (0.091)	0.013 (0.090)
					[− 0.168, 0.189]	[− 0.164, 0.191]
					0.907	0.882
vPGS						− 0.056 (0.070)
						[− 0.193, 0.082]
						0.426
Business closure × vPGS						0.258* (0.092)
						[0.077, 0.439]
						0.001
BMI (t − 2)	0.918 (0.014)	0.929 (0.017)	0.932 (0.024)	0.921 (0.026)	0.921 (0.026)	0.923 (0.025)
	[0.890, 0.946]	[0.896, 0.962]	[0.884, 0.979]	[0.871, 0.971]	[0.871, 0.971]	[0.874, 0.973]
	< 0.0001	< 0.0001	< 0.0001	< 0.0001	< 0.0001	< 0.0001
R-squared	0.82	0.832	0.868	0.869	0.869	0.870
Observations	27,593	18,795	11,941	11,934	11,934	11,934
Treated observations	895	579	376	375	375	375
Control observations	26,698	18,216	11,565	11,559	11,559	11,559

*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.

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.

Figure 2

Predicted BMI for treatment and control groups by mPGS and vPGS. The figure depicts the predicted BMI for treatment and control groups using coefficients from the empirical models reported in Table 3 (Column 5) for the mPGS interaction, and Table 3 (Column 6) for the vPGS interaction. Error bars represent 95% confidence intervals.

Event time study

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.

Figure 3

Difference in predicted BMI for treatment and control groups by the year relative to job loss for the entire analytic sample and stratified at the vPGS median. The figure plots the coefficients from an event time study model for the full sample and by samples stratified at the vPGS median (see Methods). The reference category is BMI in t − 2, or BMI at baseline before the job loss occurred for treated individuals. The dotted line depicts the approximate time point that the job loss occurred—i.e., at some point between t − 2 and the wave a job loss was reported, or time t. Error bars represent 95% confidence intervals.

Discussion

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 literature, this suggests that similar to race, gender, age, or other social moderators that exacerbate the scarring effects of job loss[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 vPGS—perhaps 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.

Strengths and limitations

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 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-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.

Conclusion

Control of phenotypic variability, both within and between individuals, is a fundamental property of biological systems that impacts how species adapt to environmental changes[76-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].

Methods

Standardized bias estimates

The standardized bias compares the distance between the marginal distributions, or the difference in sample means between the treated and matched control subsamples as a percentage of the square root of the average of the sample variances in both groups for a covariate [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.

HRS genotyping and quality control

Genotyping was conducted by the Center for Inherited Disease 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, 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 of the number of reference alleles (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 (R2 = 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 and denote the BMI and age of respondent at time point , respectively ( and ). 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., ) 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 R2 of 7.2%.

Construction of the variance polygenic score (vPGS)

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.

Validation of the variance polygenic score (vPGS)

Performance of the vPGS was assessed using 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, age2, age3, age sex, age2 sex, age3 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 ofwhere denotes the inverse normal-transformed BMI of individual , is the vPGS of individual , is the vector of covariates including sex, age, age2, age3, age sex, age2 sex, age3 sex, genotyping array, and the first 40 genetic principal components. Here, quantifies the effect of vPGS on the variability of BMI and is the 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, mean-effect 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 , 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 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 denotes an individual in the control group with the same characteristics as individual 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. 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 was in proportion to the closeness of their observables to treated individual :where is the propensity score for individual or in the treated or control group, respectively, is the kernel function, and is the bandwidth parameter. We used the program psmatch2[96] in Stata 14 to compute 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 non-treated 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 , 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 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 between treated and non-treated 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 relative to the survey year the respondent reported a job loss and their treatment status . The omitted category is the survey year prior to treatment . We also present ETS results for the full sample that includes controls for the vPGS. Each estimate of 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 should be close to zero and not statistically significant.

Supplementary Information.