Ambient Fine Particulate Matter Air Pollution and Risk of Weight Gain and Obesity in United States Veterans: An Observational Cohort Study

Background: Experimental evidence and studies of children and adolescents suggest that ambient fine particulate matter [particulate matter ≤2.5μm in aerodynamic diameter (PM2.5)] air pollution may be obesogenic, but the relationship between PM2.5 and the risk of body weight gain and obesity in adults is uncertain. Objectives: Our goal was to characterize the association between PM2.5 and the risks of weight gain and obesity. Methods: We followed 3,902,440 U.S. Veterans from 2010 to 2018 (median 8.1 y, interquartile range: 7.3–8.4) and assigned time-updated PM2.5 exposures by linking geocoded residential street addresses with satellite-based estimates of surface-level PM2.5 mass (at ∼1-km2 resolution). Associations with PM2.5 were estimated using Cox proportional hazards models for incident obesity [body mass index (BMI)≥30 kg/m2] and a 10-lb increase in weight relative to baseline and linear mixed models for associations with intra-individual changes in BMI and weight. Results: A 10-μg/m3 higher average annual PM2.5 concentration was associated with risk of incident obesity [n=2,325,769; hazard ratio (HR)=1.08 (95% CI: 1.06, 1.11)] and the risk of a 10-lb (4.54 kg) increase in weight [HR=1.07 (95% CI: 1.06, 1.08)] and with higher intra-individual changes in BMI [0.140 kg/m2 per year (95% CI: 0.139, 0.142)] and weight [0.968 lb/y (95% CI: 0.955, 0.981)]. Nonlinear exposure–response models indicated associations at PM2.5 concentrations below the national standard of 12 μg/m3. As expected, a negative exposure control (ambient air sodium) was not associated with obesity or weight gain. Associations were consistent in direction and magnitude across sensitivity analyses that included alternative outcomes and exposures assigned at different spatial resolutions. Discussion: PM2.5 air pollution was associated with the risk of obesity and weight gain in a large predominantly male cohort of U.S. Veterans. Discussions about health effects of PM2.5 should include its association with obesity, and deliberations about the epidemiology of obesity should consider its association with PM2.5. Investigation in other cohorts will deepen our understanding of the relationship between PM2.5 and weight gain and obesity. https://doi.org/10.1289/EHP7944


Introduction
Experimental evidence in animals suggests that ambient fine particulate matter [particulate matter ≤2:5 lm in aerodynamic diameter (PM 2:5 )] air pollution may be obesogenic (Madrigano et al. 2010;Sun et al. 2009;Xu et al. 2010;Zou 2010). For example, compared with mice exposed by inhalation to filtered air, mice exposed to PM 2:5 exhibited increased subcutaneous and abdominal fat mass (Sun et al. 2009;Xu et al. 2010). Reports in children and adolescents suggest that exposure to higher levels of PM 2:5 may be associated with adiposity and higher risk of weight gain (Bloemsma et al. 2019;de Bont et al. 2019;Jerrett et al. 2010Jerrett et al. , 2014McConnell et al. 2015). However, to our knowledge, the associations between PM 2:5 and the risks of body weight gain and obesity in adults have not been investigated. A greater understanding of the relationship between PM 2:5 air pollution and the risks of weight gain and obesity-upstream risk factors for several noncommunicable diseases, including cardiovascular disease and diabetes-would enhance our understanding of the health effects of air pollution.
We hypothesized that exposure to higher levels of PM 2:5 may increase the risk of weight gain and obesity. In this study, we built a large longitudinal cohort of 3,902,440 U.S. Veterans and followed them for a median of 8.1 y (27,628,688 total personyears) to investigate the relationship between PM 2:5 and the risks of weight gain and obesity.

