Distinct gene expression dynamics in developing and regenerating crustacean limbs

To investigate transcriptional dynamics and assess individual variation in developing embryonic legs we performed RNA-seq on individual fourth thoracic (T4) legs during the time course of leg development, from young limb bud stages to fully patterned and differentiated legs (21) (Fig. 1A). We collected entire T4 legs every 6 h, from 96 to 192 h post fertilization (hpf) (E_f [full] series). In order to account for the progressive regionalization of the primordium, we also collected, when possible, the distal portion of the T4 leg: the distal one-third of the leg from 120 to 138 hpf, and the carpus, propodus, and dactylus from 144 to 192 hpf (E_d [distal] series). The full and the distal leg samples were collected in pairs, from contralateral T4 legs of the same embryos, yielding a total of 70 samples covering the time course of leg development (SI Appendix, Table S1).

Principal component analysis (PCA) on the complete RNA-seq dataset (Datasets S3 and S4) showed that the principal axis of gene expression variation, explaining 14% of the variance, strongly correlates with developmental time in both the E_f and the E_d leg samples (PC1, Fig. 1B). A weaker source of variation was linked to specific samples (PC2, see SI Appendix, Fig. S1.1). To probe the strength of the temporal signal in those data, we applied RAPToR, a method that allows predicting the developmental stage of a sample from its gene expression profile, relative to a reference time series (22). We built a reference using the E_f leg samples and used this to estimate the stage of each sample. The predictions of the model match the real developmental age of each sample accurately, not only for the E_f samples (which were used to train the model) but also for the E_d leg samples (Fig. 1C and SI Appendix, Table S2). These results indicate that the dynamics of gene expression in developing legs are highly stereotypic, and that the temporal dynamics captured in the E_f and E_d series are highly coherent with each other.

Further comparing the transcriptional dynamics in the E_f and E_d samples, we found 7,963 and 1,354 genes to be differentially expressed during embryogenesis in E_f and E_d samples, respectively (DESeq2, padj < 0.001, in a total of 43,212 gene models; SI Appendix, Fig. S1.2), with an overlap of 1,121 genes (differentially expressed genes given in Datasets S5 and S6). We attribute the higher number of differentially expressed genes in the E_f samples to the fact that this dataset spans a longer developmental period and includes additional tissues. The tissue dissections to collect the E_d samples were also more challenging, possibly contributing to lower sample quality (SI Appendix, Fig. S1.1 and Table S1). Despite these differences, we noticed a high similarity in gene expression dynamics between these two datasets (Fig. 1D). Overall, this analysis shows that the temporal signal of the distal part of the leg (E_d) is largely recapitulated in the full leg series (E_f).

We performed a similar series of RNA-seq experiments to investigate the temporal dynamics of gene expression during the course of regeneration in adult T4 legs, amputated at the distal end of the carpus. Previous studies have shown that the cellular activity associated with regeneration occurs within 200 μm from the amputation site (19, 23), which in these experiments corresponds approximately to the distal half of the carpus. Samples were collected every 12 h, from the moment of amputation until 120 h post amputation (hpa), when the legs appear to be fully patterned (19). To ensure that we have sampled patterned and differentiated legs, we also collected samples at the onset of expression of a late distal leg marker (the Distal^DsRed exon trap, collected ∼150 hpa; 23, 24) and after the first molt following regeneration.

From each of these legs we collected two fragments (Fig. 2A): one consisting of the distal-most end of the leg stump (the carpus, including the blastema and newly regenerating structures; R_d [distal] series) and one from a more proximal podomere (the distal part of the merus; R_p [proximal] series). The R_d samples capture the entire region that participates actively in regeneration (19), while the R_p samples serve as controls, intended to capture transcriptional variations associated with the physiological status of each individual (e.g., molting stage, nutritional state). Overall, we collected 120 samples from 37 individuals (paired samples were collected from the left and right T4 legs in 23 individuals), spanning 13 time points, yielding a total of 60 R_d and 60 R_p samples (Fig. 2A, listed in SI Appendix, Table S1; Datasets S7 and S8).

Principal component analysis including both the R_d and R_p series reveals several distinct sources of variation in these data. PC1 captures the difference between the R_d (regenerating) and R_p (control) samples, with the notable exception of the 0 hpa (pre-amputation) and the postmolt R_d samples, which group together with the R_p series (Fig. 2B). This distribution shows that tissues undergoing active regeneration are transcriptionally distinct from the non-regenerating samples.

