A Quantitative Meta-Analysis of the Relation between Occupational Benzene Exposure and Biomarkers of Cytogenetic Damage.

BACKGROUND: The genotoxicity of benzene has been investigated in dozens of biomonitoring studies, mainly by studying (classical) chromosomal aberrations (CAs) or micronuclei (MN) as markers of DNA damage. Both have been shown to be predictive of future cancer risk in cohort studies and could, therefore, potentially be used for risk assessment of genotoxicity-mediated cancers.
OBJECTIVES: We sought to estimate an exposure-response curve (ERC) and quantify between-study heterogeneity using all available quantitative evidence on the cytogenetic effects of benzene exposure on CAs and MN respectively.
METHODS: We carried out a systematic literature review and summarized all available data of sufficient quality using meta-analyses. We assessed the heterogeneity in slope estimates between studies and conducted additional sensitivity analyses to assess how various study characteristics impacted the estimated ERC.
RESULTS: Sixteen CA (1,356 individuals) and 13 MN studies (2,097 individuals) were found to be eligible for inclusion in a meta-analysis. Studies where benzene was the primary genotoxic exposure and that had adequate assessment of both exposure and outcomes were used for the primary analysis. Estimated slope estimates were an increase of 0.27% CA [(95% CI: 0.08%, 0.47%); based on the results from 4 studies] and 0.27% MN [(95% CI: -0.23%, 0.76%); based on the results from 7 studies] per parts-per-million benzene exposure. We observed considerable between-study heterogeneity for both end points (I2>90%). DISCUSSION: Our study provides a systematic, transparent, and quantitative summary of the literature describing the strong association between benzene exposure and accepted markers of genotoxicity in humans. The derived consensus slope can be used as a best estimate of the quantitative relationship between real-life benzene exposure and genetic damage in future risk assessment. We also quantitate the large between-study heterogeneity that exists in this literature, a factor which is crucial for the interpretation of single-study or consensus slopes. https://doi.org/10.1289/EHP6404.

Introduction

Benzene is a well-known environmental contaminant that was classified as a human carcinogen (IARC Group 1) in 1974 (IARC 1974), which was recently reconfirmed (IARC 2018). Benzene and its metabolites may cause leukemia via genotoxic effects on the pluripotent hematopoietic stem cells, resulting in chromosomal changes (IARC 2009). A substantial number of cross-sectional studies have demonstrated that benzene exposure is associated with the occurrence of (classical) chromosomal aberrations (CA) and micronuclei (MN) that are considered early markers of genotoxicity [see reviews by Zhang et al. (2002) and McHale et al. (2012)]. A CA is often defined as the appearance of missing, extra, or irregular portions of chromosomal DNA (NHGRI 2016), and an MN is considered to be the small nucleus that forms whenever a chromosome fragment or whole chromosome is not incorporated into the daughter nuclei after cell division (Fenech 2002). Information on intermediate end points is increasingly used in the risk assessment for many chemicals; for example, the International Agency for Research on Cancer (IARC) integrated mechanistic evidence with evidence from other data streams to support conclusions regarding carcinogenicity (Smith et al. 2016), and recently the Dutch Health Council Committee recommended a health-based occupational exposure limit for benzene that was (primarily) based on hematological effects in humans (Health Council of the Netherlands 2014).

To date, however, there have not been many attempts to assess in a systematic and quantitative way how exposure affects some of these intermediate end points, such as the induction of genetic damage. For example, Angelini et al. (2016) recently reviewed the published data on benzene exposure and MN, but the authors restricted their findings to a nonquantitative comparison between exposed and unexposed individuals.

We hypothesize that it would be useful to investigate the exposure–response curve (ERC) of benzene-induced CA or MN. An advantage of using biomarkers of effect for a (more) quantitative risk assessment is that they may offer more precise information on the exact shape of the exposure–response relation at low exposures, especially when the clinical outcomes of interest are relatively rare at low exposure levels. For benzene this is important because there is evidence that the exposure–response relation with cancers might be nonlinear (Kim et al. 2006; Rothman et al. 1998; Vlaanderen et al. 2010). CAs and MN are additionally relevant to study because both have been found to be predictive of an increased risk of cancer in large prospective cohort studies (Bonassi et al. 2008, 2007).