Cohort Construction
We built a cohort of U.S. Veterans from the U.S. Department of Veterans Affairs (VA) databases (Al-Aly et al. 2012;Bowe et al. 2017Bowe et al. , 2018aBowe et al. , 2018bBowe et al. , 2019aBowe et al. , 2019bXie et al. 2016Xie et al. , 2017. We included Veterans who had at least one weight measurement recorded between 1 July 2010 and 31 June 2011 (n = 4,595,876) ( Figure S1) (Goodloe et al. 2017;Muthalagu et al. 2014;Noël et al. 2010). We then excluded weight measurements of <75 lb (34:0 kg) or >700 lb (317:5 kg) (resulting in 4,595,342 unique Veterans). Weight measurements that exceeded any other weight measurement within the year before to a year after its date of record by ≥50 lb (22:7 kg) were removed from measures in this time period (Goodloe et al. 2017). Of the resulting weight measurements left after application of these criteria (Veterans = 4,543,894), the data of the first weight measurement was set as time zero (T 0 ). To obtain a height value, we selected all heights for an individual Veteran that was available in VA records from 1999 to 2019. We then excluded height measurements of <48 or >84 in (<121:92 or >213:36 cm) and height measurements that varied more than 2 in (5:08 cm) from the median value of all of the individual Veteran's height measurements (Muthalagu et al. 2014); the median height was then taken, and those with no height measure were excluded (resulting in 4,528,288 Veterans remaining in the cohort).
Overall, 1.5% of Veterans were excluded from the cohort after applying data-cleaning criteria to weight and height measurements ( Figure S1). Weight and height measurements were used to calculate body mass index (BMI) {[ðweight in pounds × 703Þ divided by the square of the height in inches] (or weight in kilograms divided by the square of the height in meters)}, and Veterans with a BMI value outside the range of 10-80 were excluded (n = 4,527,985) (Noël et al. 2010). We finally selected those Veterans who had data on all individual and contextual covariates and were linkable at baseline to PM 2:5 data, excluding those where PM 2:5 was not available for their residential location in the year prior to T 0 (e.g., Alaska), yielding an analytic cohort of 3,902,440 ( Figure S1). Veterans were followed from T 0 until death or end of follow-up, 31 December 2018, and were interval censored when their residential location was not linkable to PM 2:5 data, such as in instances where they moved outside the contiguous United States. Mortality and the date of death were obtained from the VA Vital Status file (Maynard C 2017). This study was approved by the Saint Louis Veterans Affairs Health Care System institutional review board (IRB). A waiver of informed consent was granted by the IRB because it determined that this study involved no more than minimal risk, was deemed to not adversely affect the rights and welfare of studied Veterans, and could not have been practically carried out without waivers of consent. To protect the privacy and security of study information, all study members underwent Collaborative Institutional Training Initiative training for the protection of human subjects in biomedical research, Privacy and Health Insurance Portability and Accountability Act of 1996 training from the VA, and VA Privacy and Information Security Awareness and Rules of Behavior training and followed the guidelines set forth by the Veterans Health Administration (VHA) Handbook 1200.12, Use of Data and Data Repositories in VHA Research (Department of Veterans Affairs 2009). Identifying protected health information was limited to only information required to conduct the study, did not include any names or social security numbers, and will not appear in any presentations of publications of the results of this study. In addition, all study data was accessed and stored only on VA servers and maintained behind a VA firewall, where data usage is routinely audited to ensure compliance with VHA data storage guidelines.

Exposure Data
Validated satellite-based (V4.NA.04.MAPLE) estimates of ground-level PM 2:5 over the contiguous United States were derived using satellite remote sensing, chemical transport modeling, and ground-based PM 2:5 measurement data (R 2 = 0:70 with ground-based observations) (van Donkelaar 2019). Measured PM 2:5 concentrations were obtained from the U.S. Environmental Protection Agency (U.S. EPA 2020). Estimates of annual average PM 2:5 were available for 0:01 × 0:01 grids ( ∼ 1 km 2 ) for each year from 2009 to 2017. Geocoded street addresses and county and state of residence from 1 July 2009 throughout the follow-up period were available for studied Veterans from the U.S. Department of Veterans Affairs Planning Systems Support Group Enrollee File (U.S. Department of Veterans Affairs Information Resource Center 2015). Latitude and longitude coordinates for each geographic grid were matched with the latitude and longitude of each Veteran's residence throughout follow-up. In VA data, 90% of Veteran-geocoded residential locations were based on mapping to the Veteran's street address or an interpolation based on the building number and the location of their street. Otherwise, addresses were geocoded to the ZIP code tabulation area centroid. In primary analyses, each individual's residential address in the previous year was linked to the average annual PM 2:5 in that location for that year Lepeule et al. 2012). Residential addresses were updated every quarter of a year during follow-up to account for changes in PM 2:5 exposures over time and according to residential location. In all analyses, because satellitebased estimates in the very low range of PM 2:5 were sparse and in some instances resulted in null and negative exposure values, we excluded annual PM 2:5 concentrations that were <2:4 lg=m 3 (the 0.1th percentile of all PM 2:5 values during follow-up) from exposure data for an individual, thus interval censoring Veterans during these time points.

Outcomes
Height and weight, as measured and recorded for each Veteran during routine clinical encounters, were available from the VA MedSAS inpatient and outpatient data sets (Bowe et al. 2021;Murphy et al. 2002;Oddone and Eisen 2008 Xie et al. 2020). As primary outcomes, we investigated time until obesity (BMI ≥30 kg=m 2 ) among study Veterans whose baseline BMI was <30 kg=m 2 , and time until a 10-lb (4:54 kg) increase in weight from each Veteran's baseline value. Events were assumed to occur on the date when a post-baseline BMI of ≥30 or a 10-lb (4:54 kg) weight gain was first recorded in a Veteran's record. To better account for intra-personal characteristics, we also investigated the intra-individual change in BMI and weight over time. The timing of changes in weight or BMI was determined by the date on which each measurement was recorded. Similar to data-cleaning procedures employed in the cohort construction, we excluded weight measurements of <75 lb (34 kg) or >700 lb (317:5 kg) and excluded weight measurements that exceeded any other weight measurement within the year before to a year after its date of record by ≥50 lb (22:7 kg) (Goodloe et al. 2017). The value assigned as the baseline height was used as the height value for outcome BMI assessment.

Covariates
Individual-level covariates. State and county of residence, age, race, sex, and marital status at baseline were collected from the VA MedSAS (Murphy et al. 2002) inpatient and outpatient data sets. Self-reported race was identified based on encounter data and classified as White, Black or African American, or other (including American Indian or Alaska Native, Asian, Native Hawaiian or other Pacific Islander, and multiracial Veterans). If self-reported race differed across different VA encounters, the most frequently recorded classification was used. Missing self-reported race was supplemented by information from other data sources, including the Medicare and the Beneficiary Identification Records Locator Subsystem (Bowe et al. 2018b(Bowe et al. , 2019bXie et al. 2019). Ethnicity was not assessed. Data from the VA Corporate Data Warehouse Health Factors domain were used to classify baseline smoking status based on work by McGinnis et al. (2011). In brief, an algorithm validated against self-reported smoking survey data was used to convert text entries from health factors data into categories of never, former, and current smoker (for example, a Veteran whose VA record included text indicating former or previous smoking would be classified as a former smoker.) The report prior but closest to baseline was used to classify smoking; Veterans with missing smoking status data were classified as never smokers. BMI categories, used in reporting of cohort characteristics at baseline, were defined as normal or underweight (BMI <25 kg=m 2 ), overweight (BMI between 25 and <30 kg=m 2 ), and obese (BMI ≥30 kg=m 2 ); BMI was otherwise treated as continuous.
Contextual covariates (Table S1). Population density at the county level was obtained from the U.S. Census Bureau's Small Area Income and Poverty Estimates data (Bell et al. 2016). Diet (percentage of population with limited access to healthy foods), exercise (percentage of population with adequate access to exercise opportunities), excessive alcohol use (percentage of adults reporting excessive consumption), and rurality (percentage of population living in a rural area) were obtained from the County Health Rankings (Remington et al. 2015). The 2015 area deprivation index (ADI) was obtained from the University of Wisconsin (University of Wisconsin School of Medicine and Public Health 2015). The ADI allows for the rankings of census block groups by socioeconomic status disadvantage and is a composite measure of education, employment, housing quality, and poverty measures (Kind and Buckingham 2018; University of Wisconsin School of Medicine and Public Health 2015). The normalized difference vegetation index (NDVI), an indicator for green space, was obtained from the National Oceanic and Atmospheric Administration Climate Data Record of Advanced Very High-Resolution Radiometer Surface Reflectance (Vermote 2019). The NDVI for each residential location was classified for 0:05 × 0:05 geographic grids based on the mean of the top three highest daily values during July 2010 (Crouse et al. 2017;James et al. 2016;Orioli et al. 2019;Vienneau et al. 2017). The NDVI measures surface vegetation by comparing red and near-infrared spectral bands and ranges from −1 (water), to 0 (bare ground), and, maximally, to 1 (dense vegetation). Contextual covariates were updated every 3 months to reflect changes in residential location during follow-up.

Statistical Analyses
A conceptual framework for the association between PM 2:5 and obesity, including relations between exposure, covariates, and outcomes, is provided in Figure S2.
Incident obesity and weight gain. We constructed Cox proportional hazards models to estimate hazard ratios (HRs) for incident obesity (among Veterans with a BMI of <30 at baseline only) and a 10-lb (4:54 kg) increase in weight relative to baseline (among all Veterans). Models with a linear functional form of the annual PM 2:5 concentration in the previous year (averaged over all residences in the previous year) were sequentially adjusted for a) baseline BMI for the obesity outcome, and weight and height for the weight outcome to account for differences in time to outcomes that may result from differences in baseline values; b) State of residence to account for potential differences in state-level PM 2:5 composition and potential differences in state-level policies and other factors that may also affect the outcomes (Chen et al. 2020); c) age, race, and sex; d) contextual characteristics to account for the influence of the broader contextual milieu-socioeconomic deprivation, residential greenness, population density, rurality, and other measures-which may correlate with both exposure and outcomes; and e) smoking status-a putative confounder of the association under examination. All continuous covariates (baseline BMI or weight and height, age, and the contextual factors) were modeled as restricted cubic splines unless otherwise indicated.
In separate analyses we also explored nonlinearity by constructing an ensembled estimate of the risk relationship in our cohort (Nasari et al. 2016). Specifically, we constructed Cox proportional hazards models of the form: where k 0 ðtÞ is the baseline hazard at time t; and c 0 is the vector of coefficients for covariate vector x. b is the coefficient for transformed PM 2:5 , which is the product of wðzjl,sÞ, a logistic weighting function, and f(z), which is either the identity or natural log transformation of the PM 2:5 concentration, z. The logistic weighting function was as follows: where l is a location parameter representing a percentile of the PM 2:5 distribution; s is a parameter that controls the curvature of the weighting function; and r is the range of PM 2:5 concentrations. This logistic weighting function specifies a monotonic sigmoidal hazard function for the relation between PM 2:5 and the outcome, with larger values of tau resulting in less curvature, and the value of l determining the location of the inflection point or maximum slope along the PM 2:5 distribution. Multiple models were constructed by adjusting the parameters l, s, and f(z), until optimal fit was achieved. The three models that best fit the data were averaged, weighted by their model fit (log-likelihood), and used to construct an ensembled estimate. We present both the ensembled estimates and the estimates from the best-fitting model for each outcome. Ranges of the PM 2:5 concentrations covered in the analyses of obesity were from 3:308 to 14:655 lg=m 3 and from 3:324 to 14:645 lg=m 3 in the analyses of weight gain, where the difference is the value of r. Ranges were determined by the range of PM 2:5 concentrations experienced by Veterans included in each analysis (with models of incident obesity restricted to Veterans with BMI of <30 at baseline), where the distribution included one observation per time period per cohort member at risk in the cohort and, thus, differed by corresponding cohort inclusion/exclusion criteria (Nasari et al. 2016). In addition, maximum values were excluded beyond the 99th percentile and minimum values were excluded below the 1st percentile. Models of both outcomes included two possible values of τ, 0.1 and 0.2, which were selected a priori based on prior examinations by the method authors on their influence on variation in the shape of possible curves, which suggested that values beyond these two did not appreciably change examined curves' forms (Nasari et al. 2016). Finally, we modeled seven possible values of l corresponding to the -5th, 0th, 5th, 10th, 25th, 50th, and 75th percentiles of the PM 2:5 distributions for each set of study Veterans (1.671,3.327,4.983,5.823,7.159, and 9:872 lg=m 3 for models of incident obesity and 1.752, 3.369, 4.986, 5.827, 7.185, and 9:925 lg=m 3 for models of a 10-lb (4:54 kg) gain in weight). Values for the -5th percentile were derived by subtracting the difference between the 5th and 0th percentile from the value of the 0th percentile. The percentiles used as possible percentiles for l were selected a priori based on work by the authors on the modeling approach that found these values sufficiently covered possible inflection points in the logistic curve across the distribution (Nasari et al. 2016).
Intra-individual changes in BMI and weight over time. We constructed linear mixed models to analyze the relationship between PM 2:5 and intra-individual change in BMI and weight over time (Jerrett et al. 2014). Each individual's residential address in the year prior to each recorded weight measurement was linked to the average annual PM 2:5 in that location for that year (or, for Veterans who moved during a given year, a weighted average of annual PM 2:5 concentrations for each residential location, according to the time spent at each address). When two weight (i.e., BMI) measures were recorded within 90 d for the same Veteran, we retained only the earlier value. Mixed-effects models were used, with a random intercept, a random time slope at the individual-level, and a compound symmetry covariance structure. We report the association between PM 2:5 and the trajectory of BMI over the next year based on the regression coefficient of an interaction between time and PM 2:5 . As in models of incident obesity and weight gain, we used sequentially adjusted models for a) height for the weight outcome; b) state of residence; c) demographics of age, race, and sex; d) contextual characteristics; and e) smoking status. Height was included as a covariate because changes under the same conditions in weight may vary based on a person's height; baseline BMI and weight were not adjusted for because these measures were included in the set of outcome data.
We used slope estimates from the fully adjusted linear mixed models to plot average trajectories in BMI and weight according to percentiles of cumulative average PM 2:5 exposure throughout follow-up among subgroups of Veterans with cumulative exposures within ± 0:5% of the 10th, 25th, 50th, 75th, and 90th percentiles, respectively ( ∼ 39,000 Veterans per subgroup.) Cumulative average PM 2:5 exposure was defined as the average of total exposure to PM 2:5 during follow-up, calculated by taking the sum of exposures at each time point from follow-up to the time being analyzed and dividing by the difference in time from the start of follow-up to the time point in question. These trajectories represent the estimated average change in weight or BMI that Veterans in each subgroup would have experienced based on their cumulative average PM 2:5 exposure throughout the entire study follow-up period (from the start of follow-up to 8 y later). In addition to estimating linear associations between PM 2:5 concentrations and weight or BMI change within individuals, we estimated associations with PM 2:5 concentrations modeled as a restricted cubic spline, with knots at the 5th, 25th, 50th, 75th, and 95th percentiles. PM 2:5 values outside the 1st and 99th percentile (3.72 and 15:96 lg=m 3 ) were excluded to reduce the influence of outliers on these models.

