Retinoic Acid Fluctuation Activates an Uneven, Direction-Dependent Network-Wide Robustness Response in Early Embryogenesis

Introduction

Robustness is an important feature of regulatory pathways to ensure signal consistency in light of environmental changes or genetic polymorphisms. Retinoic acid (RA), a central regulatory signaling pathway active during embryogenesis and adult tissue homeostasis (le Maire and Bourguet, 2014; Metzler and Sandell, 2016), is synthesized from vitamin A (retinol, ROL) or other retinoids or carotenoids obtained from the diet (Blaner et al., 2016; Kedishvili, 2016; Fainsod and Kot-Leibovich, 2018; Ghyselinck and Duester, 2019). During embryogenesis, excessive or reduced RA signaling induces a wide array of embryonic malformations (Durston et al., 1989; Sive et al., 1990; Kessel and Gruss, 1991; Papalopulu et al., 1991; Marshall et al., 1992; Hollemann et al., 1998; Corcoran et al., 2002; Fainsod et al., 2020). Human disorders like vitamin A deficiency syndrome, DiGeorge syndrome, Fetal Alcohol Spectrum Disorder, Congenital Heart Disease, neural tube defects, and multiple types of cancer have been linked to reduced RA signaling (Kim et al., 2005; See et al., 2008; Urbizu et al., 2013; Pangilinan et al., 2014; Hartomo et al., 2015; Fainsod et al., 2020). RA levels are tightly regulated throughout life, controlling quantitative, spatiotemporal expression and activity of the RA metabolic network (Kedishvili, 2016; Blaner, 2019; Ghyselinck and Duester, 2019). The regulatory mechanisms that confer this tight regulation under constant environmental challenges, i.e., robustness, as well as conditions under which such mechanisms fail, are yet to be fully understood.

Retinoic acid biosynthesis involves two sequential oxidation steps: first, mainly, multiple alcohol dehydrogenases (ADH) or short-chain dehydrogenase/reductases (SDR) can oxidize ROL to retinaldehyde (RAL), followed by a subfamily of aldehyde dehydrogenases, retinaldehyde dehydrogenases (RALDH; ALDH1A) that catalyze the oxidation of RAL to RA (Kedishvili, 2016; Metzler and Sandell, 2016; Ghyselinck and Duester, 2019). ROL and RAL can also be produced from retinyl esters or ß-carotene from food sources (Blaner, 2019). RA availability is further regulated by several ROL, RAL, and RA binding proteins (Napoli, 2017). In gastrula vertebrate embryos, RA signaling is triggered by the initial expression of ALDH1A2 (RALDH2) completing the biosynthesis of RA (Niederreither et al., 1999; Begemann et al., 2001; Chen et al., 2001; Grandel et al., 2002; Shabtai et al., 2018). Substrate availability for the RALDH enzymes is controlled by multiple RAL reductases (Feng et al., 2010; Billings et al., 2013; Porté et al., 2013; Adams et al., 2014; Shabtai and Fainsod, 2018). Besides the temporal and spatial regulation of this wide array of RA metabolic network components, other RA network components, including the RA receptor families (RAR and RXR) and retinoid-binding proteins (Lohnes et al., 1995; Cui et al., 2003; Janesick et al., 2015) are under tight regulation. The RA signal is sensitive to the maternal nutritional status and environmental exposure to chemicals (ethanol and others) (Paganelli et al., 2010; Shabtai et al., 2018), requiring the necessity to adapt the RA biosynthetic and signaling network to nutritional changes and environmental insults. This adaptation and maintenance of normal signaling levels under changing conditions is termed robustness (Eldar et al., 2004; Nijhout et al., 2019).

Over the years, multiple reports have supported the self-regulation of relatively few RA metabolic and signaling network components by RA itself. Among the network component genes shown to be responsive to RA levels, is aldh1a2 whose expression has been shown to be down-regulated by increasing RA, in somite stage zebrafish embryos, neurula stage Xenopus embryos, and adult human epidermal cultures, or down-regulated or insensitive in mouse 8.5 dpc embryos, among others (Niederreither et al., 1997; Moss et al., 1998; Chen et al., 2001; Dobbs-McAuliffe et al., 2004; Pavez Loriè et al., 2009). The RA hydroxylase, CYP26A1, whose expression is commonly adjacent to the site of ALDH1A2 expression and RA production has also been studied (Hollemann et al., 1998; Moss et al., 1998). RA up-regulation of the cyp26a1 gene has been described in late gastrula Xenopus and zebrafish embryos, cranial up-regulation and caudal down-regulation in 8.5 and 9.5 dpc mouse embryos, human hepatocytes, and others (Fujii et al., 1997; Hollemann et al., 1998; Moss et al., 1998; Dobbs-McAuliffe et al., 2004; Topletz et al., 2015). Another important enzyme pair in RA metabolism commonly compared are DHRS3 and RDH10 which preferentially reduce RAL to ROL or oxidize ROL to RAL, respectively. RA up-regulation of dhrs3 was observed in late gastrula Xenopus embryos, while rdh10 showed no change in chick and was down-regulated in Xenopus, both in neurula stage embryos (Strate et al., 2009; Reijntjes et al., 2010; Kam et al., 2013). Few studies have analyzed the effects of RA on the expression of other enzymes in the RA metabolic network. While these observations support an RA-dependent self-regulatory network, the multitude of experimental systems, range of developmental stages, and RA concentrations employed preclude the detailed understanding of the self-regulation of the RA network.

Many of the studies showing RA regulation of components of the RA metabolic network relied mainly on very high, non-physiological and continued exposures to excess RA, and fewer studies employed a reduction of RA. These studies, although performed in different experimental systems, anatomical regions, developmental stages, and even under different experimental conditions, have led to insights on the self-regulation of the RA metabolic and signaling network. These studies provided initial understanding of the RA signaling robustness. However, the extent to which fluctuations in RA lead to network-wide feedback responses has not been investigated in-depth, including the limits of robustness possible. The common use of supraphysiological levels of RA to test the response of the embryo is largely due to the common observation that embryonic malformations are mild in response to fluctuations in RA within the physiological range (Durston et al., 1989; Sive et al., 1990; Kessel, 1992). RA knockdown studies taking advantage of inhibitors, inverse agonists, or enzymatic degradation (Hollemann et al., 1998; Kot-Leibovich and Fainsod, 2009; Janesick et al., 2014, 2013), revealed that the developmental malformations observed were usually mild compared to excess RA (Blumberg et al., 1997; Sharpe and Goldstone, 1997; Hollemann et al., 1998; Koide et al., 2001; Janesick et al., 2014; Shabtai et al., 2018). These observations suggest two possibilities to explain the mechanisms that can yield robustness to fluctuations in RA levels: insensitivity – the physiologically relevant fluctuations are too small to trigger a significant feedback response; or, active self-regulation of the RA metabolic network specifically aimed at overcoming the fluctuations and returning the RA levels to normal. To discriminate between these two alternative hypotheses, an unbiased RA network-wide analysis has to be performed to conclusively compare these scenarios. The promiscuity of RA network components for each step raises the possibility of partially overlapping robustness responses across biological replicates, and would suggest multiple alternative regulatory and enzymatic scenarios capable of returning the RA signal to normal developmental levels. The scarcity of studies performing a network-wide analysis means that the potential compositional complexity of the RA feedback response has not been truly explored.

