Single-cell RNA Sequencing Reveals Thoracolumbar Vertebra Heterogeneity and Rib-genesis in Pigs

Abstract

Development of thoracolumbar vertebra (TLV) and rib primordium (RP) is a common evolutionary feature across vertebrates, although whole-organism analysis of the expression dynamics of TLV- and RP-related genes has been lacking. Here, we investigated the single-cell transcriptome landscape of thoracic vertebra (TV), lumbar vertebra (LV), and RP cells from a pig embryo at 27 days post-fertilization (dpf) and identified six cell types with distinct gene expression signatures. In-depth dissection of the gene expression dynamics and RNA velocity revealed a coupled process of osteogenesis and angiogenesis during TLV and RP development. Further analysis of cell type-specific and strand-specific expression uncovered the extremely high level of HOXA10 3'-UTR sequence specific to osteoblasts of LV cells, which may function as anti-HOXA10-antisense by counteracting the HOXA10-antisense effect to determine TLV transition. Thus, this work provides a valuable resource for understanding embryonic osteogenesis and angiogenesis underlying vertebrate TLV and RP development at the cell type-specific resolution, which serves as a comprehensive view on the transcriptional profile of animal embryo development.

Keywords: Angiogenesis; Osteogenesis; Rib-genesis; Thoracolumbar vertebra transition; scRNA-seq.

Figure 1

Cell composition and differentiation trajectory of developing TLV A. Overview of the experimental design. Segmentation regions between TV and LV were dissected in pigs and dissociated into single cell suspensions by micromanipulation and enzyme digestion. Red dashed rectangle indicates the TLV segmentation joint. Smart-seq2 was used for scRNA-seq. B. Integrated UMAP plot of cell clusters from TV and LV. C. Gene expression patterns of each cell cluster in (B). Cell type-enriched genes are listed on right and are labeled in the same colors as corresponding cell types. D. Bifurcation of the 799 top TLV DEGs along two branches clustered hierarchically into six modules in a pseudo-temporal order. Development trajectories of TV and LV cells are shown on the right and left, respectively. Red arrow indicates the pseudo-time of cell fate 1 (MSC to HEC); blue arrow indicates the pseudo-time of of cell fate 2 (OB to CT). Representative genes are shown on the right. E. RNA velocity recapitulating the dynamics of TLV cell differentiation. The arrows indicate the position of the future state. F. Expression pattern (left), unspliced–spliced phase portrait (middle; cells colored according to E), and u residual (right) of TLV cells are shown for CD248, CD34, COL1A1, and MATN1. TLV, Thoracolumbar vertebra; TV, thoracic vertebra; LV, lumbar vertebra; scRNA-seq, single-cell RNA sequencing; UMAP, uniform manifold approximation and projection; OB, osteoblast; FB, fibroblast; SC, stroma cell; CT, cartilage; MSC, mesenchymal stem cell; HEC, hemogenic endothelial cell; DEG, differentially expressed gene.

Figure 2

Cell heterogeneity and HOXA10 expression difference between developing TV and LV A. Median scaled ln-normalized gene expression of selected DEGs for LV and TV cell clusters. B. Scatter plot comparing the average expression levels of genes in OB cell clusters from LV and TV. C. Vinplot comparing expression of HOXA10, HOXC10, and HOXD10 in each cell cluster from LV and TV. D. Average sequencing depth of HOXA10 coding sequence, HOXA10 3′-UTR, and HOXA10-AS in 137 LV and 128 TV cells. Shade rectangles indicate the region of HOXA10 coding sequence and HOXA10 3′-UTR. E. Boxplot comparing the FPKM values of HOXA10-exon3 and HOXA10-AS in 137 LV cells. P value was obtained by unpaired two-sided Welch’s t-test with correction for multiple comparisons. F. Scatter plot showing the number of reads harboring HOXA10 poly(A) tail and HOXA10-AS poly(A) tail at the HOXA10 3′-UTR locus in 137 LV cells. HOXA10-AS indicates an antisense RNA which overlaps the 3′-UTR on the opposite strand of the HOXA10 gene. FPKM, reads per kilobase of exon model per million mapped reads.

Figure 3

Cell composition and differentiation trajectory of developing RP A. Integrated UMAP plot of cell clusters from RP. B. Median scaled ln-normalized gene expression of selected DEGs for RP cell clusters from (A). Cell type-enriched genes are listed on the right and labeled in the same colors as corresponding cell types. C. Bifurcation of the 381 top RP DEGs along two branches clustered hierarchically into five modules in a pseudo-temporal order. Development trajectories of RP cells are shown on the right and left, respectively. Red arrow indicates the pseudo-time of cell fate 1 (HEC); blue arrow indicates the pseudo-time of cell fate 2 (OB). Representative genes are shown on the right. D. RNA velocity recapitulating the dynamics of the RP cell differentiation. The arrows indicate the position of the future state. E. Expression pattern (left), unspliced–spliced phase portrait (middle; cells colored according to D), and u residual (right) of the RP cells are shown for CD248, CD34, COL1A1, and MATN1. RP, rib primordium.

Figure 4

Osteogenesis network construction and cell type-specific marker immunofluorescence analysis of developing TLV and RP A. Module visualization of the network connections and associated functions using MATN1 as a hub gene. Gene-connected intra-modular degree is simultaneously indicated by spot size and color intensity. The hub gene MATN1 is indicated in yellow. B. and C. Immunofluorescence analysis of MATN1 and HMGB2 in TV and RP (B) and in LV (C). Red and green indicate fluorescence signals of MATN1 and HMGB2, respectively. White and yellow triangles indicate MATN1⁺ and HMGB2⁺ cells, respectively. D. and E. Immunofluorescence analysis of COL1A1 in TV and RP (D) and in LV (E). Red indicates fluorescence signals of COL1A1. Yellow triangles indicate COL1A1⁺ cells. White, blue, and red dashed lines in (B–E) indicate regions of TV, LV, and RP, respectively. Scale bar, 400 μm.

Figure 5

Angiogenesis network construction and cell type-specific marker immunofluorescence analysis of developing TLV and RP A. Module visualization of the network connections and associated functions using CD34 as a hub gene. Gene-connected intra-modular degree is simultaneously indicated by spot size and color intensity. The hub gene CD34 is indicated in yellow. B. and C. Immunofluorescence analysis of CD93, CD248, and CD34 in TV and RP (B) and in LV (C). Green, white, and red indicate fluorescence signals of CD93, CD248, and CD34, respectively. White, yellow, and blue triangles indicate CD93⁺, CD248⁺, and CD34⁺ cells. White, blue, and red dashed lines in (B and C) indicate the regions of TV, LV, and RP, respectively. Scale bar, 400 μm.