Coke Oven Emissions (COEs) is a broad term used to refer to mixed gases and solid particles produced during coke production, and its main component is polycyclic aromatic hydrocarbons (PAHs). Numerous studies have demonstrated that exposure to COEs/PAHs raises the risk of several diseases, including cardiovascular diseases, respiratory diseases, and lung cancer ^[1]. Workers at coke oven plants are more likely to be exposed to COEs than others, which increases their likelihood of getting associated diseases.

Numerous studies have demonstrated that COEs exposure affects the mitochondrial DNA copy number (mtDNAcn) and inhibits mitochondrial biogenesis ^[2]. Mitochondria are dual-membrane organelles found in all eukaryotic cells that regulate various cellular activities, such as adenosine triphosphate generation, cell homeostasis, apoptosis, and macromolecule biosynthesis ^[3]. Mitochondrial DNA (mtDNA) is a circular, double-stranded DNA molecule containing 16,569 base pairs without introns, and its number in the genome is called mtDNAcn. Usually, tens to thousands of mtDNAcn are found in a cell ^[4]. For mitochondrial function and cellular homeostasis, the integrity of the mitochondrial genome is of great importance. Due to the underdevelopment of the damage repair system of the mtDNA, it is susceptible to damage caused by reactive oxygen species (ROS) and thus becomes the target of many environmental pollutants ^[5]. As research progresses, the mtDNAcn is expected to become a biomarker of mitochondrial dysfunction ^[6]. The fact that COEs exposure causes a significant amount of ROS to be produced in vivo, which may inhibit mtDNA synthesis and result in mitochondrial dysfunction, suggests that the mtDNAcn may be a sensitive biomarker for diseases associated with COEs exposure ^{[7, 8]}.

It is important to establish appropriate limits for COEs concentrations in the workplace environment to protect coke oven workers from various health hazards caused by exposure to COEs-induced mitochondrial damage. There are certain limitations to extrapolating exposure limits to human health risks determined only by the size of the no-observed-adverse-effect level (NOAEL) or lowest-observed-adverse-effect level (LOAEL) metrics from animal studies ^[9]. To compensate for the limitations of the NOAEL/LOAEL, Crump et al. proposed the benchmark dose (BMD), which has been widely used in the risk assessment of environmental contaminants ^[10]. Our research group previously found a dose-response relationship between COEs exposure and telomere damage and estimated its exposure limit using the Benchmark Dose Software (BMDS) program published by the United States Environmental Protection Agency (USEPA) ^[11]. Currently, several programs have been developed to evaluate the BMD. In the present study, we selected three common methods to estimate the BMD and its lower 95% confidence limit (BMDL) of COEs-induced mitochondrial damage, including the BMDS version 3.20, Bayesian benchmark dose modeling (BBMD), and possible risk obtained from animal studies (PROAST) version 67.0.

Based on the findings of the previous study, 782 people were chosen as research subjects to investigate the dose–response relationship between COEs exposure and mtDNAcn and to assess the BMDL of COEs-induced mitochondrial damages using the BBMD, BMDS, and PROAST.

A total of 782 volunteers, including 238 participants in the control group and 544 workers in the exposure group, were enrolled in the study. The exposure group in this study comprised workers aged 18 to 60 who had at least one year of work experience. The control group members did not have any occupational COEs exposure, although living in the same city as the exposed group. After training, investigators gathered basic information from all the volunteers in a questionnaire survey. “Smoking” was defined as smoking more than one cigarette a day for over six months, while “alcohol consumption” was defined as having consumed more than two drinks per week in the past six months since the survey began. Na₂EDTA anticoagulation tubes were used to collect peripheral blood from every individual. This research protocol was approved by the Ethics Committee of Zhengzhou University (ZZUIRB 2021−153). Informed consent forms were signed by the subjects before implementation, and international and national ethical standards for biomedical research were strictly followed. The previous article published by our research team has more detailed basic information ^[11].