In the present study, we have taken an integrated dynamic transcriptomics approach to address several unresolved questions regarding the RA signaling robustness mechanisms during early embryogenesis. We designed a pulse-chase, transient RA perturbation protocol to explicitly study the kinetics of the robustness response using an unbiased, transcriptomic-based (RNAseq), RA network-wide analysis of the robustness response to increased versus decreased RA levels. The composition of the RA metabolic network during gastrulation, and the individual and network-wide responsiveness to physiological RA fluctuations were studied. We tested the limits of this robustness in both directions of RA perturbations and examined the clutch-to-clutch variations in terms of deviations from the normal response kinetics. Lastly, we combined the results from multiple clutches and different experimental techniques to propose a regulatory scheme that may underlie the robustness response.

Materials and Methods

Embryo Culture

Xenopus laevis frogs were purchased from Xenopus I or NASCO (Dexter, MI or Fort Atkinson, WI, United States). Experiments were performed after approval and under the supervision of the Institutional Animal Care and Use Committee (IACUC) of the Hebrew University (Ethics approval no. MD-17-15281-3). Embryos were obtained by in vitro fertilization, incubated in 0.1% MBSH and staged according to Nieuwkoop and Faber (1967).

Embryo Treatments

All-trans Retinoic acid (RA), Dimethyl sulfoxide (DMSO), and 4-Diethylaminobenzaldehyde (DEAB), were purchased from Sigma-Aldrich (St. Louis, MO, United States). Stock solutions of RA, and DEAB, were prepared in DMSO. Two-hour treatments of 10 nM RA, or 50 μM DEAB, were initiated during late blastula (st. 9.5) and terminated at early gastrula (st. 10.25) by three changes of 0.1% MBSH and further incubation in fresh 0.1% MBSH for the desired time.

Expression Analysis

For each sample, 5–10 staged embryos were collected and stored at −80°C. RNA purification was performed using the Bio-Rad Aurum Total RNA Mini Kit (according to the manufacturer’s instructions). RNA samples were used for cDNA synthesis using the Bio-Rad iScript^TM Reverse Transcription Supermix for RT-qPCR kit (according to the manufacturer’s instructions). Quantitative real-time RT-PCR (qPCR) was performed using the Bio-Rad CFX384 thermal cycler and the iTaq Universal SYBR Green Supermix (Bio-Rad). All samples were processed in triplicate and analyzed as described previously (Livak and Schmittgen, 2001). All experiments were repeated with six different embryo batches. qPCR primers used are listed in Table 1.

For spatial expression analysis of hoxb1 and hoxb4 we performed whole-mount in situ hybridization as previously described (Epstein et al., 1997). Probes for in situ hybridization used were: hoxb1, clone CX19/21 (Godsave et al., 1994) cut with EcoRI, transcribed with SP6 polymerase; hoxb4, clone E16 (Gont and De Robertis, personal communication) cut with BamHI, transcribed with T7 polymerase.

RNASeq Data Analysis

Sequencing was performed at the Thomas Jefferson University Genomics Core using Illumina HiSeq 4000. Reads were mapped to the genome using the Xenopus laevis 9.1 genome with STAR alignment and a modified BLAST and FASTQC in the NGS pipeline (STAR average mapped length: 142.34). Annotation of the mapped sequences using verse identified 31,535 genes. Raw counts were further filtered for non-zero variance across all samples resulting in 31,440 scaffold IDs. We analyzed the time series transcriptomic data set for differentially expressed genes using a two-way ANOVA approach [time points: 0 (wash), 1.5, 3, 4.5 h; treatments: RA, DEAB, control; n = 6 biological replicates].

High Throughput Quantitative Real-Time RT-PCR

cDNA samples were directly processed for reverse transcriptase reaction using SuperScript VILO Master Mix (Thermo Fisher Scientific, Waltham, MA, United States), followed by real-time PCR for targeted amplification and detection using the Evagreen intercalated dye-based approach to detect the PCR-amplified product. Intron-spanning PCR primers were designed for every assay using Primer3 and BLAST for 24 genes from Retinoic Acid metabolism and target pathway (Table 1). The standard BioMark protocol was used to preamplify cDNA samples for 22 cycles using TaqMan PreAmp Master Mix as per the manufacturer’s protocol (Applied Biosystems, Foster City, CA, United States). qPCR reactions were performed using 96.96 BioMark Dynamic Arrays (Fluidigm, South San Francisco, CA, United States) enabling quantitative measurement of multiple mRNAs and samples under identical reaction conditions. Each run consisted of 30 amplification cycles (15 s at 95°C, 5 s at 70°C, 60 s at 60°C). Ct values were calculated by the Real-Time PCR Analysis Software (Fluidigm). Samples were run in triplicate for the 24 genes. The primers in the first pre-amplification group selectively bind for the L homoeologs of the genes. In the second pre-amplification group, primers were selected for the S homoeologs only. Third pre-amplification group binds to all 24 genes. In order to remove the technical variability caused due to L and S gene homoeologs, Ct value for each gene in a sample was selected as the median value of the three pre-amplification runs. Relative gene expression was determined by the ΔΔCt method. gapdh was used as a housekeeping gene for normalization of the data.

Data Normalization and Annotation

Data analysis on raw gene counts was performed using the R statistical analysis tool version 3.6.0 on a 64-bit Windows platform. For the RNA-seq data, the raw gene counts were first converted into log2-transformed values using the “regularized log (rlog)” transformation from the DESeq2 package, which minimizes differences between samples for rows with small counts (Love et al., 2014). The gene expression data was then normalized across samples against the experimental variation and batch effects using the COMBAT method in R using a non-parametric adjustment (Johnson et al., 2007). Following batch correction, the gene list was filtered for a minimum expression threshold to remove genes with normalized counts less than 5 across all 72 samples. The expression data for the remaining genes was normalized using quantile normalization. RNA-Seq transcript/array probes IDs were transformed to Official Gene Symbol using merged list from 3 sources: the Xenopus laevis scaffold-gene mapping, DAVID Bioinformatics Resource 6.8 (Huang et al., 2009) or AnnotationData package “org.Xl.eg.db” maintained by Bioconductor (Carlson, 2017). The original scaffold IDs were retained along with the Official Gene Symbols for cross-reference purposes.

Differential Gene Expression Analysis

The normalized data was analyzed using an Empirical Bayes Statistical (eBayes) model that considered the following two variables and their interactions: (1) Treatment (Control, RA, DEAB) and (2) Time post-treatment-washout (t = 0, 1.5, 3, 4.5 h) with n = 6 biological replicates. Differentially expressed genes were identified based on statistically significant effects of Time, Treatment or an interaction between these two factors. P-values were adjusted for multiple testing using topTable from limma (Ritchie et al., 2015) package in R (q ≤ 0.05). The significance-filtered differential gene expression data was used in an established Principal Component Analysis (PCA) approach using the prcomp function implemented in R. The samples were annotated based on a combination of treatment and time, yielding 12 distinct sample groups. For each of the selected PCs, expression of 100 top positively loaded and 100 top negatively loaded genes was visualized as a heat map.

Dynamic Pattern Analysis and Comparative Matrix of Pattern Counts Analysis

First, for each time point, for each of the treatment conditions, the gene expression data for all six clutches (A,B,C,D,E,F) was averaged. Within treatment groups, RA, DEAB, Control, the gene expression data at time points t = 1.5, 3, 4.5 h was normalized by subtracting the corresponding “t = 0 h” group. This average differential gene expression data was then discretized to three levels (+1, 0, −1) based on a fold-change threshold [±2 (up, no or down- regulation)]. Within the three treatment groups, this discretization yielded a dynamic response pattern vector for each gene, encoded by one of 27 (3 levels^3 time-points) possible ordered sets. Counts of genes in each treatment group that follow each of the 27 ^∗ 27 (=729) possibilities were compared. Functional enrichment analysis was performed for geneset in various dynamic pattern vectors in the Control conditions, using functions enrichGO and simplify from the R package “clusterProfiler” (Yu et al., 2012).