We aimed to summarize all currently available human observational data and use this information to quantitatively describe the relation between occupational exposure to benzene and CA and MN frequencies. This would be a first step in the use of preclinical cancer-predictive end points in risk assessment. We also assessed the between-study heterogeneity effect estimates and conducted sensitivity analyses to assess how various study characteristics affect the estimated ERC. This information is important when evidence from this literature (either from a single study or from a meta-analysis) is interpreted in a risk assessment context.

Methodology

Study Identification and Selection

A literature search was conducted from until August 2018 in PubMed and Scopus, using the following search terms for “benzene” AND chromosomal aberrations: (“chromosomal aberrations” OR “chromosome*”), and the following search terms for “benzene” AND micronuclei: (“micronuclei” OR “micronucleus”). References in all identified publications were checked for additional relevant studies. Publications were selected for inclusion in the meta-analysis if they had been peer reviewed and when they fulfilled the following criteria: a) the study population had been occupationally exposed to benzene; b) exposure assessment was based on quantitative benzene exposure estimates in air or relevant benzene biomarker measurements; and c) results from cytogenetic tests were available and shown per exposure category. Studies that were excluded were those that measured general air pollution (e.g., occupational exposure to traffic air), studies with considerable co-exposure to other carcinogenic substances (IARC classification 1 or 2a), studies that reported benzene exposure at levels that could not reasonably be distinguished from general background levels (; IARC 2018). Case–control studies were also excluded from further analysis. If more than one publication was published on the same cohort, the most recent update was included.

We decided to restrict our analyses to classical CA assessments and excluded fluorescent in situ hybridization (FISH) CA studies. FISH technology is based on fluorescent DNA probes to allow for the detection aberrations in specific chromosomes (e.g., an alteration in the number of specific chromosomes or loss of particular chromosomal regions). FISH has been reported to be more sensitive in detecting the genotoxic effects of benzene, compared with classical CA analysis (Smith et al. 1998). Eleven publications (of which about half were from a single study) have employed FISH to report changes in specific chromosomes in relation to quantitatively assessed benzene exposure. However, the large heterogeneity in the various FISH protocols used (e.g., in terms of detecting specific chromosomes or measuring specific changes ranging from alterations in the number of specific chromosomes), acquisition of specific translocations or loss of particular regions of certain chromosomes, and assessments in different stages of the cell cycle (interphase vs. metaphase) prohibited the inclusion of these studies into this meta-analysis.

Study Evaluation

Each study to be included in the meta-analysis was evaluated for overall study quality and for the (expected) quality of the exposure and cytogenetic assessments. We selected a number of criteria for both end points.

Exposure assessment.

A tiered approach was used to evaluate exposure assessment quality (see Table S1). First, we assessed whether benzene exposure was the main exposure or if other co-exposures may have existed. Studies designed primarily to evaluate the effects of benzene exposure, that is, in which benzene exposure was the main exposure were classified as “A.” These were mostly studies in the shoemaking or petroleum industries. Studies that were classified as “B” were studies in jobs where significant co-exposure to other genotoxic (air) pollutants was not reported, but could not be ruled out, for example, in filling station attendants. Studies that were classified A were further divided based on the quality of their exposure assessment, with studies assessed exposure with the specific aim of exposure assessment for the respective study with at least some (personal or stationary) benzene exposure measurements or data on benzene metabolites that had been measured in blood or urine. All other A studies were classified as .

Cytogenic assessment.

The ranking of CA and MN assessments focused on the quality of the cytogenetic analyses. Studies that were classified A met the following requirements: For CA, the number of metaphases counted should be and the culture time (Carrano and Natarajan 1988). For MN, the number of bi-nucleated cells scored should be , and cells had to be incubated initially for 44 h, followed by the addition of cytochalasin B, and, subsequently, another 24–48 h of incubation (Fenech 2007). Classification as an A-class study further required that both the frequency and variability [standard error (SE) or standard deviation (SD)] of CA or MN were reported. All studies not meeting these criteria were classified as “B” (see Table S1).

