Binary matrix, associations, and biological coherenceBinary matrix

The binary matrix, depicting the presence or absence of endocrine and reproductive effects in EOGRT studies, encompasses 85 × 112 cells (equivalent to 85 investigations × 112 EOGRT studies) and catalogs 530 treatment-related effects. This matrix is illustrated in Fig. S1 of the Supplementary Information.

The heatmap and the boxplot in Fig. S2 of the Supplementary Information depict the total number (prevalence) and distribution of all 530 observed effects at the level of the GALS groups. On average, each effect was observed in 27 studies per group across the entire binary matrix, with a standard deviation of about 13, a minimum occurrence in 14 studies, and a maximum occurrence in 63 studies.

At the level of the individual investigations, the most frequent findings included sperm morphology (34 studies), sperm motility (27 studies), irregular estrus cyclicity (28 studies), pre-implantation loss (27 studies), decreased number of delivered pups (37 studies), post-natal loss (33 studies), delayed balano-preputial separation (28 studies), and delayed vaginal opening (30 studies).

Such a high number of observations (N) results in more robust statistical analyses because this large dataset provides a better representation of the underlying population, reducing sampling error and increasing the statistical power, reliability, and validity of the statistical estimates. However, substances with mode of actions underrepresented in the dataset cannot achieve statistical significance due to their low N. Therefore, the provided analysis cannot identify all existing associations between effect pairs.

Associations of effects in investigations

Associations were identified by using Fisher’s Exact Test across all the specific effects in the entire Binary Matrix. For the 530 observed effects in the Binary Matrix, a total of 193 statistically significant unique, pairwise associations based on Fisher’s Exact Test were observed at α ≤ 0.05. All statistically significant pairwise associations are presented in Table S3 of the Supplementary Information and visualized in the heatmap of Fig. 2. Refer to the Section ‘Association matrix’ for further details.

Fig. 2

Pairwise associations of the observed effects for all investigations are illustrated in this heatmap. It illustrates pairwise effect-effect associations within and across different generations and life stages: initial adult animals (aF0), the parental animals of the aF0 generation (P0), their developing offspring (dF1) and adult offspring (aF1), the parental animals of the aF1 generation (P1), and their developing offspring (dF2). Association is measured in the form of the Jaccard (J) coefficient, which identifies how often associated effects co-occur in the same study compared to all occurrences. Only J coefficients with p ≤ 0.05 are displayed, based on Fisher’s Exact test. Red data points indicate J coefficients with values closer to 1.0, representing very strong associations. Darker squares signify J coefficients with lower strengths. The heatmap is divided into segments by black horizontal and vertical lines, indicating the generations and life stages aF0, P0, dF1, aF1, P1, and dF2 (refer to the top x-axis and right y-axis). The black diagonal line divides the heatmap into two symmetrical parts, the upper and lower triangles. Data points in segments divided by the black diagonal represent effect-effect associations within the same generation, meaning aF0 × aF0, P0 × P0, dF1 × dF1, aF1 × aF1, P1 × P1, and dF2 × dF2. Data points in segments not divided by this black diagonal represent associations across different generations and life stages, for example, aF0 × P0, aF0 × dF1, aF0 × aF1, aF0 × P1, and aF0 × dF2. Data points on the green diagonal represent associations in identical investigations (same effects) between different generations of the same life stages, meaning aF0 × aF1, P0 × P1, and dF1 × dF2. C1 to C11 identify clusters of associations. The total number of unique data points is 193, with a mean (average) of 0.278, a standard deviation of 0.134, a minimum of 0.091, the 25th percentile at 0.200, a median of 0.231, the 75th percentile at 0.333, and a maximum value of 0.857. To ensure all effect labels are displayed in this heatmap, the y-axis presents every second effect label starting with the first, and the x-axis presents every second effect label starting from the second

Biological coherence of associated effects