Comparative Matrix of Pattern Counts Analysis of Retinoic Acid and 4-Diethylaminobenzaldehyde After Normalizing to Control

First, for each time point, for each of the treatment conditions, the gene expression data for all six clutches (A,B,C,D,E,F) was averaged. Within treatment groups, RA and DEAB, the gene expression data were normalized to Control at each time point by subtracting the expression of the Control group at the corresponding time point. This Control-normalized differential gene expression data for both treatment groups was then discretized to three levels (+1, 0, −1) based on a fold-change threshold [± log2(1.3) (up, no or downregulation)]. This discretization yielded a dynamic response pattern vector for each gene, encoded by one of 81 (3 levels^4 time-points) possible ordered sets. Subsequently, RA and DEAB groups were compared to count the number of genes corresponding to each of the 81 ^∗ 81 (=6561) possibilities; to create an 81 × 81 matrix representing the comparative dynamic response pattern counts (COMPACT) (Kuttippurathu et al., 2016). For a given COMPACT matrix of comparative conditions RA (vs. Control) and DEAB (vs. Control), the element at the ith row and jth column of the matrix contains the number of genes that show an “i”th pattern in DEAB and “j”th pattern in RA. For a coarse-grained version of the detailed 81×81 COMPACT, pattern-vector counts for each treatment group were further aggregated based on the first time-point, yielding 9 groups of pattern vectors per treatment group. The pair of treatment groups (RA and DEAB) was then compared to count the number of genes corresponding to each of the 9 ^∗ 9 (=81) possibilities.

Retinoic Acid Network Map and Visualization

A schematic representation for the position and function of genes involved in RA biosynthesis, metabolism, translocation, and transcription during gastrula stages was derived by data mining. The gastrula expression of each gene was confirmed by qPCR. For each time point, for each treatment, expression value for each gene was mapped to the corresponding label in the schematics using a color scale.

Gene correlation and clustering analysis – Gene expression data was analyzed for both with and without Control-normalization using Weighted Gene Coexpression Network Analysis (WGCNA) (Langfelder and Horvath, 2008) to identify modules of genes with highly correlated differential expression. We used a soft threshold value of 9 (8 for Control normalized data) to identify the initial gene coexpression modules, followed by a dissimilarity threshold of 0.25 to merge the initial modules into the final set of gene coexpression modules. Identified modules were further correlated with the traits (batch, time, treatment).

Robustness Score Calculation and Principal Curve-Based Trajectory Analysis

RNA-seq expression data of (clutches: A,B,C,D,E,F) and normalized HT-PCR expression data (clutches: G,H,I,J,K,L) were independently genewise Z-score transformed. Transformed data was then combined to result into 144 samples (12 clutches ^∗ 3 treatments ^∗ 4 time points). We further selected two subsets of genes: (1) “RA metabolism genes” [aldh1a2.L, aldh1a3.L, crabp2.L, crabp2.S, cyp26a1.L, cyp26a1.S, cyp26c1.L, cyp26c1.S, dhrs3.L, rbp1.L, rdh10.L, rdh10.S, rdh13.L, rdh14.L, sdr16c5.L, stra6.L]. This set represents the feedback regulatory mechanism utilized in response to the RA levels perturbations. (2) “HOX genes” [hoxa1.L, hoxa1.S, hoxa3.S, hoxb1.S, hoxb4.S, hoxd4.L] representing the phenotypic outcome from the treatment. For each gene set, the PCA scores for the combined data was first calculated using R function prcomp and then the scores from the first three principal components (PC1, PC2, and PC3) are used to learn a 3-dimensional principal curve using the function principal_curve from R package princurve. Two sets of points are specified for the “start” parameter for this function, which determines the origin and direction of the curve: (1) The centroid of the 0 h-Control samples from all twelve clutches, and (2) centroid of all the remaining 132 samples. For each geneset (RA metabolism genes and HOX genes), for each clutch (A-L), for each treatment (RA or DEAB), a “net absolute expression shift” can be calculated as the sum of distances along the principal curve, of the treatment samples from the corresponding control samples:

where, $\mathrm{\lambda }.ho{x}_{cl}^{RA(t)}$is the arc-distance from the beginning of the curve, for the “RA treatment samples” at time point “t,” for clutch “cl” for the principal curve learned for the geneset “Hox genes” (Returned as the parameter “lambda” from the principal curve function). The expression shift calculation allows ranking and sorting clutches based on “HOX shift” ($s.ho{x}_{cl}^{treatment}$) and “RA Network shift” ($s.metabolis{m}_{cl}^{treatment}$). A lower HOX shift for a clutch for a treatment indicates that the clutch is more robust to that particular treatment/perturbation. Furthermore, clutches were mapped onto the conceptual map across the spectrum of “HOX shift” and “RA Network shift” divided into hypothetical quadrants of Effective or Ineffective regulation, or Efficient or Inefficient regulation.

Results

Transient Physiological Perturbations to Study Retinoic Acid Signaling Robustness

To directly investigate the robustness of RA signaling, we designed a pulse-chase RA manipulation protocol to transiently perturb this signaling pathway beyond its normal, and developmentally changing levels, and study the feedback regulatory changes involved in restoring its normal signaling levels (Figure 1A). The network-wide changes of the transient RA disturbances and restoration of normal signaling levels were monitored by kinetic transcriptomic analysis (RNAseq). Parameters such as developmental stages and timing between samples were empirically optimized (not shown). This experimental approach expands on reports describing RA-driven regulation of its own metabolic and signaling network components. To reduce RA levels below the normal content in early embryos, we partially inhibited the oxidation of RAL to RA with the RALDH inhibitor 4-diethylaminobenzaldehyde (DEAB) (Shabtai et al., 2016). We chose 50 μM DEAB for RA reduction as our enzymatic (IC₅₀ and K_i) and phenotypic studies supported a partial inhibition (Shabtai et al., 2016, 2017). To increase RA levels for robustness analysis in a physiologically relevant manipulation, we treated embryos with 10 nM all-trans RA which is below the reported physiological RA content of early/mid gastrula Xenopus laevis embryos (about 100–150 nM all-trans RA) (Durston et al., 1989; Schuh et al., 1993; Chen et al., 1994; Creech Kraft et al., 1994, 1995; Kraft et al., 1994, 1995).

Transient pharmacological RA manipulations of early Xenopus embryos were performed by treating with RA or DEAB for 2 h from late blastula (st. 9.5; t = −2 h) to early gastrula (st. 10.25; t = 0 h), washed, and samples collected every 1.5 h (t = 1.5, 3, 4.5 h) after treatment termination (Figure 1A). For each biological repeat, all treatments, controls, and time points were collected from a single clutch (eggs) of a single female. Treatment efficiency and quality of the RNA samples was first monitored by quantitative RT-PCR (qPCR) for changes in cyp26a1.L, hoxa1.L, and hoxb4.S gene expression as internal readouts of the changes in RA signal levels. Only those biological repeats where both treatments exhibited the expected upregulation (RA) or downregulation (DEAB) of gene expression were selected for RNAseq analysis.

Principal Component Analysis (PCA) of the gene expression variation showed separation along the first principal component (PC1) which corresponds to transcriptomic changes of the normal progression through gastrulation (Figures 1B,C). RA-manipulated samples clustered with the control samples of the same developmental stage suggesting a robust maintenance of transcriptome composition. Normal transcriptomic changes as a result of progression through embryogenesis appear to be the dominant variable distinguishing the samples irrespective of the RA manipulation (Figures 1B,C). The second component (PC2) separated the sample groups at intermediate time points (t = 1.5 and t = 3 h) from those at t = 0 and 4.5 h, indicating a transient differential expression shift as the next most dominant pattern in the data (Figures 1B,C). These results show that overall, the dynamic transcriptome is maintained along the normal developmental trajectory despite the perturbations to RA levels.

