The Relationship between Gene Network Structure and Expression Variation among Individuals and Species

Summary

The variation generating mechanisms of development interact with the variation sorting mechanism of natural selection to produce organismal diversity. While the impacts of natural selection on existing variation have received much study, those of development on the generation of this variation remain less understood. This fundamental gap in our knowledge restricts our understanding of the key processes shaping evolution. In this study, we combine mathematical modeling, and population-level and cross-species assays of gene expression to investigate the relationship between the structure of the gene interactions regulating limb development and variation in the expression of limb genes among individuals and species. Results suggest that the way in which genes interact (i.e., development) biases the distribution of variation in gene expression among individuals, and that this in turn biases the distribution of variation among species.

Introduction

Phenotypic variation within populations is a prerequisite of evolution by natural selection, and in theory has the potential to bias the trajectory and rate of evolutionary change [1–6]. As such, identifying the processes that shape phenotypic variation has long been a fundamental pursuit of evolutionary biologists. Historically, evolutionary biologists have tended to focus on the sorting of population-level variation by selective processes, rather than on the production of that variation by developmental processes [7]. As a result, the effect of developmental processes on the distribution and magnitude of phenotypic variation among individuals and species remains unclear for most systems. In this study we use the mammalian limb as a study system to investigate two questions that address the relationship between developmental processes and phenotypic variation at the level of gene expression dynamics: (1) Does the structure of the gene network affect the distribution of variation in gene expression among individuals?, and (2) Is the distribution of variation in gene expression among individuals correlated with the evolutionary divergence in gene expression among species?

The mammalian limb is an ideal system for examining these questions because its development is well characterized, its morphology diverse, and since its form is central to many mammalian behaviors, its morphology is certainly under selection [8–12]. Many of the critical gene interactions that regulate limb outgrowth and patterning in mouse, the traditional mammal model, have been identified [9,10,13]. Initial budding of the limb from the body and limb outgrowth (embryonic day [E] 9.5 –E11) are regulated by interactions between several genes, including Bmp4, Gli3, Grem1, Shh, AER-Fgf’s (e.g., Fgf4, Fgf8), Fgf10, and Hox genes (Fig 1A). Knockouts of these genes result in pathological phenotypes ranging from severe (e.g., complete limb agenesis; AER-Fgf’s, Fgf10) to moderate (e.g., limb truncations; Bmp4) to mild (e.g., malformed digits; Shh, Gli3, Grem1) [13–18]. Most of these genes (e.g., Bmp4, Gli3, Grem1, Shh, AER-Fgf’s, and Hox genes) are also involved in later limb outgrowth and patterning (E11 –E13), but some of their interactions differ (e.g., Hox genes and Gli3, Shh, AER-Fgf’s; Fig 1B). As a result, the structure of the gene regulatory network differs for earlier (E9.5 –E11) and later (E11 –E13) limb development. This structural difference provides two opportunities to investigate the relationship between network structure and gene expression variation among individuals.

This structural difference also provides an opportunity to contrast earlier and later limb development. Research suggests that the main segments of the limb (e.g., stylopod, zeugopod, and autopod) are specified by or during the time of initial limb outgrowth [19,20]. As a result, disruption of early limb development could have potentially catastrophic effects on limb formation that are not likely to be selectively advantageous (e.g., limb agenesis). In contrast, disruptions of later limb development are less likely to have as severe an impact on the overall limb structure. While later disruptions might impact the relative size of limb segments, they are less likely to result in no limb at all. Following this logic, we might hypothesize that genes regulating early limb development generally exhibit less variation in expression among individuals than those regulating later limb development [21–33]. Additionally, it is possible that select early limb genes might vary at a level equal to or greater than that of individual later genes, but that this variation is dampened at the system level by the interactions among genes that characterize the gene network (i.e., developmental buffering) [34–37]. As population-level variation provides the raw material upon which natural selection acts, we can further hypothesize that the genes regulating early limb development also exhibit less variation in expression among species [38]. Support for these hypotheses would reinforce the importance of network structure (i.e., development) in shaping variation in mammalian limbs among individuals and over evolutionary time, while failure to support these hypotheses would suggest that network structure does not play a critical role in the generation of limb variation.

