Asymmetric robustness in the feedback control of the retinoic acid network response to environmental disturbances

Retinoic acid (RA) is a developmental signal whose perturbation is teratogenic. We show that early embryos exhibit effective RA signaling robustness to physiological, non-teratogenic, disturbances of this pathway. Transcriptomic analysis of transient physiological RA manipulations during embryogenesis supported the robustness of RA signaling by identifying mainly changes consistent with the progression of embryogenesis and not dramatic treatment-induced changes. Transcriptomic pattern comparisons revealed that RA manipulation led to a network-wide feedback regulation aimed at achieving signaling robustness and normalizing RA levels. A trajectory analysis of target gene and RA network responses identified an asymmetric robustness with a high sensitivity to reduced RA levels, and an activation threshold to increased levels. Furthemore, high robustness to increased RA inversely correlated with a low response to reduced RA. Biological replicates with similar robustness levels mounted responses whose composition likely varies based on genetic polymorphisms to achieve similar outcomes providing insights on the robustness mechanisms.


INTRODUCTION
Retinoic acid (RA) is one of the central regulatory signaling pathways active during vertebrate embryogenesis as well as in adult tissue homeostasis regulating the transcription of numerous downstream target genes (Campo-Paysaa et al., 2008;Clagett-Dame and DeLuca, 2002;le Maire and Bourguet, 2014;Marill et al., 2003;Mark et al., 2009;Metzler and Sandell, 2016). RA is synthesized from vitamin A (retinol) or other retinoids or carotenoids obtained from the diet (Ghyselinck and Duester, 2019;Kedishvili, 2016). During embryogenesis, changes in RA levels, signal timing or signal localization, result in severe developmental malformations arising from both abnormally low and increased RA signaling. Excessive RA signaling induces developmental malformations including brain defects, organ malformations and additional anatomical anomalies (Clagett-Dame and Knutson, 2011;Collins and Mao, 1999;Cunningham and Duester, 2015;Marill et al., 2003;Mark et al., 2009;Shenefelt, 1972).
Syndromes linked to reduced RA signaling include vitamin A deficiency syndrome (VAD), DiGeorge/VeloCardioFacial syndrome (DG/VCF), Fetal Alcohol Spectrum Disorder (FASD), Congenital Heart Disease (CHD), neural tube defects, and multiple types of cancer (Coberly et al., 1996;El Kares et al., 2010;Hartomo et al., 2015;Kim et al., 2005;Kot-Leibovich and Fainsod, 2009;Pangilinan et al., 2014;See et al., 2008;Timoneda et al., 2018;Urbizu et al., 2013). RA levels are tightly regulated at multiple levels throughout life to prevent aberrant gene expression as a result of diet and other environmental changes. Discrete regulatory roles of RA are usually separated temporally and spatially, taking place in different tissues, embryonic regions, or cell types, requiring the fine-tuned regulation of the source, the level, and the gene-regulatory response to this signal. This quantitative, spatial and temporal regulation relies in part on the regulated expression and activity of RA biosynthetic and metabolizing enzymes (Dobbs-McAuliffe et al., 2004;Duester et al., 2003;Hollemann et al., 1998;Sakai et al., 2001). 4 RA biosynthesis involves two sequential oxidation steps: first, mainly alcohol dehydrogenases (ADH) or short-chain dehydrogenase/reductases (SDR) oxidize vitamin A (retinol, ROL) to retinaldehyde (RAL), followed by the retinaldehyde dehydrogenase (RALDH) catalyzed oxidation of RAL to RA (Duester, 2008;Kedishvili, 2016;Parés et al., 2008). RA availability is further affected by ROL, RAL, and RA binding proteins (Kono and Arai, 2015;Napoli, 2017Napoli, , 2016Schroeder et al., 2008). ROL and RAL can also be produced from retinyl ester stores or from ß-carotene from food sources (Blaner, 2019;O'Byrne and Blaner, 2013). In vertebrate gastrula embryos, RA signaling is triggered by the activation of raldh2 (aldh1a2) transcription whose protein product completes the last enzymatic step in RA biosynthesis (Begemann et al., 2001;Chen et al., 2001;Grandel et al., 2002;Niederreither et al., 1999). Then, RALDH2 expression and availability is the earliest rate-limiting step in RA biosynthesis. During RA biosynthesis, substrate availability for the RALDH enzymes, RAL, is controlled by members of the SDR, ADH and AKR families (Adams et al., 2014;Billings et al., 2013;Feng et al., 2010;Porté et al., 2013;Shabtai et al., 2016;. Importantly, expression of many of the enzymes involved in RA biosynthesis is spatially regulated, resulting in a gradient of RA activity peaking in the caudal end of the embryo (Dubey et al., 2018;Niederreither et al., 1997;Schilling et al., 2016). Additional spatial and temporal regulation of this signaling pathway is provided by regulated expression of other components, including retinoic acid receptors (RAR and RXR) and retinoid-binding proteins (Cui et al., 2003;Janesick et al., 2015;Lohnes et al., 1995;Mendelsohn et al., 1994;Xavier-Neto et al., 2015). Besides the maternal nutritional status that can affect the levels of RA signaling in the developing embryo, environmental exposure to chemicals such as alcoholic beverages (ethanol) or other chemicals can affect the biosynthesis of RA or the status of this signaling pathway (Paganelli et al., 2010;. These observations point to the close interaction of RA signaling and the environment and the necessity to adapt the RA signaling 5 network to nutritional changes and insults (e.g., ethanol). This adaptation and maintenance of normal signaling levels under changing conditions is termed robustness (Eldar et al., 2004;Nijhout et al., 2019). Taken together, RA metabolic and gene-regulatory components are under feedback regulation by RA signaling which may provide robustness to the RA signal homeostasis.
A deeper understanding of the RA signaling pathway during embryogenesis is required to elucidate its multiple regulatory roles, and its regulation of the signaling robustness in the presence of environmental disturbances. Very commonly, the RA pathway is studied by increasing the levels of this signal from exogenous sources (Durston et al., 1989;Kessel, 1992;Sive et al., 1990). Alternatively, loss-of-function studies take advantage of RAR inhibitors, inverse agonists, inhibitors of RA biosynthesis, or degradation of the signal (Hollemann et al., 1998;Janesick et al., 2014Janesick et al., , 2013Kot-Leibovich and Fainsod, 2009). In multiple RA loss-offunction studies, the developmental malformations observed are milder than expected suggesting the presence of a compensatory mechanism conferring robustness to perturbations in RA signal (Blumberg et al., 1997;Hollemann et al., 1998;Janesick et al., 2014;Koide et al., 2001;Sharpe and Goldstone, 1997). Paradoxical teratogenic outcomes were observed in a number of RA manipulation studies, suggesting the activity of regulatory feedback mechanisms (D'Aniello et al., 2013;D'Aniello and Waxman, 2015;Lee et al., 2012;Rydeen et al., 2015). To further characterize the robustness of RA signaling during early embryogenesis we employed clearly defined and transient RA manipulations that were terminated during early gastrula, and the kinetics of recovery of the embryos was monitored by RNAseq, qPCR and phenotypic analysis for several hours thereafter. These results demonstrated a high robustness to physiological perturbations of the RA signal with relatively small transcriptome-wide changes. Further transcriptomic analysis showed that components of the RA metabolic and signaling network exhibited expression changes in a manner that is 6 dependent on the direction of RA manipulation, suggesting a mechanistic explanation for the robustness observed. Multiple embryo clutches, i.e. biological repeats, were analyzed.
Comparative analysis of the biological repeats revealed differences between clutches in their individual robustness response to enhanced or suppressed RA signaling. These results exposed potential consequences of the underlying genetic differences between clutches as each one originated from a single female, and most Xenopus laevis stocks in laboratories are outbred.
Our results suggest an asymmetric capacity for robust control of RA signal in the early embryo, likely contributing to the human developmental defects that arise due to imbalance, most often, a reduction) in Vitamin A levels during early development.

