Enhanced Loss of Retinoic Acid Network Genes in Xenopus laevis Achieves a Tighter Signal Regulation

Retinoic acid (RA) is a major regulatory signal during embryogenesis produced from vitamin A (retinol) by an extensive, autoregulating metabolic and signaling network to prevent fluctuations that result in developmental malformations. Xenopus laevis is an allotetraploid hybrid frog species whose genome includes L (long) and S (short) chromosomes from the originating species. Evolutionarily, the X. laevis subgenomes have been losing either L or S homoeologs in about 43% of genes to generate singletons. In the RA network, out of the 47 genes, about 47% have lost one of the homoeologs, like the genome average. Interestingly, RA metabolism genes from storage (retinyl esters) to retinaldehyde production exhibit enhanced gene loss with 75% singletons out of 28 genes. The effect of this gene loss on RA signaling autoregulation was studied. Employing transient RA manipulations, homoeolog gene pairs were identified in which one homoeolog exhibits enhanced responses or looser regulation than the other, while in other pairs both homoeologs exhibit similar RA responses. CRISPR/Cas9 targeting of individual homoeologs to reduce their activity supports the hypothesis where the RA metabolic network gene loss results in tighter network regulation and more efficient RA robustness responses to overcome complex regulation conditions.

For over half a century, gene duplication has been studied as one of the processes important for the evolution of species [29][30][31]. Following duplication of a gene or gene family, one of the copies can be lost, or one member retains the original function and regulation whereas the extra copy can gain novel functions, called neofunctionalization and subfunctionalization. Generally, duplications occur on a small scale involving restricted genomic regions, but in extreme cases, gene duplication of all the genes present in the genome according to Nieuwkoop and Faber [43]. All-trans retinoic acid and Dimethyl sulfoxide (DMSO) were purchased from Sigma-Aldrich (St. Louis, MO, USA). Stock solutions of RA were prepared in DMSO. For transient RA treatment, embryos were placed in 10 or 25 nM RA from late blastula (st. 9.5) and washed two hours later, at early gastrula (st. 10.25) by three changes of 0.1% MBSH and further incubated in fresh 0.1% MBSH for the desired time. Samples were collected 1, 1.5, 2, and 2.5 h after washing.

Quantitative Reverse Transcription Real-Time PCR (qPCR)
Total RNA from embryos was extracted with Aurum™ Total RNA Mini Kit (Bio-Rad, Hercules, CA, USA). cDNA was synthesized using iScript cDNA Synthesis Kit (Bio-Rad). The real-time PCR reactions were performed using the CFX384 Real Time System (Bio-Rad) and iTaq Universal SYBR Green Supermix (Bio-Rad). Each experiment was repeated at least three independent times and each time the samples were run in triplicate. slc35b1.L was used as the housekeeping reference gene. The primers used for qPCR analysis are listed in Table 1.

Generation of CRISPant Embryos
For gene-specific single guide RNA design (sgRNA), genomic DNA sequences were selected from Xenbase.org [44] for the L and S homoeologs when present and analyzed using CRISPRdirect [45] and CRISPRscan [46] for target site search. Computational estimation of the sgRNA efficiency was determined using the inDelphi software [47,48]. For the generation of F0 CRISPant embryos, we injected one-cell stage embryos with Cas9 ribonucleoprotein (RNP) complexes employing the two-RNA component (crRNA:tracrRNA) approach [49]. Briefly, chemically synthesized and modified for stability (Alt-R) RNAs (crRNA and tracrRNA; IDT, Coralville, IA, USA) ( Table 1) were annealed to generate the double guide complexes (crRNA:tracrRNA) and were incubated (10 min at 37 • C) with S. pyogenes Cas9 protein (IDT) to generate RNP complexes. Eight nanoliters of the RNP complex solution were injected into the cytoplasm of one-cell stage embryos.
To determine the efficiency of indel induction, genomic DNA was extracted from 5 individual embryos at mid-gastrula (st. 11) or later employing the GenElute Mammalian Genomic DNA Miniprep Kit (Sigma). The genomic region containing the CRISPR/Cas9 targeted region was PCR amplified using a nested PCR approach (Table 1) and the sizeselected and cleaned product was sequenced. Genome editing efficiency was analyzed by decomposition analysis [50] using the Synthego (Menlo Park, CA, USA) ICE algorithm [51].

Statistical Analysis
All statistical comparisons were carried out using the Prism software package (Graph Pad Software Inc., San Diego, CA, USA). Results are given as the mean ± standard error of the mean (SEM). Tests used were the 2-tailed t-test for two-sample comparisons, Dunnett's (ANOVA) multiple comparisons test, or Fisher test. Differences between means were considered significant at a significance level of p ≤ 0.05.