The consistency of all pairwise effect associations listed in Table S3 of the Supplementary Information was assessed against established scientific understanding. Mechanistic relationships were evaluated in accordance with the published literature, particularly the OECD Guidance Document (GD) 150 (OECD2018b). This GD summarizes the investigations in EOGRT studies diagnostic and indicative of endocrine pathways, including estrogen-mediated (both agonistic and antagonistic), androgen-mediated (both agonistic and antagonistic), and steroidogenesis-related activity (refer to Table B.1 in OECD 2018b, page 70).

Associations between identical effects across generations are inherently consistent as they are driven by the same mechanism(s) of toxicity exerted by the substance. This is illustrated by associations, such as post-implantation loss observed in both P0 and P1 generations. Associations involving effects in structurally, anatomically, and functionally related parameters or those with obvious causal relationships—for example, testis and epididymis histopathology; testis or epididymis pathologies and decreased sperm quality; and pre-implantation loss in P0 leading to fewer delivered F1 pups—are logically coherent. This coherence is documented in Table S3 of the Supplementary Information.

The analysis identified statistically significant associations that exemplify causal relationships between associated effects. For instance, decreased numbers of implants directly impact the quantity of delivered pups, highlighting an unambiguous causal link (OECD 2008). Decreased offspring numbers can be traced to diminished sperm motility, emphasizing the male factor’s role in reproductive success (OECD 2008). Decreased sperm quality can directly affect fertilization processes, leading to reduced implantation rates and fewer offspring (Janny and Menezo 1994). Reduced body temperature in sucklings, a consequence of decreased milk intake, emphasizes the importance of maternal provisioning and neonatal capability in early life stages (Henning and Romano 1982). This leads to less milk present in the stomachs of pups, directly correlating with a decline in offspring viability and illustrating the cascade of effects starting from nutritional deficits. Dystocia represents an immediate health concern linked to mortality in both dam and offspring (Baker et al. 2013), and it can be caused by prolonged gestation (Galal et al. 2012). The decrease in implantation rates is linked to structural changes in the endometrium influenced by the complex interplay of endocrine pathways (Timeva et al. 2014). Additionally, a reduction in the number of ovarian follicles may predict decreased female fertility (Jirge 2016).

This analysis shows that all the statistically significant associations can be reconciled with a biological mechanism underpinning the association, and demonstrates that the statistical analysis does not yield associations without a biological underpinning.

Analyzing response patternsAssociation matrix

For the identified associations, the J coefficients were calculated as a measure of the rate of co-occurrence. The resulting symmetrical J matrix is shown in Fig. 2 and Table S3 of the Supplementary Information lists the J coefficients.

To provide a brief overview of the Association Matrix, seven association clusters, denoted C1–C7, are visible along the main diagonal in Fig. 2. Clusters C1 and C5 represent associations for reproductive parameters in aF0 and aF1 males, respectively, with the cluster in aF1 exhibiting a higher density of data points. The cluster C2 consists of only two associated effect pairs in pregnant P0 females, whereas no statistically significant associations are observed within P1. C3 and C7 relate to perinatal effects observed in dF1 and dF2 neonates, with several associating effect pairs until weaning in dF1 compared to only two in dF2. Associations regarding female reproductive tissues are visible in cluster C6 in aF1, whereas no such associations are observed in aF0.

Additionally, association clusters extend along the green diagonal in Fig. 2. Effect–effect associations found on or in the vicinity of this line represent identical or different investigations, respectively, across different generations during equivalent life stages, namely aF0 × aF1 (the symbol × stands for across), P0 × P1, and dF1 × dF2. The dense cluster C8 relating to male reproductive parameters in aF0 × aF1 and the less dense cluster C9 for female reproductive parameters are also distinct along the green diagonal, comparable to C1 and C6, respectively. The faint clusters C10 (P0 × P1) and C11 (dF1 × dF2) contain six and four associations, respectively. Remarkably, aF0 × aF1 C9 has no equivalent in P0 × P0 and the P0 × P1 C10 no counterpart in P1 × P1.

Across different generations, life stages, and investigations, numerous pairwise associations were observed, extending beyond the areas around the black and green diagonals. For instance, pre-implantation loss in P0 animals and a decreased number of dF1 pups were associated with 13 and 15 other effects across the heatmap, respectively.