Primary study set.

Studies that scored an for exposure assessment and an A for cytogenic assessment were included in the primary analyses. The other studies were included in sensitivity analyses (i.e., the full study set). The quality of eligible studies was evaluated jointly by four reviewers for exposure assessment (R.V., J.V., L.P., and B.S.) and three reviewers for cytogenetic analysis (J.V., R.S., and B.S.).

Extraction and Preparation of Data from Selected Papers

The following information was extracted from each relevant publication that was used: the country where the study was performed; the method used for exposure assessment and cytogenetic analyses, including culture time and the types of aberrations included; and the availability of biomonitoring data. For each exposure group within a study we also extracted the number of subjects, benzene exposure level, the recorded CA and/or MN frequencies with either the SE, SD, or range of frequencies. We also extracted information on benzene exposure levels and CA or MN frequencies separately by smoking status (yes/no) when this information was available.

Several studies lacked quantitative estimates of benzene exposure but did provide information on blood or urine levels of benzene or its metabolites. For those studies, we estimated benzene exposure levels through exposure reconstruction by using the physiologically based kinetic model by Knutsen et al. (2013). When the publication provided no benzene exposure levels for controls, we assumed a background level of .

The frequency of CA and MN were expressed as the number of cells with aberrations per 100 cells (a percentage). CA or MN that were reported in other units were calculated back to 100 cells. For studies that did not report the SE for CA or MN frequencies, we calculated the SE based on the reported SD and sample size. For the studies that did not report either the SE or SD (and, hence, were not included in the primary study set), we estimated the SE under the assumption that the SD equaled the median SD for the primary study set ( for CA).

Statistical Analysis

The main goal of this paper was to summarize the available information on the exposure–response relation between benzene and markers of cell damage. Because most of the individual studies reported aberration frequencies for only two groups (one of which is either unexposed or low-exposed), there was very little opportunity to investigate nonlinear exposure–response relations using this data, unless one was willing to make strong assumptions. We therefore decided to perform a meta-analysis on summary estimates of the available exposure–response information within each study (i.e., the study-specific slope), rather than attempt to build a (hierarchical) meta-regression model (for the group-level data).

The information available consisted, for most studies, of two average exposure estimates (denoted as X0 and X1) and two estimates of the average aberration frequencies (denoted as Y0 and Y1) together with an estimate of the precision of these latter estimates (denoted as variance V0 and V1). Although the aberration frequency as a proportion could theoretically be modeled as a binomial outcome for an individual, this was not the case for the reported average aberration frequencies and would have required quasi-likelihood methods to account for overdispersion even if the individual-level data had been available. We therefore regarded the average aberration frequencies Y0 and Y1 as random variables with known variances. Standard algebra provides the rules to calculate the expected value and variance of the random variable that results from subtraction of Y0 from Y1 and dividing it by the difference between X0 and X1 (i.e., the slope factor), without reference to any specific (e.g., normal) distribution. The expected value of this random variable equals , with variance .

For example, for Bogadi (1997) the reported average proportion plus or minus the standard error (SE) of chromosomal aberrations was in the unexposed and in the exposed (with an average exposure of ). The estimated slope () is then (units %/ppm) with an estimated variance of . This slope factor and associated variance was then used in the meta-analysis.

Meta-analyses were performed using the metafor package (Viechtbauer 2019) in R (version 3.6.1; R Development Core Team). Meta-analysis using this approach failed for the full MN study set because the range of sampling variances was extremely large, making the results of the algorithm used numerically unstable. We therefore chose to analyze the MN primary study set only. We calculated the statistic to assess the heterogeneity in slope estimates between studies. Funnel plots, the trim-and-fill method, and the Egger test were used to investigate potential publication bias.