To test these hypotheses, we computationally modeled the gene networks regulating mouse limb development, and determined the sensitivity of network genes to system perturbation and the ability of network genes to perturb the system when altered. We also assessed the sensitivity of the system as a whole to perturbations in gene interactions and expression. We experimentally quantified naturally occurring variation in the expression of several network genes within a population of mouse individuals. To compare variation among species, we used transcriptomic data (RNASeq) from four mammals with divergent limb morphologies (bat, opossum, pig and mouse). We then assessed the relationship between gene and network sensitivity and gene expression variation among mouse individuals and among mammalian species. Our results suggest that the gene network that regulates early limb development is more robust than that regulating later limb development, and that this robustness buffers variation in early limb gene expression among individuals, and constrains the evolution of early limb development among species.

Results

Model construction (computational)

We assembled early (ESN) and late (LSN) stage networks for key genes regulating limb development from previously published experimental studies [9,13] (Fig 1). The ESN regulates initial limb outgrowth and the initiation of the epithelial-mesenchymal interactions that are critical to continued limb development. These events occur from embryonic days (E) 9.5 to E11 in mouse. The LSN, in contrast, regulates the limb’s differentiation along its anterior-posterior (i.e., thumb to pinky) axis and elongation along its proximal-distal (shoulder to fingertips) axis from E11 to E13 in mouse. To describe the temporal behavior of the activity (i.e., expression) levels of genes in these networks we built mathematical models (see S1 Methods).

Model simulations (computational)

After building the models, we ran a series of simulations in which we computationally interrupted interactions between genes and compared the resulting expression levels with those of the unaltered, default model (Table 1).

Within the ESN, removal of the Hox to Grem1 or Gli3R to Grem1 link affects Grem1 and Bmp4 expression levels (i.e., alters expression by 10% or more), but does not affect the expression levels of other genes. Removal of the AER-Fgf’s to Shh or Hox to Shh links only affects the expression levels of Gli3R and Shh, removal of the Fgf10 to AER-Fgf’s or Bmp4 to AER-Fgf’s links only affects the expression level of the AER-Fgf’s, and removal of the Hox to Fgf10 or AER-Fgf’s to Fgf10 links only affects the expression level of Fgf10. Removal of the Shh to Gli3R link affects the expression level of Gli3R, while removal of the Grem1 to Bmp4 affects the expression level of Bmp4. Removal of the Bmp4 to Grem1 link results in no significant change in expression levels. In total, 14 of 77 possible interactions are affected (i.e. expression levels change by 10% or more) by alterations in the ESN (18%) (Table 1).

For the LSN, removal of the Shh to Gli3R link affects Gli3R expression levels, but does not affect the expression levels of other genes. Removal of the AER-Fgf’s to Shh and Hox to Shh links affects only Shh and Gli3R expression levels. Removal of the Repressor X to Grem1 or Grem1 to Bmp4 link disrupts the expression levels of Bmp4, Grem1, and the AER-Fgf’s. Removal of the AER-Fgf’s to Repressor X link also disrupts the expression levels of Bmp4, Grem1, and the AER-Fgf’s, but also affects Repressor X expression levels. Removal of the Hox to Grem1 link affects the expression levels of Bmp4, Gli3R, Grem1, Shh, and the AER-Fgf’s. Removal of the Gli3R to Hox, AER-Fgf’s to Hox, or Gli3R to Grem1 links disrupts the expression levels of all genes save Repressor X. Finally, when the Bmp4 to Grem1, Bmp4 to AER-Fgf’s link is removed, expression levels of all genes are affected. In total, 52 of 84 possible interactions are affected (i.e. expression levels change by 10% or more) by alterations in the LSN (62%) (Table 1).