The effects of RA and DEAB on the transcriptome were not readily apparent in the first seven principal components. PC8 showed some separation of the DEAB group from the RA and Control groups (Supplementary Figures 1A,C), whereas PC10 showed some separation of the RA treatment group from Control and DEAB treatments (Supplementary Figures 1B,D). The top-ranked genes along PC8 and PC10 showed distinct dynamic expression patterns across the treatments, whereas top-ranked genes along PC1 and PC2 largely corresponded to in-common dynamic changes over time (Figure 1B; Supplementary Figure 1C). These observations suggest that RA signaling robustness dampens the gene expression changes otherwise induced by abnormal RA levels in agreement with a robust system that functions to limit the gene expression changes in the majority of the transcriptome.

Retinoic Acid Signaling Exhibits Efficient Robustness

Further support for the robustness of RA signaling was obtained by phenotypic and molecular analysis of embryos treated for RA knockdown. Xenopus embryos treated with the RALDH inhibitor, DEAB (50 μM), or overexpressing CYP26A1 to render RA inactive (Hollemann et al., 1998; Yelin et al., 2005; Topletz et al., 2015; Shabtai et al., 2017). These embryos developed mild microcephaly and slight axial shortening as previously described (Shabtai et al., 2017), which are unexpectedly mild defects for such an important signaling pathway (Figures 2A–C,E). The induction of severe phenotypes by RA loss-of-function required the combined RA knockdown with both treatments targeting the metabolism at different steps (Figures 2D,E) supporting the observed robustness of RA signaling. Analysis of changes in hoxa1 or dhrs3 expression supported the conclusion that the individual RA-knockdown treatments result in moderate gene expression changes (Figures 2F–H), and severe gene expression changes are observed when both treatments were combined (Figures 2F–H). These results further support the notion that RA robustness needs to be overcome to elicit strong gene expression and phenotypic changes.

To begin exploring the limits of the robustness response to increased RA levels we treated embryos with increasing RA concentrations in the physiological range and above (1 nM – 1 μM). Treated embryos were analyzed for morphological malformations, microcephaly, trunk shortening, and loss of tail, by early tailbud stages (st. 32–33) and revealed that the RA treatments induced developmental malformations whose severity increases in correlation to the amount of RA added (Figures 3A,B). While exposure to 10 nM RA resulted in very mild malformations involving mainly the anterior head (forebrain) region, 50 nM RA induced a clear and relatively strong phenotype, and 100 nM RA resulted in severe defects, both affecting the head (severe microcephaly) and trunk (anterior-posterior axis) regions (Figures 3A,B). In similarly treated embryos, we analyzed the expression changes of genes encoding RA metabolic enzymes and RA targets (hox genes) during early (st. 10.25) and late gastrula (st.12) (Nieuwkoop and Faber, 1967). Genes positively regulated by RA at both stages, hoxb1, cyp26a1, and dhrs3 exhibited dose-dependent responses (Figures 3C–E). The genes, aldh1a2 and rdh10, were down-regulated by most RA doses (Figures 3F,H), while aldh1a3 (raldh3) was upregulated at early gastrula and was repressed late gastrula (Figure 3G). Importantly, in many instances, treatments above 10 nM RA resulted in stronger changes in gene expression or the transition to more severe malformations suggesting an upper limit to the robustness response. We also studied the effects of graded increase in RA levels on the spatial expression patterns of hoxb1 and hoxb4 (Supplementary Figure 2). We performed a visual estimation of the expression level changes compared to control embryos, or compared to the previous (lower) RA concentration. As determined by qPCR for hoxb1 (Figure 3C), for both genes we observed the expected RA concentration-dependent in situ hybridization signal intensity increase (Supplementary Figure 2). In parallel, we observed a widening of the expression domain of both genes along the anterior-posterior axis (Supplementary Figures 2A–C,E–G).

Further support for the robustness of RA signaling came from the in-depth analysis of the kinetic transcriptomic changes. We analyzed the averaged data for all biological repeats using an unbiased dynamic pattern analysis approach along all possible discretized patterns of gene expression (Kuttippurathu et al., 2016). Following normalization to t = 0 h, up- or down-regulated genes beyond the significance threshold at each time point were scored as either +1 or −1, respectively, or 0 for no significant change (Supplementary Figure 3A). Such an exhaustive approach allows us to enumerate the dominant as well as subtle patterns in the data overcoming limitations of conventional cluster analysis that could miss or mask the smaller groups of genes with distinct expression profiles over time.

Expression of a total of 4,693 genes significantly changed more than two-fold over time (compared to t0; multiple testing corrected q-value < 0.05), and surprisingly they only grouped in 10 out of the possible 27 (3^∗3^∗3) dynamic patterns for any of the experimental groups (RA, DEAB, Control) (Supplementary Figure 3A). Only 5 patterns included 100 or more genes in any experimental group, corresponding to up- or down-regulation at later time points (3 and 4.5 h) in agreement with a return to normal expression levels after the initial RA-dependent disruption (Supplementary Figure 3A). There were no significant differences in the gene counts between RA or DEAB and control groups (two-tailed Z test, p > 0.05). Gene Ontology enrichment of the control group highlighted development, cell cycle, morphogenesis, RA biosynthesis, gastrulation, and others (Supplementary Figure 3B and Supplementary Table 1). At this level of analysis, we could not observe distinctive dynamic gene expression patterns that occur only in response to the RA or DEAB treatments. Venn diagram analysis revealed limited (about 20–30%) treatment-specific gene expression changes that do not overlap with control samples for each pattern. For instance, in the subset of genes that showed late downregulation (pattern 0, 0, −1; Supplementary Figure 3C) only 216 genes in DEAB or 187 genes in RA, did not overlap with the control sample (1,175 genes) showing that only about 18% of the genes changed as a result of the treatments beyond the normal developmental changes. The expression changes relative to the end of the RA perturbation (t = 0 h) did not vary substantially (Supplementary Figure 3), consistent with the observed robustness.

Retinoic Acid Network-Wide Responses to Restore Normal Signaling Levels

We exhaustively compared the RA and DEAB treatments adapting our unbiased COMPACT approach for analyzing time-series differential transcriptomic profiles across multiple experimental conditions (Kuttippurathu et al., 2016). For statistically significant changing genes within each treatment group (two-way ANOVA; q < 0.05), we constructed discrete patterns based on differential gene expression between treatment (RA or DEAB) and control at each time point (81 theoretically possible patterns). Also here, only 32 distinct patterns with at least one gene in either comparison were observed (Figure 4A). The effect threshold, 1.3-fold average difference, was chosen lower than 2-fold, as the differential expression analysis revealed that the RA/DEAB perturbations were leading to a smaller magnitude of changes at each time point, as compared to the larger changes occurring normally over time (Figure 1B and Supplementary Figure 3), providing additional evidence of the robustness. Intersecting the two distributions yielded a matrix where each element corresponds to a distinct pair of patterns corresponding to the gene-specific response to increased (RA vs. control) and decreased (DEAB vs. control) RA signaling (Figure 4B).