Physiological RA manipulation uncovers signaling robustness
To understand the regulatory role of RA signaling in early embryos, we routinely employ two approaches to reduce the levels of this ligand (Fig. 1A). We either partially inhibit the oxidation of retinaldehyde to RA with pharmacological RALDH inhibitors like 4diethylaminobenzaldehyde (DEAB) (Kot- Leibovich and Fainsod, 2009;Shabtai et al., 2016;, or we overexpress the RA hydroxylase, CYP26A1 (Hollemann et al., 1998;Yelin et al., 2005), that renders this signal biologically inactive (Dobbs-McAuliffe et al., 2004;Niederreither et al., 2002). Over a wide range of concentrations these treatments result in clear, but unexpectedly mild, developmental defects for such a central regulatory signaling pathway (compare Fig. 1B with Fig. 1C, D). An efficient approach to induce severe phenotypes by RA signaling loss-of-function (LOF) is to combine knock-down treatments that target the retinoid metabolism at different steps (Fig. 1A, E). These observations demonstrate 7 that RA metabolism and signaling exhibit strong robustness in the context of experimental manipulation of RA levels.
To study the physiological robustness of retinoic acid signaling we focused on moderate perturbations of this pathway within the range determined in X. laevis early embryos (about 100-150 nM all-trans RA) (Durston et al., 1989;Kraft et al., 1995;Kraft and Juchau, 1992;Sive et al., 1990). To empirically determine the RA concentrations to use, we tested the effects on gene expression of concentrations (1 nM -1 µM) spanning the physiological range and above (Fig. 2). Quantitative real-time PCR (qPCR) was performed on RNA samples of treated embryos collected during early (st. 10.25) and late gastrula (st.12) (Nieuwkoop and Faber, 1967). We analyzed genes encoding enzymes involved in RA metabolism and known RA target genes during early embryogenesis. Genes positively regulated by RA at both stages, hoxb1, cyp26a1, and dhrs3 exhibited dose-dependent responses ( Fig. 2A-C). The genes raldh2 and rdh10, are down-regulated by most RA doses (Fig. 2D, F), while raldh3 (aldh1a3) is upregulated earlier on and is repressed later on by the treatment (Fig. 2E). Importantly, RA concentrations as low as 10 nM had a significant effect on gene expression.
To specifically study the robustness of RA signaling in early Xenopus laevis embryos, we established an experimental protocol where RA levels were transiently manipulated pharmacologically within a physiological range. The perturbation was terminated by several washings to remove the treatment and then the embryos were monitored post-treatment during the recovery period. Embryos were treated with RA (10 nM) for 2 hours starting from late blastula (st. 9.5) and washed during early gastrula (st. 10.25), about 2 hours after the treatment was initiated. At different time points during the post-wash recovery period, samples were collected to perform kinetic analysis of the changes in gene expression by qPCR (Fig. 3). To monitor the changes in RA signaling levels by gene expression, two well-characterized RAregulated genes, hoxb1 and hoxb4, were studied. Expression analysis of both genes showed 8 that at the time of treatment washing (t0), both genes were upregulated when compared to control sibling embryos (Fig. 3A). Two hours into the recovery, both RA target genes were back to control expression levels when compared to stage-matched, untreated sibling embryos from the same clutch. To begin understanding the robustness of RA signal and the dynamic regulation of RA metabolism we studied the expression of dhrs3, which encodes an enzyme that preferentially reduces retinaldehyde to retinol to attenuate RA biosynthesis (Feng et al., 2010), cyp26a1, which targets RA for degradation, and raldh2 that produces RA from retinaldehyde (Shabtai et al., 2016). As expected, both negative regulators of RA signaling, dhrs3, and cyp26a1, are upregulated by the increased RA at t0 (Fig. 3B). Their return to normal expression levels is only observed in the 4-hour samples in contrast to the hox genes that returned to normal expression 2 hours earlier. As expected, raldh2 exhibited at t0 an RApromoted downregulation (Fig. 3B). It also took over 4 hours for this gene to return to normal expression levels. These results indicate that the expression of genes encoding RA metabolic enzymes, which themselves are regulated by RA (Fig. 2), is shifted to achieve normal RA levels in the face of external perturbation to RA levels, even though these genes remain abnormally expressed for a longer period.
A complementary study was performed by inhibiting RA biosynthesis taking advantage of DEAB ( Fig. 3C,D). All genes studied exhibited fluctuations during the recovery period. Also, in this case, hoxd1 and hoxb1, key targets of RA signaling, reached almost normal levels at an earlier stage than the genes encoding RA metabolic components (Fig. 3C). As expected, genes encoding anabolic enzymes, e.g., raldh2, were upregulated, and catabolic components, e.g., cyp26a1, were downregulated (Fig. 3D).
The results of the kinetic analysis of the recovery from RA manipulation by qPCR provided a novel and important support and insight into the robustness of RA signaling. We observed that while RA downstream target genes reach normal expression levels relatively 9 quickly, expression of genes encoding RA metabolic enzymes is maintained at seemingly abnormal levels for a longer period of time, such that their combined abnormal activity levels likely result in almost normal RA signaling levels. Such expression changes suggest the hypothesis that during early development RA signaling robustness is achieved, at least in part, via feedback control of network components aimed at preventing a significant differential system-wide gene expression response and the resulting teratogenic outcomes.

Non-teratogenic RA perturbations uncover signaling robustness at the transcriptomic scale
To gain a better understanding of the robustness of RA signaling and adaptation of the metabolic network to disturbances, we performed a kinetic transcriptomic analysis. In an attempt to optimize the kinetic study, parameters such as developmental stages and time between samples were empirically tested. The experimental design involves treating embryos for 2 hours from late blastula (st. 9.5) to early gastrula (st. 10.25) and collecting samples every 1.5 hours after terminating and washing the treatment (Fig. 4A). Embryos were treated with RA (10 nM) or DEAB (50 µM). For each biological repeat, all treatments, controls, and time points were collected from a single fertilization, from the eggs (clutch) of a single female. The efficiency of the treatments and quality of the RNA samples was initially ascertained by qPCR for changes in hox gene expression as a readout of the RA signal levels. Only those biological repeats where both treatments exhibited the expected up-regulation (RA) or down-regulation (DEAB) in hox gene expression were selected for RNA-seq.
We analyzed the time series transcriptomic data set for differential expressed genes using a two-way ANOVA (t=0, 1.5, 3, 4.5 hours; treatments: RA, DEAB, control; n=6 biological replicates). Principal Component Analysis (PCA) of the gene expression variation showed that all samples separated progressively along the first principal component (PC1) 10 which corresponds to the developmental stage (Fig. 4B). Surprisingly, RA manipulated samples clustered with the control samples of the same developmental stage. Normal transcriptomic changes as a result of progression through embryogenesis appear to be the dominant variable distinguishing between the samples irrespective of treatment (Fig. 4B). The second component (PC2) separated the sample groups at intermediate time points (1.5 and 3h) from those at 0 and 4.5h time points, indicating a transient differential expression shift as the next most dominant pattern in the data. The effect of the RA and DEAB treatments on the transcriptome was not readily apparent in the next seven principal components. The the eighth principal component (PC8) showed some separation of DEAB group from RA and Control groups ( Fig. 4C), whereas the tenth principal component (PC10) showed some separation of the RA treatment group from Control and DEAB treatments (Fig. 4D). The top-ranked genes along PC8 and PC10 showed distinct dynamic patterns across the treatments, whereas topranked genes along PC1 and PC2 largely corresponded to in-common dynamic changes over time ( Fig. 4D; Supplemental Fig. S1). The patterns of both opposed RA manipulations, RA and DEAB, closely resembled the control sample for the top-ranked genes along PC1 and PC2, further supporting the normal developmental changes as the dominant pattern. The overall magnitude of induced changes in gene expression in response to the RA or DEAB treatments appears to be less than the normal transcriptome changes occurring during early developmental stages, showing that RA signaling robustness possibly dampens the gene expression changes otherwise induced by abnormal RA levels. Taken together, these results suggest that the RA and DEAB treatments, applied at moderate physiologically relevant concentrations, do not alter the transcriptome extensively in the whole embryo. These results are consistent with the response of a robust system that functions to limit the gene expression changes in the majority of the transcriptome.

11
Dynamic pattern analysis reveals the molecular genetic mechanisms underlying the RA robustness response In order to characterize the molecular transcriptomic mechanism for the robustness response following the RA and DEAB treatments, we analyzed the clutch-averaged data using an unbiased dynamic pattern analysis approach to categorize the gene expression profiles along all possible discretized patterns (Kuttippurathu et al., 2016). Our analysis revealed a total of 4693 significantly differentially expressed genes with greater than 2-fold changes over time (compared to t=0 within each treatment group; multiple testing corrected q-value < 0.05). For each gene, a pattern based on the direction of up-or down-regulation above a 2-fold threshold (three possibilities: up-regulation, down-regulation, or no change) at each of the three recovery time points relative to early gastrula (t=0) was determined. This pattern analysis theoretically should yield 3*3*3 = 27 possible discretized patterns. Such an exhaustive approach allows us to enumerate the dominant as well as subtle patterns in the data and overcomes limitations of conventional cluster analysis that is likely to miss or mask the smaller groups of genes with distinct expression profiles over time. In our analysis, not all possible dynamic patterns had representative genes. Out of the 27 possible dynamic patterns, only 10 patterns were exhibited by genes in the transcriptome in at least one of the three experimental groups (RA, DEAB, Control) (Fig. 5A). Of these, only 5 patterns were exhibited by a substantial number of genes (>100) in at least one of the three experimental groups. These five patterns correspond to upor down-regulation at later time points (3h and 4.5h), and persistently induced or suppressed expression at all three recovery time points, including progressive changes in gene expression over time (Fig. 5A). At the 2-fold threshold, there were no genes that showed down-regulation first and then up-regulation and vice versa. There were no significant differences in the gene counts between RA or DEAB and control groups (two-tailed Z test, p > 0.05). Analysis of statistically enriched pathways and processes (Gene Ontology analysis) in the control group highlighted several functional annotations including development, cell cycle, morphogenesis, RA biosynthesis, gastrulation, and others. (Fig. 5B; Supplemental Table S1). These expression dynamics are consistent with the changes occurring in the early gastrula stage when the endogenous RA system is activated. The patterns with relatively fewer genes (<100) correspond to transient up-or down-regulation at 3h that is normalized by 4.5h, as well as progressive down-regulation.
The results shown in Fig. 5 highlight two aspects of the effect of perturbations in the RA signal. First, at the 2-fold threshold, there are no distinctive dynamic patterns of gene expression changes that occur only in response to RA or DEAB treatment conditions. Second, the number of genes per pattern is of a similar order of magnitude between the two treatments and control. These results suggest that, at physiologically relevant concentrations employed here, RA and DEAB treatments modulate the transcriptional and post-transcriptional regulators at a scale resembling the changes in controls, albeit with differences in potential targets. For instance, only a subset of genes that showed late down-regulation (first row in Fig. 5A) were common across all the three groups (772 genes, Fig. 5C). Several genes showed similar late down-regulation between the treatments or between one of the treatments and the Control group (Fig. 5C). In addition, many genes showed a late down-regulation pattern only in one of the treatment groups or in the Control (Fig. 5C). While the Venn diagram analysis (Fig. 5C) revealed extensive overlap between the RA and DEAB treatments and Control samples within a given differential regulation pattern, the overlap across patterns is not immediately clear from this analysis.
In order to exhaustively compare the RA and DEAB treatment groups for their effect on gene expression relative to the unperturbed temporal pattern observed in the control samples, we adapted a recently developed unbiased approach named COMPACT for analyzing time-series differential transcriptomic profiles across multiple experimental conditions 13 (Kuttippurathu et al., 2016). In this scheme, we started with the statistically significant genes within each treatment group (two-way ANOVA; q<0.05), and constructed discrete patterns based on differential gene expression between RA and control groups (81 theoretically possible patterns, with up-, down-and no-regulation at t=0, 1.5, 3 and 4.5 h). In parallel, we constructed similar discrete patterns for the DEAB versus control groups. Only the patterns with at least one gene in either comparison were included in the subsequent analysis, yielding 32 distinct patterns (Fig. 6A). Finally, we intersected the two distributions to create a matrix of comparative patterns where each element corresponds to a distinct pair of patterns corresponding to RA vs. control and DEAB vs. control (Fig. 6B). The discrete patterns were based on a 1.3-fold average difference between RA (or DEAB) and control groups at each of the four time points. The effect size threshold was chosen at a lower level 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, additional suggestive evidence of robustness.
Analysis of the COMPACT matrix (Fig. 6B, Extended Data 1) demonstrates that the majority of the genes sensitive to the RA signal manipulation exhibited sensitivity to only one direction of the RA perturbation. A total of 193 genes showed a response only to the addition of exogenous RA relative to Control (middle row), whereas 224 genes were only responsive to inhibition of RA production by DEAB (middle column). Of the genes that showed responses to both RA and DEAB perturbations, 88 genes showed up-regulation or down-regulation irrespective of RA or DEAB treatment (quadrants a and d in Fig. 6B). A set of 48 genes showed opposite transcriptional outcomes in response to exogenous RA increase versus inhibition of RA production (quadrants b and c in Fig. 6B). Interestingly, the set of genes that showed upregulation by exogenous RA addition and down-regulation by DEAB (quadrant b) contained several genes involved in the biosynthesis and metabolism of RA (Fig. 6C).
We mapped the expression changes shown in Fig. 6B to the RA signaling and metabolic network (Fig. 6D,E). Notably, genes encoding proteins involved in suppressing RA levels (cyp26a1, dhrs3) were up-regulated and the genes encoding for proteins involved in RA production (aldh1a2, aldh1a3, and rdh10) were down-regulated in response to transient increase in RA. These genes showed opposite regulation in response to exogenous inhibition of RAL to RA oxidation by DEAB (Fig. 6E). As an independent evaluation, we analyzed the differential expression time series data using another unbiased approach, weighted gene correlation network analysis (WGCNA) (Langfelder and Horvath, 2008). Manual examination of WGCNA indicated that one of the correlated gene expression modules (highest correlation with treatment) contained a similar gene set as shown in Fig. 6C, supporting the findings from our COMPACT analysis (green module in WGCNA results, Supplemental Table S2; Extended Data 2). Taken together, our results support a mechanism aimed at maintaining homeostasis of the RA signaling levels to counter the effects of exogenous perturbations, nutritional and environmental, and prevent teratogenic effects by a transcriptional feedback control system.

Clutch-wise heterogeneity demonstrates multiple alternative mechanisms to recover from RA perturbations
The analysis above describes RA as a robust signaling pathway capable of responding to environmental disturbances. Our results show that an important aspect of the regulation of robustness involves changes in RA network components, and identified the main RA network components responding during early gastrulation (Fig. 6). Interestingly, the PCA analysis suggested clutch-specific heterogeneity in gene expression changes where all samples from the same clutch tended to cluster together but the different clutches separated slightly from each other (Fig. 4B,C;Supplemental Fig. S2). For these reasons, we analyzed the transcriptional changes in the RA metabolic network within each clutch. We sought to determine whether the 15 clutch-to-clutch variation observed in the RNAseq data is evident even if we utilize a different technical approach to measure the changes in gene expression in a separate set of manipulated embryos. Hence, we generated an additional set of samples to study the dynamic expression changes in RA metabolic network genes using high-throughput real-time PCR (HT-qPCR; Fluidigm Biomark). This separate cohort of six additional clutches (biological repeats labeled G-L) treated with RA or DEAB following the same experimental design (Fig. 4A) used for RNAseq study of the original six clutches (A-F). We designed PCR primers for multiple members of the RA metabolic network and a number of Hox genes as targets of the RA signal (Supplemental Table S2). Our results show that the heterogeneity between clutches observed in the HT-qPCR data resembles the clutch-to-clutch variation evident in the RNAseq data (Supplemental Figs. S2 and S3).
We ordered the 12 biological repeats (clutches) according to the earliest time at which hoxa1 returned to the baseline levels ( Fig. 7). This ordering intermingled the clutches of the RNAseq and HT-qPCR analysis. Our analysis identified high variability in the response to the RA and DEAB treatments among the individual clutches. Clutches with significant upregulation of hoxa1 expression as a result of RA addition showed limited down-regulation in expression in response to inhibition of RA biosynthesis (DEAB treatment; clutches C, L, J, K, H, I in Fig. 7). The opposite response was also observed where DEAB induced a strong response on hoxa1 expression while the RA treatment induced a mild to very weak response (clutches E, D, F, G, A, B in Fig. 7). Examining the clutch-to-clutch variation in the extent of changes in gene expression revealed a heterogenous correspondence between hoxa1 and RA metabolic network genes. The suppressors of RA signaling, e.g., cyp26a1 and dhrs3, showed significantly altered expression in clutches (C, L, J, K, H) with larger deviations in hoxa1 expression and a strong and extended response to the addition of RA. In contrast, the extent of 16 differential expression of the RA producers (e.g., aldh1a2 and rdh10) after RA manipulation was relatively mild and did not fully align with the deviation in hoxa1 expression (Fig. 7). Such a heterogeneous correlation between extent of feedback regulation and clutch-to-clutch variation in hoxa1 expression was observed across the RA metabolic network ( Supplementary   Fig. S4).

Asymmetry of robustness in response to increase versus decrease in RA signal
We sought to quantitatively rank the robustness of each clutch in an integrated manner based on the extent of shift in multiple hox genes and the complementary feedback regulatory response of the RA metabolic network genes. To this end, we pursued a trajectory-based approach in which the temporal evolution of RA or DEAB treatment groups were compared to their controls based on the regulatory outcome of target gene expression and RA network feedback regulatory gene expression. In this analysis, 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 of RA metabolic network genes. 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 Interestingly, the scattering distribution depended on the direction of the RA manipulation. The range of clutch distribution along the x-axis in Fig. 8F suggests an immediate response to any reduction in RA levels, and an upper limit to this response that points to a constrained ability to mount a corrective feedback regulatory action for large reductions. In contrast, the clutch distribution along y-axis suggests that the response to increased RA requires a minimal change to be activated and it exhibits an upper threshold (Fig. 8F). These observations suggest a high sensitivity to RA reduction but a lower sensitivity to slight RA increase.

An efficiency-efficacy matrix organizes the clutch variability in robustness and feedback regulation
We examined if the asymmetric robustness to direction of RA change in correlates to the extent of feedback regulatory action (Fig. 9). We formulated a "efficiency-efficacy matrix" We examined if the clutches with similar levels of robustness and overall feedback regulation use similar strategies for feedback regulation, i.e., whether they show similarity in the genes that are differentially regulated in response to RA manipulation. Interestingly, we found that the identity of differentially regulated genes in the RA metabolic network can vary between clutches co-localized in the efficiency-efficacy matrix (Figs. 6E and 9C). For example, in response to addition of RA, clutches A, E and F showed similar differential regulation of dhrs3, aldh1a2, and cyp26a1, but diverged in the regulation of other genes such as, stra6, sdr16c5, adhfe1, rdh13, aldh1a3, and cyp26c1 (Fig. 9C). Similarly, in response to DEAB treatment, clutches A, B and D showed downregulation of dhrs3 and cyp26a1, and upregulation of rbp1, suggesting feedback aimed at reducing activity of biochemical processes leading to reduction in RA levels, and increasing import of retinol. However, there was much variability across these clutches in the differential expression of aldh1a2 and aldh1a3, i.e., the feedback control aimed at increasing the production of RA (Figs. 6E and 9D). Taken together, the results described in Figs. 8 and 9 provide strong evidence that the early embryo is differentially capable of responding to increase versus decrease in RA levels, with a distinct range of feedback regulatory responses elicited between the two scenarios.

Robustness of the retinoic acid signaling pathway
Signaling pathway robustness is a central characteristic of all regulatory networks to ensure signaling consistency and reliability. Variations in signaling levels can arise as a result of changing environmental conditions or genetic polymorphisms that can affect the expression or 20 activity level of network components. We describe a systems biology approach to study the robustness of RA metabolism and signaling based on transient manipulation of ligand levels.
Gastrula stage Xenopus laevis embryos have been reported to contain RA levels, around the 100 nM -150 nM range (Durston et al., 1989;Kraft et al., 1995;Kraft and Juchau, 1992). In contrast, due to the clear and severe developmental defects induced, Xenopus embryos, and many other systems, are commonly treated with RA concentrations in the 1 µM to 10 µM range to perturb this signaling pathway (Sive et al., 1990;Taira et al., 1994). To focus our robustness study of RA signaling to the physiological range, we showed that by increasing the RA levels by as little as about 10% (10 nM) we could consistently induce expression changes in RAregulated genes. Interestingly, embryo treatments with RA concentrations in the physiological range exhibit very slight developmental malformations suggesting the induction of compensatory mechanisms to control morphogen signaling and prevent abnormal gene expression (Hollemann et al., 1998;Reijntjes et al., 2005;Sive et al., 1990). We obtained similar results when we reduced the levels of RA by either blocking the biosynthesis (DEAB treatment) or by targeting this ligand for degradation (CYP26A1 overexpression). These observations are supported by multiple loss-of-function studies describing mild developmental malformations induced by RA signaling reduction (Blumberg et al., 1997;Hollemann et al., 1998;Janesick et al., 2014;Koide et al., 2001;Kot-Leibovich and Fainsod, 2009;Sharpe and Goldstone, 1997;Shukrun et al., 2019). Therefore, increased or decreased RA levels in the physiological range result in mild developmental defects, suggesting the activation of compensatory mechanisms. We show that one approach to overcome the robustness of RA signaling is to interfere with the metabolic/signaling network at two different steps as in the case of DEAB together with CYP26A1. A possible explanation for this outcome is that manipulating the RA network at more than one point hampers its ability to efficiently elicit a feedback regulatory response.

21
RA signaling robustness was further analyzed taking advantage of transient RA manipulations and the temporal kinetic, transcriptome-wide (RNAseq) analysis of the restoration of normal gene expression patterns. The large expression differences observed corresponded to normal transcriptome changes resulting from the progression through embryogenesis. The close clustering of the RA-manipulated (increased or decreased RA) and control samples at each time point further supports the robustness of this signaling pathway.
Alternatively, the treatments were inefficient, a possibility we could rule out by qPCR screening of all RNA samples for Hox expression changes prior to sequencing and subsequent computational analysis of the RNAseq data for individual gene expression shifts. Therefore, the PCA analysis suggests the activation of an efficient robustness response of the RA network to maintain normal, non-teratogenic, target gene expression levels during early gastrulation.

The RA robustness response emerges from autoregulatory changes in the RA metabolic network
Our results clearly support the robustness of RA signaling when manipulated within the physiological range. Then, how is this robustness achieved? What is the mechanism activated to achieve RA robustness? Is the RA metabolic network modified to bring about this robustness? The qPCR analysis already provided insights on the RA robustness mechanism.
RA target genes affected by the treatment revealed abnormal expression levels at the end of the manipulation (t=0), and for some target genes (hox), we observed a speedy return to normal expression levels already after 1.5 hours from the end of the treatment. Genes encoding components of the RA metabolic network also exhibited abnormal expression levels at t=0.
Interestingly, their return to normal transcript levels was delayed beyond the time required for the RA targets to reach normal expression. These observations suggest a scenario where the RA metabolic network through an RA-dependent feedback regulatory mechanism is altered in 22 an attempt to restore normal signaling and target gene expression levels. Multiple reports have described the regulation of select individual RA network component expression by RA. Most of these studies deal with RA treatments resulting in the up-regulation of enzymes involved in the suppression or reduction of RA signaling like CYP26A1, ADHFe1, and DHRS3, and downregulation of RA producers (anabolic enzymes) like RALDH2 and RDH10 (Chen et al., 2001;Dobbs-McAuliffe et al., 2004;Fujii et al., 1997;Hollemann et al., 1998;Kam et al., 2013;Sandell et al., 2012;Shabtai et al., 2017;Sonneveld et al., 1998;Strate et al., 2009).
To understand the full extent of kinetic response to transient RA manipulation, we pursued a transcriptome-wide analysis of the recovery kinetics. Out of the 27 possible kinetic patterns consisting of upregulation, downregulation or no change following treatment and washing, only 10 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 among the RA-regulated genes, direct or indirect, there is a limited number of possible regulatory outcomes either in an attempt to normalize RA levels, i.e. robustness, or as target genes. Although this classification is qualitative and there might be differences in intensity, it suggests a limited repertoire of RA regulatory responses. To further understand the regulation of RA targets we performed a comparative analysis of the patterns observed during increased and decreased RA levels using our unbiased approach, COMPACT (Kuttippurathu et al., 2016). This analysis revealed that most genes (75.40%) responding to RA manipulation and recovery exhibited a response to either increased or decreased RA levels. Only 48 genes out of 553 (8.67%) exhibited reciprocal responses to both, increased and decreased RA levels.
Interestingly, many of these genes are involved in RA metabolism or signaling, in agreement with their function in maintaining non-teratogenic RA levels.
Our results provide an integrative network-wide view that incorporates the concerted feedback regulation of multiple pathway components to achieve normal signaling under 23 changing environmental conditions. The temporal discrepancy in the return to normalcy between targets and network components and the network-wide changes in expression patterns suggests a response taking place at multiple levels. The initial response to changes in RA levels will most probably be mediated by the actual enzymes and factors already present in the cell.
In parallel, the same change in RA signaling will initiate a transcriptional response which in the case of increased RA will involve not only the up-regulation and eventual enhancement of the enzymatic RA suppressor activity, i.e. DHRS3, ADHFe1, and CYP26A1, but importantly, a complementary down-regulation and activity reduction in RA producing enzymes (RALDH2 and RDH10). This transcriptional response will necessarily exhibit a slight delay until fully functional as the network component genes will have to undergo transcription, followed by translation and post-translational modifications, when required.
We observed that the expression of some RA network components exhibits oscillatory behavior close to control expression levels, probably as a result of the fine-tuning of the RA signal. This 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 (D'Aniello et al., 2013;Lee et al., 2012;Rydeen et al., 2015). These observations identify a very dynamic feedback regulatory network continuously fine-tuning itself in response to perturbations.

Asymmetric response to increased and decreased RA levels
Our study provided new insights on the network responses to increased and decreased RA levels. One important outcome from our study is the observation that the responses to increased and decreased RA are not inverse but rather governed by different regulatory rules. These observations were reached by performing a trajectory-based analysis in an attempt to rank the 24 RA robustness based on the response of hox genes as RA targets, and the extent of transcriptomic changes of the components of the RA network. We derive the following primary conclusions regarding this asymmetry of RA signal robustness: First, the response to reduced RA signaling is activated after very slight reduction, while the response to increased RA is only activated above a threshold. Also, the response to reduced RA reaches a lower upper threshold than the response to increased RA. The alignment of the hox responses along a diagonal suggested the establishment of an RA gradient among the clutches. In contrast, the network responses showed no clear trend across clutches. Notably, the network responses show different thresholds depending on the direction of the RA manipulation. These observations suggest a high sensitivity to RA reduction but a lower sensitivity to slight RA increase, i.e., an asymmetric capacity to mount responses dependent on the direction of RA changes.
Second, very efficient responses to reduced RA are accompanied by weak responses to increased RA, within the same clutch. The inverse situation was also observed suggesting an asymmetric response to RA manipulations. The alignment of the robustness responses, hox changes, fitted a diagonal with a negative slope. 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. The inverse situation, and also more "balanced" clutches were observed. These observations suggest that not only sensitivity response thresholds are involved, but also the robustness capacity to increased and decreased RA are interconnected in an inverse fashion. These observations led us to formulate an efficiency-efficacy matrix to link the robustness response of the RA network to the outcome of their activity as detected by the hox response.
Third, different clutches with similar robustness levels to increased or reduced RA levels can mount robustness responses by incorporating different components of the RA network. The 25 transcriptomic and HT-qPCR analyses allowed us to perform a detailed determination of the network components comprising the robustness response, i.e. the mechanism of the robustness.
Clutches in proximity in the trajectory-based analysis or in the efficiency-efficacy matrix mount responses by changing different sets of genes. The RA network includes multiple components, i.e. enzymes, whose biochemical activity overlaps or is very similar. The RA producing genes, aldh1a1, aldh1a2, and aldh1a3, also known as raldh1, raldh2, and raldh3 respectively are one example. All three enzymes oxidize retinaldehyde to produce RA (Cunningham and Duester, 2015;Kedishvili, 2016Kedishvili, , 2013Shabtai et al., 2016). Some of the differences between them apparently include enzymatic efficiencies and their expression patterns including location, timing, and intensity (Blentic et al., 2003;Chen et al., 2001;Lupo et al., 2005;Romand et al., 2004;. A similar situation can be argued for enzymes with RA degrading or biosynthetic suppressing activity like CYP26A1, B1, and C1, or DHRS3, ADHFe1, RDH13 and others (Belyaeva et al., 2017(Belyaeva et al., , 2008Hollemann et al., 1998;Shabtai et al., 2017;Sonneveld et al., 1998). Some of these enzymes can reduce or prevent the production of retinaldehyde while others make the RA biologically inactive targeting it for degradation. Irrespective of the biochemical mechanism, these enzymes function to prevent excessive RA signaling. These observations suggest that in the RA network where redundant enzymes and factors performing similar activities are common, a robustness response can be established using different components to reach the same level of robustness.
Fourth, genetic polymorphisms probably explain in part the different responses between clutches to increased or decreased RA levels, and the metabolic and regulatory composition of such responses. In our experiments, we utilized commercially available Xenopus laevis laboratory stocks which are outbred and exhibit extensive genetic variability (Savova et al., 2017). Therefore, we propose that one of the main differences between clutches in their response to RA manipulation is probably the result of genetic variability. Formally, technical 26 issues could somehow account for some of the clutch-to-clutch variability uncovered in the PCA analysis. To further support the genetic basis of the RA response, we analyzed six additional clutchesby high-throughput real-time PCR. This qPCR analysis showed that the new clutches exhibit overlapping responses to RA manipulation to the clutches analyzed by RNAseq. These observations further supported the involvement of genetic variability. We mined the Savova et al. (2017) data identifying multiple polymorphisms in RA network components (Supplemental Table S4). All these observations support the potential contribution of genetic variability in the differential robustness, intensity and response composition, observed between embryo clutches. For several decades, retinoic acid signaling has been the focus of intensive study due to the severe teratogenic effects of abnormal levels and its involvement in the regulation of numerous embryonic processes, oncogenes, and other signaling pathways. While the exposure of embryos or cells to RA induces dramatic developmental changes and malformations, treatment of the same embryos with retinol, the RA precursor, requires much higher concentrations (~100X) to induce similar defects (Durston et al., 1989). Several models could 27 explain the different teratogenic potential of the precursor (retinol) or its final product (RA). In both instances, members of the CYP26 family of cytochrome P450 enzymes should partially neutralize the RA added or produced (Dobbs-McAuliffe et al., 2004;Duester et al., 2003;Hollemann et al., 1998;Sakai et al., 2001). The main difference then becomes the fact that retinol has to go through RA biosynthesis while RA is already the final active signaling ligand (Duester, 2008;Kedishvili, 2016;Parés et al., 2008). The oxidation of retinol to retinaldehyde is a reversible reaction which can reduce substrate availability for the ALDH enzymes. On the other hand, treatment with RA can only be countered through inactivation by CYP26 enzymes. Therefore, the reduced teratogenic efficacy of the retinol treatment could be the result of reduced conversion to RA as part of the feedback regulation of this network, i.e. robustness, while RA treatment is very restricted in its robustness response.

Retinoic acid signaling changes due to environmental changes and disease risk
We obtained evidence of the RA signaling robustness employing several experimental approaches. In all instances, we observed a fast response to the RA manipulation such that at the transcriptome level, the treated samples were not significantly different from controls. All our experiments were performed by either partially inhibiting the endogenous levels of RA, or by slightly increasing (about 10% increase) the physiological content of RA in the embryo.
Under these conditions, the RA robustness of the embryo efficiently regulates and normalizes the transcriptome as a whole via feedback. Based on the comparative analysis of the 12 clutches (genetic backgrounds), we can 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 polymorphisms affecting enzymatic activity or gene expression parameters. In support, 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 showed that the network response, i.e. the genes actually up-or down-regulated, is also dependent on genetic variability like promoter polymorphisms. Then, the 28 developmental malformations arising from environmental insults on RA signaling largely depend on genetic polymorphisms which will determine the efficiency and threshold of the response and the actual network components comprising such a response.

Embryo culture
Xenopus laevis frogs were purchased from Xenopus I or NASCO (Dexter, MI or Fort Atkinson, WI). 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 (Nieuwkoop and Faber, 1967).

Total RNA purification from embryos and cDNA preparation
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 iScriptTM Reverse Transcription Supermix for RT-qPCR kit (according to the manufacturer's instructions).

Expression analysis
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 (Supplemental Table S2).

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 31535 genes. Raw counts were further filtered for non-zero variance across all samples resulting in 31440 scaffold IDs.

High throughput qPCR
cDNA samples were directly processed for reverse transcriptase reaction using SuperScript VILO Master Mix (Thermo Fisher Scientific, Waltham, MA), followed by realtime 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 (Supplemental Table S2). The Software (Fluidigm). Samples were run in triplicate for the 24 genes. The primers in the first pre-amplification group selectively bind for the long isoforms of the genes. In the second preamplification group, primers were selected for the short isoforms only. Third pre-amplification group binds to all the 24 genes. In order to remove the technical variability caused due to long and short isoform of the genes, 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 genes count 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 genes 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 COMBAT method in R using a nonparametric adjustment (Johnson et al., 2007). Following batch correction, the gene list was filtered for a minimum expression threshold to remove genes with normalized count 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. 31 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.5h). 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 PC, expression of 100 top positively-loaded and 100 top negatively-loaded genes was visualized using a heat map.

Dynamic pattern analysis and COMPACT 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.5h was normalized by subtracting the corresponding 't=0h' 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 downregulation)]. 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).

RA network map and visualization
A schematic representation for the position and functioning of the genes involved in the RA biosynthesis, metabolism, translocation and transcription was formulated from the literature. Gene correlation and clustering analysis -Gene expression data was analysed 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 33 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 additional normalized qPCR 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 timepoints). We further selected two subsets of genes: 1) "RA metabolism genes" where, !. ℎ$% !" #$(%) 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 along measures 'HOX shift' (!. ℎ$% !" %+,-%.,*% ) and 'RA Network shift' (!. -01()$23!-!" %+,-%.,*% ). A higher 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.

Availability of supporting data/additional files
The raw and normalized datasets for the RNAseq and HT-qPCR data are available online as Gene Expression Omnibus datasets via SuperSeries GSE154408 containing RNAseq data: GSE154399 and HT-qPCR data: GSE154407.

ETHICS DECLARATIONS
The authors declare no competing interests.      Table S1. (C) Venn diagrams to compare the overlap between the two treatments and control, illustrated for three differential gene expression patterns. The number of genes that show a similar differential expression pattern across two treatments and control is indicated in the middle of the Venn diagrams. The number of genes that showed similar differential expression patterns only in one or two experimental groups are also indicated in the corresponding overlapping regions of the Venn diagrams. Black, control; Orange, DEAB; Blue, RA.