Conservation of the RA Network
In a recent extensive data mining effort searching the KEGG and Xenbase databases [44,52] and the literature, we assembled the putative components of the RA metabolic and signaling network assembling components described in multiple experimental systems and tissues [1,6,9]. To determine the composition of the RA network active during gastrula in Xenopus laevis embryos we assessed which components are expressed during this developmental stage based on analysis of transcriptomic data sets [27,53,54] and corroborated the gastrula expression by quantitative RT-PCR (qPCR) [27]. This view of the RA metabolic and signaling network exhibits a rather uncommon characteristic that for each enzymatic step, multiple genes encoding enzymes have been described that are capable of performing the same reaction ( Figure 1). For some of the reversible reactions, such as the oxidation of ROL to RAL and the corresponding reduction of RAL to ROL, a preferred activity has been identified, but the reverse reaction might be possible under certain conditions [55][56][57]. Because it is necessary to maintain non-teratogenic levels of RA at these stages, there appear to be several ways to control levels of RA signaling. Substrate availability for the retinaldehyde Figure 1. Evolutionary conservation of the RA metabolic and gene-regulatory network genes in Xenopus laevis. Composition of the RA metabolic and gene-regulatory network in the Xenopus laevis genome based on KEGG and Xenbase database analysis and literature searches [44,52]. Expression of the RA network components during gastrula stages was summarized from our own and published transcriptomic datasets. The relative expression shown (blue shades) during gastrula stages is based on Session et al. [32]. Dark blue, ≥10 TPM; middle blue, 0.5-10 TPM; light blue, ≤0.5 TPM; white, no data in the transcriptomic dataset. The homoeolog/singleton status of each gene is marked (L and/or S). The relative expression levels between homoeologs are summarized: =, similar expression levels; <, 3-6 fold difference; <<, more than 6 fold difference. Asterisks indicate whether temporal expression patterns of the homoeologs are similar (no asterisk), partially divergent (*), or highly divergent (**).
Analysis of the allotetraploid status of the genes encoding all the identified RA network components in the Xenopus laevis genome revealed that in 25 out of 47 genes both the L and S homoeologs [32] have been retained during evolution ( Figure 1; Table 2). Thus, 22 genes (46.8%) in the RA network have lost one of the homoeologs, bringing this metabolic and signaling network close to the genomic average (43.6%) of singletons among protein-coding genes [32]. This observation seemingly contradicts previous studies of other main signaling pathways critical for normal embryogenesis, like TGFβ, FGF, Wnt, Hh, Notch, Hippo, and of transcription factors for which the conservation of the L and S homoeologs is very high (>83.3%; Table 2) [40][41][42].  [44,52]. Expression of the RA network components during gastrula stages was summarized from our own and published transcriptomic datasets. The relative expression shown (blue shades) during gastrula stages is based on Session et al. [32]. Dark blue, ≥10 TPM; middle blue, 0.5-10 TPM; light blue, ≤0.5 TPM; white, no data in the transcriptomic dataset. The homoeolog/singleton status of each gene is marked (L and/or S). The relative expression levels between homoeologs are summarized: =, similar expression levels; <, 3-6 fold difference; <<, more than 6 fold difference. Asterisks indicate whether temporal expression patterns of the homoeologs are similar (no asterisk), partially divergent (*), or highly divergent (**).
Analysis of the allotetraploid status of the genes encoding all the identified RA network components in the Xenopus laevis genome revealed that in 25 out of 47 genes both the L and S homoeologs [32] have been retained during evolution ( Figure 1; Table 2). Thus, 22 genes (46.8%) in the RA network have lost one of the homoeologs, bringing this metabolic and signaling network close to the genomic average (43.6%) of singletons among protein-coding genes [32]. This observation seemingly contradicts previous studies of other main signaling pathways critical for normal embryogenesis, like TGFβ, FGF, Wnt, Hh, Notch, Hippo, and of transcription factors for which the conservation of the L and S homoeologs is very high (>83.3%; Table 2) [40][41][42]. Further analysis of the distribution of homoeologs and singletons within the RA network genes, however, revealed a surprising, non-random distribution of gene loss events ( Figure 1). Among the enzymes involved in the metabolic steps leading up to the production of RAL, including retinoid storage, about 75% of the genes (21 out of 28) are encoded as singletons ( Figure 1; Table 2). Interestingly, for genes involved in the oxidation of RAL to RA, hydroxylation of RA, or actual RA-driven gene regulation, almost all (94.8%; 18 out of 19) are still encoded by both L and S homoeologs ( Figure 1; Table 2).
This suggests a preferential loss of homoeologs involved in the production of RAL, regulation of RAL production, or storage of retinoids. One possible explanation for this asymmetrical gene loss in the metabolic side of the RA network could be the observation that RAL availability for oxidation by retinaldehyde dehydrogenases is like a "commitment" step ( Figure 1). The oxidation of RAL to RA by the Aldh1a1, Aldh1a2, and Aldh1a3 enzymes cannot be reversed and either promotes RA-driven gene regulation or the RA produced has to be neutralized and degraded. RA signaling is one of the major embryonic signaling pathways dependent on the maternal nutritional status and its function can be altered by environmental factors [5,[15][16][17]. Fluctuations in RA signaling, increase or decrease, can be extremely teratogenic. Therefore, we explored the possibility that the extensive gene loss observed preferentially achieves tighter regulation of the RA signal, providing an evolutionary advantage.
To assess whether this preferential gene loss in the RA metabolic and gene-regulatory pathway was restricted to RA signaling, we also analyzed the genomic evolution of two additional nuclear receptor signaling pathways closely linked to RA: vitamin D and thyroid hormone signaling [64,65]. Thyroid hormone biosynthesis and signaling in Xenopus includes 23 genes out of 52 that are already singletons (44.2%), bringing this pathway to the genomic average with no obvious distinctive distribution (Table 2). Interestingly, vitamin D biosynthesis and gene regulation exhibit a distribution resembling the RA network where part of the pathway is rich in singletons ( Figure 2). Also in this case, the whole pathway (21 genes) exhibits 38% singletons (8 genes) close to the protein-coding average ( Table 2) [32]. However, from the production of cholecalciferol (vitamin D 3 ) the pathway has 6 genes out of 9 that are already singletons ( Figure 2; Table 2). These 66.7% singletons in the vitamin D-specific part of the network show a preferential loss of genes involved in metabolism and gene regulation by this ligand.
the genomic average with no obvious distinctive distribution (Table 2). Interestingly, vitamin D biosynthesis and gene regulation exhibit a distribution resembling the RA network where part of the pathway is rich in singletons ( Figure 2). Also in this case, the whole pathway (21 genes) exhibits 38% singletons (8 genes) close to the protein-coding average ( Table 2) [32]. However, from the production of cholecalciferol (vitamin D3) the pathway has 6 genes out of 9 that are already singletons ( Figure 2; Table 2). These 66.7% singletons in the vitamin D-specific part of the network show a preferential loss of genes involved in metabolism and gene regulation by this ligand.
As both the RA and vitamin D signaling pathways involve biosynthesis of the regulatory ligand through a metabolic pathway, to assess the generality of this observation we explored the L or S gene loss in additional metabolic pathways. For the de novo purine biosynthesis pathway we scored 7 genes out of 15 that are already singletons ( Table 2). The 47% singleton status is close to the genomic average suggesting the normal rate of gene loss for coding sequences. Analysis of glycolysis + gluconeogenesis + Krebs cycle identified 9 singletons among 44 genes (20.4%), which is a low singleton proportion suggesting conservation of both homoeologs in these metabolic pathways ( Table 2). Analysis of the folic acid metabolic network indicated the reverse: a high proportion of singletons. From the information in KEGG, we identified 27 enzyme-encoding genes in the X. laevis genome out of which 19 (70.4%) are encoded by singletons (Table 2), indicating a preferential loss of one of the homoeologs during evolution. KEGG analysis of the vitamin D metabolic and signaling network identified 21 genes in the X. laevis genome. The homoeolog/singleton status of each gene is marked (L and/or S). In the metabolic part of the pathway leading to cholesterol production (above the red dotted line), the pathway runs in parallel from lanosterol and dehydrolanosterol. The relative expression levels between homoeologs are summarized: =, similar expression levels; <, 3-6 fold difference; <<, more than 6 fold difference. Asterisks indicate whether temporal expression patterns of the homoeologs are similar (no asterisk), partially divergent (*), or highly divergent (**). KEGG analysis of the vitamin D metabolic and signaling network identified 21 genes in the X. laevis genome. The homoeolog/singleton status of each gene is marked (L and/or S). In the metabolic part of the pathway leading to cholesterol production (above the red dotted line), the pathway runs in parallel from lanosterol and dehydrolanosterol. The relative expression levels between homoeologs are summarized: =, similar expression levels; <, 3-6 fold difference; <<, more than 6 fold difference. Asterisks indicate whether temporal expression patterns of the homoeologs are similar (no asterisk), partially divergent (*), or highly divergent (**).
As both the RA and vitamin D signaling pathways involve biosynthesis of the regulatory ligand through a metabolic pathway, to assess the generality of this observation we explored the L or S gene loss in additional metabolic pathways. For the de novo purine biosynthesis pathway we scored 7 genes out of 15 that are already singletons ( Table 2). The 47% singleton status is close to the genomic average suggesting the normal rate of gene loss for coding sequences. Analysis of glycolysis + gluconeogenesis + Krebs cycle identified Cells 2022, 11, 327 8 of 23 9 singletons among 44 genes (20.4%), which is a low singleton proportion suggesting conservation of both homoeologs in these metabolic pathways ( Table 2). Analysis of the folic acid metabolic network indicated the reverse: a high proportion of singletons. From the information in KEGG, we identified 27 enzyme-encoding genes in the X. laevis genome out of which 19 (70.4%) are encoded by singletons (Table 2), indicating a preferential loss of one of the homoeologs during evolution.
Based on the analysis of homoeolog loss in the Xenopus laevis genome to date [32,[40][41][42], signaling pathways and metabolic networks can be classified into three groups. There are pathways that exhibit homoeolog retention rates similar to the protein-coding gene average in the whole genome. From our analysis, we identified the thyroid hormone synthesis and signaling, vitamin D biosynthesis and signaling, de novo purine biosynthesis, and RA metabolism and signaling pathways as belonging to this group ( Table 2). The second group includes previously analyzed signaling pathways described as having very high homoeolog retention rates (low or suppressed singletons) [32,[40][41][42], and we added Krebs cycle, glycolysis, and gluconeogenesis to this group ( Table 2). Our analysis identified the third group as having a high rate of gene loss creating a high proportion of singletons (Table 2). Interestingly, apart from the folate metabolic network, for both RA metabolism and vitamin D biosynthesis and signaling, the high rate of homoeolog loss localizes to a specific region of the pathway (Figures 1 and 2).

Genomic Changes in the Loss of a Homoeolog
The high incidence of singletons in the RA metabolic network from retinyl ester storage to the production of RAL prompted us to try to understand the genomic rearrangements that resulted in this enhanced homoeolog loss that involved 21 out of 28 genes (Table 2). We focused our analysis on enzymes with alcohol dehydrogenase activity to oxidize vitamin A to RAL (Producers 1); enzymes, mainly of the short-chain dehydrogenase/reductase family, that reduce RAL to retinol (Suppressors); and proteins involved in β-carotene cleavage to RAL (Producers 3) (Figure 1; Table 3). We compared the appropriate genomic regions between the L and S chromosomes choosing the first pair of conserved genes flanking the deleted or rearranged region as boundaries (Table 3). Using these flanking genes, we could determine the distance between them in the L and S chromosomes and analyze the region between them (Table 3). This type of analysis revealed cases in which single or multiple genes were deleted. Additionally, the length of the modified region changed from 0.1 to 410 Kb (Table 3). We could group the rearrangements into three groups. The first group represents singletons in which the loss of a homoeolog involved a relatively small (0.1-36 Kb) and simple deletion ( Figure 3A; Table 3). The second group involved large (102-410 Kb) deletions ( Figure 3B; Table 3). In the third group, we identified large deletions (81 and 238 Kb) combined with extensive rearrangement of the genomic region ( Figure 3C; Table 3). These results show that most deletions leading to the loss of a homoeolog are relatively simple although the genomic region lost can be small or very large. In a few cases, the loss of a homoeolog involved complex genomic rearrangements in addition to the deletion of genes. In the locus on chromosome 1 ( Figure 3C) there are multiple adh genes suggesting the possibility that this region contained duplicated sequences that could contribute to the rearrangements in this genomic region. On chromosome 9_10 we also observed a complex deletion but did not observe gene duplications that could contribute to its creation.

Expression Overlaps and Responsiveness of RA Network Components
Our recent analysis of the RA metabolic and signaling network revealed a high degree of robustness following disruption of this pathway within the physiological range [27]. Moreover, this study showed that enzymes performing the same metabolic reaction and expressed in partially overlapping patterns might be regulated differently. These differential responses to RA fluctuation are part of the mechanism to keep this critical signal within an appropriate, non-teratogenic range [27]. One possibility for the preferential loss of homoeologs in genes encoding RA network components is the selective or non-selective reduction of gene copies to achieve tighter regulation of the signal. To begin exploring these possibilities as possible driving forces for gene loss, we searched for gene pairs expressed during gastrula stages that have the same enzymatic activity but one of them is a singleton and the other is still presented as a homoeolog pair. Based on our previous studies we chose rdh10 and sdr16c5 for genes encoding enzymes that oxidize ROL to RAL (Producers 1), dhrs3 and rdh14 for genes encoding enzymes that reduce RAL to ROL (Suppressors), and aldh1a2 and aldh1a3 for genes encoding enzymes that oxidize RAL to RA and for which no singletons are known (Producers 2) (Figure 1). The temporal pattern of expression was determined for these genes including analysis of the individual homoeologs by qPCR (Figure 4). In all three gene groups, we observed extensive overlap in the temporal expression pattern between the selected singleton and at least one of the homoeologs of the paired gene with similar enzymatic activity ( Figure 4A-C). While sdr16c5 (singleton) exhibits a pattern similar to the rdh10.S homoeolog both having significant maternal expression ( Figure 4A), the rdh10.L homoeolog is mainly expressed zygotically. Among the Suppressors tested ( Figure 4B), rdh14 exhibits what might be maternal transcripts, and its expression levels decline during gastrulation ( Figure 4B). Both dhrs3 homoeologs retain extensively overlapping temporal expression patterns that peak at late gastrulation and subsequently decline ( Figure 4B). The retinaldehyde dehydrogenase encoding genes aldh1a2 and aldh1a3 exhibit mostly zygotic expression patterns with a marked upregulation with the onset of gastrulation ( Figure 4C). These expression patterns support the partial overlap between the genes selected and are part of the RA network during gastrulation.
Cells 2022, 10, x 11 of 24 expression ( Figure 4A), the rdh10.L homoeolog is mainly expressed zygotically. Among the Suppressors tested ( Figure 4B), rdh14 exhibits what might be maternal transcripts, and its expression levels decline during gastrulation ( Figure 4B). Both dhrs3 homoeologs retain extensively overlapping temporal expression patterns that peak at late gastrulation and subsequently decline ( Figure 4B). The retinaldehyde dehydrogenase encoding genes aldh1a2 and aldh1a3 exhibit mostly zygotic expression patterns with a marked upregulation with the onset of gastrulation ( Figure 4C). These expression patterns support the partial overlap between the genes selected and are part of the RA network during gastrulation. For the four homoeolog pairs studied, we observed divergence in their temporal expression patterns. In the extreme case, rdh10.S exhibits high levels of maternal transcripts and a gradual decline during gastrula stages, whereas rdh10.L expression is activated after the midblastula transition as a zygotic gene ( Figure 4A). A similar expression pattern for this homoeolog pair was determined from transcriptomic data [32]. This divergence in temporal expression patterns suggests changes in regulatory elements and initial subfunctionalization of the homoeologs [36]. For two of the other homoeolog pairs, dhrs3 and aldh1a3, there are subtle differences in their temporal expression patterns, with extensive overlap but also new gene-specific changes ( Figure 4B,C). In contrast, the aldh1a2 homoeologs exhibit temporal expression patterns that are very similar ( Figure 4C). Interestingly, the early expression of aldh1a2 at the onset of gastrulation has been linked to the onset of RA signaling as the enzyme encoded by this gene is the last component needed to complete the biosynthesis of RA [15,[66][67][68][69]. Additionally, we showed that within the aldh1a gene family, aldh1a2 is expressed at the highest levels during early gastrula stages [15]. Thus, there appears to be selective pressure to conserve this expression pattern and the early gastrula activity of aldh1a2.
To better understand the contribution of these genes to the response of the RA network to fluctuations in ligand levels, we studied the response of these genes to subtle manipulation of RA levels. Analysis of the RA content of Xenopus laevis embryos during early gastrula estimated that they contain about 100-150 nM RA [70][71][72][73][74][75][76]. To perform physiologically relevant manipulations of RA we increased that level by about 10-25% using 10 and 20 nM treatments, respectively [27]. Embryos were treated from late blastula to early gastrula (st. 10.25) and collected for expression analysis by qPCR. This analysis revealed robust responses by the dhrs3 homoeologs (p < 0.0001) and weak (not significant) responses by the rest of the genes tested ( Figure 5). The self-regulation of the RA metabolic and signaling network to maintain or restore normal signaling levels is widely accepted, and the transcriptional responses of aldh1a2, cyp26a1, dhrs3, and rdh10 to increased RA levels are the basis of this suggestion [18][19][20][21][22][23][24]69,77]. In a recent study, we performed a detailed analysis of the RA responsiveness and requirement for the RA network genes expressed during early gastrula [27]. While some genes exhibited robust and concentrationdependent responses, others showed no significant changes in response to RA fluctuations. Additionally, the same gene was shown to exhibit different responsiveness at different developmental stages [27]. These changes in responsiveness could be explained in part by the different temporal expression patterns and the RA responsiveness of the homoeologs that were not addressed in previous studies.  Table 1. Expression changes were normalized to transcript levels in control embryos.

Homoeolog Response to Transient RA Manipulation
The enhanced homoeolog gene loss observed within the RA metabolic network leading to the production of RAL raised the question as to the possible selective pressure driving this phenomenon. Maintenance of normal RA signaling levels is central for the pre-  Table 1. Expression changes were normalized to transcript levels in control embryos.

Homoeolog Response to Transient RA Manipulation
The enhanced homoeolog gene loss observed within the RA metabolic network leading to the production of RAL raised the question as to the possible selective pressure driving this phenomenon. Maintenance of normal RA signaling levels is central for the prevention of the teratogenic effects of increased or decreased RA signaling levels [13,21,70,[78][79][80][81]. Thus, homoeolog gene loss might be a "solution" to achieve tighter signaling regulation, i.e., robustness [27]. Several models can be suggested that could drive this gene loss to achieve higher RA signaling robustness. One possibility is that one of the homoeologs in the pair has a looser regulation or enhanced sensitivity, "noisier" regulation, exhibiting enhanced responses to RA changes. It is important to note that this "noisy" gene could in theory mediate fast responses and might not necessarily be advantageous to lose. Alternatively, coordinated regulation of numerous genes performing the same enzymatic function might be more complicated to achieve, so having fewer genes would provide tighter regulation. Then, the "noisy" gene is lost preferentially to reach tighter regulation. To discriminate between these possibilities, we took advantage of an experimental protocol for transient RA manipulation and kinetic monitoring of the recovery process by qPCR ( Figure 6A) [27]. This assay allows us to monitor the robustness response as it takes place by analyzing the expression changes of RA network components and downstream, RA-regulated genes. To perform physiologically relevant RA manipulations, based on our homoeolog responsiveness analysis ( Figure 5) and our previous studies [27], embryos were treated with 10 and 25 nM all-trans RA for 2 h from late blastula to early gastrula (st. 10.25). The treatment was terminated by washing (T0) and samples were collected during the recovery period at 1.0, 1.5, 2.0, and 2.5 h post-washing (T1, T1.5, T2, and T2.5, respectively) ( Figure 6A). RNA samples were prepared for comparative expression analysis to control samples. For enzymes that oxidize ROL we analyzed rdh10.L, rdh10.S, and sdr16c5 (Producers 1), for enzymes that reduce RAL to ROL we chose dhrs3.L, dhrs3.S, and rdh14 (Suppressors), and among the genes that produce RA (Producers 2), we studied both homoeologs of aldh1a2 and aldh1a3 (Figure 1). For each homoeolog or singleton, the relative expression (fold change; ∆) was calculated at each time point relative to the expression in sibling control embryos at the same developmental stage, and the average ∆ of four biological replicates was calculated. Analysis of rdh10.L, rdh10.S, and sdr16c5 revealed very slight fluctuations of all three genes irrespective of whether 10 or 25 nM RA was used for the treatment ( Figure 6B,C). These weak responses are in agreement with the previous results of the RA responsiveness in which all three genes responded similarly ( Figure 5). Additionally, these results suggest that none of the genes exhibit heightened responses in our experimental protocol which aims to mimic physiological RA fluctuations.
Analysis of the genes encoding enzymes preferentially involved in reducing RAL to ROL, dhrs3.L, dhrs3.S, and rdh14, revealed the hypothesized situation where one of the homoeologs exhibits an enhanced response to changes in RA levels and delayed restoration of normal expression levels. Expression of dhrs3.L shows the strongest upregulation at T0 irrespective of the amount of RA employed of the three genes analyzed ( Figure 6D,E). By comparison, dhrs3.S shows a robust but weaker response at the end of the treatment, and rdh14 only exhibits a weak response. Importantly, while dhrs3.S and rdh14 reached almost normal expression levels (<2.3 fold) at the end of the recovery period (T2.5), dhrs3.L is still significantly upregulated (>2.9 fold) ( Figure 6D,E). Analysis of the aldh1a2 and aldh1a3 homoeologs revealed that by the end of the transient treatment (T0), some of the homoeologs exhibit a clear upregulation, but already one hour into the recovery period all genes are almost back to normal expression levels ( Figure 6F,G). These results show that the transient RA treatments can induce robust responses (dhrs3), but many of the genes studied exhibit mild expression changes and strong robustness responses, i.e., return to normal levels. A previous study that employed the same transient RA manipulation protocol but collected samples up to 5.5 h after washing showed efficient recovery to normal expression levels of most RA network genes [27].   To make comparisons between samples, genes, or biological replicates easier, we calculated the fold change of each gene for all time points (∆) and we summed up the values into a summary fold change score (∑∆). This tightness regulation score should be low for tightly regulated genes and high for enhanced responders with slow recovery to normal values. We calculated this regulation tightness score for all genes studied ( Table 4). The results show that for genes exhibiting moderate or limited gene responsiveness to RA and expression changes throughout the recovery period, rdh10.L, rdh10.S, sdr16c5, rdh14, aldh1a2.L, and aldh1a2.S, the average ∑∆ score ranged from 4.3 to 5.8 (Table 4). For each one of these genes, the difference in the amount of RA added had a very limited effect on the variation in their expression. For genes with higher fluctuation in their expression, aldh1a3.L, aldh1a3.S, dhrs3.L, and dhrs3.S, the ∑∆ score increased in correlation with the kinetic analysis result (Table 4; Figure 6). The score for dhrs3.L reached ∑∆ > 28 making it the least tightly regulated gene or the gene with the most extreme responses to RA fluctuations of those analyzed. Importantly, both singletons studied, rdh14 and sdr16c5 exhibited tight regulation with low responsiveness to RA fluctuations even though they have been shown to be involved in vivo in the metabolism of RA [82][83][84].

RA Responsiveness in Homoeolog CRISPant Embryos
The results of the individual homoeolog responses to transient RA manipulation identified gene pairs that represent all the possibilities initially suggested. The rdh10 and aldh1a2 genes have tightly and similarly RA-regulated homoeolog pairs. Similar but slightly enhanced responses were observed for the aldh1a3 homoeologs, whereas the dhrs3 homoeolog pair showed strong responses to RA changes and marked differences between the L and S genes. To begin to address the possible force driving the gene loss that gives rise to singletons, we took advantage of the CRISPR/Cas9 technology to create a partial, homoeolog-specific gene loss. We designed homoeolog-specific single guide RNAs (sgRNA) for the dhrs3 and rdh10 genes to knock down the expression of one homoeolog without affecting the second one. We also designed a sgRNA targeting the sdr16c5 singleton. To determine the efficiency of the sgRNAs, DNA was extracted from CRISPant embryos and the genomic region containing the sgRNA targeted sequence was PCR-amplified and sequenced. Decomposition analysis of the sequence traces [50,51] provided a quantitative assessment of the genome editing efficiency. Analysis of the sequencing traces demonstrated the creation of indels around the sequence targeted by the sgRNA and the deterioration of the sequencing quality (Supplementary Figure S1). According to the decomposition analysis, we observed a relatively robust effect of the sgRNAs inducing indels (Supplementary Figure S1).
To study the effect of losing one of the homoeologs on the RA robustness response, embryos were injected at the one-cell stage with the appropriate sgRNA/Cas9 riboprotein complex to generate CRISPant embryos, which were then subjected to the transient RA manipulation protocol using low RA concentrations (10 nM) ( Figure 6A). RNA samples were collected at T0 (wash) and at 1.0, 1.5, 2.0, and 2.5 h after treatment termination. To understand the effect of this gene knockdown on the RA signaling levels, we analyzed the expression of a panel of RA-regulated genes including cyp26a1.L, hoxd1.L/S, hoxa1.L, hoxa1.S, hoxa2.L/S, hoxb4.S, hoxb1.L, and hoxb1.S (Figure 7 and Supplemental Figure S2). Comparison at T0 of the expression levels of RA-responsive genes between RA treated CRISPants (rdh10, dhrs3, and sdr16c5) and control siblings treated with RA supported the efficiency of the sgRNAs ( Figure 7A,B and Supplemental Figure S2A,B). The RA treatment alone induced upregulation of all RA targets ranging from 3.6 to 35-fold increase, while RA treatment of the rdh10.L, rdh10.S, and sdr16c5 CRISPants resulted in a weaker RA-induced upregulation irrespective of the gene being knocked down ( Figure 7A and Supplemental Figure S2A). In agreement with the similar and limited responses to RA exposure ( Figure 6B,C), the three genes individually targeted, rdh10.L, rdh10.S, and sdr16c5, resulted in similar outcomes. This response of the RA target genes agrees with the suggestion that between similarly regulated RA network genes that encode enzymes performing the same metabolic reaction, their loss is equivalent in the early embryo. Maintenance of the singletons might reflect different spatial-temporal regulation to perform the same enzymatic reaction in different tissues in the embryo or the adult.
Kinetic analysis of the RA robustness response in the rdh10.L, rdh10.S, and sdr16c5 CRISPants was monitored by following the expression of the RA target genes during the recovery period (T0-T2.5; Figure 6A). To better understand the contribution of the RA network components studied, the CRISPant samples treated with RA were compared to siblings treated with RA only. In most instances, knockdown of each of these three genes had a mild effect on gene expression, reducing the response to the transient RA treatment (Figure 7C,E,G,I and Supplemental Figure S2C,E,G,I). It is important to note that weaker responses in the RA-treated CRISPants support a tighter regulation of the signal as a result of the gene loss phenocopy. In a few instances we observed enhanced responses to the RA exposure mainly linked to the rdh10.L CRISPant ( Figure 7E and Supplemental Figure S2C), the rest of the samples exhibited more restricted responses to RA exposure supporting a tighter regulation. To simplify the comparative analysis between CRISPants, we calculated the regulation tightness score (∑∆) of the RA target genes for all time points compared to their response in RA-treated embryos ( Figure 8A). The scores for the RA-regulated genes in the three RA-treated CRISPants showed that all of them reduced the target gene expression changes. This result suggests that in the case of rdh10.L, rdh10.S, and sdr16c5, the gene activity reduction results in tighter regulation of the RA robustness response in agreement with a gene loss model for better regulation of the signal. The low significance of the changes further supports the efficient robustness response of the RA network. In this analysis, we can also observe the enhanced responses linked to the rdh10.L CRISPant ( Figure 8A).
The dhrs3 homoeologs exhibit enhanced responses to transient RA exposure, with the dhrs3.L gene showing the strongest responses ( Figure 6D,E; Table 4). Since the Dhrs3 enzyme preferentially reduces RAL back to ROL [25,58,85], the dhrs3 CRISPants should exhibit enhanced RA signaling unless the RA network self-regulation and robustness response compensates for this loss of activity [27]. Supporting the robustness scenario, the dhrs3 CRISPants alone had almost no effect on the RA responsive genes with the exception of the two hoxb1 homoeologs exhibiting a 1.5-7.5 increase in expression compared to control samples (not shown). Analysis at T0 of both dhrs3 RA-treated CRISPants showed that these responses were not enhanced as expected, and the partial knockdown of one of the dhrs3 homoeologs resulted in reduced responses. In agreement with the loss of the loosely regulated, noisy homoeolog, these results show that the dhrs3.L CRISPants exhibit a stronger reduction in the RA response compared to knockdown of the dhrs3.S homoeolog ( Figure 7B and Supplemental Figure S2B). Then, knockdown of the loosely regulated homoeolog achieves tighter regulation of the response.

Conclusions
The RA metabolic and signaling network is strongly dependent on the nutritional status and is influenced by the environment. Fluctuations in the RA signaling levels during embryogenesis can result in severe teratogenic outcomes. The preferential gene loss of the RA network components involved in the metabolism leading to RAL production suggests a selective pressure to achieve tighter regulation of the robustness response. Eliciting a robustness response by transient RA manipulation together with knockdown of specific gene homoeologs support the suggestion that gene loss might be linked to more efficient regulation of the network. The RA robustness response efficiently overcomes the reduction of one of the homoeologs. Tighter network regulation might involve loss of homoeologs similarly regulated, or homoeologs with enhanced responses. While the allotetraploid condition of X. laevis is convenient to explore these genomic changes and their regulatory outcomes, in diploid organisms, besides gene duplications and deletions, mutation, addition, and deletion of regulatory elements might take place to achieve the similar outcomes. While the T0 analysis supports the loss of the "noisier" gene to achieve tighter regulation of the response during the RA treatment, analysis of the full recovery kinetics compared to RA-only manipulated embryos provides information as to the effects of the homoeolog knockdown on the RA robustness response. Analysis of the same panel of RA-responsive genes showed that by about 1.5 h into the recovery period (T1.5), the expression levels of the target genes analyzed in the dhrs3.L and dhrs3.S CRISPants was almost back to the same as the samples treated only with RA ( Figure 7D,F,H,J and Supplemental Figure S2D,F,H,J). We could observe slight fluctuations in expression levels but in most instances, both CRISPants gave similar variations although the responses in the dhrs3.L tend to be lower than the RA-only samples, while the dhrs3.S CRISPants gave slightly enhanced responses. For multiple genes, at T0 we observed the upregulation characteristic of the treatment before RA washing ( Figure 7F,H,J and Supplemental Figure S2D,J). Additionally, at T0 and T1, CRISPants of the more tightly regulated homoeolog, dhrs3.S, exhibit larger fluctuations in expression of the RA responsive genes. Calculation of the regulation tightness score showed the opposed outcomes of the homoeolog-specific knockdown ( Figure 8B). The RA robustness response is enhanced by knockdown of the dhrs3.S homoeolog, while knock-down of dhrs3.L results in a reduced response. In addition, in this case most changes observed were hardly significant compared to control RA-treated embryos. Analysis of the dhrs3 homoeologs identified the dhrs3.L gene as the one exhibiting enhanced responses to RA treatment and less tight regulation ( Figure 6D,E). Then, knockdown of the homoeolog exhibiting tighter regulation exposes the system to the homoeolog with the gene with the apparent looser regulation. Interestingly, within 1.5 h in the recovery, the system appears to stabilize irrespective of the homoeolog manipulated even though both dhrs3 homoeologs take longer to reach normal expression levels. These observations suggest that the RA network robustness response efficiently restores normal RA signaling levels irrespective of the homoeolog knocked down. Removing the more loosely regulated homoeolog slightly improves the robustness response.

Conclusions
The RA metabolic and signaling network is strongly dependent on the nutritional status and is influenced by the environment. Fluctuations in the RA signaling levels during embryogenesis can result in severe teratogenic outcomes. The preferential gene loss of the RA network components involved in the metabolism leading to RAL production suggests a selective pressure to achieve tighter regulation of the robustness response. Eliciting a robustness response by transient RA manipulation together with knockdown of specific gene homoeologs support the suggestion that gene loss might be linked to more efficient regulation of the network. The RA robustness response efficiently overcomes the reduction of one of the homoeologs. Tighter network regulation might involve loss of homoeologs similarly regulated, or homoeologs with enhanced responses. While the allotetraploid condition of X. laevis is convenient to explore these genomic changes and their regulatory outcomes, in diploid organisms, besides gene duplications and deletions, mutation, addition, and deletion of regulatory elements might take place to achieve the similar outcomes.