Local mechanical stimuli correlate with tissue growth in axolotl salamander joint morphogenesis

(a) Axolotl experiments

Larval animals (3–5 cm) were bilaterally amputated just proximal to the elbow joint. GSK1016790A was reconstituted in dimethyl sulfoxide (DMSO) and injected intraperitoneally at 50 μg kg⁻¹ at 22 days post amputation (dpa, n = 6). Control animals (n = 6) were injected with 50 μg kg⁻¹ DMSO. Injections were repeated at 48 h intervals. At 32 dpa, all animals were injected intraperitoneally with 5-Ethynyl-2’-deoxyuridine (EdU) and L-Azidohomoalanine (AHA). Limbs were collected 18 h later, fixed and stained.

We imaged nascent macromolecule synthesis in the regenerated forelimbs with light sheet fluorescence microscopy following the whole-mount click-it-based technique in Duerr et al. [47]. EdU is incorporated into newly synthesized DNA, which allowed for the quantification of cell proliferation through EdU-positive nuclei segmentation. AHA enabled the visualization of chondrocyte protein translation, most likely ECM, which provided a well-defined outline of the bone rudiment’s perichondria. Quantification of the three-dimensional shape was then possible through the analysis of the humerus outline (figure 1a).

• Download figure
• Open in new tab
• Download PowerPoint

Electronic supplementary material, figure S1A illustrates the timeline of the experiments and electronic supplementary material, figure S1B shows an example of the animal size used, and the location of the amputation. Injections started at 22 dpa, which is roughly when joint cavitation occurs in regenerating limbs in 3–5 cm sized animals, and continued throughout the joint morphogenesis stage of the joint formation process until 30 dpa. It is difficult to ascertain complete penetration of GSK101 via intraperitoneal delivery. Treated animals showed a systemic response to injections as they appeared debilitated and lethargic, though there was no visible difference in size of the animals or limbs.

Electronic supplementary material, figure S1C shows a central slice of a three-dimensional image stack obtained for an exemplary control elbow at 33 dpa. No ossification was observed in the fully regrown limbs; all bone rudiments were cartilaginous at this final stage. All light sheet images were acquired using a Zeiss light sheet Z.1 microscope paired with Zen software.

(b) Experimental data analysis

We segmented the regrown bone rudiments (33 dpa), identified the proximo-distal longitudinal axis of the humerus and ulna through computation of the minimum principal axis using the Fiji plugin BoneJ [48], and then aligned all limbs in three-dimensional space (figure 1b). The alignment process included mirroring of right limbs so that all limbs had the dorsal and ventral condyles in the same relative position in space. A cylinder was fitted to the aligned humerus surface using the Matlab [49] File Exchange function ‘cylinderfit’ (a regression modelling tool), and a hemispherical cap was placed on top to create the reference surface (figure 1c). These were shifted vertically upwards until the hemispherical cap was tangent to the distal end of the humerus surface. The distance from the reference surface to the humerus surface was mapped onto the reference surface and normalized using the fitted cylinder diameter (figure 1d) to account for animals of different size. We quantified dorsal and ventral condyle shapes and sizes based on the corresponding normalized areas and normalized volumes, respectively, which were extracted from the two-dimensional standardized representation of the humerus surface (figure 1e).

To quantify proliferating cells, we manually generated a small training set to train the deep learning algorithm Stardist3D [50], which was used to identify the EdU-stained cell nuclei in the three-dimensional image stack. The Fiji plugin three-dimensional Objects Counter [51] was used on the cell nuclei masks produced by Stardist3D to identify proliferating cell positions and volumes. Outliers were removed based on cell volume and we used a fixed-length cut-off to ensure quantification of cell proliferation was performed in an equivalent humerus volume across different limbs.

Electronic supplementary material, figure S2 provides a visual summary of the complete workflow, which was implemented using a combination of Fiji [52], the ZeroCostDL4Mic implementation of Stardist3D [53] and a customized code in Matlab [49]. We grouped all limb results for each measurement and ran a Shapiro–Wilk normality test. Except for the proliferating cell count, all other data measurements were normally distributed. A one-way ANOVA was used to check for statistically significant differences between the control (n = 10) and GSK101 (n = 11) groups of normally distributed data. The proliferating cell count p-value was obtained using a Kruskal–Wallis test.

(c) Experimental results