Surprisingly, the COMPACT matrix (Figure 4B and Supplementary Data 1) shows that the majority of the genes exhibited sensitivity to only one direction of the RA perturbation. A total of 193 genes only responded to the addition of RA (middle row), whereas 224 genes were only responsive to RA knockdown (middle column). Of the genes responsive to both RA and DEAB perturbations, 88 genes responded in the same direction to both treatments (quadrants a and d; Figure 4B). Only a set of 48 genes showed opposite transcriptional outcomes to RA increase versus RA decrease (quadrants b and c; Figure 4B). Importantly, the set of genes that showed upregulation by exogenous RA addition and downregulation by DEAB (quadrant b; Figure 4B) was enriched for genes involved in RA metabolism (Figures 4B,C).

Analysis of the RA signaling robustness revealed a complex, network-wide transcriptional response aimed at altering enzyme levels to redirect the metabolism of RA toward restoration of normal signaling levels. To better understand the network-wide response we proceeded to determine the components of the RA network expressed during early gastrula to specifically analyze them from the transcriptomic data. The RA network composition during early/mid gastrula was determined by data mining starting from a list of 27 possible RA metabolic and signaling network components. For 12 of the 27 genes, we found transcriptomic support in our data or additional RNAseq data sets (Xenbase.org) (Yanai et al., 2011; Savova et al., 2017; Karimi et al., 2018). All putative RA network components were validated by determining their temporal expression pattern by qPCR, analyzing both homeologs (L and S genes) when present in the X. laevis genome (Supplementary Figure 4). The basic RA metabolic network analyzed includes during early/mid gastrula 12 genes out of which 4 are singletons and 8 are present as L or S homeologs (Figure 5A; Session et al., 2016).

We analyzed the averaged expression change from all the replicates of the hox genes to monitor the change in RA signal intensity (Figure 5B) and the RA metabolic network genes to begin understanding the robustness response (Figure 5B). From this quantitative estimation of the gene expression changes we can observe the robustness of RA signaling as the expression shift decreases with recovery time for the hox targets and the RA network genes (Figure 5B). Expression changes in response to RA treatment are larger that to RA knockdown (DEAB). Also importantly, the t = 0 h time point summarizes the response during early gastrula of the RA network genes to increased and decreased RA signaling levels (Figure 5B black box). These results clearly show a network-wide response as part of the robustness of the pathway.

We mapped the expression changes in response to RA perturbation (Figures 4B,C) to the RA metabolic network (Figure 5A). Notably, based on their reported preferential enzymatic activity, not all genes exhibited their expected gene expression change (Figures 5B–D). For instance, following transient RA addition, genes encoding enzymes with a preferential retinaldehyde reductase activity (rdh14 and adhFe1) undergo at later time points downregulation, and not the expected upregulation (Figure 5C). As an independent validation, we analyzed the differential expression time series data using another unbiased approach, Weighted Gene Correlation Network Analysis (WGCNA) (Langfelder and Horvath, 2008). A WGCNA gene expression module with the highest correlation with treatment contained a similar gene set as shown in Figure 4C, supporting the COMPACT analysis (green module, Supplementary Table 2; Supplementary Data 2). Taken together, our results support a network-wide mechanism to maintain RA signaling levels to counter the effects of external perturbations and prevent teratogenic effects by a transcriptional feedback control system.

To cross-validate the observed kinetics of the RA robustness response, we performed the pulse-chase RA perturbations in other independent clutches of embryos (Figure 1A) and analyzed the expression changes by qPCR (Supplementary Figure 5). In RA-treated embryos, the two RA-regulated genes, hoxb1 and hoxb4, were upregulated at the time of treatment washing (t = 0 h) (Supplementary Figure 5A), and 2 h into the recovery (t2), both RA target genes were back to control expression levels. We studied the expression of dhrs3, which preferentially reduces RAL to ROL to attenuate RA biosynthesis (Feng et al., 2010), cyp26a1, which targets RA for degradation, and aldh1a2 that produces RA from RAL (Shabtai et al., 2016), to begin understanding the RA robustness response. Both negative RA regulators of RA biosynthesis, dhrs3, and cyp26a1 were upregulated at t = 0 by increased RA and return to normal expression levels after a 4-h recovery (Supplementary Figure 5B). As expected, aldh1a2 exhibited an RA-promoted weak downregulation at t = 0, returning to normal levels after 4 h (Supplementary Figure 5B). These results exemplify how the gene expression levels of RA metabolic enzymes are regulated to achieve normal RA levels as part of the robustness response. A complementary study performed by inhibiting RA biosynthesis with DEAB (Supplementary Figures 5C,D) revealed that all genes studied exhibited fluctuations during the recovery period. The genes hoxd1 and hoxb1 reached almost normal levels at an earlier stage than the genes encoding RA metabolic components, again showing a concerted, network-wide, response of the RA metabolic network to restore normal developmental RA levels and hox gene expression (Supplementary Figure 5C). Genes encoding anabolic enzymes, e.g., aldh1a2, were upregulated, and catabolic components, e.g., cyp26a1, were downregulated (Supplementary Figure 5D), consistent with a network-wide response model of adjusting the pathway component expression (and activity) to restore normal RA signaling levels.

Direction-Dependent Uneven Retinoic Acid Robustness Response Efficiency

To better characterize the gene responses observed in the RNAseq study to RA manipulation, we performed a detailed titration of RA levels and determined the response of the main RA network genes. Multiple embryo clutches (n = 5) were exposed to increasing concentrations of RA (1, 10, 100, and 1,000 nM) or several concentrations of DEAB (1, 5, 10, 25, and 50 μM) to induce a gradual reduction in RA levels during mid gastrula (st. 11). Gene expression changes were studied by qPCR for both homeologs (Figures 6A–F). The results show that genes encoding RA production limiting enzymes (DHRS3) or enzymes promoting RA degradation (CYP26A1) exhibit a robust, dose-dependent response to increased RA similar to well-characterized RA target genes such as hoxa1 (Figures 6A–C). In contrast, as a result of RA reduction, no clear response trend could be observed, but by 25 μM DEAB, expression of both, dhrs3 and cyp26a1 decreased by about 30% similar to hoxa1 (Figures 6A–C). For genes encoding enzymes active in RA production, the situation is more complex. rdh10 exhibited a dose-dependent reduction trend with increasing RA concentrations (Figure 6D). In contrast, reducing RA levels did not uncover a clear regulatory response. During mid gastrula, aldh1a2 exhibited a mixed response to different levels of increased or decreased RA (Figure 6E). Also, aldh1a3 is down-regulated as a result of RA addition but no clear trend is observed and exhibits a mixed response to RA reduction (Figure 6F). These results recapitulate in more detail the basic observations of the COMPACT (Figure 4B) and the pattern response analysis (Supplementary Figure 3A) where many of the genes respond to only one direction of RA perturbation.

The PCA analysis revealed that all samples (RA, DEAB, and control) from the same biological repeat (clutch) tended to cluster together, but the different clutches separated slightly from each other suggesting some clutch-dependent heterogeneity in gene expression (Figure 1B; Supplementary Figure 6). For this reason, we analyzed the transcriptional changes in the RA metabolic network within each clutch from the RNAseq data (labeled A-F). In parallel, we generated an additional set of biological replicates (n = 6) to further study the dynamic expression changes in RA metabolic network genes using high-throughput real-time PCR (HT-qPCR; Fluidigm Biomark; Supplementary Figure 7) to validate the RNAseq results. These six additional biological repeats (labeled G-L) were treated with RA or DEAB following the same experimental design used for the RNAseq study (Figure 1A). We observed a similar heterogeneity between clutches analyzed by HT-qPCR to the clutch-to-clutch variation in the RNAseq data (Supplementary Figures 6, 7).