In the present study, the exposure levels of COEs in the coking plant’s various workshops were examined. The sampling point was set up in accordance with the Code for Monitoring and Sampling of Hazardous Substances in the Air of the Workplace (GBZ159−2004). There were four workshops, each with two coke ovens. In each workshop, a total of 18 representative locations, such as furnace top, furnace side, furnace bottom, duty room, and office, were selected. The total suspended particulate medium-flow sampler (Laoying 2050, Laoshan, Shandong, China) was used to gather air samples at the location. The relative height of the sampler inlet was approximately 1.5 m from the ground. Then, the sampler was started, and the flow rate was adjusted to 100 L/min. After 5 min of sampling, the sampling flow rate, atmospheric pressure, temperature, and meteorological parameters were recorded during sampling and 5 min before the end of sampling. The flow rate was often observed during the sampling period when the sampling flow rate decreased by more than 5% due to the high pollutant concentration at 4 h. The sampling time and flow rate of each segment were recorded.

The concentration of COEs was measured in the laboratory according to the standard “fixed source emissions–determination of benzene soluble particulate matter–Soxhlet extraction method” (HJ690−2014). The cumulative exposure dose (CED) of COEs (CED-COEs) for each participant in the study, comprising the control and exposure groups, was determined using the following equation ^[12]: CED [(mg/m³) · year] = Σ C (mg/m³) × T (year), where C stands for the average daily workplace exposure concentration of COEs and T represents the hours worked at a particular site. In the exposed group, the time-weighted average (TWA) concentration was calculated according to the type of work (Supplementary Table S1, available in www.besjournal.com). In the control group, C refers to the living environment concentration, and T is the age. The coefficient of variation was 0.46%–13%, and the limit of determination was 0.004 mg/m³. Additional details of COEs collection and concentration measurements were described in previous publications of our research group ^[13].

The mtDNAcn of peripheral leukocytes was detected via real-time fluorescence-based quantitative PCR (RT-PCR). Mitochondrial NADH dehydrogenase subunit 1 gene (MT-ND1) is an internal reference gene of the mtDNA, while human β-globin is an internal reference gene of nuclear DNA. Two parallel samples were set for each sample, and each reaction system was 10 μL in total, consisting of the DNA template, Mix, DNase-free and RNase-free water, and primers. This study utilized the same measurement method as early studies ^[14].

Statistical analysis was performed using SPSS 21.0 (SPSS Inc., Chicago, USA). The number and percentage of qualitative data were described. Median and quartile were used to describe the quantitative data that did not obey a normal distribution. The quantitative data following a normal distribution were described by (Mean ± SD) deviation. The two-independent sample t-test was used to analyze the distribution of quantitative variables between the exposed group and control group, and the χ² test was used to analyze the distribution of qualitative variables between the exposed group and control group. The mtDNAcn was compared between groups by the two-independent sample t-test. The Jonckheere-Terpstra test was used for comparisons between multiple groups, such as BMI. The Jonckheere-Terpstra test is a nonparametric test used to test whether there are significant differences in the distribution of multiple populations from multiple independent samples. The classic t-test is a comparison between two groups, and the data obey the normal distribution, which needs to be based on a specific population distribution. The relationship between the demographic characteristics and mtDNAcn was analyzed by the two-independent sample t-test and variance analysis. A generalized linear model was used to analyze the relationship between the CED-COEs and mtDNAcn. Two-sided tests were used for all statistical analyses with test level α = 0.05.

The BMD estimate of the reduction in the mtDNAcn caused by CED-COEs was analyzed using the BMDS version 3.20 (USEPA), web-based BBMD ^[15], and PROAST version 67.0 ^[16]. The software fitted the dose-response relationship of the total population, male population, and female population. The benchmark response in this study was 10%, estimating the COEs exposure to the BMD and BMDL. According to the BMDS guidelines ^[17], a goodness-of-fit P-value greater than 0.1 and a minimum fit equation of the Akaike information criterion (AIC) were selected as the optimal model. According to BBMD technical guidelines, the Model-average was selected as the model with the best goodness of fit ^[15]. On the PROAST website, the model with the lowest AIC will be selected as the best model ^[18].

