Cross‐species transcriptomics uncovers genes underlying genetic accommodation of developmental plasticity in spadefoot toads

That hardcoded genomes can manifest as plastic phenotypes responding to environmental perturbations is a fascinating feature of living organisms. How such developmental plasticity is regulated at the molecular level is beginning to be uncovered aided by the development of ‐omic techniques. Here, we compare the transcriptome‐wide responses of two species of spadefoot toads with differing capacity for developmental acceleration of their larvae in the face of a shared environmental risk: pond drying. By comparing gene expression profiles over time and performing cross‐species network analyses, we identified orthologues and functional gene pathways whose environmental sensitivity in expression have diverged between species. Genes related to lipid, cholesterol and steroid biosynthesis and metabolism make up most of a module of genes environmentally responsive in one species, but canalized in the other. The evolutionary changes in the regulation of the genes identified through these analyses may have been key in the genetic accommodation of developmental plasticity in this system.


| INTRODUC TI ON
Adaptive phenotypic plasticity is a common feature of biological organisms that allows a given genotype to appropriately express alternative phenotypes under different environmental conditions (Schlichting & Pigliucci, 1998;West-Eberhard, 2003). The genetic regulation of such plasticity can diverge under selection in contrasting environments, resulting in evolved differences across lineages in their degree of responsiveness to environmental variation (Casasa & Moczek, 2018;Schlichting & Wund, 2014;West-Eberhard, 2003).
Environmental heterogeneity is the norm in nature, and developmental processes can either be responsive or resistant to such environmental fluctuations. The evolution of developmental responses to environmental factors therefore must vary along this gradient between plasticity and robustness (Lafuente & Beldade, 2019;Schwab et al., 2019). Adaptive plasticity evolves readily in organisms exposed to environmental heterogeneity given reliable cues to assess the environment, capitalizing on the modularity of development (Londe et al., 2015;Moczek, 2010) and its underlying regulatory gene networks (Projecto-Garcia et al., 2019;Schneider et al., 2014;Snell-Rood et al., 2009). During transitional phases to a novel or different environment, plasticity increases the viability and persistence of populations, contributing to the maintenance of genetic variability due to both simple demographic effects and to the shielding of genotypic variants from selection (Draghi & Whitlock, 2012;Gomez-Mestre & Jovani, 2013). If different plastic lineages (species or populations) experience consistently divergent environmental conditions, they can evolve differences in the environmental | 2221 LIEDTKE ET aL. sensitivity of the regulation of adaptive traits (i.e., experience genetic accommodation; West-Eberhard, 2003). Trait plasticity can therefore increase or decrease, and even be lost if organisms are exposed to stable environments for long enough (Kulkarni et al., 2011;Pigliucci et al., 2006;Suzuki & Nijhout, 2006). Such reduction in environmental sensitivity and increased developmental robustness may be achieved through a variety of mechanisms (Lafuente & Beldade, 2019), such as redundancy in gene enhancers (Frankel et al., 2010) and regulatory microRNAs (Brenner et al., 2010).
Environmentally induced changes in development often result from environmentally induced changes in gene expression (Beldade et al., 2011;Lafuente & Beldade, 2019). Consequently, a promising approach to the study of genetic accommodation is to compare the environmental sensitivity in gene expression of species or populations with divergent responses to standardized environmental changes i.e., their transcriptomic reaction norms. By characterizing the networks of coexpressed genes in each species and their responses to standardized environmental stimuli, we could compare the transcriptomic sensitivity to environmental conditions in the different lineages (Casasa et al., 2020), and identify modules within networks of orthologous genes associated with differences in plasticity across lineages.
Ancestral plasticity in developmental rate in spadefoot toad larvae has evolved into genetically accommodated differences among species in their developmental sensitivity to pond drying (Gomez-Mestre & Buchholz, 2006;Kulkarni et al., 2017). Pelobatoid frogs encompass the largest differences in developmental rate within anuran amphibians. At the two ends of this spectrum are the Western spadefoot toad (Pelobates cultripes) and Couch's spadefoot toad (Scaphiopus couchii). In benign conditions of high-water levels and food abundance, P. cultripes larvae develop slowly and attain large sizes at metamorphosis, although they are capable of accelerating their development upon detection of a reduction in water level, shortening their larval period by as much as 40%. This environmental sensitivity is controlled by the hypothalamic-pituitary-thyroid (HPT) and the hypothalamic-pituitary-inter-renal (HPI) axes (Denver et al., 2002), which regulate developmental acceleration by increasing thyroid hormone and corticosterone levels, expression of thyroid hormone receptors and metabolic rate Kulkarni et al., 2017). However, developmental acceleration comes at the expense of reduced size at metamorphosis, consumption of their fat bodies and increased oxidative stress Kulkarni et al., 2011Kulkarni et al., , 2017. In comparison, S. couchii tadpoles have greatly reduced the ancestral plasticity and have evolved an overall faster developmental rate (the fastest known in anurans; Buchholz & Hayes, 2002), and accordingly, show a constitutively higher metabolic rate and base levels of thyroid hormone and corticosterone compared to P. cultripes (Kulkarni et al., 2017). Still unclear however, is how these accommodated changes in developmental plasticity across species are reflected at the transcriptomic level.
Under the null hypothesis that tadpoles of both species show at least some sensitivity to ponds drying out, a reduction in water level will trigger the HPT and HPI axes to accelerate development. Gene expression in S. couchii is expected to experience greater differences over the same absolute time because development proceeds at a faster pace than in P. cultripes. Alternatively, because greater plasticity is expected to be the result of greater differential gene expression (e.g., Casasa et al., 2020), we also hypothesize that P. cultripes, with its significantly greater developmental plasticity, will experience accordingly greater changes in differential gene expression between normally developing and induced, fast developing tadpoles. Taking advantage of recently published de novo transcriptomes of P. cultripes and S. couchii (Liedtke et al., 2019), we test these hypotheses by performing RNA-Seq experiments in which we exposed tadpoles of both species to simulated pond drying at the same developmental stage and compare their gene expression profiles over multiple time points to control tadpoles reared in constant water levels.