Several significant associations were identified between various parameters within and across generations and life stages, all of which are detailed in Table S3 of the Supplementary Information. Associations with J coefficients of 0.60 or higher are discussed as examples here. The strongest association was noted for histopathological findings in the epididymis between aF0 and aF1, with a J coefficient of 0.857 (p < 0.001), suggesting a very strong linkage. A similarly significant association for post-natal loss across dF1 and dF2 had a J coefficient of 0.769 (p < 0.001). The causally related association from pre-implantation loss to decreased delivery rate from P1 to dF2 and from P0 to dF1 demonstrated J coefficients of 0.750 and 0.636 (both p < 0.001), respectively. Further, a logical association within the same generation was identified between absence of milk in dF1 pups’ stomachs and reduced body temperature, with a coefficient of 0.667 (p < 0.001). Associations between histopathological findings in the epididymis and testis within both aF0 and aF1 were notably strong (J = 0.625). A robust association (J = 0.625) was also seen for testicular histopathology across these generations. These associations are logically coherent, reflecting identical investigations across generations or effects in male reproductive organs that are anatomically and functionally related (for example, testis and epididymis for spermatid and sperm production). Remarkably, an increased anogenital distance in dF2 females correlated significantly with histopathological findings in the epididymis and testis of both aF0 and aF1 (J = 0.600, p < 0.01), potentially due to chemical insults on shared endocrine pathways (Table S3 of the Supplementary Information).

Generational and effects associations

In the initial step, the dataset of associations was categorized into groups based on (a) occurrences within the same generations versus different generations, as well as (b) for identical effects versus different effects. This categorization was designed to facilitate observation of the dataset’s general behavior, allowing for comparisons of its patterns with respect to generational dependency and the nature of the investigations.

The results of this analysis are presented in the boxplots and Venn diagrams of Fig. 3. The boxplots display the distribution of the J coefficients as a measure of the strength of the association. The graph includes five boxplots: one representing the entire dataset, and others respectively highlighting associations within the same generations, between different generations, involving the same effects, and involving different effects. The Venn diagrams complement these findings by visually representing the % overlap and % single occurrence of associations across these categories, thereby enhancing understanding of the dataset's complexity and the relationships between the different types of associations.

Fig. 3