Our results reveal significant differences in cell proliferation and bone rudiment shape between the humeri of the control group and the GSK101 group. A central slice of a representative EdU-stained humerus from the control and GSK101 groups (figure 2a) and the three-dimensional distribution of the cell nuclei identified for each (figure 2b) depict the differences in cell proliferation between groups. A central slice of each humerus analysed is provided in the electronic supplementary material, figure S3. The mean value of the proliferating cell count in the humeri of the control group is fourfold that of the GSK101 group (p-value < 0.001; figure 2c). However, the diameter of the humeri shaft is similar for both groups (figure 2d). To ensure this is not owing to an insufficiently sensitive measurement method, we computed additional metrics using an alternative methodology (electronic supplementary material, S4). All measures of humeri shaft size computed indicate there is no significant differences between the two groups.

• Download figure
• Open in new tab
• Download PowerPoint

The normalized volumes of both dorsal and ventral condyles are larger for the control group than the GSK101 group (p-value = 0.005 and <0.001, respectively). The normalized areas of the ventral condyles in the control group are also larger (p-value < 0.001), while no significant difference was found for the normalized areas of the dorsal condyles (figure 2e). These findings are visually reflected in the mean humerus surfaces computed (figure 2f), where a darker shade of red in the condyles of the mean control humerus indicates more prominent condyles for this group with respect to the GSK101 group. To compute the mean two-dimensional surface maps, we aligned the individual two-dimensional surface maps (electronic supplementary material, figure S5) based on the position of the dorsal condyle centroids.

(a) Modelling cartilage tissue growth within a poroelastic framework

The biphasic approach defines tissue as a mixture of an elastic solid skeleton with free-flowing fluid circulating within its pores. In cartilage, the fluid can be assimilated to the interstitial fluid in the tissue, i.e. water and dissolved ions, growth factors and other molecular components. The solid component represents the proteoglycans and collagen of the ECM and chondrocytes. Chondrocyte proliferation, hypertrophy, migration and/or intercalation as well as ECM production in cartilage can then be modelled together at tissue level through continuum growth of this solid phase (electronic supplementary material, figure S6). Following a common approach in the field [33–35,59], we consider growth rate to be a sum of biological and mechanical contributions.

The biological contribution represents the intrinsic morphogenetic biological factors that globally mediate tissue growth. Similar to past studies of joint morphogenesis [34,35], we assumed it is proportional to chondrocyte density in the bone rudiments. However, unlike these studies, our experimental measurements of chondrocyte density in a regenerating axolotl humerus revealed an approximately constant value throughout the bone rudiment at this stage of regrowth (electronic supplementary material, figure S7). Chondrocyte density quantification of the humeri in our experiments was not possible based on the AHA staining (electronic supplementary material, figure S9). We assumed a constant biological growth rate in time and space, within the humerus geometry and throughout the whole simulation time period in all our simulations. This implies cell density is the same in both groups.

The mechanical contribution is a function of the selected mechanical stimulus locally driving tissue growth. Mechanical loading is known to modulate the synthesis of ECM in chondrocytes. Collagen and aggrecan production, the main components of ECM in cartilage, depends on the magnitude, duration and type of loading. In particular, in vitro experiments have shown that cyclic compression promotes ECM production while static loading either has no effect on collagen and aggrecan levels, or inhibits cartilage growth [38–42]. Our poroelastic model is able to capture the differences between static and dynamic loading by defining mechanical growth proportional to a dynamic variable linked to the movement-induced fluid flow. We selected pore pressure of the fluid component, a hydrostatic measure akin to the hydrostatic stress used in past models, as the mechanical stimulus.

The discretized governing equations and continuum growth model were implemented in the open source finite element library deal.II [60]. The code used in this study is an extension of the poro-viscoelastic numerical framework in [61]. Growth was implemented following the algorithm in the electronic supplementary material, figure S8. Further details of the poroelastic formulation, the growth model and their numerical implementation are provided in the electronic supplementary material, S6.

(b) A finite element model of joint morphogenesis

We generated a finite element model of a generic humerus bone rudiment after cavitation, i.e. at the start of the experiments, with the goal of predicting the grown humerus shape at the end of the joint morphogenesis stage. Given that our model is a tool to probe potential mechanisms of load mechanotransduction in joint morphogenesis, we strove to keep its parameters as generic as possible.

The geometry and loading conditions (figure 3a) were informed by experimental data. A normally regenerating forelimb at 17 dpa in a 3 cm sized animal, which corresponds to the time point just after joint cavitation, was used (electronic supplementary material, figure S10A). We segmented the bone rudiment shapes from the three-dimensional image stack (electronic supplementary material, figure S10B).