In previous papers published by the research group, the demographic characteristics of the research subjects have been fully described ^[13]. The results show that there were statistically significant differences in age, gender, smoking, and drinking (P < 0.05), whereas there were no statistically significant differences in the BMI between the two groups (P > 0.05). The comparison between the exposed group and control group showed that the CED-COEs in the exposed group were 1.12 (0.34, 2.14) higher than that in the control group 0.07 (0.06, 0.09), with a statistically significant difference (P < 0.001). The mtDNAcn in the exposed group was lower than that in the control group (0.60 ± 0.29 vs. 1.03 ± 0.31, P < 0.001) (Table 1).

The intergroup analysis results (Table 2) show that there were statistically significant differences in the mtDNAcn between the control group and exposed group in terms of age, gender, smoking, drinking, and BMI, and the mtDNAcn of the exposed group was lower than that of the control group (P < 0.001). The results of stratified analysis within the group showed that in the stratified gender of the control group, the mtDNAcn of the females was higher than that of the males (1.08 ± 0.34 vs. 1.00 ± 0.28; P = 0.039), and there was no statistical significance between the mtDNAcn and other demographic characteristics of the control group or exposed group (P > 0.05).

The CED-COEs by the quartile were divided into four groups as the independent variable, mtDNAcn as the dependent variable, and age, gender, smoking, drinking, and BMI as the covariates. The low-dose group was used as the control group, and the relationship between CED-COEs and the mtDNAcn was analyzed using the generalized linear model. The trend test showed that the mtDNAcn decreased as the CED-COEs exposure increased (P < 0.001) after adjusting for all covariates (Table 3).

In this study, according to the EPA (USEPA) technical guidelines for the BMD and taking into account the mtDNAcn damage, the cut-off point of the mtDNAcn was defined as the 5th percentile of the control group. That is, ≤ 0.57 was defined as the impaired group, and > 0.57 was defined as the normal group. Table 4 shows that the mtDNAcn damage rate increased with increasing CED-COEs (P < 0.001), further stratified by gender, and a trend for males and females can also be seen.

Based on the BMDS user guide, a goodness-of-fit P-value greater than 0.1 and a minimum fit equation of AIC were selected as the optimal model. The model types are Dichotomous-Hill, Gamma-Model, Log-Logistic, Multistage, Weibull, Logistic, Log-Probit, Probit, and Quantal-Linear. The results of the BMDS show that the total population and female population models were fitted unsuccessfully, and the male population model was fitted successfully as the Dichotomous-Hill model, as shown in Table 5. The Dichotomous-Hill dose-response model formula is P [mtDNAcn damaged] = g + (v − v × g)/{1+exp[-a − b × Log(dose)]}. The BMD and BMDL of males were 0.087 and 0.076 mg/m³ per year, respectively.

According to the BBMD technical guidelines, the Model-average was selected as the model with the best goodness of fit. As shown in Table 6, the model of the total population, male population, and female population all chose the Model-average. The BMD and BMDL in the total population were 0.100 and 0.068 mg/m³ per year, respectively. The BMD and BMDL in the males were 0.411 and 0.062 mg/m³ per year, respectively. In the females, the BMD and BMDL were 0.439 and 0.061 mg/m³ per year, respectively (Table 6, Supplementary Figures S1–S3, available in www.besjournal.com).

In the PROAST method, the optimal model is selected when the AIC value is the lowest. The model types are Gamma, Log-Logistic, Weibull, Log-Probit, Two-State, LVM-Expon.m5-, and LVM-Hill m5-. BMD of the total population, males, and females were 0.0787, 0.0791, and 0.0751 mg/m³ per year, respectively. The BMDL in the total population, males, and females were 0.0735, 0.0713, and 0.0671 mg/m³ per year, respectively (Table 7).

