Associations between Long-Term Air Pollution Exposure and Risk of Osteoporosis-Related Fracture in a Nationwide Cohort Study in South Korea

Bone health is a major concern for aging populations globally. Osteoporosis and bone mineral density are associated with air pollution, but less is known about the impacts of air pollution on osteoporotic fracture. We aimed to assess the associations between long-term air pollution exposure and risk of osteoporotic fracture in seven large Korean cities. We used Cox proportional hazard models to estimate hazard rations (HRs) of time-varying moving window of past exposures of particulate matter (PM10), sulfur dioxide (SO2), carbon monoxide (CO), nitrogen dioxide (NO2), and ozone (O3) for osteoporotic fracture in Korean adults (age ≥ 50 years) in the National Health Insurance Service-National Sample Cohort data, followed 2002 to 2015. HRs were calculated for an interquartile range (IQR) increase. Comorbidity and prescription associated with osteoporosis, age, sex, body mass index, health behaviors, and income were adjusted in the models. Effect modification by age, sex, exercise, and income was examined. We assessed 56,467 participants over 535,481 person-years of follow up. Linear and positive exposure-response associations were found for SO2, while PM10 and NO2 showed nonlinear associations. SO2 was associated with osteoporosis-related fracture with marginal significance (HR for an IQR [2 ppb] increase = 1.04, 95% CI: 1.00, 1.09). The SO2 HR estimates were robust in analyses applying various moving windows of exposure (from one to three years of past exposure) and two-pollutant models. The central HR estimate of O3 implied positive associations but was not significant (HR for 0.007 ppm increase = 1.01, 95% CI: 0.97, 1.06). PM10, CO, and NO2 did not show associations. Vulnerable groups by sex, age, exercise, and income varied across air pollutants and there was no evidence of effect modifications. Long-term exposure to SO2, but not PM10, CO, NO2 and O3, was associated with increased osteoporotic fracture risks in Korean adults.


Introduction
One of the major public health concerns of aging is bone-related health. The global prevalence of osteoporosis is more than 200 million people, which constitutes 1% of all disability and causes 8.9 million osteoporosis-related fractures every year globally [1]. Osteoporosis can be associated with bone fracture, which can lead to disability or mortality [2]. Studies projected a 2.28-fold increase in the number of osteoporosis-related fractures in the next few decades (1124,060 in 2018 to 2563,488 in 2050) in Asia and a 1.59-fold increase in health care costs (from 9.5 billion US$ in 2018 to 15 billion US$ in 2050) [3,4]. Given that ethnicity is a vulnerability factor for osteoporosis [21]. Calcium is a critical source for maintaining bone mass in adults, and studies have demonstrated that calcium intake from food was significantly lower in the Korean population compared to the American population based on the national nutrition health study [22]. Information on multiple time-varying measures of air pollution and individual-level confounders of osteoporosisrelated fracture in the nationwide large cohort offers a novel opportunity to examine the associations in long-term exposure to air pollution and incidence of osteoporosis-related fracture. This research is unique in its focus on various exposure windows, improved ascertainment of osteoporotic bone fracture, and examination of effect modification by individual-and community-level factors.