Several sensitivity analyses were performed to explore the potential bias/variance tradeoff of including all studies (the full study set) in the analysis and to evaluate the impact of individual studies on the estimated meta-slope, or consensus ERC. In the first sensitivity analysis, we allowed the full set of studies to contribute to the estimation of the consensus ERC. The second (set of) sensitivity analyses was performed by leaving out one study at a time from the sample set (i.e., jackknifing) for both the primary and full study set.

To evaluate potential effect modification by smoking, we restricted the analyses to data from studies that reported exposure and genotoxic effects by smoking status, and we evaluated the interaction between smoking and benzene exposure on the frequency of CA and MN by including smoking as a moderator variable.

Results

Chromosomal Aberrations

Of the 745 identified studies on CA and benzene exposure, 16 (1,356 individuals) were considered eligible for inclusion in the meta-analysis. An overview of all eligible studies and a flow-chart indicating how studies were selected is provided in Figure S1. Four studies (477 individuals) were included in the primary study set. Of the 16 selected studies, the vast majority (), across a range of occupations, showed higher CA frequencies in jobs involving benzene exposure when compared with unexposed jobs (see Table S2). Only the benzene-exposed fuel tanker drivers in Lovreglio et al. (2014) had lower CA frequencies in comparison with their control group. The individual ERCs for all 16 studies (the full study set) are presented in Figure 1, which clearly shows the large variation across studies in baseline CA frequencies and the derived benzene ERCs.

Figure 1.

Chromosome aberrations (number of aberrations per 100 cells) and benzene exposure (ppm). The primary study set contained studies where benzene was the primary genotoxic exposure and that had adequate assessment of both exposure and outcomes. The term other studies refers to studies not meeting these criteria. The numbers shown with each line indicate the study name. Note: CA, chromosomal aberration.

A meta-analysis of the slopes on the primary study set confirmed that benzene exposure was positively and significantly associated with an increase in CA [ (95% CI: 0.08%, 0.47%)] but with an of 94%, indicating very strong heterogeneity between studies (Figure 2). Expanding the analysis to include the full study set resulted in a similar positive but heterogenous association [ (95% CI: 0.16%, 0.42%), ].

Figure 2.

Meta-analysis of slope estimates (in units%/ppm) of benzene-induced chromosomal aberrations. The primary study set contained studies where benzene was the primary genotoxic exposure and that had adequate assessment of both exposure and outcomes. The term other studies refers to studies not meeting these criteria. The of the primary study set was 94%; the of the full study set was 91.1%. Note: CI, confidence interval.

To assess the sensitivity of the outcome of the meta-analysis, one study at the time was excluded from both the primary and the full study set using a jackknife approach. The results are presented in Table S3. These analyses showed that excluding a single study had a considerable impact on the point estimate of the meta-slope (range: 0.22–0.35) when limiting analysis to the primary study set. When considering all studies, Berlin et al. (1985) was found to be most influential on the meta-slope estimate [i.e., the average slope was 0.29, whereas excluding Berlin et al. (1985) would result in a slope estimate of 0.22, ].

Of the selected CA studies from the full study set, four (Jablonická et al. 1987; Major et al. 1994; Fracasso et al. 2010; Sram et al. 2004), contained specific information on CA aberrations for smokers and nonsmokers in the control and benzene-exposed group (see Figure S2). There was no statistical interaction effect of smoking ().

There was evidence of publication bias for the CA studies (). We attempted to adjust the results for this bias using the trim-and-fill method (see Figure S3).

Micronuclei

Of the studies that were found, 13 (2,097 individuals) were considered eligible for inclusion in the meta-analysis (see Table S4). Nine of these studies (1,672 individuals) were included in the primary study set (see Figure S4). Figure 3 indicates that most studies showed higher MN frequencies in the exposed group (), but the slopes varied considerably.

Figure 3.

Micronuclei and benzene exposure. The primary study set contained studies where benzene was the primary genotoxic exposure and that had adequate assessment of both exposure and outcomes. The term other studies refers to studies not meeting these criteria. The numbers shown with each line indicate the study name. The study by Liu et al. (1996) is not displayed in this figure because it was the only study that measured very high exposure levels (up to ). Note: MN, micronuclei.

