Introduction

During critical developmental phases following conception, the human embryo undergoes key transitions including zygotic genome activation and the differentiation of embryonic and placental trophoblast stem cells (aka ESC and TSC) from bipotential stem cells5,6. Access to human embryos is limited and ethically constrained for stress testing, necessitating the use of pluripotent embryonic ESC in high throughput screens (HTS) to assess the impact of maternal and environmental stressors on miscarriage and birth defects. A key concern is the progress of the ESC lineage through distinct Oct4-dependent pluripotent states—naïve, formative, and primed—before gastrulation, offering a model to evaluate stress risks. Previous HTS studies utilizing ESC expressing Rex-1 promoter red fluorescent protein (RFP), fluorescence ubiquitinated cell cycle indicator (FUCCI) system, and Pdgfra promoter- green fluorescent reporter (GFP) have identified risk levels associated with environmental toxicants7,8,9,10,11,12, yet comprehensive transcriptomic analyses linking risk levels, that affect ESC proliferation, to shared stress markers are lacking.

Here we discovered pan-stress transcriptomic markers that were categorized by three statistical means as linked to normal development as weak in strength and distinguished from two strong stressors. The four, diverse, weak stressors shared with normal differentiation (ND), completely concordant directional changes in expression of ninety-four gene markers which had an inferred, highly non-random co-regulation. These ninety-four weak pan-stress markers were compared with the strong stressors to produce fifty-eight markers unique to weak stimuli, twenty-one pan-stress markers concordant for all strong and weak stimuli, and fifteen stress markers discordant for strong and weak stressors. Overall, these three marker sets, once this predictive method is validated, might offer new opportunities to develop biomarkers for miscarriage and for optimizing naive pluripotent stem cells for regenerative medicine, and identify as yet uncategorized environmental stressors that are capable of disrupting early developmental processes for pluripotent stem cells that are precursors to all adult cells.

Results and discussion

Four weak and two strong stressors were identified

It was previously established that control RA is like hyperosmotic sorbitol stress in the induction of extra-embryonic endoderm (XEN) using Pdgfra-promoter-GFP extra-embryonic endoderm (XEN, 1st ESC differentiated lineage) reporter ESC and XEN marker mRNA9,13,14. Here we confirm these two stimuli are similar in the strong induction of large numbers of significant differentially expressed genes (DEG) by pValue and false discovery rate (FDR), compared with initial normal stemness (NS, AKA naïve pluripotency) (Fig. 1).

Fig. 1
figure 1

Four weak embryotoxic stimuli and normal differentiation (ND) induced fewer numbers of pValue or FDR significant DEG versus normal stemness (NS) compared to the two strong embryotoxic stimuli.

Full size image