• Download figure
• Open in new tab
• Download PowerPoint

A mesh was generated based on the smoothed-out surface of the segmented humerus with a total of 512 hexahedral elements. We scaled the geometry size to achieve a cross-sectional humerus size closer to the values identified in our experiments. Meshing of the geometry inevitably entails a slight loss of surface detail. We computed and visually compared the two-dimensional surface maps of both the segmented geometry and the meshed geometry (electronic supplementary material, figure S10C) following a procedure analogous to the one used in the humerus three-dimensional shape analysis. Comparison of two-dimensional surface maps confirmed that the meshed surface retained the main characteristics of the original humerus.

Free-flow boundary conditions across all external surfaces except the bottom (proximal) one were set in the finite element model. Vertical displacements of the bottom surface were fixed, and lateral displacements of nodes in the bottom surface were fixed (electronic supplementary material, figure S10D). These boundary conditions allowed for outward growth of the humerus shaft while avoiding spurious translations as well as the rotation of the whole bone rudiment.

The loading conditions applied (figure 3a; electronic supplementary material, figure S10E), modelled a 1 s flexion-extension cycle of the elbow. The growth resulting from a single cycle was extrapolated for multiple cycles. Loading was applied as a pressure over a roughly circular surface representing the contact areas between the radius/ulna and the humerus. A sine-like loading profile over this area was considered, with the loading area sweeping over the humerus surface. The sweep path was estimated based on anatomical observations of the axolotl elbow joint. The value of the load profile changed throughout the cycle to mimic the effect of muscle contractions, reaching the maximum value for the peak flexion position. Load step increments of 0.01 s were applied. We studied the effect of varying loading and boundary conditions on our computational results (electronic supplementary material, figure S11). In this way, we ensured the robustness of our computational set-up to produce results from which to extract meaningful insights.

The material properties were either estimated from the literature or based on an educated guess, except for the initial intrinsic permeability of the biphasic material. Preliminary simulations identified this parameter as having a considerable impact on the predicted patterns, while the rest of the material properties did not substantially alter the predicted growth patterns. Hence, we adjusted the value of the initial intrinsic permeability based on experimental stress-relaxation data obtained through nanoindentation tests on an axolotl forelimb (electronic supplementary material, figure S12). Electronic supplementary material, S9 provides further details of all model parameters.

(c) Computational results

A regenerating humerus model based on local changes in fluid pressure (figure 3b) induced by an elbow flexion-extension loading cycle predicted a final humerus morphology (figure 3d, left) that resembled our experimental observations of the control group. When the mechanically driven growth component (figure 3c) was removed, shape prediction (figure 3d, right) resulted from constant volumetric biological growth only and was in accordance with the experimental observations of the GSK101 group. When the mechanically driven growth component was included, local mechanical growth occurred in regions of high compressive pressure, which were observed underneath the surface load representing the radius contact area throughout the cycle. However, we did not observe an analogous pressure below the load representing the ulna contact area. Pressure was most pronounced at the posterior proximal part of the humerus shaft. Complete predicted patterns for pressure as well as other mechanical stimuli that were initially considered as potential drivers of the mechanical growth model are provided in the electronic supplementary material, S10.

To quantify the differences between the healthy and mechanosensitively impaired cases, we computed at each surface node the magnitude of the distance between the original surface and the grown surface, and normalized this measure with the maximum value of the two cases. We then mapped the resulting patterns onto a reference surface, fitted to the original surface mesh, and flattened it to obtain a two-dimensional representation of normalized growth (figure 3e). The mapping procedure was analogous to the one used to obtain the two-dimensional surface maps of the experimental humeri (figure 1c–e). In both predictions, humerus surface growth increased towards the distal portion of the bone rudiment, but the healthy growth case resulted in larger values as well as a notably asymmetrical pattern. A larger surface growth was predicted in the area corresponding to the future ventral condyle for the healthy growth case (figure 3e, left). The contour of the condyles from the corresponding mean experimental surfaces in figure 2f is shown on the two-dimensional maps.

(a) Transient receptor potential vanilloid 4 desensitization during joint morphogenesis altered final humerus shape

Our analysis of the regenerating axolotl limbs revealed an altered humerus morphology for the GSK101 group (figure 2c–e). The mean two-dimensional surface maps computed for each group (figure 2f) illustrate the main findings: the condyles of the control animals have larger normalized volumes than the GSK101 group (darker shade of red in contour map). The shape of the ventral condyle, as measured based on the normalized area, is more affected by TRPV4 desensitization than the dorsal condyle.