For each gene in each model, we then determined the number of genes whose removal alters expression of the gene in question (i.e. expression levels change by 10% or more), and the number of genes that exhibit expression changes when the gene in question is removed. We used the resulting values to generate a simulation space that was used to evaluate the ability of genes to affect other genes, and to be affected themselves by network perturbations (Fig 2C and 2D). Simulation spaces for the ESN and LSN were generated using the same scales to allow comparisons. For the ESN (Fig 2C), all genes fall in the lower left quadrant of the space, suggesting that they do not greatly affect expression of other genes, and are not greatly affected by others. In contrast, most LSN genes (Fig 2D) fall in the right upper and lower quadrants of simulation space. Genes in both the upper (e.g., Bmp4, Gli3R) and lower (e.g., AER-Fgf’s, Grem1, Shh) right quadrants and their boundaries are affected by perturbations in other genes, but genes in the upper right quadrant also affect the expression of other genes while genes in the lower right quadrant do not. Hox A/D falls in near the middle of the simulation space, suggesting that it moderately affects others and is affected by them. Only Repressor X falls in the lower left quadrant for the LSN, suggesting that it does not affect others and is not affected itself.

Sensitivity analyses (computational)

We next varied the parameter values used in the models and compared the resulting gene expression levels to those of the unaltered, default model (S1 Table). Results indicate that the ESN is most sensitive to changes in Hox A/D parameters (39% of Hox A/D parameter changes result in a ≥10% change in the expression level of another gene), followed by Shh (20%) and Bmp4 (18%), Gli3R (18%), and Grem1 (17%). The ESN is less sensitive to changes in AER-Fgf (10%) and Fgf10 (9%) parameters. Bmp4 is the most sensitive of the ESN genes to changes in parameters of other genes (32% of parameter changes result in a ≥10% change in the expression level of Bmp4), followed by Gli3R (26%), Grem1 (19%), Fgf10 (17%), and Shh (13%). Hox A/D (9%) and Fgf8 (7%) expression levels are less sensitive to ESN parameter changes.

The LSN is more sensitive to changes in Repressor X (73%), Bmp4 (60%), and the AER-Fgf’s (53%), and less sensitive to changes in Hox A/D (39%), Grem1 (39%), Shh (25%) and Gli3R (18%). Within the LSN, the AER-Fgf’s (61%) and Gli3R (63%) are the most sensitive to parameter changes in other genes, followed by Bmp4 (51%), Grem1 (47%), and Shh (40%), while Hox A/D (30%) and Repressor X (17%) are less sensitive.

The percentages listed above were used to generate a sensitivity space, similar to the simulation space described above (Fig 2A and 2B). ESN and LSN sensitivity spaces were generated using the same scales to facilitate comparisons. Similar to the simulation results, all ESN genes group within or on the boundary of the lower left quadrant of the space, suggesting that alteration of the values of their related-parameters does not greatly affect expression of other genes, and that their expression levels are not greatly affected by alterations in the values of the related-parameters of other genes. LSN genes are more distributed in the sensitivity space. Bmp4 and the AER-Fgf’s lie in the upper right quadrant, similar to their location in the simulation space (Fig 2C). Grem1, Shh, Gli3R and fall in or on the boundary of the lower right quadrant, indicating that they do not affect others but are affected themselves. Of these, Grem1 and Shh also fall within the lower right quadrant of the simulation space (Fig 2C). Repressor X lies within the upper left quadrant, suggesting that alteration of the values of its related-parameters affects the expression of other genes but that its expression is not greatly affected alterations in the values of the related-parameters of other genes. Repressor X also falls on the left side of the simulation space, but in the lower quadrant. Similar to its location in the simulation space, Hox A/D falls near the center of the plot.

Among individual variation in gene expression (qPCR)

We performed a series of real-time quantitative PCR (qPCR) assays to quantify the expression levels of genes that appear in both the ESN and LSN models (Bmp4, Gli3, Grem1, Shh, and the AER-Fgf Fgf8) in mouse embryos. For the early developmental stages (ES), the averaged, scaled expression level was 2.02 for Bmp4, 2.58 for Fgf8, 2.19 for Grem1, 1.68 for Shh, and 0.02 for Gli3. For the later developmental stages (LS), the average, scaled expression level was 2.29 for Bmp4, 0.83 for Fgf8, 1.92 for Grem1, 2.32 for Shh, and 0.02 for Gli3.