Boxplots and corresponding Venn diagrams illustrate the distribution of Jaccard coefficients (J) to show the degree of co-occurrence (DCO). These visualizations span the entire heatmap (“All”), covering comparisons within the same generations (“Same generations”), across different generations (“Different generations”), and for associations involving positive findings in the same or different investigations (“Same effects” and “Different effects”, respectively). The y-axis of the boxplots depicts the J distribution. Statistical parameters are summarized in Table S4 of the Supplementary Information. The × symbol in the boxplots denotes the mean (average). Comparisons of distributions of J within the same generations versus across different generations did not reveal any significant differences at a significance level α ≤ 0.05 (#, Mann–Whitney U test). However, the distributions of correlations for the same versus different effects were found to be significantly different at α < 0.001 (***, Mann–Whitney U test). These findings indicate that, although there is no significant difference in distributions of J coefficients between the same and different generations, a significant difference exists between associations of the same and different effects. The Venn diagrams demonstrate the frequency in percent of associated effects occurring in different studies (represented by blue and orange circles) versus within the same studies (indicated by grey overlap) across the categories

(a) Associations within and across generations

For the entire dataset of statistically significant associations between two investigations which have an effect, the median J coefficient was 0.231, with a mean of 0.278 ± 0.134. Associations within the same generation had a median J of 0.226, and the mean was 0.252 ± 0.107. For associations across different generations, the median J was 0.239, and the mean was 0.290 ± 0.144, as shown in Fig. 3. The Mann–Whitney U test did not reveal any significant differences in the J coefficient distribution between associations within and across generations, at a significance level of α ≤ 0.05. The Venn diagrams in Fig. 3 illustrate that the DCO for effects co-occurring within the same generations is 24%, while it is slightly higher, at 27%, across different generations.

Thus, it was concluded that, for all generations in the entire dataset, the co-occurrence of associated effects, whether within the same or across different generations, does not significantly differ, displaying a relatively weak DCO of approximately one-fourth.

(b) Associations for Effects on the Same and Different Investigations

Associations involving an effect on the same investigation demonstrated a median J coefficient of 0.400 and a mean of 0.433 ± 0.183. In contrast, associations including effects on different investigations exhibited a median J of 0.222, with a mean of 0.264 ± 0.120, as depicted in Fig. 3. The Mann–Whitney U test identified statistically significant differences between the associations of the same effects versus different effects, at a significance level of α < 0.001. The Venn diagrams in Fig. 3 show that the DCO for associated same effects is 42%, compared to 25% for associated different effects, indicating that associations for the same effects are 1.7 times more likely to co-occur within a given study than different effects.

For all generations in the entire dataset, this suggests that effects of the same nature significantly co-occur more frequently than those of different natures.

The strength of associations between GALS groups is heterogenous

In the second step, the analysis was deepened by examining associations for individual generations, beyond the aggregate view of the entire dataset. This approach aimed to uncover more nuanced patterns of behavior within the data.

It was investigated whether the GALS groups had different J coefficients. Figure 4 presents boxplots of the J coefficients for associations between investigations categorized into ‘animal groups’, sorted by descending medians. These categories form groups based on the magnitude of their J coefficients. Specifically, they divided into two distinct groups: one where associations have a median J coefficient > 0.30, and another with a median J coefficient < 0.24. A notable difference between these groups is that categories with medians above 0.30 include associations with effects in P1 and/or dF2, whereas categories with medians below 0.24 contain associations with all other generations, excluding those with P1 or dF2.

Fig. 4

Boxplots of all Jaccard coefficients (J) for associations at the animal group level (categories), sorted by descending medians. Upon visual inspection, the categories can be divided into two groups: the first group has associations with a median > 0.30 and the second group with a median < 0.24. Notably, the distinct feature between the first and second group is that the categories with medians > 0.30 encompass associations with effects in P1 and/or dF2. Conversely, the categories with medians < 0.24 comprise all other associations (no associations with P1 or dF2). Statistical parameters are summarized in Table S4 of the Supplementary Information

For further analysis, the categories were aggregated into two groups based on their median J coefficients—greater than 0.30 and less than 0.24. The results of this comparison are visualized in the boxplots and Venn diagrams in Fig. 5.

Fig. 5

Boxplots and their respective Venn diagrams across the categories: ‘All,’ encompassing the entire set of associations; ‘J > 0.3,’ which includes categories with medians greater than 0.3 (indicating associations across different generations with at least one effect in P1 or dF2); and ‘J < 0.24,’ covering associations with medians less than 0.24 (encompassing all other associations). The Mann–Whitney U test, comparing the distribution of J coefficients between the ‘J > 0.3’ and ‘J < 0.24’ categories, yielded a p-value below the significance level α of 0.001 (***). This very low p-value indicates a significant difference in the distributions of J coefficients between the two groups, supporting the hypothesis that these groups are distinct in terms of their medians. Specifically, it suggests that the J coefficients for associations between different generations including effects in P1 and dF2, demonstrate stronger associations compared to all other associations. This implies that associations across different generations including at least one effect in P1 and dF2 co-occur within the same study more frequently than all other associated effects. The Venn diagrams illustrate the percentage of correlated effects occurring in different studies (represented by blue and orange circles) versus those within the same studies (indicated by the overlap) across the three categories. Here, blue and orange denote earlier and later generations, respectively, in the order of aF0, P0, dF1, aF1, P1, and dF2. Statistical parameters for these categories are provided in Table S4 of the Supplementary Information

A Mann–Whitney U test comparing the J coefficient distributions between these two groups revealed a p-value < 0.001, indicating a significant difference in the distributions. The DCO for the group with associations to P1 and/or dF2 was 44%, compared to 20% for the group with associations to all other generations, showcasing a 2.2-fold different DCO.

This comparative analysis concludes that effects associated with P1 and/or dF2 significantly co-occur more frequently within the same study, albeit at a moderate rate of around 44%, than effects associated with other generations, at a low rate of around 20%.

Comparing associations of analogous generations and life stages

Building on the previously mentioned findings that the effects associated with P1 and/or dF2 occur significantly more frequently than others, a direct comparison was conducted. This comparison examined the effect associations between (a) P0 × dF1 and P1 × dF2, (b) aF0 × P0 and aF1 × P1, and (c) aF0 × aF1, P0 × P1, and dF1 × dF2. The objective is to determine if there is a biologically plausible connection for association between GALS groups.

(a) Associations between effects measured in parents and their offspring

Associations for P0 × dF1 displayed a median J coefficient of 0.200 and a mean of 0.244 ± 0.128. The associations for the equivalent generations P1 × dF2 had a median of 0.375, and the mean was 0.418 ± 0.185 (Fig. 6). The distribution of J coefficients significantly differed between P0 × dF1 and P1 × dF2, with a significance level α < 0.01, according to the Mann–Whitney U test. The Venn diagrams in Fig. 6 illustrate that the DCO for P0 × dF1 associations is 25%, whereas it is 44% for P1 × dF2 associations, indicating approximately a 1.8-fold difference. Based on this comparison of associations in equivalent generations, it is concluded that associated effects co-occur significantly more frequently in P1 × dF2 compared to P0 × dF1 in the same study.

Fig. 6

Distribution of Jaccard (J) coefficients displaying the associations between effects in equivalent animal groups across different generations, specifically between P0 × dF1 and P1 × dF2, as well as aF0 × P0 and aF1 × P1. The strength of association, measured by the degree overlap, which is the co-occurrence of associated effects within the same studies compared to their occurrence across all studies, are higher for associations that include P1 and dF2 effects. The degree of co-occurrence (overlap) based on the number of associated effects for P0 × dF1 is 25%, for P1 × dF2 43%, for aF0 × P0 17%, and aF1 × P1 39%, as shown in the respective Venn diagrams. The distribution of J coefficients is significantly differed between P0 × dF1 and P1 × dF2, as well as between aF0 × P0 and aF1 × P1, with a significance level of p ≤ 0.01 (**) and p ≤ 0.05 (*), respectively, using the Mann–Whitney U test. The boxplot for all associations is given as reference. Statistical parameters for these categories are provided in Table S4 of the Supplementary Information

(b) Associations between effects measured in adults and the corresponding parents

Associations for aF0 × P0 displayed a median J coefficient of 0.174 and a mean of 0.169 ± 0.022. The associations for the equivalent generations aF1 × P1 had a median of 0.333, and the mean was 0.367 ± 0.164 (Fig. 6). The distribution of J coefficients significantly differed between aF0 × P0 and aF1 × P1, with a significance level of p ≤ 0.05, using the Mann–Whitney U test. The Venn diagrams in Fig. 6 illustrate that the DCO for aF0 × P0 associations is 17%, and for aF1 × P1 associations, it is 39%, meaning it is approximately 2.3-fold higher for aF1 × P1. Based on this comparison of associations, it is concluded that associated effects co-occur significantly more frequently in the same study for aF1 × P1 compared to the equivalent generation aF0 × P0.

(c) Associations between effects measured in different generations for adult, parental, and offspring investigations

Associations for the consecutive equivalent life stages and generations aF0 × aF1, P0 × P1, and dF1 × dF2 exhibited median J coefficients of 0.231, 0.343, and 0.430, with means of 0.282 ± 0.140, 0.364 ± 0.054, and 0.465 ± 0.237, respectively (Fig. 7A). These three associations correspond to data points displayed in the sections intersected by the green diagonal in the heatmap shown in Fig. 2. Like the previous results, the strength of associations is higher for associations that include P1 and/or dF2 effects. The distribution of J coefficients significantly differed between aF0 × aF1 and P0 × P1, as well as between aF0 × aF1 and dF1 × dF2, with p ≤ 0.05 in the Mann–Whitney U test. However, no significant difference was observed between P0 × P1 and dF1 × dF2, suggesting that P1 and dF2 may enhance the strength of association in a comparable manner. The Venn diagrams in Fig. 7A and B show that the DCO for aF0 × aF1 associations is 26%, compared to 36% for P0 × P1 and 49% for dF1 × dF2, which are approximately 1.4 and 1.9 times higher compared to aF0 × aF1, respectively. Based on this comparison of associations across consecutive equivalent life stages and generations, it is concluded that associated effects occur significantly more frequently together in the same study for associations involving P1 and/or dF2 effects.

Fig. 7

Panel A: Distribution of Jaccard (J) coefficients displaying the associations between effects in consecutive equivalent animal groups across different generations, specifically between aF0 and aF1 (aF0 × aF1), P0 and P1 (P0 × P1), and dF1 and dF2 (dF1 × dF2). This analysis corresponds to the associations displayed on and around the green diagonal in the heatmap shown in Fig. 2. Like the previous results, the strength of association is higher for associations that include P1 and dF2 effects. The degree of co-occurrence based on the number of associated effects for aF0 × aF1 is 26%, for P0 × P1 36%, and for dF1 × dF2 49%, as can be seen from the respective Venn diagrams. The distribution of J coefficients significantly differed between aF0 × aF1 and P0 × P1, as well as between aF0 × aF1 and dF1 × dF2, with a significance level of p ≤ 0.05 (*, using the Mann–Whitney U test). Statistical parameters for these categories are provided in Table S4 of the Supplementary Information. There is no statistically significant different distribution of J between P0 × P1 and dF1 × dF2 (#, using the Mann/Whitney U test). Panel B: Illustrating the occurrence of associated effects in the rat life cycle of the EOGRT study

Summary of results

The binary matrix contains 530 effects, among which 193 pairwise associations have been identified at a significance level of α ≤ 0.05 (Fig. 2). Across the entire matrix, the strength of these associations is relatively low, with a J coefficient of approximately 0.23 and a DCO of 26%. This analysis reveals no statistical difference between associations within the same generations and those across different generations. However, median J coefficients are statistically significantly higher for associations involving the effects in the same investigations (median J = 0.400, DCO = 42%) compared to those involving effects in different investigations (median J = 0.222, DCO = 25%) as illustrated in Fig. 3.

Further analysis of the data reveals that associations involving dF2 and/or P1 exhibit statistically significantly higher median J coefficients, ranging from 0.300 to 0.430, in contrast to all other associations, which range from 0.182 to 0.231 (Fig. 4). Consequently, two distinct groups of associations have been identified: one with median J coefficients > 0.30 and another with coefficients < 0.24. Notably, all associations with median J coefficients greater than 0.30 involve effects observed in P1 and/or dF2: the DCO for associations involving P1 and/or dF2 is approximately 44%, compared to 20% for all other associations (Fig. 5).

An analysis of equivalent generations and life stages yields similar findings: associations involving P1 and/or dF2 are associated with statistically higher strengths compared to those excluding P1 and/or dF2, with median J coefficients for P1 × dF2 (0.375) and aF1 × P1 (0.333) being higher than those for P0 × dF1 (0.200) and aF0 × P0 (0.174), respectively. Moreover, median J coefficients increase in the following order of consecutive equivalent associations: aF0 × aF1 (0.231), P0 × P1 (0.343), and dF1 × dF2 (0.430) again supporting the finding that associations with effects observed in P1 and dF2 are stronger than others.

This analysis uncovers a consistent pattern: associations across generations involving effects observed in P1 and/or dF2 are significantly stronger and exhibit higher DCOs compared to associations among other generations, indicating intergenerational differences in response to endocrine and reproductive challenges. Associations for aF0, P0, dF1, and aF1 demonstrate low median co-occurrence rates of approximately 17 to 25%, whereas associations with P1 and/or dF2 show moderate median rates of around 40 to 44%.