We also analysed the shapes and sizes of the anterior and posterior concavities (blue regions in the two-dimensional surface maps; electronic supplementary material, figure S5) following an analogous procedure to the condyle measurements and did not observe significant differences between groups for any measurement. It could signify that these shape characteristics of the humerus were already present at the onset of the experiment; starting treatment with GSK101 sooner after amputation may result in more severe changes. We analysed a regenerating limb after cavitation (start of our experiments) for the purposes of developing the initial computational model. The two-dimensional surface map obtained (electronic supplementary material, figure S10C, top) supports the notion that the basic humerus shape could be present already at this stage. The concavities seem to already be present and the dorsal condyle is clearly defined, similar in shape to those of the fully regenerated limbs in the experiments (electronic supplementary material, figure S5). However, the ventral condyle is barely discernible after cavitation. This implies that the concavities and dorsal condyle may form in the earlier stages of the joint formation process, which is probably why we found little change in their shapes.

Taken together, this data indicates that TRPV4 desensitization during joint morphogenesis in regrowing forelimbs alters the final humerus shape. During regeneration, blastema cells dedifferentiate and may assume a different role in the regrowing limb. These cells have been shown to retain distinct roles in axolotls [62], which could affect the way chondrocytes in regenerating limbs respond to mechanical stimuli. However, joint morphogenesis occurs at a much later time point than the dedifferentation process within the limb regrowth timeline. It seems reasonable to assume that mechanotransduction pathways in joint formation of regenerating axolotl limbs are probably the same as those in developing limbs.

Numerous studies have shown that chondrocytes have several separate but overlapping mechanotransduction pathways [28,30]. Other channels of the TRP family have been suggested to have load-associated effects in cartilage [63], but TRPV4 is undoubtedly the major regulator of mechanical and osmotic signal transduction in this family. The Piezo1 and Piezo2 channels have also been identified as key stretch-induced mechanotransducers in chondrocytes [64]. It would be interesting to see whether altering these other channels has effects on morphology similar to those seen in this study, to further explore the interrelated roles of each channel in cartilage mechanotransduction. Furthermore, axolotls have much larger cell sizes and longer cell cycles than the vast majority of vertebrates, which probably influences their mechanosensitive response, and would also be a fascinating topic for further study.

Alternative ways of blocking mechanics in developing joints have been used in the past to study the effect of mechanical stimuli on joint morphogenesis, namely muscle paralysis in chicks [6–10] and genetically modified altered-muscle mice [11–13]. These studies also revealed morphological differences. Here, we used a TRPV4 agonist, which represents the clinical genetic deficits associated with abnormal skeletal development [65,66]. Our three-dimensional analysis of the humerus surface allows the assessment of shape changes that are not evident in more simple measures used in the past, such as cross-sectional outlines or linear anatomic measurements like humeral head width.

(b) More prominent condyles and increased chondrocyte proliferation were not associated with larger humeri

The substantial reduction in cell proliferation of the GSK101 group (figure 2a–c; electronic supplementary material, figure S3) did not result in smaller humeri sizes (figure 2d). Axolotls have long cell cycles, which have been recorded to be up to 88 hours in regenerating tissues [67,68]. Throughout the 10-day experimental treatment, few complete cycles would have occurred. Also, proliferating cells were only a relatively small percentage of the total chondrocytes in the bone rudiment. Therefore, the total amount of cell proliferation may not have been sufficient to produce actual changes in bone rudiment size. In addition, our quantification of cell proliferation corresponded to an 18 h window at the end of the experiment, which may not be representative of the complete treatment period of 10 days.

The decrease in condyle normalized volumes and in the ventral condyle normalized area for the GSK101 group may be owing to matrix production instead of cell proliferation. Proliferation was not localized to the condyles, rather it was homogeneously distributed. Our data seem to indicate that TRPV4-mediated proliferation is unlikely to be a major contributor to growth during axolotl joint morphogenesis in regenerating forelimbs.

(c) Local fluid pressure may promote tissue growth during joint morphogenesis

To link the experimentally observed changes in humerus morphology owing to TRPV4 desensitization with impaired mechanosensitivity in the growing tissue, we built a computational model of joint morphogenesis. Through hypotheses and simplifying assumptions, we have isolated a potential contributor to the mechanotransduction of mechanical loading into local tissue growth and subsequent shaping of the joint.

