A comprehensive search strategy was used to identify articles in the EMBASE, PubMed, Web of Science, and Cochrance Library databases in December, 2021. The search strategy was provided in the supplementary material (Supplementary Table A.1). In addition, the reference lists of the included articles were screened manually in order to identify the potentially relevant new articles.

There were two screening phases in the study selection, namely initial screening based on the title and abstract, which was then followed by full-text screening of the eligible articles for final inclusion. In each phase, two observers (WH and YW) independently assessed each article. Discrepancies were resolved through discussion, or by consulting a third investigator (FQ).

The studies were included for full text screening when they met all the inclusion criteria: (1) original full paper presented unique data; (2) all female rodents (with all species and all ages) were exposed to DEHP/MEHP; (3) at least one of the following outcome measures related to folliculogenesis and ovarian steroidogenesis was examined (including ovary weight, serum estradiol levels, serum progesterone levels, mean estrous cyclicity, postnatal day (PND) of puberty onset, age of vaginal opening (VO)); (4) controlled studies with a separate control group. In case of doubt, a full text analysis was performed. During the full text screening, the studies were excluded if they met one of the following criteria: (1) duplicates; (2) not primary study; (3) the intervention/exposure of DEHP/MEHP not one of gavage, oral or inhalation; (4) without outcomes of interest (folliculogenesis and ovarian steroidogenesis related outcomes); (5) not a rodent study; (6) unhealthy rodents like animal model of disease.

Two review observers (WH and YW) independently extracted the bibliographic data from the included studies, including the characteristics of animal model, the exposure, the study design and the relevant outcome measures (i.e., ovary weight, serum estradiol levels, serum progesterone levels, follicle counting, oocyte counting, mean estrous cyclicity, PND of puberty onset, age of VO). If disagreement occurred, two authors discussed and arrived at a consensus with the investigator (FQ). We considered each analysis of a specific outcome measure with a specific dose as a separate comparison, however, when a study was carried two or more experiments in the same condition, we would merge the data with formulae. In addition, analyses in different ovaries (even derived from the same animal) and analyses with different time windows of exposure were considered as separate comparisons. Therefore, multiple comparisons from one study were included.

If the outcomes were repeatedly measured at different timepoint(s) in the same animals, we selected the timepoint at which the measured effect was greatest (i.e., the timepoint with the strongest association with the outcome). The timepoint of greatest efficacy was selected when the absolute difference between the mean of the exposure group and the mean of the control group divided by the sum of standard deviations (SDs) was highest. The direction of the effect was not considered when selecting the timepoint of greatest efficacy by using the absolute difference [24].

All data was collected as mean and SD. Standard error of the mean was recalculated to SD. In case the number of animals was unclear, a conservative estimate would be made, usually the mean. When the data were reported as median and interquartile range, the authors have contacted for raw data. The studies were excluded from the analyses when there existed missing data, unavailable author contact details, or no response from authors or remainders within 6 weeks.