PC2 reveals marked differences between the groups of samples collected from each individual (particularly in individuals marked in bold, Fig. 2C). This variation reflects real biological differences between individuals, as we find a much higher correlation of gene expression in samples collected from the same individual than among different individuals (Fig. 2C and SI Appendix, Fig. S2.1A). These samples were processed in a randomized order, so these correlations could not arise from postprocessing batch effects. We return to the source of this inter-individual variation in the next section.

PC3 and PC4 capture a temporal signal corresponding to the progress of regeneration in the R_d samples (Fig. 2D and SI Appendix, Fig. S2.2A). On PC3, the transcriptional profile of pre-amputation samples (0 hpa) matches the profile of samples collected at the end of the regenerative time course, consistent with our expectation that regeneration is largely completed by ∼120 hpa and after the following molt.

Principal component analysis on the R_d samples alone captures the temporal signal in PC1 (Fig. 2E) and both temporal and inter-individual variation in PC2 (SI Appendix, Fig. S2.2B). Overall, our analysis reveals that regeneration and inter-individual variation are the major sources of variation in the Rd samples.

Using a reference timeline based on the R_d samples, the stage of regenerating samples (R_d) can be partly predicted based on their transcriptome (linear regression r² = 0.29) (Fig. 2F and SI Appendix, Table S3). In contrast, most of the control samples (R_p) are assigned to the fully differentiated state, confirming that these samples do not carry a substantial regenerative signal (linear regression r² ∼ 0; Fig. 2F).

We hypothesized that the observed “individual signal” (PC2 in Fig. 2C) is linked to the physiological state of each animal, as it is shared by all the samples collected from each individual. Since molting is a major physiological variable in adult crustaceans, we decided to test directly the impact that the molting cycle might have on the leg transcriptome.

Selecting animals of the same age/size as in the regeneration RNA-seq experiments, we monitored the molting status of 66 animals over two successive molts; we observed that this cohort molted with a mean period of 27 d (SD 7.2 d). We then collected entire T4 legs from 20 of these animals at different stages of the molting cycle (Fig. 3A and SI Appendix, Table S1) and performed RNA-seq on these samples (Datasets S9 and S10). Principal component analysis on these 20 samples shows that different stages of the molt cycle are well separated on PC1 and PC2, representing almost half of the total variation (large circles in Fig. 3B). Major transcriptional changes can be observed in the 5 d that precede molting (orange, brown, and yellow circles in Fig. 3B), followed by more stable transcriptional profiles post molting (blue and purple circles in Fig. 3B).

On the principal components describing the molting cycle, we projected the expression data of our regeneration time series in order to assess the molting status of each sample. We observed that most samples are associated with the intermolt phase, except those highlighted in bold in Fig. 2C, which are associated with near-molt stages (Fig. 3B). Molt-associated genes are a major driver of the interindividual variation (seen in PC2, Fig. 2C) in the regenerating leg samples (SI Appendix, Fig. S3.1).

Applying a soft clustering approach (Mfuzz) on the molting cycle dataset, we defined eight distinct sets of coexpressed genes (SI Appendix, Fig. S3.2 for clustering parameters, SI Appendix, Fig. S3.3 for cluster content, and Datasets S11 and S12 for cluster data). The samples collected shortly before molting show the largest changes in gene expression (orange, brown, and yellow phases in Fig. 3C). We identified 131 transcription factors whose expression changes during the molt cycle (SI Appendix, Fig. S3.4 and Dataset S30). These factors, which include the ecdysone receptor and other known mediators of molt responses in arthropods, are prime candidates for future studies on the interplay between molting and regeneration in Parhyale.

These analyses confirm our hypothesis that molting status has a strong transcriptional influence on the regenerating leg transcriptomes. In particular, the imminence of molting deeply modifies the transcriptional state of an adult leg (Fig. 2C).

