Data of infertile women who visited a single fertility center for fertility evaluation was used. The center is the largest single fertility center of the country located in Seoul, South Korea, and a half of patient population are from outside of the capital area (9 provinces of the country). The result of ovarian reserve test and residential address of those who newly visited between January 2016 and September 2018 were obtained. Because our exposure assessment is based on air quality monitoring data, our analysis was restricted to the subjects living within 6 km from a monitoring station to minimize measurement error [18,19,20]. Having further excluded women previously diagnosed with chromosomal abnormality, having a history of unilateral or bilateral oophorectomy, and aged < 20 or > 49 years, the final study population included 2276 women.
Estimation of ambient air pollution exposure
Hourly concentrations of particulate matter less than or equal to 10 or 2.5 μm in diameter (PM10 and PM2.5, respectively), nitrogen dioxide (NO2), carbon monoxide (CO), sulfur dioxide (SO2), and ozone (O3) measured at the 269 air quality monitoring sites located throughout the country for 2016 and 2018 were used. PM10 and PM2.5 are inhalable particles containing chemical compounds which can reach deep into the lungs and even into the bloodstream . These data were from the National Institute of Environmental Research (NIER, https://www.nier.go.kr/). The following daily representative concentrations of the six pollutants were determined for each woman: 24-h averages for PM10, PM2.5, NO2, and SO2, and maximum of seventeen 8-h moving averages for CO and O3. In order to obtain consistent measurements, the daily averages calculated only for days in which > 75% of hourly measurements (18 h) were used at each site. The maximum concentrations of CO and O3 were used because the majority of their production is affected by commuter traffic and sunlight, respectively . Using these daily representative concentrations, average concentrations for the four periods of 1, 3, 6, and 12 months before the ovarian reserve test were computed. These four periods of exposure allow us to explore the critical period of exposure, given the immediate and long-term change in ovarian reserve was observed when ovaries are affected by chemical injury [23, 24]. The exposure estimate of each period at the nearest monitoring sites was assigned to the women based on their geocoded home addresses at the time of the test assuming women’s addresses remained the same within a year.
Outcome variable: measurement of ovarian reserve
Anti-Müllerian hormone (AMH) is a widely used indicator of ovarian reserve in women of reproductive age . Those with serum AMH as low as 0.5–1.1 ng/mL are likely to respond poorly to ovarian stimulation and thereby show low pregnancy rates in assisted reproductive technology treatments . In addition, DOR manifested as low AMH has predictive values for the risk of cardiovascular disease [27, 28]. The ovarian reserve of each woman was determined by measuring levels of AMH. Serum obtained on menstrual day 2 or 3 was separated from peripheral blood by centrifugation and stored at − 80 °C until analysis. AMH was measured using the Elecsys® AMH assay (Roche Diagnostics), which is a sandwich assay based on electrochemiluminescence technology. The total duration of the assay is 18 min; the sample volume is 50 μL. The assay is calibrated against the Beckman Coulter AMH Gen II ELISA assay with a measuring range of 0.01–23 ng/mL. Considering the age-dependent change in AMH, the AMH ratio, defined as observed AMH divided by age-specific AMH level, was used as an age-adjusted measure of ovarian reserve . The age-specific reference level of AMH was calculated using a previously described quadratic model (logAMH = −1.442 + 0.225 × age − 0.004 × age2) [30, 31]. Given the lower limit of AMH in the minimal criteria for DOR is 0.5 ng/mL , further examination was conducted for the risk of AMH < 0.5 ng/mL (“Low AMH”) [33, 34].
Information of women’s age, bodyweight, height, previous smoking history, working status, and residential address was retrieved from medical records. Bodyweight and height were measured at the time of initial visit. Body mass index (BMI) was calculated by dividing person’s weight in kilograms with their height in meters squared and categorized it into one of three groups (low, normal, overweight, and obese) based on the recommended BMI cut-off points for determining overweight and obesity in Asian populations . History of smoking and working status were recorded as binary variables (yes or no).
Descriptive statistics were calculated for the total population over the country. Given women living in capital area may be less deprived socioeconomically and exposed to higher air pollution compared to those living outside, a secondary analysis was conducted restricting to those living in Seoul. The characteristics between Seoul residents and the others were compared. Pairwise correlation structure between air pollutants for different averaging periods was examined with Spearman correlation test. We conducted generalized linear regression analyses for the AMH ratio and logistic regression analyses for low AMH controlling for woman’s age. Effect estimates are presented as regression coefficients for AMH ratio and odds ratios (ORs) for presence of low AMH with their 95% confidence intervals (95% CIs) per an interquartile range (IQR) increment in each pollutant concentration. The IQRs were 8.0 μg/m3 for PM10, 2.4 μg/m3 for PM2.5, 8.8 ppb for NO2, 1.2 ppb for SO2, 870.0 ppb for CO, and 11.0 ppm for O3. Risk estimates were adjusted for age (excluded in models for AMH ratio which already adjusted for age), BMI, season at the time of testing, previous smoking history, and district of residence in all models. Age in years (< 35, 35–40, or ≥40 years), BMI, working status, history of smoking (yes or no), and season (March to May, June to August, September to November, or December to February) were recorded as categorical variables. The clustering effect of 177 administrative districts was also included in the model using the cluster function of the R package “survival.” The AIC and residual deviance of the generalized linear model (which includes all covariates) for AMH ratio were 4498.6 and 865.23. The AIC and residual deviance were 4448.8 and 856.09 when 1-month average PM10 is added to the model. We fitted generalized additive models with non-parametric smoothing splines to further assess a non-linear exposure–response relationship between air pollutant concentration and AMH ratio. The model was fitted with gam function of package “gam.” We conducted all the analyses in R (R Version 3.2.1).