Statistical tests reveal that the mean-standardized variances of expression levels significantly differ among genes in the earlier (ES, E10-E11; Bartlett’s Test, F-ratio = 8.614, DF = 4, P < 0.001*) and later (LS, E11-E13; Bartlett’s Test, F-ratio = 5.823, DF = 4, P < 0.001*) stages of development. In the ES, Shh displays the highest average mean-standardized variance (coefficient of variation, CoV) (0.847), followed by Bmp4 (0.567), Grem1 (0.531), Fgf8 (0.474), and Gli3 (0.380). Fgf8 displays the highest average CoV in the LS (1.701), followed by Shh (1.523), Bmp4 (1.037), Gli3 (0.910), and Grem1 (0.703).

Litter membership also has the power to significantly explain the variance in expression levels in a given gene (e.g., Bmp4) that are observed among individuals (ANOVA; Bmp4 F-ratio = 2.379, DF = 8, P = 0.026*; Gli3, F-ratio = 5.742, DF = 8, P = < 0.001*; Grem1, F-ratio = 4.412, DF = 8, P = < 0.001*; Fgf8, F-ratio = 7.097, DF = 8, P < 0.001*; Shh F-ratio = 2.162, DF = 8, P = 0.043*).

Relationship between model predictions and gene expression variation among individuals

We next compared the among-individual, standardized variation in the expression level of a gene in vivo (CoV, from qPCR) with the: (1) number of genes whose removal alters expression of the gene in question (i.e., alters expression level by 10% or more), and (2) number of genes that exhibit expression changes when the gene in question is removed (from the simulation analyses).