The computational results show that compressive fluid pressure can predict humerus morphology during joint morphogenesis. In the predicted normalized surface growth map for the healthy growth case (figure 3e, left) the ventral condyle exhibited a considerably larger amount of growth than the dorsal condyle, while the mechanosensitively impaired case (figure 3e, right) showed similar (smaller) growth values for both condyles. This agrees with the larger normalized area observed in the ventral condyle of the experimental control group with respect to the GSK101 group (figure 2f). The predictions for the healthy growth case (figure 3e, left) exhibited more growth towards the distal area than the mechanosensitively impaired one (figure 3e, right), which only had a slight gradient in the proximo-distal direction. Experiments also showed more growth (larger normalized volume) in both condyles of the control group with respect to the GSK101 group (figure 2e, left).

Certainly, our model points to a relationship between the fluid pressure distribution and the shaping of the joint. Chondrocytes might not be sensing interstitial hydrostatic pressure directly, but rather a different biophysical factor related to it. Osmotic stresses have been repeatedly identified as the stimuli triggering a series of signalling events in relation to the TRPV4 channel, that are propagated into changes in gene expression and ECM synthesis. Yet, studies have shown that osmotic loading as well as mechanical loading elicit responses of the TRPV4 channel [24,28,29,31]. Recent publications suggest TRPV4 is a cell volume sensor and is activated regardless of the molecular mechanism underlying said volume change [69]. Furthermore, hydrostatic and osmotic pressures have similar effects on cartilage formation [70], and they both affect intracellular ion signalling in chondrocytes [71,72]. It is not within the scope of this study to determine the complex interrelations between the osmotic and hydrostatic pressures induced by mechanical loading on cartilage. Many studies have shown that hydrostatic pressure increases cartilage gene expression and extracellular matrix formation (see review in [73]). Our computational results indicate that fluid pressure can predict local tissue growth in the experimentally informed model of joint morphogenesis developed in this study.

(d) Poroelasticity can be used to explore how dynamic loading dictates bone rudiment morphology

Owing to the nature of the poroelastic tissue, compressive dynamic loading generates the non-homogeneous fluid pressure pattern within the humerus that dictates tissue growth in our computational model (figure 3b). By contrast, static loading generates an initial pressure distribution that quickly dissipates as fluid seeps out of the bone rudiment (electronic supplementary material, figure S12C). Such behaviour is in agreement with experimental studies showing that cartilage growth is promoted by repetitive compressive loading while static loading inhibits cartilage growth [38,–42]. Unlike our previous models of joint morphogenesis [35,36], we are now able to inherently capture the effect because of the type of loading imposed owing to the biphasic approach that incorporates the fluid flow component into the modelling. An earlier computational study [8] used poroelasticity to relate local patterns of biophysical stimuli to the emergence of joint shape in a model of a chick knee, but could not predict growth morphologies. Through the solid component growth, our model goes a step further and can more confidently relate local tissue growth to final bone rudiment morphology based on cyclic loading-induced mechanical stimuli.

We explored alternatives to the compressive pore pressure as mechanical stimuli for our growth model (electronic supplementary material, S10), including measures of solid compression and pore fluid velocity. The positive divergence of the seepage velocity (electronic supplementary material, figure S14A, top row) stood out because, unlike the other measures, its distribution within the humerus is quite different from the fluid pressure pattern. Hence, we implemented this measure of the rate of solid compression as an alternative mechanical growth stimulus in our formulation. The resulting local tissue growth owing to the mechanical contribution was distributed more evenly towards the distal part of the humerus (electronic supplementary material, figure S15A), instead of being localized below the radius contact loading (figure 3c). In addition, less growth was observed in the proximal part of the humerus for the alternative model. Interestingly, this produced an apparent rotation of the humerus grown surface (electronic supplementary material, figure S15B) rather than the slight bending and outward growth observed broadly around the ventral condyle region for the pressure-based mechanical growth (figure 3d, left). Further study would be required to ensure artefacts owing to inadequate loading or boundary conditions are not at play here before discarding the rate of tissue compression as a potential biophysical stimuli within the joint morphogenesis process.

These exploratory simulations demonstrate the potential of the proposed model as a tool to unravel the mechanisms at play in the shaping of the joint. Through the computational study of how different measures of pressure, compression and fluid flow evolve in response to loading set-ups representative of in vivo conditions, we could identify potential biophysical stimuli for further study in experiments.