We ordered the 12 clutches according to the earliest time point at which hoxa1 returned to control levels (Figure 6G and Supplementary Figure 8) intermingling clutches of both studies. We observed high variability in the response to the RA manipulation among the individual clutches. Interestingly, clutches with significant hoxa1 upregulation following RA addition showed limited downregulation in response to RA biosynthesis inhibition (clutches C, L, J, K, H, I in Figure 6G). Clutches E, D, F, G, A, B (Figure 6G) showed the opposite response, i.e., robust hoxa1 change following RA knockdown, while RA increase induced a milder response. The RA suppressors, e.g., cyp26a1 and dhrs3, showed significantly altered expression in clutches with larger deviations in hoxa1 expression (C, L, J, K, H) and a strong and extended response to the addition of RA. In contrast, the differential expression of RA producers (e.g., aldh1a2 and rdh10) after RA manipulation was relatively mild and did not fully align with the deviation in hoxa1 expression (Figure 6G). Such a heterogeneous correlation between the extent of feedback regulation and clutch-to-clutch variation in hoxa1 expression was observed across the RA metabolic network (Supplementary Figure 8).

Uneven Robustness in Response to Increased Versus Decreased Retinoic Acid Signaling

To better understand the different clutch responses to RA manipulation, we sought to quantitatively rank the robustness of each clutch in an integrated manner based on the expression shift kinetics of the hox genes and to analyze the feedback response of the RA metabolic network genes. We implemented a trajectory-based approach in which the samples of all biological repeats (clutches) at all time points were projected onto the first three principal components based on the expression of hox genes (hoxa1, hoxa3, hoxb1, hoxb4, hoxd4) or the RA metabolic network genes (Figure 7A and Supplementary Figure 9A). A principal curve was fit to the projected data, representing the trajectory in which the system evolves over time for all clutches combined. Each clutch was visualized separately along the principal trajectory, allowing us to compare the temporal evolution of deviations between RA manipulations and controls (Figure 7A and Supplementary Figure 9A). As a multi-gene measure of the robustness of each clutch, the net absolute distance between samples and controls was computed at all time points along the principal curve. The distance between treatments and controls decreased over time, with significant differences between clutches in the time taken to close the gap (Figure 7B and Supplementary Figure 9B). The wide range of clutch-to-clutch robustness variability assessed by analysis of multiple hox genes showed that clutches with low robustness to increased RA showed larger distances between samples and controls at any time point (Figure 7A, clutch C), compared to robust clutches with decreasing distances (Figure 7A, clutch E). In contrast, DEAB treatment resulted in an opposite robustness pattern (Figure 7A), with clutch C showing the least deviation and clutch E the widest deviation (Figure 7A bottom row and Figure 6B). Clutch J exhibited intermediate robustness responses to both RA manipulations (Figures 7A,B). The anti-correlated pattern of RA robustness suggests an uneven tradeoff in the gastrula embryo to counter changes in RA levels (Figure 7C). The hox responses to altered RA levels revealed that all biological repeats aligned along a negative-slope diagonal and covered the whole range of responses (Figure 7C). This result suggested the establishment of an RA robustness gradient among the clutches and a trade-off in the intrinsic capacity to counteract increased or reduced RA levels, i.e., uneven robustness.

A similar analysis of the genes encoding RA network components showed a scattered pattern among the biological repeats with no clear trend, but a characteristic pattern of reduced deviation over time in several clutches (Supplementary Figures 9A,B). There was no significant correlation in the net shift in feedback regulatory action between RA addition versus reduction (Supplementary Figure 9C). Importantly, the scattered clutch distribution depended on the direction of the RA manipulation. The range of clutch distribution for DEAB treatment suggests a response to very small RA reductions and an upper response limit (Supplementary Figure 9C). In contrast, the clutch distribution in response to RA addition suggests that the response requires a minimal threshold change to be activated and it also exhibits an upper threshold. These observations suggest a high sensitivity to any amount of RA reduction but a lower sensitivity to small increases in RA levels.

Independent validation of these observations was performed by transiently treating multiple additional embryo clutches (n = 5; M-Q) to disturb RA levels (RA or DEAB). Individual gene responses were analyzed by qPCR at the end of the treatment (t0) and 2 h into the recovery from the treatment (t2). Following transient RA treatment, the expression changes of RA-responsive genes revealed that in the majority of cases we could observe the expected initial up-regulation and the subsequent relative down-regulation (Figure 7D). As we observed previously in the RNAseq and HT-qPCR individual clutch analysis where some clutches (L, H, G) exhibited apparently reversed responses (Figure 6G and Supplementary Figure 8), also in this case, clutch Q exhibited a relative up-regulation during the recovery period (Figure 7D). Following the transient reduction of RA levels (DEAB), we also observed clutches (N and P) that did not exhibit the characteristic up-regulation during the recovery period (Figure 7D), but rather exhibited a relative down-regulation like clutch K (Figure 6G and Supplementary Figure 8). Importantly, the same clutches with abnormal responses (Q, N, P), exhibited normal responses to the complementary RA manipulation supporting the uneven response to increased versus decreased RA levels.

Robustness Variability Impacts the Efficiency-Efficacy of the Response to Retinoic Acid Manipulation

Plotting the robustness response outcome (trajectory analysis) of the RA target genes (hox genes) (Figure 6A) as a function of the change in the RA metabolic network required to achieve the recovery from the RA manipulation (Supplementary Figure 9A) can provide insights on the effort invested by each clutch to achieve its own specific level of RA signaling robustness (Figures 8A,B). The response to increased RA levels showed some correlation between the RA network and hox shifts (Figure 8A). Interestingly, increased RA prompted a strong feedback response in 9 out of 12 of the clutches, but only half of them achieved a high robustness outcome (Figure 8A). Also, two clutches exhibited a highly efficient robustness response with a mild feedback response (Figure 8A). In agreement with the proposed uneven robustness response, for reduced RA levels, the outcome (hox shift) as a function of the network shift revealed no clear correlation (Figure 8B). While we observe hox responses, i.e., robustness outcomes, throughout most of the range, the network response exhibits an upper limit (Figure 8B). This clutch response distribution suggests that although the RA network exhibits a mild feedback response for all the clutches, this is sufficient to achieve a sensitive, high robustness to reduced RA levels in at least 8 out of 12 clutches (Figure 8B; lower left quadrant). The comparison of the robustness response (hox) and the effort to achieve it (network response) yielded an efficiency-efficacy matrix, within which clutches whose network response was low and the hox shift was also minimal were identified as establishing an efficient robustness response (Figure 8C). At the other extreme, some clutches exhibited extensive shifts in the RA network response that were unable to minimize the hox shifts resulting in an ineffective robustness response (Figure 8C). The different robustness responses of the clutches analyzed suggested that for each clutch, the robustness response is the result of a dynamic adaptation of the components of the network available. We examined whether neighboring clutches in the efficiency-efficacy matrix (Figure 8C) employ similar strategies for feedback regulation. Interestingly, we found that the identity of the differentially regulated genes in the RA metabolic network can vary between clutches co-localized in the matrix (Figures 8D,E). For example, in response to increased RA, clutches A, E and F showed similar differential regulation of dhrs3, aldh1a2, and cyp26a1, but diverged in stra6, sdr16c5, adhfe1, rdh13, aldh1a3, and cyp26c1 (Figures 8A,D). Similarly, in response to reduced RA, clutches A, B, and D showed downregulation of dhrs3 and cyp26a1, and upregulation of rbp1, suggesting feedback aimed at reducing enzymes active in RA reduction, and increasing import of retinol. However, there was high variability in the differential expression of aldh1a2 and aldh1a3 to increase RA production (Figures 8B,E).