| Experimental design and transcriptome assembly
Sample collection, transcriptome sequencing and de novo assembly are detailed in Liedtke et al. (2019) and published as part of the NCBI BioProject PRJNA490256. In brief, we individually raised tadpoles from single clutches (i.e., siblings) of P. cultripes and S. couchii in 3 L aquaria (water column = 18.5 cm) under laboratory conditions until they reached a stage at which developmental responses to water level reduction are maximum in these species (Gosner stage 35;Gosner, 1960;Kulkarni et al., 2011). At this point, we dropped the water level of half the tadpoles from each species to 675 ml (4 cm) to trigger developmental acceleration. We very briefly removed tadpoles from their containers to extract the water, also manipulating in the same way tadpoles remaining in the constant water treatment. After that, tadpoles in the constant water treatment of both species were maintained for 24 h. We maintained P. cultripes tadpoles in reduced water for a period of 24, 48 or 72 h, and S. couchii tadpoles for 24 or 48 h. Scaphiopus couchii develops so quickly that we did not include a 72 h time point in our RNA-Seq analysis because development may have progressed within that timeframe. Within the timeframe used, however, we did not observe advancement in developmental stage for either species, and the comparison was therefore stage-matched (at Gonser 35). This experimental procedure follows  and elicits accelerated development in these species, even if the response of S. couchii is much less pronounced than that of P. cultripes (Kulkarni et al., 2011(Kulkarni et al., , 2017. At the end of each time period, and after 24 h for the controls, we euthanized three tadpoles per species, eviscerated and homogenized them for total RNA extraction. We extracted RNA using a Trizol extraction method following the manufacturer's protocol (Invitrogen), and sequenced 21 TruSeq libraries using a HiSeq2000 (Illumina) sequencer in paired end mode with the read length 2 × 76 bp, generating on average 38 million paired-end reads per library. We assembled de novo transcriptomes for both species, based on 12 libraries for each species, using the Trinity pipeline (Haas et al., 2013), annotating them with Trinotate (https://trino tate.github.io/).

| Transcript abundance and differential gene expression
A schematic representation of the bioinformatic workflow and analytical procedures followed in this study is shown in Supporting Information 1. We estimated transcript abundance using kallisto v04.3.1 (Bray et al., 2016) at the Trinity "gene" level and we explored clustering of individuals by treatment using a principal component analysis on the standardized counts (log 2 counts per million). We used edger v3. 26.6 (McCarthy et al., 2012) to identify differentially expressed transcripts for across all possible combinations of pairwise comparisons of treatment group per species (a total of six comparisons for P. cultripes and three for S. couchii). Each group consisted of three biological replicates, and we deemed genes with a false discovery rate of q < 0.05 and absolute log fold change >2 as significantly differentially expressed.
To identify genes that show correlated expression patterns over time, we singled out genes that were differentially expressed for at least one pairwise comparison and performed a clustering analysis on their standardized expression values across all treatments. We normalised the count data for these genes using the trimmed mean of M-values (TMM) method in edger and averaged across replicates to obtain mean expression per treatment. Subsequently, we standardized expression values per treatment to have mean = 0, SD = 1 and clustered using the soft clustering in mfuzz v2.44.0 (Kumar & Futschik, 2007). We set the fuzzification parameter to 1.8 and determined the number of clusters by graphically exploring the minimum centroid distance between a range of cluster numbers, choosing as cutoff the number of clusters above which this distance was no longer substantially reduced. We then assigned genes to clusters based on their maximum membership probability.
We performed functional enrichment analyses using g:Profiler

| Coexpression analysis -orthologue assignment
We identified orthologues between P. cultripes and S. couchii using the HaMStR pipeline (Ebersberger et al., 2009). HaMStR compares query sequences (based on open reading frames: ORFs) with reference proteins from curated databases and profile hidden Markov models (pHMMs) derived from alignments of reference proteins. Orthologues are identified when there is reciprocity between query sequence, reference protein and pHMM. We downloaded our reference proteins from orthodb v.9 (Kriventseva et al., 2007(Kriventseva et al., accessed 28/03/2018. As only one amphibian species was available from orthodb (the western clawed frog, X. tropicalis), we completed our reference protein database with two reptile species (the Chinese softshell turtle Pelodiscus sinensis and the American Alligator Alligator mississippiensis) and the West Indian Ocean coelacanth Latimeria chalumnae. We aligned reference protein sequences using prank v.170427 (Löytynoja, 2014), and created the pHMMs using hmmer v.3.1b2 (Zhang & Wood, 2003). HaMStR uses hmmsearch to match query sequences with pHMMs and then BLASTP to compare query sequences to the original reference proteins from X. tropicalis. Query sequences for which the top scoring X. tropicalis BLASTP hit also contributed to the pHMM hit passed the reciprocity test and were added to a list of target species orthologues.
When multiple orthologues matched a single ORF, we removed the hits based on: (i) the taxonomic specificity of the reference protein (tetrapod proteins took precedence over metazoan proteins); (ii) whether the ORF was the most representative hit (if an ORF was the best hit for one orthologue, but not for one or more other orthologues, we removed it in these other cases); and (iii) E-values, in cases where an ORF was the best hit for more than one orthologue, or in cases where it was never the best hit.

| Coexpression analysis -network analysis
We normalised transcript count data (Trinity "isoform" level) for each species separately using the trimmed mean of M-values (TMM) method in edger and extracted counts per million (cpm) values.
When a transcript contained more than one ORF for which there was a matching orthologue, we divided cpm values by the number of ORFs contained within that transcript. We then matched these adjusted cpm values to their respective orthologue. In many cases an orthologue was represented by multiple isoforms. In these cases, we averaged cpm values across isoforms to determine a single representative orthologue count value. We then combined results for both species and only took forward for network analysis orthologues that were present in both species.
We performed a weighted gene correlation network analysis (WGCNA; Langfelder & Horvath, 2008) for the combined data set containing count values for all 21 libraries. Such combination of transcriptomes from multiple species in network analyses can highlight evolutionarily shared and divergent patterns of gene expression (Morandin et al., 2016). We checked data to ensure that orthologues did not contain too many missing values, and assessed the scale-free topology to select an appropriate soft thresholding power by which to raise pairwise correlations of expression. We transformed orthologue coexpression similarity into a signed adjacency matrix based on the selected soft threshold power, and hierarchical clustering identified modules of highly interconnected genes (with a minimum module size of 30). We merged highly similar modules, i.e., those for which module eigengenes were highly correlated (cutoff set at 0.2). The eigengene is defined as the first principal component of a module and provides a single representative profile of orthologue expression for a given module.

| Relation of coexpression modules to experimental factors
We assessed the relationships between coexpression modules and species, developmental time and water level treatment by fitting linear models with eigengene values for each module as the response variable. We considered species and water level as factors, and we included time (24 and 48 h only) as a covariate. We assessed the effects of the factors and covariate using the ANOVA function from the car package in r (Fox & Weisberg, 2011). We also initially included the interaction terms species:time, and species:water, but removed them when not significant (p > .05), refitting the model.

| Coexpression analysis -enrichment analysis and visualization
We carried out functional enrichment analysis of gene ontology terms, KEGG and Reactome gene pathways for modules of interest, i.e., those significantly associated with species, water level, time or one of the interactions. We carried out the enrichment analysis in the r package gprofiler2, using as background domain the previously identified orthologues from the X. tropicalis proteome.
To further explore the orthologues involved in response to water level treatment, we calculated gene significance (GS) values in WGCNA (defined as the absolute value of the correlation between orthologue and a given variable, factor or trait ;Langfelder & Horvath, 2008). For all modules, we extracted orthologues with a GS p-value (p.GS) <.05 and visualised them using cytoscape v 3.8.0 (Shannon et al., 2003), with nodes (orthologues) linked by edges (expression correlations) determined by WGCNA. We only visualised edges with edge values >0.025, removing nodes without edges.

| Transcript abundance and differential gene expression
The principal component analysis on the standardized count data (counts per million) for P. cultripes samples resulted in overlapping clusters of treatment groups on the first two axes, together explaining 32.76% of the variance (Figure 1a). Out of a total of 428,406 Trinity "genes", 875 were differentially expressed in at least one pairwise treatment comparison (0.20%), which mfuzz soft clustering grouped into seven distinct expression profiles (Figure 1b).
The majority of genes (446 genes; clusters 1 and 2) were upregulated after 48 or 72 h of exposure to low water levels, whereas 236 genes (clusters 4, 6 and 7) were downregulated either after 24 or    cultripes was therefore experiencing a delay of up to 48 h in changes in activity of some pathways and processes whose response to pond drying is conserved across the two species. Among these is apoptosis   Module eigengenes were significantly associated with one or more explanatory variables or interactions for five of the eight modules (Table 1; Figure 3). The "blue", "green", "turquoise" and "greenyellow" modules were the four largest modules (together containing 11,368 out of 11,939 orthologues), and each of these was highly significantly associated with species (see Table 1 for F values and p-values). The "blue" (2700 orthologues) and "green" (473 orthologues) modules showed higher expression in P. cultripes than S. couchii (Figure 3b,c), while the "turquoise" (7731 orthologues) and "greenyellow" (464 orthologues) modules showed higher expression in S. couchii relative to P. cultripes (Figure 3d,e). The "cyan"

| Coexpression analysis -network analysis
module was relatively small (58 orthologues), but showed significant species-by-water interaction (F = 8.15, p = .014; Table 1): P. cultripes exhibited a steep decrease in eigengene expression in response to water drop, while S. couchii eigengene expression did not respond to water level (Figure 3a).   Table 2. Here, we briefly describe some of the main functionally enriched pathways from the Xenopus Reactome. For the "cyan" module, associated with species differences in developmental plasticity to decreased water level, the five most significantly enriched Reactome pathways (Table 2; Supporting Information 7) were associated with general metabolism The "green" module, which showed higher eigengene expression in P. cultripes than S. couchii was associated with immune system and complement pathways (REAC:R-XTR-168256, REAC:R-XTR-166658

| Coexpression analysis -functional enrichment analysis and network visualisation
and REAC:R-XTR-166663), while the "greenyellow" module, which showed higher expression in S. couchii than P. cultripes, was associated with the cell cycle and mitosis (R-XTR-69278, R-XTR-1640170, REAC:R-XTR-68882, R-XTR-68886, and R-XTR-2555396). The "blue" and "turquoise" modules were not functionally enriched for any Reactome pathways, probably due to their large sizes and unspecialised natures, although various GO biological processes and cellular components associated with muscular development were enriched in the "blue" module (higher eigengene expression in P. cultripes) and the GO molecular function "transferase activity" was F I G U R E 3 Comparison of P. cultripes and S. couchii eigengene expression for five modules (a-e) in which either species as a main effect or the species:water interaction was a significant factor (as indicated  Across all modules, 149 orthologues were found to have a significant (p.GS < .05) association with water (Supporting Information 8). The network visualization (Figure 4) shows 52 of these waterassociated orthologues which had at least one interaction (adjacency >0.025) with another water-associated orthologue. Orthologues from the "cyan" and "turquoise" modules showed greater numbers of intramodule connections. The strongest pattern of coexpression (many adjacency values >0.1) was observed in six orthologues in the "cyan" module: mevalonate kinase, mevalonate (diphospho) decarboxylase, NAD(P) dependent steroid dehydrogenase-like, phosphatidylserine synthase 2 and cystathionine gamma-lyase. All six of these orthologues contributed to at least one of the significantly enriched Reactome pathways (general, lipid, and steroid metabolism, cholesterol biosynthesis and cysteine formation from homocysteine).
The orthologues with the strongest associations with water are not necessarily shown in the network due to low adjacency values with other nodes, although stomatin (EPB72)-like 2 and TP53 regulating kinase from the "cyan" module, and uroporphyrinogen III synthase and pancreas specific transcription factor 1a (PTF1) from the "turquoise" module were among the top 10 most significant (p < .004) and are included in Figure 4.

| DISCUSS ION
Transcriptomic approaches permit a more comprehensive screening for environmentally sensitive gene regulatory networks than candidate gene approaches, and are thus advantageous for studying complex evolutionary processes such as developmental plasticity and genetic accommodation (Aubin-Horth & Renn, 2009; F I G U R E 4 Network of coexpressed orthologues with significant water effects (p.GS < .05) across all modules. Nodes are coloured based on module of origin (see Figure 3). Edges with adjacency values <0.025 are not included, and nodes without at least one edge were excluded. Edge thickness corresponds to adjacency value and node size corresponds to gene significance for water (the largest nodes having p.GS values < .001). Those orthologues in the water-sensitive cyan module which contributed towards enriched Reactome functions are highlighted with node borders (thicker borders indicating contribution to more functions) and enlarged labels. The Reactome functions and full names of abbreviated orthologues can be found in Supporting Information 8 Edge adjacency Beldade et al., 2011;Sommer, 2020;Young, 2013). The prevailing signal in the transcriptomic data was attributable to species differences, regardless of experimental manipulation of water level. This is not surprising given more than 120 million years of independent evolution between these species and considerable divergences of their genomes (Liedtke et al., 2018;Zeng et al., 2014). As has been shown for other systems (e.g., Casasa et al., 2020), we initially hypothesized that P. cultripes would manifest greater regulatory response to a reduction in water level than S. couchii, due to its more pronounced developmental plasticity. However, the developmentally more canalized S. couchii showed around eight times the number of differentially expressed genes in response to reduced water level compared to P. cultripes. The development of S. couchii is also known to be responsive to pond drying (Kulkarni et al., 2011;Newman, 1989), but this result is nonetheless unexpected, given that in comparison to P. cultripes, its developmental rate is very robust to external perturbations including variation in water level, the addition of exogenous corticosterone, or the inhibition of corticosterone synthesis (Kulkarni et al., 2017). However, physiological systems consist of complex networks that arise from interacting units and feedback loops (Burggren & Monticino, 2005;Sturmberg et al., 2015) and physiological processes in organisms in their prime state are tightly regulated so that every small deviation is immediately corrected (Goldberger et al., 2002;Vargas et al., 2015). Maintaining the course of development in S. couchii in the face of environmental fluctuations may require a buffering mechanism that constantly adjusts the underlying developmental and physiological processes, thus also causing frequent adjustments in gene expression.
The developmental and metabolic rates of S. couchii tadpoles in benign conditions are substantially higher than those of P. cultripes, and thus, the observed results also fit our null hypothesis; only slight acceleration (i.e., lower plasticity) of an already accelerated system (S. couchii) may be enough to result in greater gene expression change in very short periods of time. This is supported by the network analysis which shows that the majority of orthologues are assigned to modules with consistently higher expression levels in S. couchii regardless of water level (the S. couchii associated modules "turquoise" and "greenyellow" contained more than 2.5 times as many orthologues as the P. cultripes associated modules of "blue" and "green"). These highly expressed orthologues make up pathways related to cell cycle, mitosis and DNA replication, possibly reflecting the faster absolute development in this species. In addition, P. cultripes shows a lagged or slower response, which may further explain the lower number of differential gene expression in the first 72 h.
The single largest mfuzz cluster in S. couchii characterizes expression changes in the first 24 h (cluster 1), but for P. cultripes, the biggest response is seen only after 72 h (cluster 1), representing a clear asynchrony in the peaks of expression changes in important processes such as apoptosis. Apoptosis is a major process in amphibian metamorphosis, responsible for the cell death in larval-specific tissue, but also in larval-to-adult remodelling of organs (Ishizuya-Oka et al., 2010; Kerr et al., 1974).
Despite the large differences in the transcriptomic profiles of the two species, the network analysis identified a module of genes whose expression in response to pond drying showed a divergence in transcriptomic reaction norms between species (i.e., a significant species-by-treatment effect). More precisely, genes of this module experienced significant differences in expression in response to pond drying in P. cultripes, but not S. couchii. We interpret this pattern as a signal of genetic accommodation (Sikkink et al., 2019) (Mani et al., 2015), showed moderate coexpression with these orthologues and also occupied an important position between the environmentally sensitive "cyan" module and the S. couchii associated "turquoise" module. These findings suggest that synthesis and regulation of cholesterol is a key process underlying developmental plasticity in P. cultripes, and that its regulation has diverged in S. couchii. Differences in pathways associated with lipid metabolism reflect the fact that P. cultripes tadpoles experience a >2-fold increase in metabolic rate and deplete their fat reserves during developmental acceleration. In turn, S. couchii tadpoles maintain a higher and constant metabolic rate regardless of water level, and produce almost no fat bodies even if metamorphosis is chemically blocked and larvae are fed ad libitum under benign conditions (Kulkarni et al., 2011(Kulkarni et al., , 2017. Another two orthologues in the divergent plasticity module which had particularly high gene significance were TP53 regulating kinase and transmembrane Bax inhibitor motif containing 6 (TMBIM6). TP53 is an important regulator of apoptosis (Yogosawa & Yoshida, 2018), whereas Bax is known to induce apoptosis during amphibian metamorphosis (Sachs et al., 2004). Reduced expression of TP53 regulating kinase and the Bax inhibitor TMBIM6 are probably involved in the programmed cell death of larval tissues during accelerated development.
The timing of amphibian larval development is primarily controlled by thyroid hormone acting together with corticosterone (Denver, 1997;Sachs & Buchholz, 2019). Although the role of corticosteroids in accelerating development is less well understood than that of thyroid hormone, it is clear that CORT accelerates thyroid hormone-induced metamorphosis with peak blood levels occurring at the climax of natural metamorphosis (Kulkarni & Buchholz, 2014). Moreover, endogenous CORT levels are known to be facultatively increased in P. cultripes, but consistently high in S. couchii in the face of pond drying (Kulkarni et al., 2017). In this study, we find that the genes coding for deiodinase 3 (dio3), an enzyme important for the activation of thyroid hormone and the transcription factor Krüppel-like factor 9 (KLF9), downstream target genes of corticosterone in amphibians (Bonett et al., 2009;Kulkarni & Buchholz, 2014), both experienced a significant positive fold change at 72 h of low water exposure in P. cultripes, but were not significantly differentially expressed in S. couchii. We interpret this to be confirmatory that genes directly related to the HPI axis were differentially expressed in P. cultripes in response to reduced water level, but were not as environmentally sensitive in

S. couchii.
This first genome-wide scan of the molecular mechanisms of plasticity and canalization of developmental acceleration in two spadefoot toad species is a significant step in our understanding of the evolution of plasticity. Despite the great divergence due to the long period of independent evolution between these species and the fundamentally different rates of development, we could identify a relatively small set of genes that are involved in their divergent responsiveness to the risk of pond drying, i.e., genetic accommoda- The role of developmental plasticity in evolution has gained prominence in recent years, but how plasticity itself evolves at the molecular level is only recently being elucidated (Casasa et al., 2020;Levis et al., 2020;Sikkink et al., 2019;Sommer, 2020). By studying two species of spadefoot toad with differing degrees of plasticity in response to pond drying, we identified orthologues and functional gene pathways whose environmental sensitivity in expression (i.e., their transcriptomic reaction norms) have diverged. In doing so, we discerned transcriptomic modules whose evolutionary changes in regulation may be linked to genetic accommodation of developmental plasticity in this system.

ACK N OWLED G EM ENTS
This project was funded by the Ministerio de Economía, Industria y Competitividad (MINECO) through grant CGL2017-83407-P. We also thank D. Buchholz for important and insightful discussions.

AUTH O R CO NTR I B UTI O N S
Ivan Gomez-Mestre and Hans Christoph Liedtke designed the study and carried out all laboratory work. Hans Christoph Liedtke, Ewan Harney and Ivan Gomez-Mestre conducted analyses and wrote the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The raw RNAseq ad the transcriptome assemblies have been deposited on the NCBI's Sequence Read Archive (SRP161446) and Transcriptome Shotgun Assembly database (P. cultripes: GHBH01000000; S. couchii: GHBO01000000), under BioProject PRJNA490256. All processed data sets used in this article are available as Supporting Information (1-8) deposited on dryad (https:// doi.org/10.5061/dryad.cnp5h qc3z).