In triplicate bulk RNAseq experiments performed after 72 h exposure using doses based on previously established no observable adverse effects levels (NOAEL) and IC75 levels for four weak embryotoxic stressors7,8, the number of significant differentially expressed genes (DEG) versus normal stemness (NS) identified with and without false discovery rate (FDR) correction was 64% to 98% lower comparing the four weak embryotoxic stimuli and normal differentiation (ND) to the two strong embryotoxic stimuli (Supplement Table 1, (Poisson regression model, pValue < 0.05). The black and red histogram bars and labels describe the numbers of significant differentially expressed genes before and after false discovery rate correction, respectively.

The other controls were ESC maintained in the initial normal, naïve pluripotent stemness (NS conditions), due to culture with leukemia inhibitory factor (LIF). ND is caused by taking away LIF while RA, without LIF is used as a positive control of XEN differentiation. All other stressors are added with LIF present, as an override of NS.

The number of significant DEG versus NS without and with false discovery rate correction (FDR) was significantly reduced by 64–99% comparing the four weak embryotoxic and ND to the strong embryotoxic stimuli (p < 0.0001, Fig. 1, Supplement Table 1). Thus, ND was similar to the four weak stimuli (BP, Cort, DEP, PFOA).

The lack of FDR significance has been reported for a similar set of environmental toxicants in non-human primate ESC15. This is a property of stress resistance in ESC as reported16,17.

The baseline for normal differentiation without stress is NS to ND, and the range for fold-change (FC) DEG was -63FC downregulated to + 78FC upregulated compared with NS (Fig. 2). In contrast RA and 300 mM sorbitol had a 10.5 × and 10.9 × higher range of FC change than ND, and the four weak stressors (at low, NOAEL, and high IC75 doses) had a 0.8–1.6 × FC higher range than ND. Thus, the override of stemness of ESC by the weak stressors (LIF present), was similar in FC range with ND (LIF removal). Similar to the three stimulus categories in Fig. 1, the FC versus NS were statistically significantly lower comparing four weak and ND to strong embryotoxic stimuli (Fig. 2; exact Wilcoxon Two-Sample Test p-value = 0.036; exp(β) 0.11, 95% confidence interval (CI) 0.05–0.23)18.

Fig. 2
figure 2

Strong stressors, RA and 300 mM sorbitol, induce ~ 10FC higher upregulated and downregulated FC than ND (B), and four weak, environmental stressors induced a similar FC range as ND (A). Increase expression in green and decreased expression in red. For four environmental stressors, the changes in gene expression are not very different (0.8–1.6FC) from ND (1 FC), but for RA and Sorbitol, the changes in gene expression are huge (10.5–10.9 FC). The Blue boxes in (A) and (B) show the identical range of up- and down-regulation of unstressed ND compared with NS. The dotted green line is the highest FC upregulated ND and red dotted line is the lowest FC downregulated gene (A, B). The solid green line is the highest FC upregulated from strong and weak stimuli and red solid line is lowest FC downregulated gene upregulated from strong and weak stimuli (A, B). Far right on (B) S300 enabled expression of 15,544 genes, of these 37 had > highest ND upregulated DEG, and 25 < than lowest ND downregulated DEG. Thus, few DEG constitute the up- and down-regulated DEG for S300 > ND.

Full size image

Analysis of volcano plots also showed a substantial increase in FDR significance in RA and sorbitol (Supplemental Fig. 1) versus the weak stressors (Supplemental Fig. 2). Analysis of the top ten significant GO groups suggests a strong developmental component to the response of the weak stressors. Many transcriptomic programs associated with naïve pluripotent ESC development in culture partitioned between weak and strong embryotoxic stress. Strong and weak stressors defined by DEG here align well with those defined by dosimetry in a previous ESC HTS19,20. In in a previous report14, it was noted that fewer transcripts/cell and genes/cell are expressed with increasing sorbitol (i.e., less with stress), but here more p-value and FDR significant genes were induced by S300 (Fig. 1) from fewer genes (i.e., more from less), and some of the increase in FC range of S300 > ND is accomplished by relatively few DEG (Fig. 2B) (i.e., the thin veneer of stress). The high FC DEG above and below the ND range populate many GO groups which are also much higher than ND, despite lower transcript number14. Thus, it appears that stress requires a more focused, economized gene usage, to accomplish simultaneous stress homeostasis and emulate a normal or imbalanced differentiation that sustain a chance of survival, until stress subsides.

The balance to 1st and 2nd lineage from 0th naive pluripotent lineage was different for weaker ND compared with strong stimuli RA and S30014. We tested whether the four environmental stressors were like ND or strong stimuli for exiting Naïve pluripotency Supplemental Fig. 3, Part 1), entering 1st (Part 2), or 2nd lineage (Part 3), inducing HOX genes (Part 4), equilibrating between ESC and 2CEL subpopulation transcriptomic markers (Part 5).

Fig. 3
figure 3