The total population was subdivided according to the total length of service (40 years for males and 35 years for females). The BMDL of the male population was 0.076 mg/m³ per year. The occupational exposure limits (OELs) of the male COEs was 0.00190 mg/m³. Based on the application guidelines of the BBMD, our results determined that the BMDL of the mtDNAcn damage was 0.068 mg/m³ per year in the total population, 0.063 mg/m³ per year in the males, and 0.061 mg/m³ per year in the females. The OELs is 0.00170 mg/m³ for the total population, 0.00158 mg/m³ for the males, and 0.00174 mg/m³ for the females. Based on the same calculation, in PROAST, the OELs of the total population, males, and females are 0.00184, 0.00178, and 0.00192 mg/m³, respectively (Table 8).

The sources of COEs in our environment include coke production, asphalt mining, incomplete combustion of organic materials, and automotive exhaust. When working in a coke oven facility, various particles and volatile organic substances, especially PAHs, can enter the body through the skin or respiratory system and cause some harm to the mitochondria in the organism ^{[19, 20]}. Mounting evidence shows that long-term exposure to environmental contaminants may result in mitochondrial malfunction, which may eventually lead to human diseases ^{[21, 22]}. An abnormal mitochondrial function is associated with various complex human diseases, such as respiratory diseases ^[23], Alzheimer’s disease ^[24], cardiovascular diseases ^[25], and even cancer ^[26]. The mtDNAcn, a biomarker representing changes in the tissue metabolic state brought on by exposure to external detrimental stimuli, will indicate changes in the mitochondrial content and mitochondrial integrity ^[27].

Based on the study’s findings, the mtDNAcn was lower in the COEs exposure group (0.60 ± 0.29 vs. 1.03 ± 0.31), and this difference is statistically significant (P < 0.001). A dose-response connection between the CED-COEs and mtDNAcn was observed. The fact that the mtDNAcn decreased as the CED-COEs increased (P < 0.001) indicates that the mtDNAcn is a biomarker of COEs exposure. Pieters N’s study ^[28] found that when indoor non-volatile COEs levels doubled, the mitochondrial damage was negatively correlated with COEs exposure levels. Du et al. ^[29] examined the urine of 697 coke oven workers who had worked at the Shanxi coking plant for more than a year and found that PAHs urinary metabolites were negatively correlated with the peripheral blood mtDNAcn. The above results are consistent with the results of this study. Some studies have found that short-term exposure to environmental toxins may increase the mtDNAcn and increase the number of mitochondria to compensate for the damage, and long-term continuous exposure to pollutants may lead to a decrease in the mtDNAcn ^[30]. According to the analysis, CED-COEs in the total population, males, and females showed a dose-response relationship with the mtDNAcn damage rate. The mitochondrial damage rate increased with the increase in CED-COEs. Based on such interesting results, we further investigated to fully understand their impact on health risk assessment. This is the exposure limit that can be used as a new COEs exposure.

With the development of our society and economy, the importance we attach to occupational safety increases every year. Different countries have established OELs. China’s regulations published in 2019 limit the time-weighted allowable concentration of COEs to 0.1 mg/m³ (benzene soluble part) ^[31]. The TWA of COEs recommended by the National Institute for Occupational Safety & Health of the United States is 0.20 mg/m³ (benzene soluble part) for coke oven workers ^[32]. In addition, the permissible exposure limit set by the Occupational Safety and Health Administration in 1998 was 0.15 mg/m³ (benzene soluble portion) ^[33].

At present, the traditional level of unobserved adverse events has been replaced by the BMD to assess human health risks. In this study, we used three different calculation methods. The BBMD and PROAST use frequency-based statistical methods, but the BMDS is probabilistic ^[15]. The BMDS recommended by the USEPA can fit various types of data, among which the two most commonly used are dichotomy data and continuous data. The BMDS mainly conducts dose-response model fitting and parameter estimation based on the maximum likelihood method ^[34]. The BBMD is based on Bayesian analysis, and MCMC sampling is used for model parameter estimation and reference dose estimation. The BBMD can not only estimate the BMD of binary and continuous data but can also provide a method to incorporate prior information through the prior distribution of model parameters, which can improve the reliability of the dose-response model of inferior data. In addition, the BBMD can calculate the weight of the fitting model and further calculate the average reference dose of the model ^[35]. The web application for PROAST was created in the R programming language to make it simple to use. Its ability to compare dosage responses across several subgroups is an essential feature.