Study Population and Study Design
The Korea University Institutional Review Board approved this study. This cohort study used the National Health Insurance Service-National Sample Cohort (NHIS-NSC) data, in which information is available for age, district of address of that year, month of death if occurred during the follow-up period, income level, all prescription drugs, and treatment claim records for enrollees for the period between 2002 and 2015. Based on the Act on Social Security, all Koreans are eligible to become members of this national health insurance system. Enrollment of the national insurance system covers about 97% of the entire population in South Korea [18]. The NHIS-NSC dataset provides the cohort data on 2% of enrollees sampled from the total national insurance enrollees by age, sex, region, insurance type (e.g., employer-sponsored insurance beneficiaries, self-employed insurance beneficiaries, dependent of insurance beneficiaries), and income. This data are an open cohort where participants alive on January 1 of the year following their enrollment in the NHIS remain in the cohort and each participant has varying enrollment year. The NHIS is legally obligated to provide the national health checkup for insured adults. Fully-insured persons through their employment and self-insured head of household can take health checkup regardless of their age, while others in the household (e.g., dependents) can undergo the checkup from age 40 years old or older [23]. Employees working in non-office environments are required to take the health examination every year, and those in other occupations can take it every 2 years [23]. We limited our study regions to the target regions of the capital city (Seoul) and 6 metropolitan cities (Daejeon, Daegu, Incheon, Gwangju, Busan, Ulsan) in South Korea (n = 541,253) and targeted persons ≥50 years during the first 3 years of their enrollment in the NHIS-NCS dataset in these target regions (n = 84,544). The NHIS-NSC database provides the district of address of the participant for every year during the follow-up. The spatial basis for residential address of the enrollees in the database is district (administrative unit equivalent to the U.S. borough). Follow-up periods were calendar years.
To conduct this longitudinal study, we targeted those who had the health checkup data and lived in our target study regions in the baseline period (i.e., the first 3 years from entry) among adults ≥ 50 years from the NHIS-NSC data. We defined the baseline period as the first 3 years from entry and the follow-up start year as the 4th year from the entry to examine long-term past exposure to air pollution up to the 3 preceding years. We excluded persons who already had diagnosis of osteoporosis-related fracture at the baseline period (n = 1526), persons who were lost before the analysis start year (n = 16,641), and persons who did not have data of health checkups (n = 9910). As a result, about 66.8% (n = 56,467) of the persons age ≥ 50 years at the entry year in the target regions were included in the final analysis (n = 84,544) ( Figure 1). As we excluded participants lost for follow-up during the baseline period (i.e., death, move-out, loss of enrollment for the insurance system), higher percentages of older age groups (i.e., ≥65 y) were found among the excluded participants, but there were no significant differences for income or sex between the included participants and the excluded participants (Supplemental Table S1).  [18,20]. Other pathological fractures (ICD-9: 733.1 or ICD-10 M84.4, M84.5, M84.6, M84.7) that were not osteoporotic fractures were excluded. Participants who had diagnosis of osteoporosis-related fracture during the baseline period were excluded.