Positive and Negative Controls
We estimated associations between PM 2:5 exposures and allcause mortality as a positive outcome control given consistent evidence of associations in other populations (Burnett et al. 2018;Cohen et al. 2017; GBD 2019 Risk Factors Collaborators 2020). All-cause mortality was ascertained using mortality data available from the VA Vital Status Mini file (Maynard 2017). We evaluated ambient sodium concentration as a negative exposure control for incident outcomes and weight gain over time among Veterans living ≤30 mi (48:3 km) from a U.S. EPA air monitoring station. Because there is no prior evidence or biologic plausibility for an association between atmospheric sodium and weight or BMI (Lipsitch et al. 2010), associations with ambient sodium concentrations would suggest potential bias in associations with PM 2:5 related to the use of air monitoring station data to estimate exposures in our study population. Ambient sodium exposures were estimated based on measured concentrations from the nearest U.S. EPA air monitoring station.

Sensitivity Analyses
We performed multiple sensitivity analyses. First, to enhance the spatial resolution of PM 2:5 exposure assessment, we estimated exposure based on measurements from the air monitoring station closest to each residence in analyses restricted to Veterans living within 30 mi (48:3 km) and 10 mi (16:1 km) of a U.S. EPA air monitoring station, respectively (Miller et al. 2007). We used the Haversine formula (Robusto 1957) to link the latitude and longitude of the ZIP code tabulation area centroid for each residential address to the latitude and longitude of the nearest PM 2:5 monitoring station. In addition, we repeated analyses restricted to Veterans whose residence latitude and longitude did not change by >1 km during follow-up. We also additionally adjusted for the number of hospitalizations and number of weight measurements in the year before baseline to assess potential confounding related to the frequency of interactions with the health care system, and we adjusted for marital status at baseline (Umberson 1992).
For models of the risk of obesity and a 10-lb (4:54 kg) increase in weight we also constructed within-city models [where city was defined as a metropolitan statistical area (MSA), which can include multiple 0:01 × 0:01 PM 2:5 grids] to further account for shared regional similarities (Miller et al. 2007). Models included both city-average exposure levels (defined for each year as the population-weighted average PM 2:5 exposure for all residents of a given MSA) as estimates of between-city effects, and the difference between each individual's exposure level from the city-level exposure as an estimate of within-city effects. Veterans were interval censored during time periods where they were not residing in an MSA. We report HRs for within-city effects corresponding to a 10-lg=m 3 higher individual-level PM 2:5 relative to the city-level average.
We also repeated Cox proportional hazards models after stratifying by baseline age (≤50, >50 to ≤60, >60 to ≤70, and >70 y), sex, race (Black, White, or other), and BMI (<25, >25 to <30, and ≥30 kg=m 2 ) to allow for differences in baseline hazards. In addition, we considered alternative outcome definitions including 1-and 3-kg=m 2 increases in BMI, where we adjusted for baseline BMI, and a 20-lb (9:07 kg) increase in weight relative to baseline, where we adjusted for baseline weight and height, 5% increases in BMI and weight relative to baseline, and the risk of becoming overweight or obese (BMI ≥25) among Veterans with a BMI of <25 at baseline. Finally, to assess potential bias due to differential loss to follow-up across PM 2:5 levels, we repeated linear mixed models restricted to Veterans who had at least one weight measurement ≥6:5 y after baseline, the median time from baseline to the last recorded weight measurement during follow-up of each Veteran.
Data were complete for all Veterans included in the analyses because Veterans with missing exposure, outcome, or covariate data were excluded from the study cohort. Effect estimates were considered statistically significant if the 95% confidence interval (CI) did not include the null value. All analyses were conducted using SAS Enterprise guide 7.1 (SAS Institute Inc.) and R (version 3.5.3; R Development Core Team).

