INTRODUCTION
Congenital vertebral malformation (CVM) is a 3-dimensional (3D) deformity of the vertebral body that can lead to congenital scoliosis, congenital kyphosis, and congenital kyphoscoliosis. Severe cases of CVM may even result in cardiopulmonary dysfunction [1]. The incidence of CVM is estimated to range from 0.13 to 0.50 per 1,000 live births [2]. The etiology of CVM involves both genetic and environmental factors [3].
From a genetic perspective, various genes play critical roles in skeletal system development. For example, hemivertebrae, a type of CVM, is characteristic of TBX6-associated congenital scoliosis [4]. Additionally, Mesp2 and Tbx6 are involved in the periodic patterning and segmentation clock mechanisms during early skeletal development in mice [5]. Runx2 also plays a key role in osteoblast maturation by regulating the expression of bone matrix protein genes and bone gamma carboxyglutamate protein, which may be associated with CVM [6]. However, genetic factors contributing to CVM exhibit significant heterogeneity. From an environmental perspective, vitamin A deficiency (VAD) has been linked to congenital spinal deformities in rats [7]. VAD is a severe and prevalent public health issue in developing countries, particularly among pregnant women and children [8]. As a vital nutrient during embryonic development, maternal VAD during pregnancy can lead to congenital abnormalities or even fetal death [9]. Furthermore, VAD has been shown to affect the development of other organs, including the kidneys, genitourinary tract, diaphragm, aortic arch, lungs, and heart [10]. In addition to VAD, other environmental factors such as hypoxia, cigarette smoking, alcohol consumption, valproic acid, boric acid, and hyperthermia have also been reported to be associated with CVM [11].
The influence of vitamin A on embryo development is mediated by enzymes that convert vitamin A first to retinaldehyde and then to retinoic acid (RA). RA functions as a ligand for 2 families of nuclear receptors, the retinoic acid receptors (RARs) and the retinoid X receptors (RXRs), which bind DNA and directly regulate transcription [12] When RA binds to the RAR partner of RAR/RXR heterodimers, which are attached to a regulatory DNA element, it triggers a cascade of events that recruit transcriptional coactivators and initiate transcription [13]. Kawakami et al. [14] demonstrated that RA plays a crucial role in buffering the influence of laterality information flow on the left-right progression of somite formation, thereby ensuring the bilateral symmetry of somitogenesis. Further research suggests that RA antagonizes Fgf8 expression in the ectoderm, and a failure in this mechanism leads to excessive FGF8 signaling to adjacent mesoderm, initially resulting in smaller somite and subsequent left-right asymmetry [15]. These findings highlight the vital role of RA signaling in somitogenesis and suggest that VAD might disrupt RA synthesis. Since somite differentiate into sclerotomes, which eventually develop into vertebral bodies, we hypothesize that VAD may lead to CVM through the disturbance of RA signaling.
Although Li et al. [7] identified that VAD induces congenital spinal deformities in rats using x-ray imaging, the images were relatively vague, making it difficult to distinguish the specific type of spinal deformity. To address this, we combined microcomputed tomography (CT) and x-ray imaging to construct a 3D model of VAD-induced spinal deformities in rats. We then applied whole mount in situ hybridization (WMISH) to examine the expression patterns of genes related to RA signaling and sclerotome differentiation in embryos at various gestational stages. Following this, we utilized laser capture microdissection (LCM) in combination with RNA sequencing (RNA-seq) to analyze the transcriptional abnormalities in the sclerotome of gestational day (GD) 12.5 rat embryos caused by VAD. Finally, we validated the RNA-seq results using real-time quantitative polymerase chain reaction (RT-qPCR).
MATERIALS AND METHODS
1. Sex as a Biological Variable
Our study examined male and female animals, and sex-dimorphic effects are reported (Supplementary Fig. 1).
2. Laboratory Animal Use
The animal work was taken place in Laboratory Animal Research Facility, National Infrastructures for Translational Medicine, Institute of Clinical Medicine, Peking Union Medical College Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College. We obtained ethics approval from the Animal Welfare and Ethics Committee of Peking Union Medical College Hospital (approval number: No. XHDW-2023-003). Anesthesia of Sprague-Dawley (SD) rats were induced by isoflurane gas. SD rats were euthanized through carbon dioxide inhalation.
3. Animal Model
We obtained 24 ten-week-old virgin female SD rats and 24 ten-week-old male SD rats from SPF Biotechnology Co., Ltd (Beijing, China). We randomly divided the female rats into 2 groups: the VAD group and the control (CON) group. The VAD group (n=12) received an AIN-93G growing rodent diet without added vitamin A (D13110GC, Research Diets, USA) (Supplementary Table 1). The CON group (n=12) was fed an AIN-93G growing rodent diet (D10012GM, Research Diets) (Supplementary Table 2). After 1 week of acclimatization, we fed the rats their respective diets for 2 weeks. We then mated the rats with normal males between 6 PM and 10 PM. During gestation, the VAD group continued to receive the AIN-93G growing rodent diet without added vitamin A, while the CON group received the AIN-93G growing rodent diet. Each maternal rat has 10–15 embryos or neonatal rats. We fed the neonatal rats from 4 randomly selected gestation rats in the VAD group and 4 randomly selected gestation rats in the CON group different diets until they were 2 weeks old for x-ray and micro-CT scanning. We collected GD10.5 and GD12.5 embryos from 4 randomly selected gestation rats in the VAD group and 4 randomly selected gestation rats in the CON group for WMISH. Additionally, we collected GD12.5 embryos from 2 randomly selected gestation rats in the VAD group and 2 randomly selected gestation rats in the CON group for LCM combined RNA-seq and RT-qPCR.
4. Serum Vitamin A Measurement
We obtained 12 ten-week-old female SD rats from SPF Biotechnology Co., Ltd (Beijing, China). We randomly assigned the rats to either the VAD group or the CON group. The VAD group (n=6) received an AIN-93G diet with no added vitamin A, while the CON group (n=6) was fed an AIN-93G growing rodent diet. After 1 week of acclimatization, we started feeding the rats their respective diets. 14 days after the diet change, we collected orbital blood from the rats to measure serum retinol levels. We collected blood samples in EDTA-coated tubes. All serum vitamin A measurements were conducted within 2 weeks of sample collection. We obtained the serum by centrifuging the blood at 1,200×g for 15 minutes at 4°C. To minimize photoisomerization of vitamin A, we handled the plasma under reduced yellow light and stored it in the dark at -80°C until we determined the vitamin A concentrations. We measured serum vitamin A concentration using Liquid Chromatography Mass Spectrometry (ACQUITY UPLC BEH C18 Column, 1.7 μm, 2.1 mm×100 mm).
5. X-Ray Assessment
We anesthetized SD rats for x-ray spinal alignment evaluation. We positioned each rat in a standard imaging setup. We obtained both posteroanterior and lateral radiographs using a digital x-ray photography system (uDR588i, United Imaging, China). The x-ray settings were 2 seconds exposure, 800 mA, and 103 kV.
6. Micro-CT Assessment
We sacrificed 2-week-old rats and placed them in an upright posture, securing them in the scanner unit with their head, bilateral shoulders, and bilateral lower extremities in a neutral position. We scanned spinal alignment using a μCT 45 desktop micro-CT scanner (SCANCO Medical, Switzerland). The scanning parameters were as follows: source voltage, 70 kV; source current, 114 μA; voxel size, 34 μm. We performed 3D reconstruction of the spine using RadiAnt DICOM Viewer (64-bit).
7. Whole Mount In Situ Hybridization
We performed WMISH of rat embryos following the protocol outlined in the “Whole mount in situ hybridization protocol for mRNA detection” (http://geisha.arizona.edu/geisha/protocols.jsp). We transcribed digoxigenin-labeled antisense RNA probes in vitro using T7 or Sp6 RNA polymerases, as described in the manufacturer’s manuals. The probes used targeted rat Pax1 (490-bp fragment) and rat Raldh-2 (447-bp fragment). Sequences of the primers used for in vitro transcription of Pax1 and Raldh-2 probes are listed in the supplementary material (Supplementary Table 3). We examined stained embryos using a Leica M205 FA Fluorescence stereo microscope.
8. LCM-associated RNA-seq
We performed LCM-associated RNA-seq following a published protocol [16]. We collected GD12.5 SD rat embryos, embedded them in OCT (SAKURA Tissue-Tek O.C.T.), and cryosectioned them serially from the distal to the proximal region at a thickness of 20 μm. We mounted the sections on MMI MembraneSlides, fixed them immediately with 100% ethanol, 95% ethanol, 70% ethanol, 90% ethanol, and 100% ethanol (Sigma-Aldrich, USA). We captured sclerotome samples with a 150-μm diameter circle using LCM (MMI Cellcut Plus system) on transverse sections (Supplementary Fig. 2). We collected the samples on IsolationCaps and lysed them using Guanidine isothiocyanate (GuSCN, Invitrogen, cat. no. 15577-018). After extracting and denaturing the RNA, we obtained cDNA through reverse transcription. We performed cDNA preamplification and bead purification to create a purified cDNA library. We constructed the DNA library using the TruePrep DNA Library Prep Kit V2 for Illumina (TD503, Vazyme). Sequencing was conducted using a NovaSeq 6000 sequencer (Illumina, USA).
9. Data Analysis of RNA-seq
1) Quality control
We initially processed the raw data (raw reads) in FASTQ format using in-house Perl scripts. In this step, we obtained clean data (clean reads) by removing reads containing adapters, reads with poly-N sequences, and low-quality reads from the raw data. Concurrently, we calculated the Q20, Q30, and GC content of the clean data. All subsequent analyses were based on this highquality clean data.
2) Reads mapping to the reference genome
We downloaded the reference genome and gene model annotation files directly from the genome website. We built an index of the reference genome using Hisat2 v2.0.5 and aligned the paired-end clean reads to the reference genome using Hisat2 v2.0.5. We chose Hisat2 for mapping because it can generate a database of splice junctions based on the gene model annotation file, providing superior mapping results compared to nonsplice mapping tools.
3) Quantification of gene expression level
We used FeatureCounts v1.5.0-p3 to count the number of reads mapped to each gene. We then calculated the Fragments Per Kilobase of transcript per Million (FPKM) base pairs sequenced for each gene based on the gene length and the number of reads mapped to the gene. FPKM accounts for sequencing depth and gene length, making it the most widely used method for estimating gene expression levels.
4) Differential expression analysis
We performed differential expression analysis between the VAD and CON groups using the DESeq2 R package (1.20.0). DESeq2 provides statistical routines for identifying differential expressions in digital gene expression data, utilizing a model based on the negative binomial distribution. We adjusted the resulting p-values using the Benjamini and Hochberg approach to control the false discovery rate. Genes with a p-value <0.05, as identified by DESeq2, were classified as differentially expressed.
5) Enrichment analysis of differentially expressed genes
We conducted gene ontology (GO) enrichment analysis of the differentially expressed genes using the clusterProfiler R package, which corrects for gene length bias. GO terms with a p-value less than 0.05 were considered significantly enriched in the differentially expressed genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) provides a database resource for understanding high-level functions and utilities of biological systems, such as cells, organisms, and ecosystems, based on molecular-level information. This includes large-scale molecular datasets generated by genome sequencing and other high throughput experimental technologies (http://www.genome.jp/kegg/). We used the clusterProfiler R package to assess the statistical enrichment of differentially expressed genes in KEGG pathways. KEGG terms with a p-value <0.05 were considered significantly enriched.
6) Gene set enrichment analysis
Gene set enrichment analysis (GSEA) is a computational method used to determine whether a predefined gene set shows a significant, consistent difference between 2 biological states. We ranked genes based on the degree of differential expression between the 2 groups. We then tested whether these predefined gene sets were enriched at the top or bottom of the ranked list. GSEA can detect subtle changes in gene expression. We used the local version of the GSEA analysis tool available at http://www.broadinstitute.org/gsea/index.jsp. We performed the GSEA independently using GO and KEGG datasets.
7) RNA extraction and RT‐qPCR
To isolate the somite from GD12.5 rat embryos, we excised the head, tail, forelimbs, hindlimbs, and internal organs (Supplementary Fig. 3A and B). We extracted total RNA from the somites of GD12.5 embryos using TRIzol reagent (Invitrogen, USA) following the manufacturer’s instructions. We reversetranscribed cDNAs using PrimeScript RTase (TaKaRa, China). RT‐qPCR was performed with SsoFast EvaGreen Supermix (Bio-Rad, USA). We used GAPDH as the reference gene. Melting curves were analyzed to ensure the specificity of the amplification. The primer sequences used in this study are listed in the supplementary materials (Supplementary Table 4). We determined relative expression levels using the 2−ΔΔCt method.
10. Statistical Analysis
We analyzed statistical data using GraphPad Prism 8 (Graph-Pad Software, USA). We utilized Student’s t-test to compare differences between the VAD group and the CON group. We considered differences with p-values <0.05 to be statistically significant. Data represent mean±standard error of the mean. The sample sizes and group names are specified in the figures.
RESULTS
1. VAD Induces CVM in Rats
To investigate whether VAD during embryonic development leads to CVM in rats, we randomly assigned 10-week-old virgin female SD rats into 2 groups. One group was fed a VAD diet (VAD group), and the other was fed a normal rodent diet (CON group) for 2 weeks before mating with male SD rats. The pregnant rats continued to receive either the VAD or CON diet throughout conception and lactation. We subjected the 2-week-old rat pups to x-ray and micro-CT scanning. We observed a butterfly vertebra at the T11 vertebral section in the VAD group. The posteroanterior view of the x-ray images showed that the T11 vertebra was divided into 2 parts, whereas the T11 vertebra in the CON group appeared normal (Fig. 1A and C). The lateral view of the x-ray images indicated that the degree of thoracolumbar kyphosis was greater in the VAD group compared to the CON group (Fig. 1B and D). The anteroposterior, posteroanterior, and lateral views from micro-CT scans also showed that the T11 vertebra in the VAD group was divided into 2 parts, confirming that this CVM was a butterfly vertebra. In contrast, the T11 vertebra in the CON group appeared normal (Fig. 1E–J). The 3D reconstruction of the T11 vertebra provided additional detail on this CVM. Anteroposterior, posteroanterior, top-bottom, and bottom-top views further confirmed the butterfly vertebra phenotype (Fig. 1K–R). The left-right and right-left views demonstrated that the anterior edge of the butterfly vertebra was thinner than the posterior edge, a characteristic known as anterior wedging of thoracic vertebrae. In contrast, the anterior and posterior edges of the T11 vertebra in the CON group showed similar thickness (Fig. 1S–V). This phenotype was also observed in another VAD neonatal rat (Supplementary Fig. 4A–V).
In total, we physically examined 49 neonates from the VAD group and 41 neonates from the CON group using x-ray and micro-CT at 2 weeks of age. Compared to the CON group, all vertebral malformations in the VAD group were butterfly vertebrae, predominantly located at the T10–13 vertebrae. The butterfly vertebrae appeared either continuously or intermittently at these locations (Supplementary Fig. 5A–Q). We found that 32.65% (16 of 49) of neonates in the VAD group exhibited CVMs, while no vertebral deformities were observed in the CON group (0 of 41) (Fig. 1W). To further confirm the impact of the VAD diet, we measured the serum vitamin A concentration in maternal rats from both the VAD and CON groups. After 2 weeks of feeding, maternal rats on the VAD diet had significantly lower serum vitamin A concentrations compared to those on the vitamin A sufficient diet (p<0.05) (Fig. 1X). Overall, these results indicate that VAD induces butterfly vertebrae, a specific type of CVM.
2. Genes Related to RA Signaling and Sclerotome Development Are Downregulated at GD10.5 in VAD Group
To further investigate the mechanism by which VAD induces CVM, we performed WMISH on GD10.5 rat embryos using probes targeting Raldh-2. We selected GD10.5 because Raldh-2 is expressed in the somite of GD10.5 rat embryos but not in GD12.5 embryos, and GD10.5 is a critical time point for sclerotome development [17]. In the CON group, Raldh-2 was expressed in the somite of GD10.5 rat embryos, while no signal was observed in the somite of GD10.5 rat embryos from the VAD group (Fig. 2A and B). This finding indicates that Raldh-2 expression is downregulated in the somite of GD10.5 embryos in the VAD group. We repeated the WMISH on GD10.5 rat embryos using Raldh-2 probes 3 times (Supplementary Fig. 6A–F).
To investigate whether sclerotome differentiation is suppressed in the VAD group, we performed WMISH on GD10.5 and GD12.5 rat embryos using probes targeting Pax1. We observed Pax1 expression in the sclerotome at GD10.5 and GD12.5 in the CON group, and at GD12.5 in the VAD group. However, no Pax1 signal was detected in the sclerotome of GD10.5 rat embryos in the VAD group (Fig. 2C–F). These results indicate that Pax1 expression is downregulated in the sclerotome of GD10.5 embryos in the VAD group, while the expression pattern at GD12.5 is similar between the CON and VAD groups. We repeated the WMISH for Pax1 on GD10.5 and GD12.5 rat embryos 3 times (Supplementary Figs. 7A–F and 8A–F).
Together, these results demonstrate that Raldh-2 is downregulated in the somite, and Pax1 is downregulated in the sclerotome at GD10.5 in the VAD group. Raldh-2 is a dehydrogenase responsible for converting retinaldehyde to RA [18]. GD10.5 is a critical time point for sclerotome development, and Pax1 is an essential transcription factor for this process [19]. Therefore, we hypothesize that downregulation of the RA signaling pathway may impair sclerotome development.
3. LCM Combined RNA-seq Revealed Transcriptome Abnormalities Related to Sclerotome Development and Retinal Metabolism
To investigate the transcriptional mechanisms underlying VAD-induced CVM, we performed LCM and RNA-seq on the sclerotomes of 22 GD12.5 rat embryos (11 embryos from the CON group and 11 embryos from the VAD group). In the transverse section model of a GD12.5 rat embryo, the sclerotome is surrounded by the neural tube, notochord, dermomyotome, and myogenic cells (Fig. 3A). Based on this anatomical context, we isolated the sclerotome using LCM (Fig. 3B and C). Figures of the LCM process for all 22 embryos are provided in the supplementary material (Supplementary Fig. 9A–V).
We conducted principal components analysis (PCA) to determine whether gene expression profiles could distinguish between the CON and VAD groups. PCA based on the expression levels of all genes showed some overlap between the CON and VAD groups (Supplementary Fig. 10A and B). In contrast, PCA based on differentially expressed genes (DEGs) revealed a clear separation between the CON and VAD groups (Supplementary Fig. 10C and D). These results indicate that the CON and VAD groups can be distinguished based on the expression levels of all genes and DEGs obtained from RNA-seq. Compared to the CON group, the VAD group exhibited upregulation of 659 genes and downregulation of 3,733 genes (p≤0.05, |log2FoldChange|≥1) (Fig. 3D). A heat map illustrates the clustering of DEGs across samples from the CON and VAD groups (Fig. 3E). GO enrichment analysis for the biological processes (BPs) associated with upregulated genes in the VAD group revealed enrichment in terms such as “negative regulation of cell development,” “negative regulation of locomotion,” “embryonic skeletal system morphogenesis,” “cartilage development,” “osteoblast differentiation,” and “somite development” (Fig. 4A; Supplementary Fig. 11A). This suggests that the upregulated genes in the VAD group are primarily involved in regulating skeletal developmental processes. For downregulated genes, GO enrichment analysis for BPs indicated significant associations with “vitamin transmembrane transport,” “calcium ion transmembrane transport,” “cell maturation,” “cell fate commitment involved in pattern specification,” “osteoblast development,” “paraxial mesoderm development,” “cell fate specification,” “embryonic axis specification,” and “bone cell development” (Supplementary Fig. 11B and C). This suggests that the downregulated genes are linked to skeletal development, vitamin transport, calcium ion transport, and cell fate determination. Additional GO enrichment analyses for cellular components and molecular functions were also conducted (Supplementary Fig. 12A–H). KEGG enrichment analysis revealed that upregulated genes in the VAD group are enriched in pathways such as the “transforming growth factor (TGF)-beta signaling pathway,” “Notch signaling pathway,” (Supplementary Fig. 13A and B). This finding suggests that the upregulated genes in the VAD group are linked to signaling pathways associated with the clock-and-wavefront model of somitogenesis [20]. Additionally, previous research has identified a signaling crosstalk between the TGF-beta and bone morphogenetic protein (BMP) signaling pathways, which may contribute to vertebral forma-tion failure [21]. The KEGG enrichment analysis also showed that downregulated genes in the VAD group are enriched in pathways related to “osteoclast differentiation,” “retinol metabolism,” “vitamin digestion and absorption,” “PI3K-Akt signaling pathway,” and “Wnt signaling pathway” (Fig. 4B; Supplementary Fig. 13C). These findings suggest that the downregulated genes may be associated with impaired retinol metabolism caused by VAD, the dynamic balance between osteoclasts and osteoblasts, and the clock-and-wavefront model of somitogenesis [20].
GSEA analysis indicated that upregulated genes are enriched in the “POSITIVE_REGULATION_OF_OSTEOBLAST_DIFFERENTIATION(GO:0045669),” reflecting an upregulated trend, and “RETINOL_METABOLISM(RNO00830),” showing a downregulated trend (Fig. 4C and D). We present the gene expression heatmap and random enrichment score distribution of both terms (Supplementary Fig. 14A–D). Core enrichment genes for “POSITIVE_REGULATION_OF_OSTEOBLAST_DIFFERENTIATION( GO:0045669)” include BMP2, BMP7, Smad1, and Smad5 (Supplementary Table 5). Core enrichment genes for “RETINOL_METABOLISM(RNO00830)” include Aldh1a1, Aldh1a7, Adh1, and Cyp4a8 (Supplementary Table 6). These results suggest that Smad-dependent BMP signaling may be upregulated in the sclerotome of GD12.5 rat embryos in the VAD group. Furthermore, the downregulation of retinol metabolism in the VAD group indicates that the RA signaling pathway may be downregulated in the sclerotome of GD12.5 rat embryos in this group.
4. Validation of Genes Related to RA Signaling and Osteoblast Differentiation at GD12.5 in Sclerotome
To validate the expression levels of genes related to the RA signaling pathway, we isolated somite from GD12.5 rat embryos and performed RT-qPCR. RNA-seq results indicate that Cdx1, Drd2, Olig2, Hoxa5, Hnf1b, and Pitx2 were downregulated in the VAD group (Fig. 5A). Similarly, RT-qPCR results revealed that the relative expression levels of Cdx1, Drd2, Olig2, Hoxa5, Hnf1b, and Pitx2 were also downregulated in the VAD group (Fig. 5B). Previous reports have shown that RA stimulation upregulates the expression of these genes [22]. Thus, our results suggest that RA signaling is downregulated in the sclerotome at GD12.5.Given that both GO enrichment and GSEA results suggest a relationship between osteoblast differentiation and the DEGs in the VAD group, and that BMP2, BMP7, Smad1, and Smad5 are core enrichment genes in GSEA, we conducted RT-qPCR to validate these findings by analyzing genes involved in the BMP signaling pathway and osteoblast differentiation. RNA-seq results showed that Nog, Chrdl1, Bmp2, and Smad1 were upregulated in the VAD group, while Mkk3, Dlx5, Osterix, Runx2, and Smad6 were downregulated in the VAD group (Fig. 5C). Similarly, RT-qPCR results revealed that the relative expression levels of Nog, Chrdl1, Bmp2, and Smad1 were upregulated in the VAD group, while the relative expression levels of Mkk3, Dlx5, Osterix, Runx2, and Smad6 were downregulated in the VAD group (Fig. 5D). Key DEGs by pathway were summarized (Supplementary Table 7).
Previous studies have shown that Nog and Chrdl1 encode proteins that antagonize BMP activity, Mkk3 is involved in non-Smad-dependent BMP signaling, Smad1 participates in Smaddependent BMP signaling, Smad6 encodes proteins that negatively regulate Smad signaling, and Dlx5, Osterix, and Runx2 encode transcription factors that regulate genes related to osteoblast differentiation [23]. Thus, our results suggest that the non-Smad-dependent BMP signaling pathway is suppressed by the upregulation of antagonists such as Noggin, which inhibits the expression of osteoblast differentiation-related transcription factors. Concurrently, the Smad-dependent pathway attempts to counteract this suppression by upregulating Smad1 and downregulating Smad6. However, this attempt is insufficient, as the expression levels of osteoblast differentiation-related genes remain downregulated in the sclerotome of the VAD group.
DISCUSSION
Based on our findings, we constructed a working model to illustrate the underlying mechanism, VAD leads to the downregulation of the non-Smad-dependent BMP signaling pathway by suppressing the RA signaling pathway, this suppression inhibits osteoblast differentiation, resulting in sclerotome dysplasia and ultimately suppresses vertebral formation, leading to CVM (Fig. 5E).
Butterfly vertebrae represent a type of CVM resulting from the failure of ventral ossification of the vertebral body during gestation. They typically occur in the thoracic and lumbar regions of the spine, while cervical butterfly vertebrae are rare [24]. Micro-CT results showed that all CVMs in our study were identified as butterfly vertebrae, specifically located at the T10–13 vertebrae in rats. This pattern closely mimics the condition observed in actual CVM patients [25]. Additionally, most butterfly vertebrae in the VAD group exhibited anterior wedging, characterized by a thinner anterior edge of the vertebral body compared to the normal vertebrae in the CON group. This phenomenon closely resembles a case report of a 46-year-old woman with butterfly vertebrae at T6, who also exhibited anterior wedging of the T6 vertebra [26]. Anterior wedging of thoracic vertebrae may contribute to an increased degree of kyphosis in the thoracic curve, which was also observed in our model (Fig. 1G and J). These similarities show a strong connection between our animal model and clinical practice. Thus, potential translational implications of this animal model should be noticed, for instance, vitamin A supplement of maternal nutrition, vitamin A related public health strategies, etc. In addition, the incidence rate of CVMs in our study was 32.65%, which is higher than the previously reported incidence rate of spinal anomalies in the thoracic region (13.8%) [7]. This discrepancy may be attributed to the enhanced accuracy of micro-CT, which allows for the detection of more detailed skeletal abnormalities.
Raldh-2, also known as Aldh1a2, is a dehydrogenase responsible for converting retinaldehyde to RA [18]. The early phase of body axis extension is guided by trunk neuromesodermal progenitors (NMPs), which generate trunk somite, while the later phase is directed by tail NMPs that produce tail somites [27]. In Raldh-2−/− embryos, trunk somites are smaller than normal, indicating that the RA-generating function of Raldh-2 is crucial for trunk somitogenesis. However, Raldh-2 is not necessary for tail somitogenesis [28]. Our WMISH results showed that Raldh-2 expression was downregulated in the somite of GD10.5 rat embryos in the VAD group. This finding suggests that the GD10.5 rat embryos in the VAD group lack RA generated by Raldh-2, leading to downregulation of the RA signaling pathway in the somite. This deficiency may disrupt trunk somitogenesis.
BMP signaling plays a crucial role in osteoblast differentiation. BMPs promote osteoblast differentiation through 2 pathways: the Smad-dependent and non-Smad-dependent pathways [29]. In the Smad-dependent pathway, BMPs activate Smad1/5/8 as their R-Smad. The phosphorylated R-Smad complex then associates with Smad4 and translocates into the nucleus. In the nucleus, this complex would recruit co-factors and Runx2 to regulate the expression of osteogenic genes, including Runx2, Dlx5, and Osterix [30,31]. In the non-Smad-dependent pathway, phosphorylated TAK1 recruits TAB1 to initiate the MKK-p38 MAPK or MKK-ERK1/2 signaling cascade. MAPK then phosphorylates Runx2, Dlx5, and Osterix to enhance their transcriptional activity [23,29]. Additionally, BMP signaling is regulated by extracellular antagonists such as Noggin, Chordin, Grem1, and Grem2. Noggin expression is also regulated by TGF-β signaling, indicating crosstalk between BMP and TGF-β signaling pathways [32,33].
Our RNA-seq and RT-qPCR results showed that the expression levels of Noggin, Chrdl1, Bmp2, and Smad1 were upregulated in the VAD group, while the expression levels of Mkk3, Dlx5, Osterix, Runx2, and Smad6 were downregulated. These findings suggest that the non-Smad-dependent signaling pathway is suppressed, whereas the Smad-dependent pathway is upregulated. The suppression of non-Smad-dependent signaling pathway would downregulate osteoblast differentiation-related genes, and the upregulation of Smad-dependent signaling pathway would upregulate osteoblast differentiation-related genes. However, this upregulation did not rescue the expression levels of osteoblast differentiation-related genes. Furthermore, KEGG enrichment analysis of upregulated genes indicated that TGF-β signaling is enriched, suggesting that its upregulation may regulate Noggin expression, thereby affecting BMP signaling.
The interaction between the RA signaling pathway and the BMP signaling pathway plays a crucial role in ossification. RA signaling can activate non-Smad-dependent BMP signaling and promote bone formation by engaging various MAPK kinase cascades [34]. Additionally, RA signaling enhances ectopic bone formation induced by BMP signaling [35]. Conversely, RA signaling suppresses Smad-dependent BMP signaling and reduces ossification. Reports indicate that RA signaling inhibits chondrogenesis and heterotopic ossification by diminishing Smaddependent BMP signaling [36]. RA also represses Smad-dependent BMP signaling by lowering the level of phosphorylated Smad1 [37]. These findings suggest that the inhibition of RA signaling may lead to a decrease in non-Smad-dependent BMP signaling and an increase in Smad-dependent BMP signaling.
Pax1, a transcription factor crucial for sclerotome development, is expressed in the developing spine of rat embryos [19]. Pax1 encodes a paired-domain transcription factor that is present in early sclerotome [38]. Heterozygous Pax1 knockout mice exhibit various spinal deformities, including a kinky tail phenotype [39]. Pax1 regulates extracellular matrix genes, such as collagen and aggrecan, and is essential for mesenchyme condensation and intervertebral disc development [40]. Although numerous studies have elucidated the biological functions of Pax1 in spine development, the mechanism by which it contributes to CVM pathogenesis remains largely unknown.
In our WMISH results, we observed Pax1 expression in the sclerotome of GD10.5 rat embryos in the CON group, but not in the VAD group. By GD12.5, Pax1 expression levels appeared similar between the CON and VAD groups in the sclerotome. These findings indicate that the VAD environment downregulates Pax1 expression in the sclerotome at GD10.5, but this downregulation is no longer evident at GD12.5, reflecting a dynamic suppression effect of VAD. Furthermore, since Pax1 encodes a transcription factor involved in sclerotome development, our results suggest that sclerotome development is suppressed at GD10.5 in the VAD group.
However, this study has some limitations. WMISH results show that Raldh-2 expression is downregulated in the somite at GD10.5, and Pax1 expression is also downregulated in the sclerotome at GD10.5. Consequently, we initially attempted to perform laser-captured microdissection and RNA-seq using GD10.5 embryos. However, the GD10.5 embryos were too small to control accurately during OCT embedding. Therefore, we used GD12.5 embryos as a replacement. Other limitations should also be aware, such as absence of functional rescue experiments, lack of behavioral outcomes to assess functional significance of the deformities, and the limited applicability of rodent embryogenesis to human development.
CONCLUSION
In summary, we established a VAD-induced CVMs rat model. All vertebral malformations were butterfly vertebrae located at the 10th to 13th thoracic vertebrae. Non-Smad-dependent BMP pathway and RA signaling pathway may play a critical role in the pathogenesis of CVMs. Our findings demonstrate that VAD may be one of the causes of CVMs, which is hypothesized to serve as a novel therapeutic target for the nonsurgical treatment of CVMs in the future.
PDF Links
PubReader
ePub Link
XML Download
Full text via DOI
Full text via PMC
Download Citation
Supplement
Print