Air Pollution Assessment
We obtained the monitoring data of hourly PM10, SO2, CO, NO2, and O3 concentrations from the National Ambient air Monitoring Information System (NAMIS) for 2001-2015 [28]. PM2.5 data were not available because monitoring started in 2013 in our target cities. The system included 133 stations across 74 districts in the 7 target cities and provided verified observation data for 2001-2020. The average number of monitoring stations in each district was 2.4 (minimum = 1, maximum = 7, Q1 = 2, Q3 = 4). The average size of study districts was 73.1 km 2 (minimum 2.8, maximum 753.3, Q1 = 17.9, Q3 = 68.2). The year 2002 was selected as the start of exposure assessment in this study because that is the first year that the NHIS database provides residential location (district) of the cohort participants. We calculated daily means for each of the 5 air pollutants from hourly observations from each monitoring station. Then, annual mean of air pollution concentrations was calculated by averaging the daily mean values from monitoring stations within each district for each air pollutant. These district-specific annual mean values were used to assess We identified the initial (first) occurrence of bone fracture at osteoporosis-related sites (i.e., hip, spine, forearm, humerus [upper arm]) using the operational diagnosis [24][25][26][27] based on the primary diagnosis according to the 10th version of International Classification of Disease (ICD-10) codes that occurred during the period between 2002 and 2015 and based on procedure codes (e.g., open reduction and internal fixation, percutaneous/close pinning, external fixation, cast, skeletal traction, kyphoplastry, etc.  [18,20]. Other pathological fractures (ICD-9: 733.1 or ICD-10 M84.4, M84.5, M84.6, M84.7) that were not osteoporotic fractures were excluded. Participants who had diagnosis of osteoporosis-related fracture during the baseline period were excluded.

Air Pollution Assessment
We obtained the monitoring data of hourly PM 10 , SO 2 , CO, NO 2 , and O 3 concentrations from the National Ambient air Monitoring Information System (NAMIS) for 2001-2015 [28]. PM 2.5 data were not available because monitoring started in 2013 in our target cities. The system included 133 stations across 74 districts in the 7 target cities and provided verified observation data for 2001-2020. The average number of monitoring stations in each district was 2.4 (minimum = 1, maximum = 7, Q1 = 2, Q3 = 4). The average size of study districts was 73.1 km 2 (minimum 2.8, maximum 753.3, Q1 = 17.9, Q3 = 68.2).
The year 2002 was selected as the start of exposure assessment in this study because that is the first year that the NHIS database provides residential location (district) of the cohort participants. We calculated daily means for each of the 5 air pollutants from hourly observations from each monitoring station. Then, annual mean of air pollution concentrations was calculated by averaging the daily mean values from monitoring stations within each district for each air pollutant. These district-specific annual mean values were used to assess long-term exposure levels for each participant in that district as described in the statistical analysis section.

Covariates
We obtained information from the NHIS-NSC dataset for those who conducted health checkups. Variables included body mass index (BMI), smoking status, alcohol consumption (glasses per week), and frequency of physical exercise per week (0, 1-2, 3-4, 5-6, 7 times). We generated a binary variable for high alcohol intake using the alcohol consumption variable. A unit of alcohol varies from 8 to 10 g alcohol among different countries [19]. This is equivalent to a standard glass of beer (285 mL), a single measure of sprits (30 mL), a medium-sized glass of wine (120 mL), or 1 measure of an aperitif (60 mL) [19]. In this study, total daily alcohol intake was measured by the number of glasses of 'Soju', and 1 glass of Soju is considered to have 1 unit of ethanol. Given that the recommended amount suggested by the WHO is <5 units per day for men and <2.5 units per day for women, we defined high alcohol intake as 5 or more units for men and 3 or more units for women [19].
We considered the risk factors in the WHO Fracture Risk Assessment Tool (FRAX) algorithm. The FRAX algorithm provides 10-year absolute fracture risk of a person based on BMI, smoking status, excessive alcohol intake, use of oral glucocorticoids, diagnosis of rheumatoid arthritis, and diagnosis of secondary causes of osteoporosis [19]. We defined exposure to long-term oral glucocorticoids in each person and year if the participant was prescribed oral glucocorticoids (i.e., hydrocortisone, cortisone acetate, prednisone, prednisolone, methylprednisolone, triamcinolone) for more than 30 days in that the prior year [19]. Using all claim data based on ICD-10 codes for each participant during the follow-up period, we examined if a participant had confirmed diagnosis (both hospital visit and hospitalization) of rheumatoid arthritis during the follow-up period. Secondary causes of osteoporosis include disorders strongly associated with osteoporosis. These disorders include type 1 diabetes, hyperparathyroidism, hyperprolactinemia, hypopituitarism, Cushing's syndrome, chronic obstructive pulmonary disease, chronic renal failure, severe malnutrition, multiple myeloma, and idiopathic hypercalciuria [29,30]. Three secondary causes of osteoporosis including osteogenesis imperfecta, premature menopause, and hypogonadism were not considered in the analysis. The NHIS-NSC did not provide medical claim data on these causes for privacy and data protection issues. We generated a timevarying binary variable for osteoporosis-related conditions and assigned each participant a 'yes' for the year they were diagnosed with any of these symptoms and the following years during the follow-up period; otherwise, a 'no' was assigned. We also calculated the Charlson Comorbidity Index (CCI) for each participant and each follow-up year by reviewing the claim data in the NHIS cohort database to identify comorbid disease. Details of the included comorbid diseases and calculation methods were adopted from a previous study [31]. We also generated a time-varying binary variable for use of anti-osteoporosis agents (i.e., sodium alendronate, raloxifene hydrochloride, risedronate sodium, raloxifene hydrochloride, teriparatide, teriparatide acetate).

Statistical Analysis
Person-years of follow-up were considered for each participant in our statistical analysis. Follow-up time in the analysis was defined for each participant as the period from 1 January in the analysis start year (i.e., the first year after the baseline period) to the date of whichever occurred first: incidence of osteoporosis-related fracture, the end of follow-up (31 December 2015), death, moving out of the seven cities, or the loss to follow-up. Multivariate Cox proportional hazards models stratified by age (5-year interval) and sex were used to estimate hazard ratios (HRs) of bone fracture and 95% confidence intervals (CIs) in relation to long-term exposure to air pollution for an interquartile range (IQR) increase. For the participants included in the analysis, a right censoring occurred for those known to have moved out of our target cities, deceased, or lost enrollment for the NHIS between 2005 and 2015.
The district-specific average air pollution concentration values were assigned to each participant based on their residential district. To reflect changing levels of PM 2.5 over time, we considered time-varying variables for air pollution. Moreover, we considered various exposure windows for air pollution in our analyses because the appropriate exposure period for developing osteoporosis-related fracture is unknown. We used 1-, 2-, and 3-year moving annual average concentrations of PM 10 , SO 2 , NO 2 , O 3 , and CO in the prior year as a time-varying variable to estimate bone fracture risk in that year. We also applied a 2-year moving average of air pollution with a 1-year lag for each year (e.g., exposure in 2002-2003 for the fracture risk in 2005). We examined up to a 3-year moving window of past exposure considering the bone remodeling cycle. In normal bone structure, the median duration of the remodeling cycle in cortical bone is 120 days, and the total surface of cancellous bone is completely remodeled over a period of 2 years [32]. Hormonal disorders or some treatments can shorten the remodeling duration to 100 days or increase it to 1000 days [32].
For each air pollutant, we applied two-pollutant models additionally adjusting for air pollutants as the continuous variable to account for mutual correlation of pollutants.
In our main time-dependent models, we controlled for the potential confounding effects of exposure to oral glucocorticoids in the past year (yes, no), diagnosis of rheumatoid arthritis (yes, no), diagnosis of secondary causes of osteoporosis (yes, no), income-based insurance fee (<40th percentile, 40th-79th percentiles, ≥80th percentile; higher percentile indicating higher income), smoking status (current, former, and never), high alcohol consumption (yes, no), BMI, and frequency of exercise per week (0, 1-2, 3-4, 5-6, 7 times). All variables were included as time-varying variables. We verified the proportional hazard assumption of the time-fixed variables (i.e., sex) by plotting log-log Kaplan-Meier survival estimates [33]. We examined potential effect modification by age (50-69 y, 70+ y), sex (men, women), sex and age (men 50-69 y, men 70+ y, women 50-69 y, women 70+ y), frequency of exercise (low: never, moderate: 1-4 times per week, high: 5-7 times per week), and income (low: 0-39th, moderate: 40-79th, high: 80-100th percentile). Effect modification was examined by adding a multiplicative interaction term between air pollution and each effect modifier. Significance of the interaction term was examined by the p-value at a significance level of 0.1.
To examine the shape of the exposure-response relationship in this cohort, natural spline functions with 4 degrees of freedom (df) were applied for each air pollutant using Cox models, adjusted for all covariates. In estimating HRs, linear exposure-response relationships were applied in Cox models.

Results
Characteristics in the entry year or during the follow-up years of the included participants are shown in Table 1. Participants with and without outcomes differed in all characteristics except for diagnosis of secondary causes of osteoporosis during follow-up, which justified adjustment of the covariates. The average time of follow-up was 9.5 years and was higher in persons who did not experience osteoporosis-related fracture (6.4 years). Among the total 56,498 participants, 5398 participants (9.6%) experienced osteoporosisrelated fracture. Incidence of these fractures was significantly higher in women; 77.7% of the participants who experienced a fracture were women. About half of participants (48.6%) were age 50-64 years in their cohort entry year. Participants were exposed to PM 10 levels exceeding the national annual standard (50 µg/m 3 ). Long-term exposure concentrations of air pollution for the entire follow-up years for each participant were significantly different between the group experiencing the outcome and the group not experiencing the outcome for all air pollutants.
The levels of PM 10 , SO 2 , NO 2 , and CO were positively correlated, while the levels of O 3 showed negative correlations with PM 10 , CO, and NO 2 (Supplemental Table S2). The degree of negative correlations of PM 10 , CO, and NO 2 with O 3 were moderate (r < −0.5). Linear Cox models showed marginal significance for elevated risk for osteoporosisrelated fracture at higher levels of ambient SO 2 ( Table 2). The HR for an IQR increase in SO 2 concentration (2 ppb) was 1.04 (95% CI: 1.0, 1.08). The HRs and 95% CIs for osteoporosisrelated fracture associated with SO 2 were robust in all examined exposure windows. The HRs for an IQR increase were 1.00 (95% CI: 0.93, 1.07) for PM 10 (IQR = 12.7 mg/m 3 ) and 0.99 (95% CI: 0.94, 1.04) for CO (IQR = 0.192 ppm) in the models using 3-year moving annual average exposure levels. The HR for an IQR increase (0.012 ppm) was 1.01 (95% CI: 0.97, 1.05) for NO 2 in the model using a two-year moving annual average with a one-year lag. The HR for O 3 was 1.01 (95% CI: 0.97, 1.06, IQR = 0.007 ppm) in the models using oneand two-year moving annual average levels. Table 3 shows HRs of osteoporosis-related fracture associated with air pollution by age, sex, frequency of exercise, and income. The impacts of PM 10 , CO, NO 2 , and O 3 were not significant in any subgroups, while the HRs were significant in some subgroups for SO 2 . The effect modification by frequency of exercise per week was significant at a significance level of 0.1 (p-value = 0.064) for PM 10 , but no HRs were significant for the subgroups of exercise. The interaction terms of each effect modifier indicated no difference in the risks. As a result, there were no significant differences in the HRs among subgroups for all air pollutants. Figure 2 demonstrates that the three-year moving annual average SO 2 showed nearly linear and positive relationships with osteoporosis-related fracture. Decreasing log HR was found above the national standard of annual mean PM 10 (i.e., 50 µg/m 3 ). The associations for CO showed monotonic curves.
In models including simultaneous exposure to two air pollutants, the HR estimates of air pollution remained robust. The effect estimates of SO 2 slightly attenuated and remained statistically significant after adjustment for PM 10 and O 3 (Supplemental Table S3).  In models including simultaneous exposure to two air pollutants, the HR estimates of air pollution remained robust. The effect estimates of SO2 slightly attenuated and remained statistically significant after adjustment for PM10 and O3 (Supplemental Table S3).  The curves were estimated with Cox proportional hazards models using natural splines with 4 df for PM 10 , SO 2 , CO, and NO 2 , and 3 df for O 3 . Models adjusted for age, sex, history of rheumatoid arthritis or secondary causes of osteoporosis during the follow-up period, exposure to oral glucocorticoids in the prior year, use of anti-osteoporosis agents, the Charlson Comorbidity Index, household income-based insurance fee, BMI, smoking status, high alcohol intake, and frequency of exercise per week.

Discussion
Our study found significant and positive associations between long-term exposure to SO 2 and osteoporosis-related fracture, which is similar to the finding from a previous time-series analysis focusing on short-term air pollution exposure in Spain [17]. We did not find statistically significant associations with PM 10 exposure. This is similar to the findings from a cohort study in Norway examining PM and self-reported forearm fracture [16]. Similarly, a previous Korean study using the same cohort sample datasets, focusing on women aged >50 y, applying a time-independent PM exposure around the middle point of the follow-up period (2008-2009), and adjusting for age, income, CCI, BMI, smoking, drinking, exercise, and total cholesterol, did not find associations of PM 10 and osteoporotic fracture [18]. According to a recent systematic review [34], O 3 and hip fracture has only been studied in a time-series analysis conducted for Spain [17], and this study showed no associations between O 3 and hip fracture. While our study was based on retrospective cohort design and examined long-term air pollution exposure, we also found no associations between O 3 and osteoporotic fracture. The finding of no associations between NO 2 and osteoporotic fracture in our study was similar to results in the Norwegian cohort study showing the impact of long-term NO 2 exposure on self-reported forearm fracture [16], but it differed from the significant association between short-term NO 2 exposure and hip fracture in Spain [17]. While these varying study results warrant future research, our cohort study adds evidence on the associations between long-term air pollution exposure and osteoporosis-related fracture in adults.
Research gaps also remain regarding what is the appropriate exposure period for developing osteoporosis-related fracture. Currently, there is no consensus for windows of long-term exposure to air pollution in relation to osteoporosis-related fracture risk. The exposure windows used in the studies for bone fracture also varied. Alver et al. applied a time-independent variable for air pollution exposure (PM, NO 2 ) from 1992 to 2001 in their cohort study for self-reported forearm fracture [16]. Chiu et al. used average PM 2.5 concentration for the one year before the index date in their case-control study of osteoporotic fracture [1]. Prada et al. assessed the association of one-year PM 2.5 average and the count of hospital admissions from osteoporotic fracture for each year in 2003-2010 using generalized linear mixed models [14]. Although we examined past exposure up to three years, we were unable to assess exposure from the more distant past. It is unknown if the impacts from air pollutants may persist for a longer period, which warrants future research.
We found significant associations between SO 2 exposure and osteoporosis-related fracture in Korean adults. This finding can be reinforced by the results from a recent study in Taiwan reporting that SO 2 exposure was associated with increased osteoporosis risks in adults [8]. However, less is known regarding mechanisms for SO 2 compared to other air pollutants such as PM, O 3 , or NO 2 , which warrants more work on the physiological mechanisms for air pollution including SO 2 . Large quantities of SO 2 are emitted from fossil-fuel burning (e.g., coal, oil, diesel) and other industrial activities containing sulfur. Moreover, SO 2 and resulting pollutants can travel large distances far from the emission sources affecting a large number of residents [35]. SO 2 is generally highly correlated with other air pollutants, especially NO 2 and PM, which were positively correlated with osteoporosis in previous studies [36]. During the study period, high SO 2 emission and SO 2 concentrations were consistently observed in four of our seven target cities (e.g., Daejeon, Ulsan, Incheon, Busan) in our study, which contain many fossil-fuel-burning plants located at seaports and industrial complexes [37,38]. Although focusing on different health outcomes (e.g., mortality), a previous study using the Korean national health and nutritional examination surveys with mortality follow-up also found a similar pattern that associations between SO 2 and all-cause or cardiovascular mortality were larger than for the associations for NO 2 or PM 10 in seven metropolitan cities in South Korea [37]. Another Korean study based on mortality in these cities found that the association between short-term PM 10 exposure and all-cause or cardiovascular mortality was higher in cities with a higher ratio of SO 2 level to PM 10 level [39]. These results suggest that SO 2 may be a proxy for other air pollutants (e.g., finer particles, sulfate) that are not fully captured by PM 10 . We also identified an association between SO 2 exposure and osteoporotic fracture that persisted with adjustment for air pollution in the two-exposure models.
We did not find associations between O 3 and osteoporosis-related fracture in our analyses, which is aligned with a result by a time-series study in Spain by Mazzucchelli showing no associations between short-term O 3 exposure and osteoporotic fracture [17]. O 3 showed high correlation with PM 10 , CO, and NO 2 in our study cities indicating the difficulty in separating the influence of O 3 from that of the other air pollutants. The reduction of O 3 HR estimates in the two-pollutant models supports this as well. On the other hand, our study using the fixed-location measurements did not have information of spatial patterns of O 3 within a district. Therefore, high-resolution O 3 modeling data would be useful for assessing the association between long-term O 3 exposure and osteoporotic fracture with decreased exposure errors. Misclassification for NO 2 exposure could also be relevant for the null associations between NO 2 and osteoporotic fracture in this research. Moreover, NO 2 reacts and decays more quickly compared with PM 2.5 [40]. Uncertainties remain regarding how the exposure errors from these sources would affect the associations between NO 2 and osteoporotic fracture in our study regions. Like O 3 , high-resolution data of NO 2 would be useful in future studies.
Our study showed significant associations between long-term SO 2 exposure and osteoporosis-related fracture in women, whereas no significant associations were not found for men. Higher risks of osteoporosis in female adults than male adults have been suggested in previous studies [41]. Estrogen deficiency in postmenopausal females is a known risk factor for osteoporosis and osteoporotic bone fracture [1]. We also found significant associations between SO 2 and osteoporosis-related fracture in the younger age group (50-69 years), high income persons, and persons with higher frequency of exercise. However, a caution is needed in determining the effect modifications for these subgroups. The percentage of persons aged 50-69 years in the entry year was 83.9%, and the low percentage of persons ≥70 years would have affected the statistical power of the subgroup analysis. The percentage of persons who exercised ≥5 times per week was 71.6%, and the percentage of the high-income group (80-100th percentile of the income distribution) was 46.5%. Other possible reasons for the higher association for persons with higher frequency of exercise include potentially higher personal exposure to ambient air pollution during exercise outdoors and higher ventilation rates while exercising, although the data does not address the type of exercise or location of exercise.
Our study has some strengths. To our best knowledge, this is the first study to assess the impact of long-term SO 2 and CO exposure in relation to osteoporotic fracture. This is also the first longitudinal study to adjust for secondary causes highly associated with osteoporosis (i.e., factors considered in FRAX) in assessing the associations between air pollution exposure and osteoporosis-related fracture. Various moving windows of exposure to air pollution were examined. We identified osteoporosis-related bone fracture using operational ascertainment based on procedure codes as well as disease codes, (ICD) although we did not have information of osteoporosis based on BMD or laboratory examinations. Previous studies mainly identified osteoporotic fracture using ICD codes [14,17,18], but one study identified patients using ICD codes and surgical codes [1]. Moreover, the cohort in this study is a sample representing a large portion of the entire population. To the best of our knowledge this study is the first to examine effect modification by income and exercise levels for the air pollution impacts on osteoporosis-related fracture. Lastly, we contribute evidence on the air pollution impacts for Asian persons, and this race is known to be more vulnerable to risk of bone fracture [21].
We also have limitations. First, we did not consider temporal changes in chemical compositions of particulate matter, while different chemicals may be differently associated with BMD [7]. Second, use of station-based monitoring data of air pollution and lack of information on individual's time-activity patterns limit investigation of the spatial heterogeneity of air pollution, the differences between indoor and outdoor concentrations, and varying activity patterns that affect personal exposure. The spatial resolution of this research could obscure within-city heterogeneity of pollution levels as the Korean NHIS-NSC dataset only provides the address information at the district level. Third, we examined PM 10 rather than PM 2.5 although this smaller size of particulate matter could also be associated with these health outcomes. Air pollution monitoring data on PM 2.5 is not available for our full study period and study locations. Future work should investigate this pollutant including comparison of effect estimates from PM 2.5 and PM 10 along with consideration of different chemical structures of particulate matter. Fourth, we did not consider the timing of bone fracture within a year due to higher computational burden of person-day or person-month follow-up estimations. Incidence of bone fracture, particularly hip fracture, in the elderly is generally higher in cold seasons, and this seasonality has been suggested to be associated with deficiency in vitamin D or meteorological factors such as cold temperature, rain, ice, and snow [43][44][45], and although our study used annual averages, issues of seasonality warrant further study. Fifth, we also did not have information on BMD, so we did not consider ascertainment of osteoporotic fracture based on BMD. Dietary habits of study populations were not available as well. Sixth, we focused on adults living in urban regions to keep compatibility among participants for characteristics of air pollution exposure and unmeasured contextual confounders. Lastly, we limited our analyses to the cohort participants who had the health examination data. Limiting analysis to subjects with the health checkup data in the analysis has the advantage that various health-related individual-level variables (e.g., smoking, alcohol drinking, exercise) could be controlled, but statistical power would decrease due to lower sample size [23]. Furthermore, the national health examinations are available for fully-insured employees, self-insured heads of household, and dependents >40 y, so the included participants in our study may have different characteristics from the participants who did not have the health examination data in the NHIS-NSC dataset. According to the existing scientific literature, the factors increasing the likelihood of receiving health examinations include younger age, male sex, higher education and income, smoking, higher frequency of exercise, and lower depressive symptoms [46]. Given these patterns, our results may not be generalizable for the entire adult population >50 years in South Korea.

Conclusions
Long-term exposure to SO 2 was associated with a significant increase in osteoporosisrelated fracture incidence in adults (age ≥ 50 y) in a large cohort with detailed information on medical records and health behaviors. These results were robust to various modeling approaches for exposure window and adjustment for multi-pollutant exposure. The associations varied by frequency of exercise and income, but there were no significant effect modifications by these factors across the studied air pollutants. Our results suggest that long-term exposure to air pollution may be an emerging risk factor of osteoporosis-related