Results
We assembled a cohort of 3,902,440 Veterans followed for a median 8.1 [interquartile range (IQR): 7.3-8.4 y], corresponding to 27,628,688 person-years. The demographic and health characteristics of the overall study cohort and by PM 2:5 quartile (at cohort entry) are described in Table 1. The cohort was 94% male, 76% White, and 14% Black, with a median age at baseline of 64 y (IQR: 56-75) ( Table 1). During the course of follow-up, 28% of Veterans died. At baseline, most (58%) were married, either current (29%) or former smokers (35%), and overweight (38%) or obese (40%). At baseline, cohort members lived in counties where 16% of residents were in rural areas, 16% of county adult residents reported excessive alcohol consumption, and 6% had limited access to healthy food. Veterans had a median of 10 weight measurements recorded during follow-up (IQR: 7-15). In general, Veteran characteristics were similar across PM 2:5 Table 1.  ADI is a measure of socioeconomic status disadvantage, with a range from low to high disadvantage of 0 to 100.
c NDVI is a measure of greenspace, and ranges from low to high from −1 (water), to 0 (bare ground), and, maximally, to 1 (dense vegetation).

PM 2:5 and the Risk of Weight Gain and Obesity
Among the subgroup of Veterans who were not obese at baseline (n = 2,325,768, 60% of the total cohort), 446,113 (19%) became obese (BMI ≥30) during follow-up ( Table 2). The risk of obesity was positively associated with a 10-lg=m 3 higher average annual PM 2:5 concentration, with an HR from the fully adjusted model of 1.08 (95% CI: 1.06, 1.11). In the entire cohort, 42% gained 10 lb (4:54 kg) during follow-up relative to their baseline weight. As for risk of obesity, the risk of a 10-lb (4:54 kg) increase in weight was also positively associated with 10-lg=m 3 higher average annual PM 2:5 concentration [fully adjusted HR = 1:07 (95% CI: 1.06, 1.08)].
To characterize the shape of the relationship between PM 2:5 concentrations and risk of obesity or a 10-lb (4:54 kg) weight gain, we used an ensemble modeling approach that accommodates a variety of nonlinear relations. The resulting estimates indicated positive associations for both outcomes ( Figure 1A,B). Parameters for the top three best-fitting models for both outcomes may be found in Table S3. Resulting curves had a single inflection point, with steeper slopes for lower PM 2:5 concentrations. For both the risk of Table 2. Association of PM 2:5 with risk of obesity and gain in body mass index or body weight in a national cohort of U.S. Veterans selected from 1 July 2010 through 31 June 2011 and followed until 31 December 2018 (n = 3,902,440 Baseline measurement was BMI. c Baseline measurement was height and body weight. , and state of residence, age, race, sex, smoking status, area deprivation index, normalized difference vegetation index, county-level percentage rural residency, population density, percentage limited access to healthy food, percentage access to exercise opportunities, and percentage of adults reporting excessive alcohol consumption. Lines represent the estimated difference in risk associated with a given PM 2:5 concentration compared with the reference concentration of 1 lg=m 3 (in consideration of the log-linear nature of the response). Bands represent the 95% confidence interval. 2:205 lb = 1 kg. Model parameters of the optimal model and the second and third best-fitting models used to derive the ensemble estimates are reported in Table S3. Note: BMI, body mass index; PM 2:5 , ambient fine particulate matter (particulate matter ≤2:5 lm in aerodynamic diameter). obesity ( Figure S4A) and the risk of weight gain ( Figure S4B), the curve based on the ensemble estimate was similar in shape to that for the optimal model. For both outcomes, the ensembled estimates displayed less precision than those from the optimal models.

PM 2:5 and Intra-Individual Changes in Weight
Linear mixed models of changes in BMI and weight indicated positive associations with PM 2:5 concentrations (Table 2), and spline analyses also suggested positive associations across the range of exposure levels experienced by the cohort ( Figure S5). We then estimated intra-individual changes in weight and BMI according to the distribution of cumulative average PM 2:5 exposure. On average, Veterans with cumulative average PM 2:5 levels at the 90th percentile gained an additional 3:73 lb (1:69 kg) or 0:54 kg=m 2 BMI by the end of follow-up compared with Veterans who had cumulative exposures at the 10th percentile ( Figure 2; Tables S4-S5).

Positive and Negative Controls
As expected, a 10-lg=m 3 higher average annual PM 2:5 was associated with all-cause mortality (

Sensitivity Analyses
Results were generally consistent with the primary analyses when we assigned PM 2:5 exposure levels based on mean annual measures from the nearest air monitoring station within 30 mi of each residence (n = 1,980,554 for incident obesity and n = 3,314,883 for all other outcomes) and 10 mi (16:1 km) from each residence (n = 1,320,183 for incident obesity and n = 2,190,712 otherwise); and when analyses were restricted to nonmovers (n = 1,628,448 for incident obesity, n = 2,705,653 otherwise) (Table 3). Results were also consistent with the primary analyses after additional adjustment for the number of  inpatient hospital stays and weight measurements in the year before baseline (as indications of the frequency of interactions with the health care system), and marital status at baseline (Table 3). Models of associations between incident obesity and a 10-lb (4:54 kg) weight gain also indicated positive associations when restricted to Veterans living within an MSA throughout followup, with HRs for a 10-lg=m 3 increase in within-city PM 2:5 of 1.20 (95% CI: 1.16, 1.24), n = 1,835,720 and 1.18 (95% CI: 1.15, 1.20), n = 3,066,141 for incident obesity and a 10-lb (4:54 kg) weight gain, respectively (Table S6). Estimates were similar to primary models when Cox proportional hazard models were stratified to allow for differences in baseline hazards by age, race, sex, and BMI [HR = 1:05 (95% CI: 1.03, 1.07) and HR = 1:07 (95% CI: 1.06, 1.08) for incident obesity and a 10-lb (4:54 kg) weight gain, respectively] ( Table S6). The HR for the incidence of overweight or obesity among 835,212 Veterans with a BMI of <25 at baseline was positive but closer to the null [HR = 1:04 (95% CI: 1.01, 1.06)] than the corresponding HR for incident obesity, whereas the association was stronger for a 20-lb (9:07 kg) increase in weight [HR = 1:18 (95% CI: 1.16, 1.20)] compared with the HR for a 10-lb (4:54 kg) increase (Table S7). Associations were also positive for a 10-lg=m 3 increase in PM 2:5 and the incidence of 1-and 3-kg=m 3 increases in BMI, a 5% increase in BMI, and a 5% increase in weight relative to baseline (Table S7). Finally, estimates for intra-individual changes in BMI and weight with a 10-lg=m 3 increase in PM 2:5 were similar to primary analyses when restricted to 1,955,448 Veterans who had at least one recorded weight measurement 6.5-7.5 y after baseline, with estimated mean increases in BMI and weight of 0:139 kg=m 2 (95% CI: 0.137, 0.142) and 0:962 lb (95% CI: 0.947, 0.977) (0:436 kg (95% CI: 0.430, 0.443)), respectively (Table S8).

Discussion
In this cohort of 3,902,440 U.S. Veterans followed for a median 8.1 (IQR: 7.3 to 8.4) y, who corresponded to more than 27 million person-years of follow-up, exposure to higher levels of PM 2:5 was associated with an increased risk of weight gain and obesity, and with intra-individual increases in weight and BMI during followup. Estimated exposure-response functions for incident obesity and a 10-lb (4:54 kg) weight gain suggested positive associations at all PM 2:5 concentrations above the minimum (2:4 lg=m 3 ), including concentrations below the 2012 National Ambient Air Quality Standard of 12 lg=m 3 (U.S. EPA 2013). Findings were consistent across multiple sensitivity analyses, including models using alternate exposure definitions, outcome definitions, and covariate adjustments. As expected, we also estimated a positive association with all-cause mortality (as a positive outcome control), whereas results were null for associations with ambient air sodium concentrations (as a negative exposure control). The constellation of evidence suggests that PM 2:5 may be obesogenic.
Our results build on the seminal discoveries that inhaled nanoparticles when sufficiently small, may permeate through the alveolar space and enter the bloodstream where they may interact with extrapulmonary organs (Miller et al. 2017). Animal studies also suggest that exposure to particulate matter activates genes associated with lipogenesis in adipose tissue leading to increased adipocyte size, increased adiposity, and increased visceral fat mass (Hamanaka and Mutlu 2018;Mendez et al. 2013;Sun et al. 2009). Animals exposed by inhalation to particulate matter also exhibit increased macrophage infiltration into adipose tissue, elevated concentrations of pro-inflammatory cytokines, impaired adipose mitochondrial function, increased leptin and adiponectin, and increased insulin resistance (Hamanaka and Mutlu 2018;Mendez et al. 2013;Sun et al. 2009). These mechanistic findings are complemented by evidence from human studies indicating that PM 2:5 is a metabolic risk factor (Pope et al. 2015) that has been associated with diabetes mellitus (Bowe et al. 2018b), chronic kidney disease, and other cardiometabolic conditions (Al-Aly and Bowe 2020; Bowe et al. 2020aBowe et al. , 2020b. Furthermore, epidemiologic literature has suggested that traffic density and traffic-related pollutant levels were associated with an increased risk of weight gain in children (Jerrett et al. 2010(Jerrett et al. , 2014. Taken together, these prior studies provide biologic and epidemiologic plausibility for the findings reported in this study. This study has several limitations. The cohort included U.S. Veterans who were mostly older, white, and male, which limits the generalizability of the study results. Models of incident obesity (or overweight) were necessarily restricted to Veterans who were not obese (or overweight) at baseline, and who may have therefore differed in their susceptibility to PM 2:5 and other factors that may further limit generalizability. Although we adjusted for several known confounders, we cannot completely rule out the possibility of residual confounding. Individual information on diet and exercise was not available, and we adjusted for smoking and marital status at baseline only. Although we used high-resolution exposure data and linked it to residential addresses, exposure misclassification may have been present. Many contextual covariates were assessed at the county level, and we did not investigate interactions between contextual characteristics. The use of electronic health records, which contain only information recorded at receipt of care, may have resulted in misclassification of outcome measurements and the timing of incident outcomes. Our analyses did not account for indoor exposure to air pollutants, did not consider lags in the relation between exposure and outcome, and did not consider potential heterogeneity in effect across different populations and by differences in composition and toxic content of PM 2:5 .
Strengths of the study include the large cohort size with wellcharacterized health characteristics and long follow-up time. PM 2:5 concentrations spanned the spectrum of levels experienced by people throughout the United States. We leveraged the opportunity that the outcomes of interest (weight and BMI) are continuous variables recurrently measured in the same individual-which enabled us to develop an analytic approach to estimate intraindividual change in weight associated with PM 2:5 . We used an ensemble modeling approach to characterize the shape of the association between PM 2:5 and obesity in our cohort, which allowed for any nonlinearity in the relationship. Finally, we developed and tested a negative exposure and a positive outcome control, to lessen concerns about spurious associations.
In summary, in a cohort of 3:9 million Veterans followed for more than 8 y (corresponding to more than 27 million personyears of follow-up), PM 2:5 air pollution was associated with a higher risk of obesity and weight gain. The association was evident for the incidence of obesity and a 10-lb (4:54 kg) gain in weight and for intra-individual changes in weight and BMI during follow-up. The shape of the estimated exposure-response function suggests that risk was evident at PM 2:5 concentrations <12 lg=m 3 -the current U.S. EPA air quality standard for average annual ambient PM 2:5 concentrations. Strategies aimed at addressing obesity may need to consider air pollution and other environmental exposures as potential causes.