The observations of the uneven and only partially overlapping transcriptomic changes between clutches led us to begin exploring the alternative RA-regulated expression shifts of individual network components to achieve signaling robustness (Figure 9). Based on the analysis of clutch-specific responses, we propose a regulatory scheme to help explain how different alternative strategies of changes to the RA network activity can lead to similar robustness levels. This decision tree exemplifies how each clutch follows a decision trajectory based on the direction of RA fluctuation and clutch-specific composition of RA network variants (Figure 9). In this decision tree scheme, RA network-wide adjustments in response to increase versus decrease in RA follow distinct non-overlapping trajectories. From our results on the gradient of clutch-wise robustness (Figure 6G), it follows that each clutch is likely optimized to follow a trajectory in the decision tree with a trade-off for following other trajectories, such that high robustness to one direction of RA fluctuation may be achieved at the cost of low robustness to opposite direction of perturbation in RA levels, i.e., uneven robustness. In addition, clutches could also be optimized to follow multiple trajectories with reduced robustness capacity in each scenario, yielding moderate robustness to either direction of RA fluctuation. Taken together, these results provide strong evidence that the early embryo is capable of robust maintenance of RA levels by mounting a range of network-wide feedback regulatory responses, which are constrained by an apparent trade-off that unevenly limits robustness for one of the directions of RA fluctuation.

Discussion

Robustness of the Retinoic Acid Signaling Pathway

Signaling pathway robustness is aimed at ensuring signaling consistency and reliability overcoming changing environmental conditions or genetic polymorphisms resulting in gene expression changes or protein activity variation. Commonly, strong developmental defects by RA manipulation in Xenopus embryos, and other vertebrate models, are induced by treating with high RA concentrations in the 1 to 10 μM range (Sive et al., 1990; Taira et al., 1994). This is in contrast to the estimation that gastrula stage Xenopus laevis embryos contain RA levels around the 100–150 nM range (Durston et al., 1989; Schuh et al., 1993; Chen et al., 1994; Creech Kraft et al., 1994, 1995; Kraft et al., 1994, 1995). To study the physiological robustness of RA signaling, we determined that by treating embryos with as little as 10 nM RA (about 10% of the estimated endogenous content) we could consistently induce measurable, even if less than dramatic, expression changes in RA-regulated genes. Our results and other studies have shown that embryo treatments mimicking physiological fluctuations in RA concentrations exhibit relatively mild developmental malformations suggesting the activation of compensatory mechanisms to prevent abnormal gene expression and developmental malformations (Sive et al., 1990; Hollemann et al., 1998; Reijntjes et al., 2005). Similar mild defects were observed when we reduced RA levels by either blocking the biosynthesis (DEAB treatment) or by targeting this ligand for degradation (CYP26A1 overexpression), and these observations are supported by multiple loss-of-function studies describing mild developmental malformations induced by RA signaling reduction (Blumberg et al., 1997; Sharpe and Goldstone, 1997; Hollemann et al., 1998; Koide et al., 2001; Janesick et al., 2014; Shabtai et al., 2017, 2018). Therefore, the mild developmental defects observed by inducing physiological fluctuations, increase or decrease in RA levels can be explained by a compensatory transcriptional response, i.e., robustness. An alternative explanation could be an intrinsic insensitivity of RA-mediated gene regulation to physiological fluctuations. Our results reject the latter explanation and support the former model by showing activation of a compensatory transcriptional robustness response. We show that one approach to overwhelm this RA signaling robustness is to interfere with the metabolic/signaling network at multiple steps (e.g., DEAB and CYP26A1 induction) to hamper its ability to efficiently elicit a feedback regulatory response.

The Retinoic Acid Robustness Response Involves Autoregulatory Changes in the Retinoic Acid Metabolic Network

To study the RA signaling robustness directly, we developed an experimental protocol to induce transient, physiological disturbances in RA levels. Our study supports the conclusion that the RA metabolic and signaling network exhibits highly efficient robustness to physiological fluctuations. Transcriptome-wide analysis of the restoration of normal gene expression patterns following transient RA manipulation showed that the transcriptomes of RA-manipulated (increased or decreased) and control embryos are highly similar (PCA analysis) at each time point, supporting the robustness of RA signaling to maintain normal, non-teratogenic, target gene expression levels during early gastrulation. A deeper analysis revealed only weak, but measurable, transcriptomic changes that could be attributed to the RA manipulation. Very mild morphological malformations could be observed even when significant changes in RA target genes were observed.

Then, what is the mechanism activated to achieve RA robustness within the physiological range as observed in our study? qPCR analysis provided insights on the RA robustness mechanism showing that target genes exhibited abnormal expression levels at the end of the manipulation (t = 0), and some target genes (hox) exhibited a speedy return to normal expression after only 1.5 h from the end of the treatment. Genes encoding RA metabolic network components also exhibited abnormal expression levels at t = 0, but their return to normal expression was delayed beyond the time required for the RA targets to reach normal expression. These observations suggest that the RA metabolic network is altered through an RA-dependent feedback regulatory mechanism to restore normal RA signaling and target gene expression.

Multiple reports have described the RA-dependent regulation of individual RA metabolic components. Most studies show up-regulation of enzymes suppressing or reducing RA signaling like CYP26A1, ADHFe1, and DHRS3, or down-regulation of RA producers (anabolic enzymes) like ALDH1A2 and RDH10 as a result of RA treatment (Fujii et al., 1997; Hollemann et al., 1998; Sonneveld et al., 1998; Chen et al., 2001; Dobbs-McAuliffe et al., 2004; Strate et al., 2009; Sandell et al., 2012; Kam et al., 2013; Shabtai et al., 2017). However, many of these studies utilized non-physiological amounts of RA mostly in continuous treatments, they were performed in multiple organisms, cell types, and developmental stages, and some even obtained contradictory results (Fujii et al., 1997; Niederreither et al., 1997; Hollemann et al., 1998; Moss et al., 1998; Chen et al., 2001; Dobbs-McAuliffe et al., 2004; Pavez Loriè et al., 2009; Strate et al., 2009; Reijntjes et al., 2010; Kam et al., 2013; Topletz et al., 2015). Our experimental design based on transcriptomic analysis together with our description of the RA network during early/mid gastrula allowed us to define the RA network-wide response to signaling changes, in one organism, at a defined developmental stage, in an unbiased manner by studying all the network components expressed during gastrula Xenopus embryos.

We observed that some RA network components exhibit an oscillatory behavior close to control expression levels, probably as a result of the fine-tuning of the RA signal. Such fine-tuning of the RA signal levels coupled to the inherent delay in the transcriptional response could transiently result in the inversion of the overall signaling direction, and an oscillatory transcriptional behavior. In a few instances such paradoxical observations of RA signaling outcome following RA manipulation, i.e., overcompensation, have been reported (Lee et al., 2012; D’Aniello et al., 2013; D’Aniello and Waxman, 2015; Rydeen et al., 2015). These observations identify a very dynamic feedback regulatory network continuously fine-tuning itself to overcome perturbations.

The recovery kinetics analysis after conversion of the RNAseq data to discretized patterns revealed that only 10 out of the 27 possible expression patterns were represented in any of the treatment or control samples, and only five patterns were exhibited by a substantial number of genes (n > 100). These observations suggest that RA-regulated genes, direct or indirect, can only exhibit a limited number of possible regulatory responses although intensity differences are possible. Surprisingly, a comparative analysis of the discretized patterns using our unbiased COMPACT approach (Kuttippurathu et al., 2016) revealed that most genes (75.40%) responding to RA manipulation and recovery exhibited a response to only either increased or decreased RA levels. Only 48 genes out of 553 (8.67%) exhibited the commonly expected reciprocal responses to increased versus decreased RA levels. Many of this small group of genes encode RA metabolism or signaling components as part of the robustness response to maintain non-teratogenic RA levels.