To investigate the transcriptional dynamics driven by the regenerative process independently of the physiological/molting status of each animal, we developed a Bayesian modeling approach (using JAGS, ref. 25; model outlined in Fig. 4A) to dissect the contributions of regeneration (R, bold red in Fig. 4A) and the individual’s physiology on the the transcriptome of R_d samples (gray circle in Fig. 4A). Based on the results presented earlier, we assumed that the variation due to an individual’s physiology would be shared by all the samples collected from each individual (R_d and R_p from contralateral T4 legs). In contrast, the variation in gene expression driven by the regenerative process should be specific to each R_d sample. Previous observations suggested that individual limbs can regenerate at different speeds (19); we therefore modeled the regenerative signal separately in each R_d sample, even when we collected them from the left and right T4 legs of the same individual.

The regenerative signal R was modeled as an enrichment value, similar to a fold change between R_d and the same individual’s control/physiological signal measured in the R_p samples, taking into account sampling errors/variation (Fig. 4A; R values given in Dataset S13). An R value of 1 conveys that there is no difference in the expression of a given gene between R_d and R_p, R > 1 means that the gene is up-regulated in the regenerating sample, and R < 1 that the gene is down-regulated. Comparing the temporal profiles of R and R_d shows that R values preserve the temporal signal of regeneration but largely reduce the inter-individual variation associated with molting (SI Appendix, Fig. S4.1 A and B).

This modeling approach is successful in extracting the expression dynamics of regeneration from the overall transcriptional variation, without introducing unintended distortions or artifacts in the data (Fig. 4B and C). The principal axes of variation correlate better with regeneration time in the modeled R data compared with the R_d data (SI Appendix, Fig. S4.2), and predictions on the regenerative stage of each sample using RAPToR are also more accurate using R instead of R_d (SI Appendix, Fig. S4.3 A and C; average distance 21 versus 30 h, respectively).

A more targeted approach for removing molt-related variation is to exclude from the R_d dataset the five samples collected close to molting and the genes whose expression is significantly affected by molting (>6,000 genes). RAPToR predictions made using this alternative approach have a similar accuracy with the predictions made using Bayesian modeling on the entire dataset (SI Appendix, Fig. S4.3 B and C). Given that the targeted approach excludes five samples and >6,000 genes from the analysis, we decided to pursue our study using Bayesian modeling.

The regenerative stage of each sample is predicted more accurately in early phases of regeneration (0 to 36 hpa) than in later phases (48 to 120 hpa). This mirrors observations in live imaging experiments, in which wound closure reliably takes place in the first 1 to 2 d after amputation, but the onset of later events varies (19). Given these variations, instead of sample collection time, we decided to use the predictions made by RAPToR to place the regenerating samples on a temporal scale (pseudotime) reflecting each leg’s progress in regeneration, based on its transcriptional profile.

Having captured the transcriptional profiles of leg development and regeneration, we turned our attention to comparing these profiles, to determine the extent to which the dynamics of leg regeneration mirror those of leg development. To render the embryonic data (E_f and E_d) comparable with the modeled data from regenerating legs (R), raw counts in the embryonic datasets were transformed into enrichment values (fold changes) by applying a similar Bayesian modeling as for the R values (see Methods; SI Appendix, Fig. S4.1 D–G, and Datasets S14 and S15); we refer to these transformed data as the E and D series (based on E_f and E_d, respectively). We find that a large proportion of genes showing temporal variation during regeneration also show dynamic expression during leg development (SI Appendix, Fig. S4.4 and Dataset S16).

A combined principal component analysis on the E and the R datasets reveals that, overall, gene expression varies more during development than regeneration: the two major axes of variation (PC1 and PC2) capture well the transcriptional dynamics of leg development, but not the dynamics of leg regeneration (Fig. 5A). There is no obvious alignment between the variation seen in the embryonic and the regenerating time series.

Next, we tried to temporally align these datasets using RAPToR. We built a reference time series based on the embryonic leg (E) data and tested whether the regenerating leg (R) samples can be aligned to this reference. Pre-amputation and postmolt samples are consistently assigned to the latest stages of the leg development series, as expected of fully differentiated leg tissues (Fig. 5B and SI Appendix, Table S4). The other samples are inconsistently placed on the developmental series, some matching early and some later phases of development with no obvious pattern, suggesting that there is no straightforward relation between the phases of leg regeneration and leg development.