For both the ES and the LS, neither the relationships between the number of genes whose removal alters expression of the gene in question (#1) and the CoV (ES—Least-Squares Regression, R² = 0.202, P = 0.448; LS—R² = 0.214, P = 0.433), nor the relationships between the number of genes that exhibit expression changes when the gene in question is removed (#2) and CoV (ES—R² = 0.503, P = 0.180; LS—R² = 0.142, P = 0.532) are significant (Fig 3).

For both the ES and LS, Shh is among the genes that are the least sensitive to perturbations in other genes, has a relatively low impact on the expression of other genes when altered, and displays a relatively high CoV. The opposite is true for Gli3R during the ES. Fgf8 displays the highest CoV during the LS and is highly sensitive to perturbations in other genes, like all LS genes.

Given the large difference in the percentage of possible interactions that are affected by computational alterations in the ESN and LSN (18% and 62%, respectively), we next compared the level of variation in measured gene expression during early (ES) and later (LS) development (from qPCR). For every model gene (5 of 5), the average CoV is greater later than earlier in development (P = 0.031*) (Fig 3). When the variation around the averages is taken into account using a resampling technique, the average CoV remains significantly greater later than earlier in development for four of the five genes (Bmp4, P = 0.028*; Gli3, P < 0.001*; Grem1, P = 0.024*; Fgf8, P < 0.001*). Only for Shh, a gene with among the highest CoV in both the early and later stages, does the average CoV does not remain significantly greater later than earlier in development (P = 0.092). The average CoV of the housekeeping gene β-actin does not significantly differ during earlier (0.832) and later (0.929) development (P < 0.217).

Gene expression variation among species (RNASeq)

To calculate gene expression variation among species, we first generated transcriptomic libraries for bat (Carollia perspicillata), opossum (Monodelphis domestica), pig (Sus scrofa) and mouse (Mus musculus) forelimbs for early (ES; early limb bud) and late (LS; paddle) limb stages. We then used a set of 6,583 genes orthologous to all four species (see S1 Methods) to calculate the among-species conservation of gene expression at each developmental stage, using the mean of all species pairwise Spearman coefficients. All resulting pairwise Spearman coefficients are positive and > 0.50, suggesting that the orthologous genes might perform similar functions between species (Fig 4A and 4B). However, the degree of gene expression conservation decreases from 0.5667 at ES to 0.5612 at LS of forelimb development. To test the robustness of this difference with respect to the selection of orthologous genes, we randomly sub-sampled 500 sets of orthologous genes at early and late stages at intensities ranging from 50 to 100% of all orthologous genes (Fig 4C). For each intensity, the distributions of gene expression conservation levels between early and late stages were significantly different (T-test, P-value < 0.05*). These results suggest that gene expression patterns vary more among species during the LS than ES of limb development, consistent with patterns of variation among individuals.

We also calculated the among-species conservation of mean-standardized expression at each developmental stage for the five genes that appear in both the ESN and LSN models (Bmp4, Gli3, Grem1, Fgf8, and Shh). At the ES, Gli3 (standard deviation [SD] = 0.177) falls within the top 25% of conserved genes while Shh (SD = 0.975) and Fgf8 (SD = 0.475) fall within the top 25% of divergent genes. Bmp4 (SD = 0.257) and Grem1 (SD = 0.260) fall near the middle of the range of genes. Fgf8 (SD = 0.971) and Shh (SD = 1.196) are also among within the top 25% of divergent genes at the LS. However, Grem1 (SD = 0.123) is among the 25% most conserved genes at this stage, and Bmp4 (SD = 0.253) and Gli3 (SD = 0.235) fall near the middle of the range of genes.

Average divergence level and average variation in expression levels among mouse individuals (CoV, as measured with qPCR) are positively correlated for the ES (R = 0.906) and LS (R = 0.854) (Fig 5), and the correlation between divergence level and CoV is significant for the ES (P = 0.019*) and LS (P = 0.016*) after bootstrapping.

Discussion

The results reported here suggest that the structure of the early stage network (ESN) renders it more robust to perturbation than the later stage network (LSN). Findings also suggest that among individual variation in expression levels is lower for genes regulating early (ES) than late (LS) limb development, and that gene expression levels are heritable. Results of this study also suggest that the expression levels of genes at early stages generally vary less among species. Additionally, results suggest that among individual and among species variation in the expression levels of several model genes are significantly correlated for the early and later stages of limb development. Taken together, these findings suggest a scenario in which a robust ESN buffers among individual variation in gene expression early in limb development, and, as variation is a prerequisite of evolution by natural selection, limits the evolution of early limb development among species. The findings of this study are therefore consistent with the hypotheses that: (1) the structure of the early limb gene network influences the distribution of variation in mammalian limbs among individuals and over evolutionary time [21–23,38–43], and, more generally, (2) the process of early limb development generally precludes the random accumulation of variation in gene expression across the network [21,44]. Results of this study are also consistent with a scenario in which species-specific differences accumulate as development progresses. However, it is important to note that the results of this study are based on a limited number of RNASeq samples per species. Analyses of additional samples are needed to determine the degree to which the RNASeq-driven results of this study are robust to experimental variation among samples. Furthermore, while the network models used in this study are based on solid experimental data [9, 13], the results of this study are only as accurate as the network models being used.

This study did not find a significant correlation between individual gene’s sensitivity to network perturbation or impact on the network when perturbed and variation in expression among individuals. This result could stem from the lack of a relationship between these variables, or from the incompleteness of the limb gene network used in this study. However, this study did find evidence that variation in expression level significantly differs among genes, with some genes being more variable among individuals in a population, and others less so. Shh displays among the greatest variation in expression levels during both the early and later stages of limb development among individuals and species, and has among the least impact on the system when altered. During the early stages of limb development, Gli3 displays the least variation in expression levels and has among the greatest effect on the system when altered.

As population-level variation provides a necessary prerequisite for evolution by natural selection, we might expect genes with the greatest expression variation among individuals to also display the most variation in expression among species. In this study this would be late stage Shh, Bmp4, and Fgf8. This study did find a significant correlation between the variation in the expression of these and the other key model genes among individuals and species during late limb development. Furthermore, of the few studies that have compared gene expression among mammalian limbs, a disproportionally high number have found differences in later stage Shh, Bmp4, and Fgf8 expression among species. Evolutionary changes in Bmp4 expression contribute to digit reduction in horses and jerboa, while evolutionary changes in Shh signaling contribute to digit reduction in pigs [45]. A broader initial range and secondary redeployment of Shh signaling helps generate the unique phenotype of the bat wing [46], and hind limb loss in dolphins is initiated by a disruption in Shh signaling [47]. Shh signaling is also activated exceptionally early during the rapid outgrowth of opossum forelimbs [48,49]. Fgf8 expression is higher in the AER of bat wings than mouse limbs and is also secondarily redeployed later in bat wing development [50]. However, genes beyond Shh, Bmp4, and Fgf8 also display expression differences in mammalian limbs. For example, the expression of 5’ Hox A/D genes differs in bat and kangaroo limbs [51,52], compared to mouse [53]. Clearly more studies of limb development in diverse mammalian species are needed to resolve this issue, but results to date are consistent with alterations in the expression of genes acting during late limb development (E11 to E13 in mouse) including Shh, Bmp4, and Fgf8 frequently contributing to mammalian limb evolution.

In line with the proposed predominance of changes in late limb development in limb evolution, this study found no evidence for the existence of significant variation in the expression of early limb genes that is masked by systems-level processes. However, the expression-based findings of this study do not rule out the existence of cryptic genetic variation in genes with roles in early limb development. Cryptic genetic variation can provide a source of evolutionary potential when uncovered by environmental or genetic perturbations [54,55], and thereby can expedite evolutionary change. Thus, if the genes regulating early limb development possess significant cryptic genetic variation that has been uncovered over evolutionary time, we would expect early limb development to vary among species. Early limb development, as defined in this paper, encompasses establishment of the limb field and the initial outgrowth of the limb. The primary limb segments (e.g., stylopod, zeugopod and autopod) are also likely specified during this time [19,20]. Initial limb outgrowth appears to be generally conserved in tetrapods across large phylogenetic distances [56,57]. Initial limb outgrowth is even conserved in some tetrapods that do not possess limbs in their adult form (e.g., boas, dolphins) [47,58]. The primary segments of the limb are also broadly conserved across limbed tetrapods [8,59]. These observations together with the findings of this study suggest that the genes regulating early limb development either do not possess significant cryptic genetic variation, or that the robustness of the ESN has inhibited the ability of environmental or genetic perturbations to uncover this variation. Whichever the case, the early development of the limb appears to have been relatively conserved over the evolutionary history of tetrapods.

Materials and Methods

Ethics statement

All animal work was conducted according to relevant national and international guidelines. Animals were euthanized using CO2 inhalation followed by cervical dislocation. The University of Illinois IACUC approved this research (protocols #13128, 14159, 14199, 14209).

Model construction

The starting point for the mathematical models used in this work was the seminal paper by Bénazet et al. 2009 [13]. These models were designed to pertain to the entire limb. Two interconnected feedback loops were incorporated into their model: a fast loop between Grem1 and Bmp4 and a slower loop between Shh, Grem1 and the AER-Fgf’s. Following the findings reported in Sheth et al. 2013 [9], this network was divided into an Early (ESN) and Late (LSN) Stage Network and augmented to also include Hox genes (specifically, Hox A and D genes), Gli3R, Fgf10, and an as yet unidentified repressor, Repressor X. It is important to note that the specific genes, gene interactions, and equations that we include in our models match those presented in Bénazet et al. 2009 and Sheth et al. 2013, which are well supported by experimental, biological evidence.

Mathematical models to describe gene interactions were constructed in MATLAB and were based on ordinary differential equations, which are outlined in S1 Methods, along with the model parameters.

Model simulations

We ran a series of simulations in MATLAB in which we removed interactions between genes. Only one interaction (i.e. link between two genes) was removed in each simulation, and removal simulations were performed for all interactions.

Sensitivity analyses

We ran a series of simulations in MATLAB in which we varied the parameter values used in the models. Only one parameter was modified in each analysis.

Among individual variation in gene expression

We performed a series of real-time quantitative PCR (qPCR) assays to quantify the expression levels of the genes that appear in both the ESN and LSN, namely Bmp4, Gli3, Grem1, Shh, and the AER-Fgf Fgf8, in the forelimbs of 71 mouse embryos (outbred ICR strain, Taconic) from 9 litters.

These litters ranged in age from E10 to E13. The limbs of the animals in these litters were staged according to Wanek’s staging guide, which divides mouse limb development into 15 stages [60]. Embryonic limbs from limb ridge (Wanek stage 1) through bud formation (Wanek stage 4; 4 stages total; E10 –E11) were grouped into an “early stage” (ES) for analyses, while limbs in the paddle stages of development (Wanek stages 5–8; 4 stages total; >E11 –E13) were grouped into a “late stage” (LS). Limb samples were evenly spread over all stages. The Coefficient of Variation (CoV), which is the standard deviation divided by the mean, was used to quantify variation in expression level for a given gene for early and later limb development [61]. Additional details for the qPCR analyses are in S1 Methods, and standard and dissociation curves for each gene are in S2 Methods.

Bartlett’s Test was used to compare the variance of expression levels in the ES and LS [61]. ANOVA was used to examine the contribution of litter membership (i.e., heredity) to observed patterns of gene expression [61].

Relationship between model predictions and gene expression variation among individuals

The relationship between among individual variation in gene expression level with the number of genes whose removal alters expression of the gene in question by 10% or more from the default value (from simulation analyses) and whose values deviate by 10% or more from the default values when the gene in question is removed (from simulation analyses) was statistically assessed using Least-Squares Regression [61]. The average levels of variation in measured gene expression during early and later development (CoV) were statistically compared [61]. To determine the significance of the observed differences in CoV, we used a Monte Carlo approach, in which we shuffled the observed CoV from early and later development, to generate a null distribution of CoV differences between them. Specifically, we pooled all replicate CoV irrespective of developmental stage, then randomly drew, with replacement, two samples equal in size to the measured early and later samples, respectively. We used as a measure of significance the proportion of 10,000 replicates in which the difference between CoV’s of randomly shuffled samples was greater than or equal to the observed difference.

Gene expression variation among species (RNASeq)

Embryonic mice, opossums, bats and pigs with early (ES) or late (LS) stage forelimbs were obtained from a variety of sources (see S1 Methods). Forelimbs for the ES were harvested at Stage 14 for bat, E11 for mouse, Stage 28 for opossum and E22 for pig [62–65]. For the LS, forelimbs were harvested at Stage 15 for bat, E12 for mouse, Stage 29 for opossum and E26 for pig.

Limbs were removed from embryos and stored in RNALater in -20°C until further processing. RNA was extracted from tissues using E.Z.N.A. Total RNA Kit I (OMEGA bio-tek #R6834), and converted into RNASeq libraries with the Illumina TruSeq RNA Sample Preparation Kit (Illumina RS-122-2001). Libraries were sequenced on an Illumina HiSeq 2500 housed in the Roy G. Carver Biotechnology Center at the University of Illinois. Resulting reads were processed, aligned to published genomes or de novo assemblies, and gene expression levels assessed (see S1 Methods). All data from this study have been deposited in the Gene Expression Omnibus (GEO) with the accession number GSE71390.

We analyzed the conservation of the gene expression profiles of bat, mouse, opossum, and pig across embryonic limb development, using the mean of all species pairwise Spearman coefficients (see S1 Methods). The relationship between these average species pairwise coefficients and among individual gene expression variation was assessed using Pearson Product-Moment Correlation [61]. To account for the variation among samples we used a bootstrap approach. Specifically, for the comparisons of the average species pairwise coefficients and among individual gene expression variation, we resampled, with replacement, the CoV’s (among individuals) and mean-standardized expression levels (among species) and recalculated the Pearson product-moment correlation coefficient (R) between the resampled CoV’s and mean-standardized expression levels. We used as a measure of significance the proportion of 10,000 replicates in which the calculated R was greater than or equal to zero. We also ordered the species pairwise Spearman coefficients for each individual gene from highest (most conserved) to lowest (most divergent).

Supporting Information

Attachment 1

Attachment 2

Attachment 3