The temporal delay in the return to normalcy between RA targets and network components and the network-wide response taking place at multiple levels provide an integrative view of the feedback regulation and robustness response. Initially, the response to RA changes will most probably be managed by the actual enzymes and factors already present in the cell. In parallel, the same RA changes will elicit a self-regulatory transcriptional response which for increased RA should involve the upregulation of RA suppressor activities (DHRS3, ADHFe1, and CYP26A1), and the complementary downregulation of RA producing enzymes (ALDH1A2 and RDH10). This transcriptional response will exhibit a slight delay due to the transcription, translation, and post-translational modifications required.

Uneven Response to Increased Versus Decreased Retinoic Acid Levels

Our study provides new insights on the network responses to increased and decreased RA levels. One important conclusion is that these responses are not simple opposites, but are rather under different regulatory rules. As determined from the unbiased discretized pattern analysis, a small proportion of genes respond inversely to increased and decreased RA levels. Otherwise, most of the RA-responsive genes respond to one direction of RA change. This observation was further supported by our relative ranking of the robustness of the different embryo clutches based on the trajectory analysis. We conclude that the RA auto-regulatory robustness response is uneven based on: First, the limited number of RA target genes that respond inversely to increased and decreased RA. Second, the response to reduced signaling is activated by even a slight RA reduction, while the response to increased RA is only activated above a threshold. Also, the response to reduced RA reaches an upper response limit, while the response to increased RA appears to have a wider range. Notably, the network responses show different thresholds depending on the direction of the RA manipulation suggesting a high sensitivity to RA reduction but a lower sensitivity to slight RA increase. Third, the alignment of the robustness responses, hox changes, fitted a diagonal with a negative slope with a significant correlation. This negative slope suggested that while a clutch might very efficiently deal with increased RA, i.e., high robustness, the same clutch struggles to compensate for a reduction in RA, i.e., low robustness, suggesting an intrinsic trade-off along the lines of a Pareto multiobjective optimization (Miettinen, 1998; Schuetz et al., 2012; Shoval et al., 2012; Tendler et al., 2015). In the present case, the Pareto optimality suggests that improving robustness toward one objective (e.g., decreased RA) can only be achieved by degrading robustness toward another (increased RA). Our results showing that some clutches achieve only moderate robustness to positive and negative fluctuations in RA, whereas others show uneven robustness, are consistent with the notion of Pareto optimality of robustness to RA fluctuations in the early embryo. Fourth, different clutches with similar robustness can mount robustness responses by incorporating different components of the RA network. The transcriptomic and HT-qPCR analyses allowed us to perform a detailed determination of the network components comprising a robustness response across different clutches to uncover its mechanism. The RA metabolic network includes multiple enzymes with overlapping biochemical activities, e.g., ALDH1A1, ALDH1A2, and ALDH1A3 enzymes that oxidize retinaldehyde to produce RA (Kedishvili, 2016, 2013; Shabtai et al., 2016; Ghyselinck and Duester, 2019). The enzymatic efficiencies and expression patterns including location, timing, and intensity may be different across the enzymes (Chen et al., 2001; Blentic et al., 2003; Romand et al., 2004; Lupo et al., 2005; Shabtai et al., 2016). The enzymes CYP26A1, B1, and C1, or DHRS3, ADHFe1, and RDH14 function to prevent excessive RA signaling (Hollemann et al., 1998; Sonneveld et al., 1998; Belyaeva et al., 2008, 2017; Shabtai et al., 2017). These redundant enzymes and factors can help establish a robustness response using different components (Figures 8D,E, 9). Finally, genetic polymorphisms can probably partially explain the different responses between clutches of our outbred experimental model, Xenopus laevis (Savova et al., 2017). Formally, technical issues could account for some variability. The six clutches analyzed by HT-qPCR, six clutches analyzed by RNAseq, and the five clutches from the qPCR analysis showed a similar distribution of the same variability supporting the possible involvement of genetic polymorphisms. We mined the Savova et al. (2017) data identifying multiple polymorphisms in RA network components (Supplementary Table 3), supporting the potential contribution of genetic variability to the differential and uneven robustness across clutches.

Retinoic Acid Signaling Changes Due to Environmental Changes and Disease Risk

While RA exposure induces dramatic developmental malformations, treatment of the same embryos with ROL or RAL, the RA precursors, requires higher concentrations (70–100× and 5×, respectively) to induce similar defects to the RA treatment (Durston et al., 1989; Yelin et al., 2005). The important biochemical difference between the compounds is that ROL or RAL have to go through RA biosynthesis to become the active ligand, while RA is already the final product. The oxidation of ROL to RAL is a reversible reaction that can reduce substrate availability for the ALDH1A enzymes, while excess RA can only be countered through inactivation by retinoic acid hydroxylating enzymes, e.g., CYP26 (Hollemann et al., 1998; Sakai et al., 2001; Dobbs-McAuliffe et al., 2004). Therefore, the reduced teratogenicity of the precursors could be the result of reduced conversion to RA as part of the feedback regulation of this network, i.e., robustness, while the robustness response is more restricted in its regulatory feedback management of RA addition.

All our experiments were performed by either partially inhibiting the endogenous levels of RA, or by slightly increasing (∼10% increase) the physiological RA content. Under these conditions, the RA robustness of the embryo efficiently normalizes the transcriptome via regulatory feedback. Based on the comparative analysis of 12 clutches, we suggest that robustness efficiency will have a threshold beyond which it will become ineffective in restoring normal RA signaling. Our clutch analysis suggests that this threshold might be strongly dependent on genetic variability affecting enzymatic activity or gene expression parameters. Consistent with these observations, a threshold or toxicological tipping point for RA signaling was recently described in a cell-based model (Saili et al., 2019). The clutch analysis also suggested that the RA network response might also be dependent on genetic variability like promoter polymorphisms. Then, the developmental malformations arising from environmental insults on RA signaling largely depend on genetic polymorphisms which will determine the composition, efficiency and threshold of the response and the actual network components comprising such a response.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/geo/, GSE154408, GSE154399, and GSE154407.

Ethics Statement

The animal study was reviewed and approved by the Institutional Animal Care and Use Committee (IACUC) of the Hebrew University.

Author Contributions

RV and AF conceived and supervised the study, designed the experiments, and the analysis methodology. LB-K, MG, TA, KK, and AF performed embryo experiments and real-time PCR assessment and developed the figures. AB performed the initial analysis of the RNAseq data. SA conducted the high-throughput PCR validation of RNAseq results. MP conducted the analysis of transcriptomics and HT-qPCR data, performed the network and trajectory analyses, and developed the figures. MP, RV, and AF interpreted the results and drafted the manuscript. All authors contributed to the article and approved the submitted version.

Funding

RV acknowledged the financial support for this project from the National Institute of Biomedical Imaging and Bioengineering grant U01 EB023224 and from the Department of Pathology, Anatomy, and Cell Biology, Thomas Jefferson University. AF acknowledged the financial support from the Israel Science Foundation (grant 668/17), the Manitoba Liquor and Lotteries (grant RG-003-21), and the Wolfson Family Chair in Genetics. RV and AF acknowledged the Pilot Funding grant from Thomas Jefferson University and The Hebrew University of Jerusalem collaborative research program.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.747969/full#supplementary-material

Abbreviations

COMPACT, Comparative Matrix of Pattern Counts; DEAB, 4-diethylaminobenzaldehyde; PCA, Principal Component Analysis; PCR, Polymerase Chain Reaction; RA, retinoic acid; WGCNA, Weighted Gene Coexpression Network Analysis.

References