The risk of bias of the included studies were assessed independently by two observers (WH and YW) using the SYRCLE’s Risk of Bias (RoB) tool (version 2, SYRCLE, http://WWW.SYRCLE.NL), which is specifically designed for animal studies [23]. Briefly, the ten items in the RoB tool were scored with “yes” indicating low risk of bias, “no” indicating high risk of bias or “? / unclear” indicating that the item was not reported and therefore the risk of bias was unknown.

Any discrepancies were resolved among the authors to reach a consensus. If, upon further discussion, the evaluation team couldn’t reach an agreement on a risk of bias determination for a particular domain, the more conservative judgment would be selected (i.e., if one observer made a judgment of “yes” and the other made a judgment of “? / unclear”, the “? / unclear” judgment would be used). The figure of risk of bias was performed using R (version 4.1.2, R Core Team, Vienna, Austria).

Meta-analyses were performed using Review Manager 5.4.1 (The Cochrane Collaboration) [25]. For subgroup analysis, a minimum of 2 studies per subgroup was required. If meta-analysis was not possible, data were reported through a descriptive summary. We extracted the means and the SDs from each study. We combined the means and the SDs when the difference of experiments in the study were negligible. The following formula was applied to compute the pooled Mean: ${\mathrm{Mean}}_{\text{pooled}}=\left(N_1{M}_1+N_2{M}_2\right)/\left(N_1+N_2\right)$, and the pooled SD was calculated as: ${\mathrm{SD}}_{\text{pooled}}=\sqrt{\frac{(N_1-1)S{D}_1{}^2+(N_2-1)S{D}_2{}^2+\frac{N_1{N}_2}{N_1+N_2}\left(M_1{}^2+M_2{}^2-2M_1{M}_2\right)}{N_1+N_2-1}}$, where the N1, N2 is the sample size of each experiment. If there were data from more than two experiments that needed to be merged, we followed the above formula to merge the data from two of the experiments, and then merged the resulting data with the third experiment, and so on.

Heterogeneity was assessed using I-square (I${}^2$) and the significance level of the meta-analyses was set at p heterogeneity 0.05 or I${}^2$$>$50%. When the significant heterogeneity was observed, the random-effects model was used. Otherwise, the fixed-effects model was adopted. Subgroup analysis and sensitivity analysis were applied to explore the origin of heterogeneity. For the outcome measure which multiple subgroup analyses could be conducted, the p-value for statistical significance was adjusted to p 0.01, to account for the potential false positive results.

The assessed subgroup variables included: animal species (rats or mice), strains, route of exposure (gavage, oral, diet, inhalation), dosage of treatment (below the tolerable daily intake (TDI) of 50 ug/kg/d [26], between 50 ug/kg/d and 5 mg/kg/d, between 5 mg/kg/d and 500 mg/kg/d, above 500 mg/kg/day, time window of exposure (perinatal exposure: during gestation and lactation period; prenatal exposure: during gestation only; postnatal exposure: during lactation period only), time of outcome measurement (before or after PND21, the normal weaning period), offspring of measure (F0: female mice before mated; F1: mice born form the F0 generation; F2: the F1 generation was mated with nontreated males to generate the F2 generation; F3: the F2 generation was mated with nontreated males to generate the F3 generation).

A sensitivity analysis method was conducted to evaluate the results. It is the latest timepoint was selected instead of the greatest timepoint [24].

The quality of the included studies was rated according to the confidence rating methodology of the Office of Health Assessment and Translation (OHAT) [27]. The available studies on a particular outcome were firstly grouped by the key study design features then each grouping of the studies was given an initial confidence rating according to those features [28]. According to the OHAT, an initial high confidence could be downgraded to a moderate, low or very low confidence based on assessment of five factors (i.e., risk of bias, unexplained inconsistency, indirectness, imprecision and publication bias). In addition, four factors could help to increase confidence in the results. As the factor “plausible confounding” primarily applied to the observational studies, the factor was not accessed [28]. Furthermore, the factor “large magnitude of effect” was also not assessed. According to the Wassenaar et al. [24], it was difficult to define “large effect”, as a relatively low dose of exposure led to major public health problems. Funnel plots were used to examine the publication bias when at least 5 studies could be included.

A flow chart of the study was provided in Fig. 1. Using the comprehensive electronic database search strategy, a total of 981 unique articles were identified from EMBASE, PubMed, Web of Science, and Cochrance Library databases. Finally, 15 articles were included in this systematic review after screening, of which, 8 reported the effects on ovary weight [29, 30, 31, 32, 33, 34, 35, 36], 6 on the serum estradiol levels [19, 29, 32, 37, 38, 39], 5 on the serum progesterone levels [19, 29, 32, 37, 38], 3 on the mean estrous cyclicity [29, 30, 32] and age of VO [39, 40, 41], 2 on the oocyte counting [39, 42] and PND at puberty onset [39, 40]. The detailed characteristics of all the included studies were provided in Supplementary Table A.2.

Fig. 1.

The flow chart of study selection process.

Using the SYRCLE’s RoB tool (version 2, http://WWW.SYRCLE.NL), the main observation from the risk of bias was the existence of many “? / unclear” scores, which means that most items were not sufficiently reported, resulting in an unknown risk of bias. The individual scores of the RoB tool of the included studies were provided in Fig. 2 and Table 1 (Ref. [19, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42]). Compared to other bias, the selection bias was lower: the baseline characteristics and allocation concealment had 73% and 53% of low risk of bias respectively. The items more frequently recorded included random outcomes assessment (seven studies had a low risk of bias and one study had a high risk of bias), and incomplete outcome data (four studies had a low risk of bias and two had a high risk of bias).

Fig. 2.

Results of the risk of bias and methodological quality indicators for all included studies.

The overall effect of the exposure on ovary weight was not changed with the sensitivity analyses using the data from the latest measured timepoint (MD = –0.01; 95% CI = 0.00, 0.02; Supplementary Table A.2). However, for subgroup analyses (i.e., species, strain, route of exposure, dose of exposure, time of exposure measure, time window of exposure, offspring of measure), the effect for ovary weight was changed from a significantly negative effect to a non-significantly negative effect, so we used the latest measured timepoint, instead of the greatest measured timepoint, for the subgroup analyses. Sensitivity analyses were not conducted for other measure results because they didn’t contain the latest measure time studies [34, 36].

An initial confidence rating for all the seven outcomes were “high confidence” (Table 3, Ref. [27]), because nearly all the studies had four features: controlled exposure, exposure prior to outcome, individual outcome data, and comparison group used. However, the quality of evidence for all the outcome measures were downgraded because there were many “unclear” scores for the risk of bias (Table 1). In addition, the confidence of all the outcome measures was downgraded with respect to unexplained inconsistency because of the substantial heterogeneity (i.e., I${}^2$ $>$75%). As for the ovary weight, the different timepoints and a minimal CI overlap were observed.

Furthermore, for ovary weight, the publication bias was judged as undetected, according to the funnel plot (Supplementary Fig. A.9). The confidence for the serum levels of estradiol and progesterone were downgraded for publication bias, based on the funnel plots (Supplementary Fig. B.6 and Supplementary Fig. C.4). Consequently, the quality of evidence of ovary weight, the serum levels of estradiol and progesterone were “very low”. For other four outcome measures, the quality of evidence were “low”, which means the true effect may be different from the apparent relationship.

DEHP was identified as an index compound, since it has the most robust underlying toxicological dataset, therefore, the current TDI of DEHP at 0.05 mg/kg bw per day was reaffirmed in 2019 by the EFSA Food Contact Materials, Enzymes and Processing Aids (CEP) Panel’s updated opinion [26]. The “dose of exposure” subgroup analysis for the ovary weight showed that the dose lower than 0.05 mg/kg/day significantly increased on ovary weight and the dose higher than 500 mg/kg/day significantly decreased the ovary weight. In addition, the exposure dose between 0.05 mg/kg/day and 5 mg/kg/day significantly prolonged the mean estrus cyclicity, while the exposure dose between 5 mg/kg/day and 500 mg/kg/day significantly decreased the serum levels of estradiol and progesterone. Additional data is necessary to enhance the credibility of the results and whether the current TDI is safe enough.

Other subgroup analyses for the ovary weight showed that the exposure of DEHP/MEHP significantly decreased the ovary weight in F2 and F3, indicating that the exposure has transgenerational effects [33, 36]. However, the subgroup analysis of offspring could not be carried out for other outcome measures because of the insufficient number of studies (Table 2). According to E. Brehm et al. [33] and S. Rattan et al. [36], this transgenerational effects also impacted the follicle numbers, estrus cyclicity, and ovarian steroidogenesis in the offspring.

However, due to the limitations of the available data and low evidence of quality, the results of subgroup analyses should be interpreted with caution.

This systematic review applied SYRCLE approach as the main methodology, which is specifically designed for animal studies. In addition, some formulas were used to calculate the data from the repeated experiments in the same study. Furthermore, we analyzed the transgenerational effects of DEHP/MEHP exposure.

The main limitation of this systematic review is that the amount of the included studies is relatively small and the quality of the available data is low, therefore, other outcome measure could not be included in this study. In addition, there was significant heterogeneity across the included studies (i.e., I${}^2$ $>$75%; Table 2). Despite the utility of sensitivity analysis and subgroup analysis, the origin of heterogeneity could not be fully traced. In order to account for the heterogeneity, we applied the random effect model.

According to the funnel plots of the serum levels of estradiol and progesterone (Supplementary Fig. B.6 and Supplementary Fig. C.4), a strongly suspected publication bias was observed. The reasons for the potential publication bias may be partly attributed to the publication tendency for positive results. Meanwhile, the relatively strict inclusion criteria maybe also contribute to it.

This systematic review only focused on the rodents and excluded other species, although the rodents are specifically relevant for risk assessment purposes. Other species [43], even human-cell studies [44], would also be valuable addition to the conclusion. Out of the included studies, only the route of exposure of one study [40] is inhalation. Other routes, such as subcutaneous injection [45] and intraperitoneal injection [46] should also be fully explored.

The reason why the follicle counting and other outcome measures of folliculogenesis were not be included lied in the uncertainties, such as the time taken from each experiment, the freezing method, the slice angle, the person observed, and the data without uniform units. We excluded the mixture exposure of DEHP/MEHP and other reproduction toxicants (i.e. dibutyl phthalate, diisobutyl phthalate, butyl benzyl phthalate), because the composition and dosage of substances in different experimental groups are different [47, 48].