In the current investigation, the mtDNAcn was employed as a biomarker of exposure to occupational COEs for BMD analyses. According to the BMDS user guide, the fitting of the total population and female population model was not successful, but the fitting of the male population was successful. The above results show that the BMDS failed to fit the data of the total population and females. However, the BBMD successfully fitted the models of the total population, males, and females, which also confirmed that the BBMD procedure had less failure rate ^[15]. Through a variety of calculation methods, the exposure limit of PAHs was conservatively estimated to be 0.002 mg/m³. Although the difference between males and females was small, it still could not be ignored. Compared with current OELs, the BMDL estimated by the mtDNAcn as a biomarker of health damage caused by COEs was lower in this study. These findings suggest that the current recommended OELs may not be sufficient to protect coke oven workers and that the appropriate reduction of OELs may be considered to better protect the health of coke oven workers. Policymakers can appropriately increase the current emphasis on precision health risk assessment.

COEs exposure was measured using the BBMD at the telomere injury threshold in earlier research, resulting in COEs exposure limit of 0.024 mg/m³ for males and 0.018 mg/m³ for females (benzene soluble fraction) ^[11]. COEs exposure limits were calculated using oxidative damage values, and the reference COEs exposure limits were 0.037 mg/m³ for males and 0.021 mg/m³ for females ^[36]. These results are significantly higher than our findings, indicating that the mtDNAcn may be a more sensitive biomarker in COEs-exposed employees than telomere length and oxidative damage. These markers may not always indicate disease, but they can indicate a metabolic malfunction or accelerated aging, both of which are substantial alterations. Over the past 10 years, numerous studies have found links between telomere length and a variety of age-related illnesses, including osteoporosis, dementia, cancer, and cardiovascular disease ^[37].

However, the mtDNAcn may not be a disease endpoint and may be affected by certain conditions. Therefore, more population studies or in vitro experiments are needed to prove the results of this study.

In conclusion, this study shows that there is a dose-response relationship between COEs and mtDNAcn injury rate of peripheral leukocytes of coke oven workers; that is, the mtDNAcn decreased with the increase in CED-COEs. The results of the BDML were estimated by three methods, suggesting that lower exposure levels can better protect the susceptible population, which can provide a reference for the formulation of COEs exposure threshold and provide a good working environment for workers.

This research protocol was approved by the Ethics Committee of Zhengzhou University (ZZUIRB 2021−153). International and national ethical requirements for biomedical research were strictly observed during implementation, and informed consent forms were signed by the subjects.

Figure S1.  Dose-response relationship between CED-COEs and mitochondrial DNA copy number damage rate in the total population. (A) Logistic model (B) Log-Logistic model (C) Probit model (D) Log-Probit model (E) Quantal-linear model (F) Multistage model (G) Weibull model (H) Dichotomous-Hill model. The X-axis described the dose of CED-COEs and the Y-axis described the rates of mtDNAcn damage. This plot includes original dose-response data and a fitted curve with its 90th percentile interval shaded in blue.

Figure S2.  Dose-response relationship between CED-COEs and mitochondrial DNA copy number damage rate in male. (A) Logistic model (B) Log-Logistic model (C) Probit model (D) Log-Probit model (E) Quantal-linear model (F) Multistage model (G) Weibull model (H) Dichotomous-Hill model. The X-axis described the dose of CED-COEs and the Y-axis described the rates of mtDNAcn damage. This plot includes original dose-response data and a fitted curve with its 90th percentile interval shaded in blue.

Figure S3.  Dose-response relationship between CED-COEs and mitochondrial DNA copy number damage rate in female. (A) Logistic model (B) Log-Logistic model (C) Probit model (D) Log-Probit model (E) Quantal-linear model (F) Multistage model (G) Weibull model (H) Dichotomous-Hill model. The X-axis described the dose of CED-COEs and the Y-axis described the rates of mtDNAcn damage. This plot includes original dose-response data and a fitted curve with its 90th percentile interval shaded in blue.