The meta-slope of the primary study set of MN studies is 0.27% (95% CI: , 0.76%) with an of 99.5% (Figure 4). We were unable to derive a consensus slope for the full data set owing to the large range in sampling variances for some of the slope estimates.

Figure 4.

Meta-analysis of slope estimates (in units %/ppm) of benzene-induced micronuclei, based on the primary study set. The primary study set contained studies where benzene was the primary genotoxic exposure and that had adequate assessment of both exposure and outcomes. The of the study set was 99.5%. Note: CI, confidence interval; MN, micronuclei.

When applying jackknifing on the primary study set, the range of slopes varied between 0.17 [when Basso et al. (2011) was excluded] to 1.03 [when Ren et al. (2018) was excluded] (see Table S5). Three MN studies from the full study set (Bukvic et al. 1998; Basso et al. 2011; Sha et al. 2014) contained information on smoking status and corresponding MN frequencies in both the benzene exposed and unexposed group. There was no statistical interaction effect of smoking () (see Figure S5).

There was evidence of publication bias for the MN studies (). The trim-and-fill method pooled estimate was 0.24 with a corresponding SE of 0.98, whereas the meta-slope of the primary study set was (see Figure S6).

Discussion

A substantial number of occupational studies has been undertaken in the past decades to study the effect of benzene on genotoxicity, covering a range of occupations and entailing various co-exposures and different levels of benzene exposure. Although in the vast majority, CA and MN frequencies were associated with jobs involving benzene exposure (IARC 2018), a consensus quantitative relationship between exposure levels and CA or MN frequency is not yet defined.

To quantitatively summarize data on occupational benzene exposure and markers of genetic damage, we conducted a systematic review, followed by a meta-analysis and a set of sensitivity analyses. We found a positive slope for both CA [ (95% CI: 0.08%, 0.47%)] and MN [ (95% CI: , 0.76%)] in relation to benzene exposure for the high quality studies (primary study set), but the between-study heterogeneity in slope estimates was large, with an estimated of over 90% for both end points.

Including also studies with less stringent quality criteria (the full study set) the slope was comparable for CA [ for all studies (95% CI: 0.16%, 0.42%), ] in comparison with the primary study set. We were unable to derive a meta-slope for the full data set for MN owing to the large range in variances for the study-specific slopes.

There are several possible explanations for the observed differences in individual slopes, including the small number of individuals per study (on average about 40), the limited dose groups (in most cases there was only one exposure estimate for the exposed group), differences in exposure assessment methods, and uncertainty in both exposure levels and CA or MN counts. This latter issue is reflected by the results of the HUman MicroNucleus project, which was established in early 2000 to gather MN data from 25 labs representative of many countries and populations (Bonassi et al. 2001). An overall median MN frequency in nonexposed individuals of 6.5 per thousand was found with an interquartile range of between 3 and 12 per thousand. The authors concluded that most of the observed total variance could be explained by laboratory methods. The issue of uncertainty in measuring DNA damage rates might be even more critical for classical CAs given that no harmonized protocol has been developed for this biomarker and that CAs are known to be more difficult to score in comparison with MN (Fenech 2002). Indeed, we found that analyzing and comparing all CA studies included in this meta-analysis was difficult because the interpretation of (total) chromosomal aberrations varied among studies. A range of structural aberrations were included by some studies, whereas others reported only a specific type of aberration, thereby making comparisons more difficult. Still, all authors aimed to collect the same information (i.e., total classical CAs), whereas this was not the case for the FISH studies, where different end points in different chromosomes were selected in different division stages of the cells (e.g., interphase or metaphase). Hence, we decided to exclude studies based on FISH technology because they could not be sensibly combined. However, potentially, the observed analytical variance in CA measurements could be less for FISH-based studies because scoring is easier and can be automated. In addition, a chromosome-wide aneuploidy study described heterogeneity in monosomy and trisomy rates for specific chromosomes (Zhang et al. 2011). Focusing on the more sensitive chromosomes may then provide a more accurate estimate of benzene’s genotoxic effect that maybe be partly masked in CA and MN assays where all chromosomes are considered.