Pan-stress gene expression was discovered with the sole criterion that the four weak environmental stressors shared genes with p value significant DEG, but there was a surprising complete concordance for directional change (previous page). Step 1 (A) identified 104 genes with common pValue significance for ESC stimulated with benzo(a)pyrene (BaP), Cortisol (Cort), diethylphthalate (DEP) and perfluoro-octanoic acid (PFOA) at NOAEL levels for growth decreased by stress (defined previously in7,8 and 104 were concordant in direction. When compared to ND DEG, there were 94 shared genes and 100% were concordant. In (B) 2472 genes shared pValue significance for strong embryotoxic stimuli RA and S300. In (C) and (D) concordance of up-regulated (green) and down-regulated (red) genes are 100% for 94 genes for ND + four weak environmental stressors, or 100% concordant for 21 genes for ND + four weak + two strong embryotoxic stressors, respectively. In (C, D, E) blue numbers at top show the fraction of genes where all stimulus caused DEG concordance, yellow at the bottom show the fraction of DEG for all five stimuli (C, E) or seven stimuli (D) within one log2FC of each other, and in (F) the green number showed fraction of all DEG with shared significance. (E) Thus, at the IC75 dose, 48% of DEG are concordant and the purple boxes in diagram show lack of concordance for each stimulus and gene. In (C, D, E) an assessment of similar FC magnitude was done by testing whether all five or all seven stimuli caused an FC range within 2FC of each other. Remarkably, where there is 100% directional concordance (C), there is also 62% concordance in magnitude.

Full size image

For Part 1 (Supplemental Fig. 3) , all stimuli cause exit from Naïve pluripotency, but ND has the highest fraction of pValue significant, downregulated naïve markers14,21,22, at highest magnitude. Weak stressors have a high fraction but fewer genes with less magnitude, in contrast to RA and sorbitol which have fewer upregulated genes and less magnitude. For Part 2, RA and to a lesser extent S300 have higher fractions of 1st lineage markers14,21,22 induced to higher magnitudes than ND. For part 3, ND has the highest fraction of highest upregulated 2nd lineage markers14,21,22, four weak stressors a lower fraction of upregulated markers and RA and S300 largely suppress these markers. For Part 4, RA activated HOX genes14,23 to highest levels, S300 activated fewer HOX genes to lower levels, and ND suppressed the most HOX genes with the greatest magnitude. The four weak stressors suppressing fewer HOX genes at lower magnitudes than ND. For Part 5, 2-cell embryo like (2CEL) subpopulations of ESC, RA induces highest fraction and magnitude upregulated markers14,24 with S300 a lower fraction and magnitude, with ND highest fraction and magnitude downregulated with the four weak stressors causing downregulation but fewer and lower magnitude than ND. Thus, as in Fig. 1, DEG for five sets of markers (Supplemental Fig. 3, Part 1–5) sets ND as a weak intermediate with weak environmental stressors concordant with ND and less powerful in gene induction number and magnitude, with RA and S300 strong and largely discordant and imbalanced compared with ND. Also, ND more efficiently downregulates naive pluripotency DEG, whereas strong stressors may induce a bipotential state of naive pluripotency and XEN markers as previously suggested by scRNAseq analysis of S30036.

Hypothesized weak and strong pan-stress markers were identified, but concordant directional changes not required by the hypothesis, were also identified

The discovery of pan-stress markers for stressed naïve pluripotent mouse ESC began with the hypothesis that a sufficient, small set of mechanistically diverse stressors would share common p value significant DEG in overriding initial naïve pluripotency, and these genes would predict the risk of not-yet-tested stressors that would affect miscarriage.

There were several surprises in testing this hypothesis. First, there was a surprisingly high 104 shared p value significant genes for the four weak environmental stressors (Fig. 3A). Secondly, serendipitously, 104 significant genes were concordant in directional changes, induced by stress from initial NS. Thirdly, adding ND to the comparison, a surprisingly high 94/104 genes shared by pValue significant DEG and serendipitously all were concordant in directional DEG. Surprisingly, shared concordant DEG for four diverse weak stressors was phenocopied by ND. Perhaps the best stress response is the one that emulates ND as a homeostatic target. This is one difference between somatic and developmental toxicology. For somatic cell homeostasis, the established parenchymal function is the retrospective pre-stress homeostatic baseline. But, in the developmental toxicology of pluripotent stem cells without parenchymal function, the homeostatic baseline is necessary future parenchymal function obtained by emulating ND. It is likely that pan-stress genes regulating future homeostasis are a class of “conditional stress lethal” genes, where the knockout has little effect in a normal vivarium or normal, unstressed culture conditions but becomes essential after adding diverse stressors4,25.

To discover an overall pan-stress DEG we contrasted the four weak stressors/NS and ND/NS, with two strong stimuli RA and 300 mM sorbitol/NS (Figs. 1, 2). The 94 DEG markers from weak stressors were compared to the 2472 shared pValue significant genes for the two strong stimuli (Fig. 3B). Of the 94 genes defined by shared p value and concordance of weak stimuli, 21/36 stimuli induced DEG that were shared in significance and concordance, and 15/36 DEG had shared significance but were discordant with the two strong embryotoxic stressors. The other 58/94 concordant weak stress DEG were either not significant for DEG of one or both strong stressors.

Figure 3C shows the surprising concordance of 100% of DEG for each of 94 genes for ND and four stressors, and also a 62% similarity of magnitude for DEG (i.e., 5 stimuli within 2FC). Including the two strong stressors (Fig. 3D), there are 21 genes with 100% concordance and 29% similarity in magnitude. For the higher IC75 dose, 46% of the genes lost significance for DEG of at least one stimulus, but there was still 52% concordance of 94 genes for five weak stimuli and 48% similarity in magnitude (Fig. 3E). This suggests that the number of concordant genes decreases with stronger stimuli, or higher doses, and it will be important to test more strong stimuli and lower < NOAEL doses for robustness of pan-stress markers to obtain the set with fewest false positive or negative predictiveness.

Are there genes whose DEG identify weak and strong stress stimuli? Venn diagrams illustrate the process of selection of the 104 concordant DEG for the weak stressors, 94 DEG shared with ND, and 36 DEG shared with the two embryotoxic stimuli, for which 21 are concordant and 15 discordant (Supplemental Fig. 4). Unique strong and weak marker genes appear in the 5/15 discordant throughout all seven stimuli, and there are 15 unique weak stressor genes not shared with RA nor S300, and 2436 genes shared by RA and S300 and not the weak stressors. Potentially the 5/15/94 discordant weak and strong markers could distinguish weak and strong stressors. ENPP2 and IGF2 are examples of this. Six of the 58 weak stress concordant genes and two of the 15/94 discordant genes (HOXA5, HOXB5) could distinguish weak from strong, but more stimuli are needed to confirm this.  Hypothetically, these eight genes which may predict weak stress. Interestingly, although the predominant DEG (Supplemental Fig. 5C) direction is upregulated for the 10/104, 94 and 21/94, all eight of the potential weak stress identifiers are downregulated, which may be a second identifying criterion. Weak stress is also more permissive in function that phenocopies ND, as strong stress induces higher significance and higher FC DEG for genes that induce XEN or suppress formative pluripotency. Strong stress is more instructive of imbalanced lineage and function.

Fig. 4
figure 4

WNT5B and FZD8 are a probable ligand-receptor pair induced significantly and highly by all seven stimuli38,39. Projected protein–protein interaction for agonistic and antagonistic interactions of WNT5B and FZD8 (yellow circles) are shown for weak embryotoxic stressors (A–D), ND (E) and strong embryotoxic stressors (F, G). See red (upregulated) and blue (downregulated) gene names and (A) (activation) and (B) (binding) interaction lines. Figure modified from in silico cross-referencing of RNAseq analysis by Advaita iPathway40.

Full size image
Fig. 5
figure 5

Pan-stress transcriptomes are defined kinetically in four stages from before stress through to cellular programmed death. Stage 1 is pre-stress preparation stress response (PPSR51,55). Stage 2 is immediately after stress initiation; the integrated stress response (ISR)  is aimed at returning the stressed cell to homeostasis regarding cell autonomous energy usage and non-autonomous parenchymal function (51,52 (https://www.ebi.ac.uk/QuickGO/term/GO:0140467). Stage 3 is the novel 94 gene and 21 gene Developmentally-associated stress response (DASR) reported here, and Stage 4 is the Senescence associated secretory phenotype (SASP)53 that may have two stages56, and occurs if stage 2 and 3 do not resolve or adequately adapt cells to the stress.

Full size image

It appears that the four stressors emulate ND in several ways; 1) complete directional concordance (Figs. 2, 3) similar FC range compared with NS (Fig. 2), and a similar biological tendency to increase more formative pluripotency markers and decrease HOX and 2CEL markers (Supplemental Fig. 3, Part 3–5) during stress override of stemness.

Pan-stress markers for weak and strong stressors arise non-randomly

The largest surprise occurs when testing the null hypothesis; that the 94 concordant DEG markers of weak stress, and 21 concordant markers of weak and strong stressors together occur by random chance. By permutation analysis, the probability of observing 104 consecutive up- or down-regulated significant genes out of the 15,404 candidates for one toxicant is 3.7 × 10−28 (i.e., highly improbable). We applied a statistical calculation given that up- or down-regulation of each DEG has ~ 1/2 probability if random. This calculation indicated that once one of the five weak stimuli are set as a baseline for up- and down-regulation, the probability of the second stimulus to match that of the first was 1/294 resulting in a probability, for each stimulus, of ~ 5 × 10−29. Identical probabilities for complete concordance for all five stimuli is 5 ×  ~ 5−29, and thus cumulative probability of complete concordance for all stimuli is ~ 6 × 10−114. Similarly for the 21 pan-stress genes across 7 stimuli, the chance of this as a random event is 1/221 = 4.8 × 10−7, and thus cumulative all stimuli is 7 × 4.8 × 10−7 which is ~3.4 × 10−48. It is highly improbable that shared marker sets for weak stressors or weak and strong stressors are random.

Candidate genes that may regulate the serendipitous, concordant non-random pan-stress response

The serendipitous concordance and the non-random DEG suggests that molecular mechanisms coordinate the pan-stress response. Candidate genes for regulation of non-randomness were searched for in two functional groups, transcription factor (TF) DEG that would regulate promoters, and secreted factors (SF) that would regulate intercellular stress responses. For the TFs, 1385 mouse TF genes were downloaded from Fantom (Supplemental Methods)26 and tested for intersection with the 94 pan-stress gene subsets (Supplemental Fig. 5A). Candidate gernes for cell–cell communication were identified by intersecting 94 pan-stress subsets with 1317 mouse SFs by the same method (Supplemental Fig. 5B)27.

For TFs, three subset identifying seven TF in the pan-stress set (Supplemental Fig. 5A). Zinc finger proteins (ZFP)160 and ZFP879 are shared by all seven weak and strong stimuli in the pan-stress 21/36/94 concordant DEG, whereas homeobox (HOX)B5 and HOXA5 were in the discordant 15/36/94 DEG. HOXA3, T-cell leukemia homeobox protein (TLX)1 and Homeobox protein (MEIS)2 were specific to the subset of 58/94 weak stress-only concordant genes.

For eight SFs, and outer plasmalemma proteins, WNT5B and FZD8 are shared, and upregulated, by all seven stimuli in the pan-stress 21/36/94 concordant DEG, whereas insulin growth factor (IGF)2, ectonucleotide pyrophosphatase/phosphodiesterase (ENPP)2, and collagen (COL)1A1 DEG were in the discordant 15/36/94 DEG (Supplemental Fig. 5B). Metalloproteinase (MMP)9, tenascin (TNC) and folate receptor (FOLR)1 were specific to the subset of 58 weak stress-only significant and concordant genes.

A heat map of the seven TFs and eight SFs showed that most were downregulated and four were upregulated from NS to ND, and these included four pan-stress genes upregulated for all seven stimuli from the 21/96 pan-stress genes (Supplemental Fig. 5C). This suggests that WNT5B, FZD8 SFs and ZFP160 and ZFP879 TFs are candidates for new coordination of mostly up-regulated Pan-stress response in Supplemental Fig. 6. ZFP879 and ZFP160 are expressed in oocytes and early mouse embryos from cleavage stage to blastocyst and tends to decrease after fertilization28, transcription is active29, but both ZFP879 and ZFP160 are regulated by SAPK inhibitors through the blastocyst and are part of similar regulation (i.e., group 6, which also includes XEN markers Col4a1/a2, Agp8, Lama1, Sparc Pdgfra, and Sox17 (30 and its derived GEO data profiles). Thus, ZFP879 and ZFP160 are associated genes important in XEN differentiation in cultured embryos under stress. The importance of these seven TF and eight SF genes is suggested by their broad participation in 12 major biological process GO groups, mostly developmental in nature, with 32 sub-categorized functions (data not shown).

A heat map was performed to identify the nature of the pan-stress of the three concordant stress groups, 10/104 unique weak, 94/104 weak + ND, and 21 weak and strong pan-stress (Supplemental Fig. 6, Supplemental Table 1). The major tendency is up-regulation; 70% for the 10/104, 67% for the 94/104, and 90% for the pan-stress 21. This suggests that ESC stress adaptation is sufficiently mature to mount new transcription.  Early phases (0.5–5 h) the TSC hyperosmotic stress response has been characterized previously31,32,33,34,35 and are wholly or primarily by DEG downregulation and only by 24 h does stress adaptation enable predominant upregulated DEG3,36. Developmental programs supporting increased 1st or 2nd lineage with strong or weak stress, are key new functions. For the 94 concordant pan-stress DEG, ND was the dominant up-regulator but Cortisol and PFOA had a high level of regulation in override of NS, and DEP induced the least upregulation. For the pan-stress concordant 21, RA was the dominant u-regulator followed by ND, with S300 surprisingly low. But Cortisol and PFOA had the highest level of override of stemness to force new DEG programs similar to ND.

We next investigated the DEG for the most significant functions. The top 25 most significant DEG for each of the three controls (i.e., two strong) and four weak environmental stimuli had DEG compared by heat map (Supplemental Fig. 7). These data show that RA is the strongest inducer of TFs and SFs from supplemental Fig. 6, but developmental genes controlling initial naïve pluripotency are regulated mostly by S300 and weak stressors.

The non-randomness suggests a coordination of gene expression and cell–cell communication transcription factors and several genes are suggested as mediators of the highly coordinated, shared pan-stress response. Two cell–cell communication-mediating genes were highest for FC and significant for all seven strong and weak stimuli as suggested by several lines of evidence (Supplemental Fig. 7, Figs. 3, 4). Frizzled (FZD)8 receptor and its probable ligand wingless-related integration site (WNT5B). WNT5B is required for cell migration, proliferation, or differentiation in many cells37. WNT5B predominantly signals through the non-canonical β-catenin-independent signaling pathway and often is an antagonist of canonical β-catenin-dependent Wnt signaling. Although x-crystallography models ligand-receptor binding for WNT5B- FZD838,39, the in-silico analysis used here indicates biological activity but not direct binding through the ligand-receptor (Fig. 4)40. Like WNT5B, FZD8 also tends to mediate non-canonical WNT signaling41,42,43. Normal42 and cancer stem cells41 under stress require non-canonical WNT functions41,44 suggesting that stress in both models may prepare cells to exit stemness. Since WNT signaling is necessary to maintain naïve pluripotency or to reprogram to pluripotency45,46, hypothetical blocking of canonical signaling by non-canonical signaling might be a common pathway to override naïve ESC stemness and enable differentiation permissibly.

A UMAP cluster analysis shows that for two of the five genes with highest significance discordance for both strong and weak stressors, Igf2 , with strong XEN tendency, is strongly expressed in clusters 6 and 2, the latter of which tends to have trophoblast markers14 (Supplemental Fig. 8). In contrast ENPP2 in expressed at lower levels uniquely in clusters 0 and 1. Igf2 is mostly upregulated in hyperosmotic stressed ESC. Whereas ENPP2 is mostly upregulated in ND ESC, thus Igf2 appears to be a better stress marker for strong stress.

The pan-stress genes support differentiation, with non-canonical WNT signaling possibly decreasing initial 0th lineage ESC Naïve pluripotency. But other factors would then direct differentiation to one alternate pathway or another; 1st lineage XEN or 2nd lineage formative pluripotency. In other words, the pan-stress response here supports balanced or imbalanced differentiation but is not strongly instructive for it. The highest induction of 1st lineage XEN and suppression of 2nd lineage naive pluripotency occurs in DEG shared by RA and S300. The final support of FZD8 and WNT5B as unique enablers of stressed override of naïve pluripotency is that they are not necessary for normal development47. WNT5B regulates bone density48. FZD8 is important in stress-induced fibrosis49. Neither is necessary for embryonic survival during the embryogenic period emulated by this ESC culture method47. This is important; whereas wildtype FZD8 and WNT5B ESC respond to stress, we hypothesize that knockout of either gene will lead to impaired stress responses to all or most stressors, thus creating a “conditional stress lethal” knockout ESC stress assay discussed previously4.

The pan-stress gene set is a developmentally associated transcriptomic marker set unique from three other kinetic markers sets

Three other pan-stress transcriptomic marker sets have been reported previously. A pre-stress “health insurance policy”50 is Oct1-dependent as the Oct1 knockout made cells unable to respond to many diverse stressors51 (Fig. 5, far left, stage 1). An early post-stress response, the Integrated Stress Response (ISR), should allow the return of differentiated somatic cells or stem cells to pre-stress cellular homeostasis and features components of endoplasmic reticulum stress (Fig. 5, stage 2). This has been reviewed elsewhere52, and the ISR is established enough to have the GO group GO:0140467; Integrated stress response signaling (https://www.ebi.ac.uk/QuickGO/term/GO:0140467. Stage 4 occurs when the stress persists; the cells exhaust and activate a senescence-associated secretomic phenotype (SASP)53. The transition by somatic cells from ISR to SASP has also been studied54. To date most of these pan-stress marker sets are not well-defined for a substantial number of cell types and diverse stressors. However, the set of 94 and subset 21 pan-stress genes here were screened against the 72 genes of stage 1, the 36 and 79 genes from two versions of the ISR in stage 2, and the 125 genes from the SASP in stage 4. There were almost no shared 94 genes that were DEG here (Supplemental Fig. 9). Therefore, we hypothesize that our gene sets relate to a novel, distinct Developmentally Associated Stress Response (DASR). These categories may change with more stresses and cell types studied, but the current set is informative and indicates somatic and stem cells have homeostatic stages of preparing and adapting to ongoing stress, which enables developmentally specific homeostatic response reported here. If either stage 2 or 3 fail to resolve stress effects, then stage 4 leads to cellular and organismal death.

Discussion

We used a strategy of testing mechanistically diverse weak and strong stimuli to isolate pan-stress transcriptomic markers for weak, or both weak and strong stressors. One major surprise was the nonrandom emulation of ND after removing stemness-maintaining growth factor, by the four weak stressors overriding stemness. This suggests the most successful stress responses to weak stressors most closely resembles ND. The non-randomness of the discovered DEG markers suggests regulatory mechanisms of coordinating the response. WNT5B and FZD8 intercellular communicating factors and ZFP160 and ZFP879 TFs are immediate candidates for coordinating the pan-stress response of naïve pluripotent ESC. A lack of knockout lethality show that these two TFs and SFs are not necessary for normal development and could be used in “conditional lethal” stress assays proposed previously4,25. The non-canonical WNT signaling of WNT5B-FZD8 in stem cells under stress, would permissively block naïve pluripotency sustained by canonical WNT signaling and enable each strong or weak stressor to direct cell fate to either 1st or 2nd lineages, enabling the imbalanced differentiation previously reported during the stem cell response to strong stressors.

Biological significance of weak and strong stressors; the amazing organismal homeostasis

To characterize strong stressors; note that S200-S300 cause cell death in the first 12 h but little death through 72 h14,57. Thus, surviving cells after 12 h of strong stress adapt and survive, but at S300 do not grow. However, an analysis of dose-dependent DEG with continuously increasing or decreasing FDR-significant changing DEG through S200, S250 and S300 doses. These FDR significant DEG showed that zero growth at S300 may be due to aerobic glycolysis downregulation by GO group and KEGG pathway analysis. But these DEG at S300 induce highly significant GO groups that enable positive regulation of cell proliferation and negative regulation of apoptosis at 72 h14. In contrast a highly decreasing GO group produced by FDR significant DEG at S300 was negative regulation of cell growth. Thus, at the highest survivable stress dose at S300, the cells are not growing at 72 h, but are planning to proliferate and not planning apoptosis at 72 h, after early cell death in the first 12. Moreover, ESC have adapted to stress and are preparing, at Tfinal at 72 h, for stress subsidence so that growth and differentiation can resume with the possibility of catching up to a normal developmental trajectory14.

In general, mouse and human ESC/iPSC are more stress-resistant than their differentiated progeny58 to strong doses of reactive oxidative species (ROS)17,59, but undergo apoptosis60, without undergoing apparent senescence and retain self-renewal under stress. However, these studies emphasize strong stress effects.

This emphasis led to a deficit of analysis of cell stress response functions at stress doses not affecting growth58,61,62,63. For example, in one study, dose-dependent effects on ROS stress on cell cycle regulatory molecules occurred only at doses where there were fewer cells after Tzero59. In contrast, a feature of the study here is that weak stress effects were established at NOAEL doses from the ESC HTS7,8, where there were no observable adverse effects on cell growth for the low doses of the four weak stressors (Fig. 3A). This suggests that significant changes in GO groups for developmental functions occurred without growth effects reported by the HTS. Thus, weak stressors can have significant single DEG and the functions corresponding to significantly changing GO groups which are enriched for these genes, even at NOAEL doses without diminished growth. For both biological examples of highest survivable stress doses with no growth, and lower stress dose stress without significant growth decrease, transcriptomic follow-up reveals that stressed ESC are mounting homeostatic programs aimed at organismal, developmental survival.

Limitations of this study

This study has several limitations. The predictive nature of weak and strong pan-stress markers are unknown. Culture methods producing less initial ESC heterogeneity than those used in this study, would lead to faster transcriptomic responses of higher significance and fold change of developmental DEG and a clearer view of organismal homeostasis. How well mouse naïve pluripotent ESC predict human pluripotency is not known nor how well cultured ESC predict the equivalent ESC lineage stage in the human embryo. Further testing is needed also for female and male mouse ESC of different species to define predictivity across these variables. Although the focus here was on the shared pan-stress response, other statistical comparisons may be useful. Unique stress-specific genes and GO group functions may provide more unique markers and mechanisms for each stress type. Weak or strong DEG markers that are not shared by all weak or strong stressors, respectively, may be useful in certain circumstances. Strong and weak stress DEG that are not shared with NS to ND progression may be useful. Dose-dependence of the pan-stress was studied in Fig. 3 in little detail. The amazingly high concordance of directional change, and appreciable concordance of magnitude, were peak at the low NOAEL dose for growth, meaning that ESC grow normally at this dose. The higher IC75 dose had lower concordance and fewer significant DEG, and doses lower than the NOAEL dose were not tested by RNAseq analysis. As noted above, ESC are stress resistant compared to differentiated progeny, but the NOAEL dose is still above the relevant exposures of many drugs, diet supplements, and pollutants. Thus, it is important to test lower doses for the robustness of shared significance and concordance of the pan-stress DEG response. Clearly stressors are manifold and nearly constantly impacting cellular development and function, so simultaneous different stressors summed on the ESC lineage cell is greater than that of a single stress4,64. Thus although single NOAEL doses here are higher than environmentally relevant doses they may actually be in range with the relevance of environmental and internal milieu stressors (i.e., cortisol) when exposed simultaneously4. Finally, not all pan-stress markers are necessarily translated, and pan-stress proteomes and phosphoproteomes are needed.

Methods

Materials

FUCCI muESC were a kind gift from Dr Pierre Savatier65. ESC reagents for maintenance culture and viable fluorescent assay of ESC were described previously7,8 and in Supplemental methods. Creation and use of Rex1-RFP ESC have been described previously11,14 and used here, to supplement FUCCI muESC bulk RNAseq analysis, with scRNAseq analysis.

Methods

Embryonic stem cell culture

FUCCI mouse embryonic stem cells were cultured as described previously9,10,11. All ESC were optimized at passage for exponential growth during the stimulus period, which began 18 h after passage at 25% confluence. Other details and rationale for culture are in Supplemental methods.

RNA isolation, cDNA library prep, and RNAseq

NOAEL and IC75 demonstration doses were determined for FUCCI ESC after HTS toxicant and control exposures (see papers for methods)7,8. The pipeline for sequencing and analysis of RNAseq are from previous reports7,14 and in Supplemental methods66,67,68,69,70.

To ascertain whether other stem cell lines besides FUCCI ESC had similar transcriptomic responses, additional data was analyzed from an exposure of Rex1-RFP ESC for 72 h using NS, ND, and three doses of sorbitol for dose-dependent effects: 200 mM, 250 mM, and 300 mM. The experimental design, methods and bioinformatics for these studies are from a previous report7,14 and in Supplemental methods.

scRNA-seq data preprocessing and quality control

This methodology, from a previous paper14, is discussed in supplemental methods71.

scRNA-seq data clustering analysis

scRNAseq clustering analysis, as done in a previous paper14, is discussed in supplemental methods.

Probability of the null hypothesis

After the directional concordance of the shared significant DEG by the stressors was established, two methods were used to assess the randomness of the serendipitous concordance: a “coin flipping” test for randomness and a Permutations analysis71,72. Both are discussed in supplemental methods.

DEG analysis of normal differentiate transition from naive pluripotency to differentiated XEN or formative pluripotency with associated mRNA programs

To quantify programs of DEG markers from naïve pluripotency to XEN or formative pluripotency, we used the same reports of these in vitro and in vivo from scRNAseq data as we have previously done3,14. For associated programs we used other reports of scRNAseq analysis or markers studies for each program: HOX73, Primitive streak22,43,74,75,76, 2cell embryo-like (2CEL)24, placenta GO groups14. Stress transcriptomic data significantly impacted pathways, biological processes, and molecular interactions. mRNA as analyzed using Advaita Bio’s iPathwayGuide40.

Analysis of conservation of gene expression levels from PPSR to DASR

We examined 86 genes from a study of genes necessary for cell survival during response to diverse stressors51, utilizing bulk RNAseq data from FUCCI muESC, and found 72 genes to be shared. For each stimulus, we compared gene counts (mRNA expression) to NS. Genes with counts below 10 were excluded from the analysis. SASP gene list was from a consensus of several tissues53.

Misc analysis and figures

Basic data analysis and graphing was done in MS Excel, and figures were formatted in MS PowerPoint and Adobe Photoshop Elements13 Editor. Advaita iPathwayGuide40 was used for Fig. 4. Volcano Plots were produced by GraphPad Prism (9.2.0). Gene Ontology (GO) and Kegg Pathway analysis was done using NIH’s DAVID 202177. Heatmaps were made using ClustVis at https://biit.cs.ut.ee/clustvis/40,78. Venn diagrams are made using molbiotools’ list comparator at https://molbiotools.com/listcompare.php.and by Panther at https://geneontology.org/. More detail is in supplemental methods79.