When using PCA and RAPToR, global expression profiles could be dominated by specific groups of genes that show strong differential expression (e.g., terminal differentiation genes), concealing relevant expression dynamics that occur on a smaller scale (e.g., in genes involved in patterning, the control of cell proliferation, and morphogenesis). To dissect the temporal dynamics of genes associated with different phases of development and regeneration, we turned to a clustering approach, which classified genes into four major coexpression clusters in the E series and eight major coexpression clusters in the R series (the same analysis was also performed on the D series; SI Appendix, Fig. S5.1 and Datasets S17–S28; cluster sizes in SI Appendix, Table S5). All the E clusters and five of the eight R clusters appear to be associated with specific phases of development or regeneration (Fig. 5C). In the embryonic leg data, genes in cluster E2 are expressed predominantly in the early phases of leg development, genes in E4 and E1 in mid phases, and genes in E3 in the late phase (Fig. 5C, Left). In the regenerating leg data, genes in cluster R4 are expressed early, R1 in early–mid phases, R8 in mid–late phases, and R6 and R2 are associated with differentiated legs (both pre- and post-regeneration). The remaining three R clusters (R3, R5, and R7) do not show a consistent temporal signal (Fig. 5C, Right).

Having identified clusters of coexpressed genes in the embryonic and regenerative time series, we systematically compared the gene content of these clusters to determine whether similar sets of genes are coexpressed during development and regeneration. We measured their overlap in terms of enrichment (fold change) relative to random expectation (see Methods). The largest overlaps are observed between the early developmental cluster E2 and the mid regenerative clusters R1 and R8, and between the late developmental and regenerative clusters (E3, R2, and R6), which are associated with differentiated legs (Fig. 5E and SI Appendix, Fig. S5.4D and S5.5). Despite this enrichment, we find that these clusters show significant differences in gene content (Fig. 5F).

Examining the expression profiles of R cluster genes during leg development and E cluster genes during regeneration (SI Appendix, Fig. S5.3) did not yield any additional insights.

A finer classification of coexpressed genes yielded 32 and 37 clusters for E and R, respectively (SI Appendix, Fig. S5.9; Datasets S18, S20, S24 and S26). Similar overlaps in gene content were observed at this finer resolution (Fig. 5G). Despite these similarities, we did not detect a well-aligned temporal sequence of expression during development and regeneration (alternative strategies for temporally aligning the gene clusters were tested; see Fig. 5G and SI Appendix, Fig. S5.9 and S5.10).

To associate the identified clusters with biological functions, we performed a gene ontology (GO) enrichment analysis on each cluster (SI Appendix, Fig. S5.7 and S5.8) and we categorized the enriched GO terms based on processes that contribute to leg development and regeneration (see Methods; Fig. 5D, and SI Appendix, Fig. S5.7). In the embryonic leg data, the most noticeable feature is the strong enrichment of the early gene cluster with GO terms associated with cell proliferation (cluster E2), followed by phases enriched in patterning and morphogenesis (cluster E4), and cell differentiation (cluster E3) (Fig. 5D, Left). In the regenerating leg, we observe an initial phase associated with stress, wounding, immune responses, and cell death (cluster R4), followed by a phase associated with cell dedifferentiation, cell proliferation, patterning, and morphogenesis (cluster R8), and then a phase associated with cell differentiation (cluster R6); a gene cluster that is associated with TOR signaling and growth is expressed throughout the process (cluster R5) (Fig. 5D, Right). This temporal sequence is in agreement with the regenerative phases identified by live imaging (19).

Overall, we observe that enriched GO terms have distinct temporal distributions in these two datasets. As expected, regeneration starts with a phase of wound healing (cluster R4) that is not represented in embryonic leg development. But notably, embryonic legs express genes associated with cell proliferation and patterning/morphogenesis in distinct phases (clusters E2 and E4), whereas in leg regeneration these processes occur simultaneously (cluster R8).

Taking a more targeted approach, we also examined the temporal profiles of specific genes that are likely to be involved in immunity, cell proliferation, leg patterning, and cell differentiation (Fig. 6; gene list in SI Appendix, Table S6). The genes associated with immunity, cell proliferation, and patterning were selected based on published information (particularly on Drosophila orthologs, see SI Appendix, Supplementary Methods) and on the GO term analysis; genes associated with differentiated neurons and muscles were identified from Parhyale single-cell transcriptomic data (17).

Immunity-associated genes are markedly up-regulated during the early phases of regeneration, following wounding (Fig. 6A and SI Appendix, Fig. S6.4B). The same genes are expressed only in the later phases of embryonic development, possibly connected to the differentiation of circulating hemocytes (Fig. 6A and SI Appendix, Fig. S6.1A and S6.3).