Exposure assessments applied in the different studies were not equivalent. For example, some studies described the measurement of personal exposure in the breathing zone of workers, and others estimated exposure based on stationary samples. Although we did take uncertainty in measured CA or MN frequencies into account in this meta-analysis, this was not possible for the exposure estimates because information on the exposure assessment error was often missing (as commonly seen in meta-analyses). This could potentially add to the observed heterogeneity between studies. Note that it has been reported that has a substantial bias when the number of studies is small because an alternative confidence intervals could replace the estimates (von Hippel 2015). We provide both estimates in this paper.

The studies in the meta-analysis were weighed on the corresponding inverse variance. We did so because we assumed the SD described the variability of counted CAs or MN between individuals; this was confirmed by a number of authors of the studies used for meta-analysis in this paper (R. Sram and P. Lovreglio, personal communication). Due to the limited number of studies, the limited number of estimates per study, and the high heterogeneity between studies, we refrained from explorations on nonlinearity of the ERC. However, some evidence is available that the deleterious effects of benzene might, in fact, be nonlinear (Kim et al. 2006; Rothman et al. 1998).

We also investigated the potential interaction effect between smoking and benzene exposure. The results were not significant for either CA or MN; however, the number of studies available with information on CA or MN frequencies specifically for smokers and nonsmokers was limited. There was some evidence of publication bias, most notably for the MN studies where the trim-and-fill method indicated a (relatively small) change in pooled estimate but a very large increase in the corresponding SE (i.e., the slope estimated decreased from 0.27 to 0.24), and the corresponding SE became much larger (from 0.19 to 0.98). This effect was mainly driven by the study of Basso et al. (2011), where the observed outcome was much larger for the reported SE in comparison with other studies (see Figure S6). Publication bias thus influenced the MN slope estimate, adding to the uncertainty of the meta-slope. Evidence for publication bias for MN was also found by Angelini et al. (2016). Finally, in the past few years, several studies explored possibilities beyond measuring classical CA or MN frequencies, that is, studying other biomarkers [e.g., telomere length assays (Bassig et al. 2014)], or other cell types [e.g., cells closer to acute myeloid leukemia/myelodysplastic syndrome-related stem cells (Zhang et al. 2012)]. Potentially, these studies could also be used to derive ERCs, but we chose to focus on CA and MN assays because the evidence base for these studies is larger and their association to benzene exposure is widely accepted.

We quantitively estimated the relationship between benzene exposure and induced genetic damage to inform future risk assessments; for example, because both CAs and MN have been found to be predictive of an increase in future cancer risk (although not for specific hematopoietic malignancies, given the lack of power to investigate these associations), our derived ERC might be used to inform the shape of the benzene–leukemia curve in exposure ranges where epidemiological data is limited. One strategy to include our ERC in a benzene–leukemia risk assessment would involve a Bayesian approach that uses information on the meta-slope of CA and MN as prior in estimating the benzene–leukemia ERC. We believe the inclusion of biomarkers could potentially enhance risk assessment because these studies can be informative for early biological effects (which might not be included in cohort studies on, e.g., cancer occurrence) or relevant biological pathways, and in addition they can inform the ERC at lower exposure levels because biomonitoring studies have generally more power than outcome studies at lower exposure levels.

In conclusion, we summarized all available data on CA and MN frequencies in benzene-exposed workers and provide a quantitative estimate of the exposure–response relationship. As far as the authors are aware, no other study applied the quantitative approach used in this paper for other chemicals. We believe that further insight into the quantitative relationship between real-life benzene exposure and genetic damage could provide a starting point for further inclusion of this data in benzene risk assessments; however, the considerable observed heterogeneity between studies should be carefully considered.

Click here for additional data file.

Click here for additional data file.