The expression profiles of genes expressed in the G2/M phase of the cell cycle indicate that cell proliferation occurs mainly during the early phase of leg embryogenesis, and the mid–late phases of regeneration (Fig. 6B and SI Appendix, Fig. S6.1B, S6.3, S6.4). This is consistent with data from live imaging (19) (SI Appendix, Fig. S6.5) and with the GO enrichment analysis (Fig. 5D).

Leg patterning genes were predominantly expressed in the mid–late phases of leg development (from 120 hpf onward), after the down-regulation of genes associated with cell proliferation (Fig. 6C and SI Appendix, Fig. S6.1, S6.2, and S6.3). During leg regeneration, this set of genes is predominantly expressed during later phases, overlapping with the expression of genes associated with cell proliferation (Fig. 6C and SI Appendix, Fig. S6.2 and S6.4).

Finally, genes associated with differentiated muscles and nerve cells are strongly expressed during the late phases of leg development in the embryo and in fully differentiated adult legs before and after regeneration (Fig. 6D and E and SI Appendix, Fig. S6.1, S6.3 and S6.4). They are down-regulated during the course of regeneration, possibly reflecting neuron and muscle dedifferentiation during these stages.

Computations were performed using JAGS via the R package rjags (25). Details on modeling are provided in SI Appendix, Methods. The resulting R, E, and D transformed values are provided in Datasets S13–S15 (values have been further log transformed).

In order to infer the progression of regeneration or development based on transcriptomic data, we used the R package RAPToR v1.1.4 (Real Age Prediction from Transcriptome staging on Reference), a recently developed tool to accurately predict individual samples’ developmental age from their gene expression profiles (22).

For building RAPToR references, we used the function ge_im with the formula “X ∼ s(hpa, bs = ‘ts’)” and parameter dim_red=“pca”. For the RAPToR reference based on E_f values for Fig. 1C, we used the normalized and log-transformed tpm values and selected the genes variable in the E_f samples (the top 20,000 variable genes as calculated by DESeq2), excluding the genes with a low expression (75th percentile count above 10, with a final set of 16,199 genes) and 30 principal components (values in SI Appendix, Table S2). To build the RAPToR reference based on R_d values for Fig. 2F, we used the top 20,000 most variable genes in samples R_d as calculated by DESeq2, excluding the genes with a low expression (75th percentile count above 10, with a final set of 12,438 genes). As regeneration is less synchronous than embryonic development (19), we were concerned that the sampling timing would not faithfully reflect regeneration progression and that our samples were an imperfect reference. To avoid overfitting, we made two changes to the standard RAPToR protocol for building a reference: We used only three PCs, and for each sample we built a separate reference excluding the sample being tested. Pre-amputation samples were given an arbitrary timing of 300 hpa that would place them at the end of the regenerative sequence, and postmolt samples were excluded from the reference (values are provided in SI Appendix, Table S3). To build the RAPToR references based on R values for Figs. 4C and 5B and SI Appendix, Figs. S5.2–4, S5.9, S5.10, S6.1, and S6.2, we used the 20,000 top genes with the most variable R values, further excluding the genes with a low expression (75th percentile count above 10, for a final set of 12,434 genes). To avoid overfitting, we followed the same procedure described above (values in SI Appendix, Table S7 for R samples, SI Appendix, Table S8 for E samples, and SI Appendix, Table S9 for D samples). For the RAPToR reference based on E values for Fig. 5B, we used a gene list that was as exhaustive as possible, consisting of the union of the 20,000 top genes with the most variable R, E, or D values, and excluding the genes with a low expression (75th percentile count above 10, for a final set of 16,759 genes) and 10 PCs (values in SI Appendix, Table S4). For the RAPToR reference based on E values for Figs. 5C, 5.9C, and 6 and SI Appendix, Table S6.2A, we used the 20,000 top genes with the most variable E values, excluding the genes with a low expression (75th percentile count above 10, yielding a final set of 14,632 genes), and 3 PCs. Correlation between predicted time of amputation and real time of amputation (Figs. 2F and 4C) was computed excluding postmolt and t₀ samples, the real timing of the former being too uncertain and the later being considered as an end point to regeneration.