TBX3 regulates the transcription of VEGFA to promote osteoblasts proliferation and microvascular regeneration

Data quality control and pre-processing

Single-cell RNA sequencing data were downloaded from sample GSM4119369 (from the GSE138791 dataset) of human periosteal-derived cells (hPDCs) and treated in a serum-free growth medium for six days (Bolander et al., 2020). Low-quality cells were defined as follows: (1) cells with less than 200 genes detected; (2) cells with >4,000 genes detected; and (3) cells with more than 5% mitochondrial genes. Low-quality cells are usually considered non-viable cells and therefore removed from the dataset. Finally, we screened 11,119 cells and 16,884 genes.

Gene expression data were normalized and scaled using the Seurat (Version 4.1.0) package (Hao et al., 2021). The data were normalized using the “NormalizeData” function, the normalization method was set to “LogNormalize”, and the scaling factor was 100000. Then, the “FindVariableFeatures” function was used to filter the first 2000 high-variable genes, and the subsequent downscaling to these genes was restricted. Further, the “RunPCA” function was used to perform principal component analysis (PCA) on the single-cell expression matrix restricted to the highly variable genes. The signal-to-noise ratio was improved by selecting partially significant principal components. The significance of each principal component was calculated using the “JackStraw” method provided in the Seurat package, and the top 40 principal components were selected for subsequent analysis.

Cell clustering analyses were done using the “FindClusters” function, with the resolution set to 1.5, and by embedding cells into the graph structure in PCA space. Subsequently, RunUMAP and RunTSNE functions were used to perform UMAP and tSNE dimensionality reduction analysis.

Cell type annotation

Cell subgroups were annotated according to the expression levels and ratios of marker genes in the clustered subgroups. The marker genes corresponding to the cell types were: mitotic cell - MKI67, AURKA, CDC20, CDK1; pro-hematopoietic osteochondral progenitor- ANGPTL7; BMPR2+ osteochondral progenitor—POSTN, BMPR2; CD34+ osteochondral progenitor—POSTN, CD34; Osteochondral progenitor—POSTN, PDGFRB, THY1.

Gene pooling enrichment analysis

Gene Ontology (GO) provides a framework to describe the function of gene products in all organisms. GO annotations include three major categories: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). The KEGG knowledge base collects a large amount of manually proofread information on biological pathway networks. To enrich the set of the gene of interest into GO terms and KEGG terms, the “enrichGO” and “enrichKEGG” functions in the clusterProfiler package were used, respectively (Version 4.2.2; Yu et al., 2012; Wu et al., 2021). In addition, the “GSEA” function was used to determine whether the pathway gene sets were significantly different between the two biological states, and an FDR < 0.05 was set as significant enrichment.

Transcriptional regulation analysis

To infer the gene regulatory networks of VEGF+ cell subpopulations, the SCENIC (Version 1.2.4, Aibar et al., 2017) toolkit was used. SCENIC consisted of three steps: (1) Identification of co-expression modules between TFs and potential target genes; (2) for each co-expression module, identification of the direct target of the TF by TF-motif enrichment analysis genes, and grouping of each TF and its corresponding direct target gene into a regulon; (3) scoring each regulon activity for every cell using the AUCell algorithm (AUC).

To quantify the cell type specificity of regulon, a pre-described algorithm was used (Suo et al., 2018). A vector to represent the distribution of the activity fraction of a given regulon in the cell population was used. Then its activity was normalized to a sum of 1. Moreover, to indicate whether each cell belongs to a particular type, a vector was set with its elements taking values of 0 or 1, with 1 indicating that the cell belongs to a specific cell type and 0 suggesting that it did not. We also normalized this vector to a sum of 1. Next, we calculated the Jensen–Shannon scatter (JSD) between these two vectors, often used to quantify the difference between two probability distributions. The JSD values range from 0 to 1, 0 indicating the same distribution and 1 indicating a significant difference. Finally, the modulator specificity score (RSS) was defined by converting the JSD to a similarity score.

The smaller the JSD, the larger the RSS. Significant regulators were defined as those with the highest cell type specificity score for each cell type. Cytoscape (Version 3.8.0) was to visualize the regulatory-related networks of target cells.

Regulation module analysis

Transcriptional regulatory modules were identified using the Conjugation Specificity Index (CSI). The calculation of CSI consisted of two steps. First, the Pearson correlation coefficient (PCC) between the activity scores of each pair of regulators was calculated. Next, the Euclidean distance hierarchical clustering method based on the CSI matrix was used to identify different regulatory modules. For each regulatory module, the activity score for a cell type was defined as the average of the activity scores of its regulatory members across all cells of that type. The top-ranked cell types were then identified for each module.

Cell culture and transfection

The human osteoblast cell line was purchased from the Institute of Cell Research, Chinese Academy of Sciences, Shanghai, China, and cultured in DMEM high sugar medium, 10% FBS, and 1% double-antibody in a 5% CO₂ constant temperature and humidity incubator. After the cells were treated with TBX3 expression intervention, the cells were transfected with 100 pmol TBX3 siRNA, and total cellular RNA was extracted after 48 h of treatment.

Total cellular RNA extraction and qRT-PCR analysis

Total cellular RNA was extracted with Trizol (Beyotime, Shanghai, China) reagent, and the experiments were performed according to the commercial instructions. The extracted RNA was quantified by NanoDrop and reverse transcribed into cDNA by GoScript™ Reverse Transcription Kit (Promega, Wisconsin, USA), followed by FastFire qPCR PreMix (SYBR Green, TIANGEN, China) to determine the expression level of mRNA. The mRNA expression levels were calculated as 2^−ΔΔCt, and GAPDH was used as an internal reference.

Protein extraction and western blot

Total protein was extracted by RIPA according to the manufacturers’ protocol, containing protease inhibitors. Equal amounts of protein were separated by 10% SDS-PAGE gel and transferred to PVDF membranes under certain conditions. The primary antibody was then incubated overnight under 4 °C, the next day, the secondary antibody was incubated for 1–2 h at room temperature, followed by exposure to developer and storage of the images through Tanon 5200 (Tanon, China). The primary antibodies used above include anti-VEGFA (1:1000, abcam, USA) and anti-GAPDH (1:2000, Servicebio, Wuhan, China). Among them, GAPDH was used as an internal reference.

CCK-8 cell proliferation assay

We used the CCK-8 assay to determine the proliferative capacity of the cells. We divided the cells of TBX3 expression intervention for 48 h into 96-well plates equally, with 1*10³ cells per well. And we applied CCK-8 for cell proliferation assay every other day for a total of 5 days. For CCK-8 assay, 10 µL of CCK-8 reagent was added to each well, and OD450 value was measured after 3 h.

Statistical analysis

All experimental samples, from at least three independent experiments, were expressed as mean ± standard deviation. The Chi-square test was used for the statistical analysis of frequencies between two groups. Statistical differences between two independent samples were analyzed by t-test. P <0. 05 was considered statistically significant.

Identification of VEGF+ cell subpopulations from human periosteal-derived cells (hPDCs)

Single-cell RNA sequencing data were acquired and processed (See ‘Methods’). After quality control and discarding low-quality cells, 11,119 hPDCs were used for further analysis (Figs. S1A–S1C). Based on the expression levels of marker genes in each clustered subgroup, the clustered subgroups were annotated as CD34+ osteochondral progenitor (n = 4563), osteochondral progenitor (n = 5316), pro −hematopoietic osteochondral progenitor (n = 600), BMPR2+ osteochondral progenitor (n = 534) and mitotic cells (n = 106, Fig. S1D–S1F).

We classified all cells into VEGF+ cells (expressing at least two genes) and VEGF- cells (expressing at most one gene) based on the expression of VEGFA, VEGFB, and VEGFC (Figs. 1A–1B). VEGF+ cells were considered angiogenesis-related cells. We obtained 8697 VEGF- cells and 2422 VEGF+ cells. We next compared the significantly different GO and BP pathways between VEGF- and VEGF+ cells by gene pooling enrichment analysis. The results showed that pathways such as vascular vein development, VEGF-stimulated cellular responses, histogenesis, and VEGF pathways were significantly activated in VEGF+ cells (Fig. 1C), indicating that VEGF+ cells were indeed associated with an angiogenic activity.

Reconstruction of transcriptional regulatory networks

To comprehensively reconstruct the gene regulatory network in a subpopulation of VEGF+ cells, we analyzed VEGF+ cell expression data using the SCENIC process. SCENIC linked cis-regulatory sequence information to single-cell RNA sequencing (scRNA-seq) data. SCENIC consists of three main steps: co-expression analysis, target gene motif enrichment analysis, and scoring of regulatory activity. The main results included a series of regulators (each representing a TF and a set of co-expressed and motif significantly enriched target genes) and a regulatory activity score (RAS) for each cell. The SCENIC analysis showed that out of the 274 regulators screened, 145 were activated in more than 1% of the cells (Table S1). We focused on those 145 transcriptional regulators. The number of target genes per regulator ranged from 10 to 918, involving 6,968 genes.

We calculated the specificity scores between the 145 regulons and different cell types separately to understand if regulons are specific to a particular cell type. We found that: TCF7L2 was most associated with CD34+ osteochondral progenitor; TBX3 was predominantly associated with osteochondral progenitor; FLI1 was most associated with pro-hematopoietic osteochondral progenitor; NFKB2 was most associated with BMPR2+ osteochondral progenitor and; EZH2 was most associated with Mitotic cells (Figs. S2A–S2E). We showed the distribution of different cell types with the corresponding most relevant regulon activated cells in TSNE space (Figs. S2A–S2E). TCF7L2 is mainly enriched in CD34+ osteochondral progenitor; TBX3 is mainly enriched in osteochondral progenitor; FLI1 is mainly enriched in pro-hematopoietic osteochondral progenitor; NFKB2 is mainly enriched in BMPR2+ osteochondral progenitor; and EZH2 is mainly enriched in mitotic cells.

To systematically characterize the combination patterns among the 145 regulons, we compared the similarity between each regulon based on the CSI matrix and divided them into eight modules (M1- 21 regulons; M2- 59 regulons; M3- 15 regulons; M4- 11 regulons; M5- 24 regulons; M6- 5 regulons; M7- 8 regulons; M8- 2 regulons) (Fig. 2A). For each module, we identified several representative cell types by their average activity scores and the activity distribution of the modules in VEGF+ cells (Figs. S3A–S3B). These results demonstrate the presence of M1-M8 transcription factors in different cell types and individuals. In addition, we also list the TFs and their binding motifs that represented the cell type specificity within each module (Fig. 2B).

Finally, we used Cytoscape to visualize the regulatory network of important TFs (Fig. S4A). To explore which biological functions these TFs were mainly involved in, we performed GO and KEGG enrichment analysis on 145 TFs, showing the enrichment results of the top 10 KEGG pathways, BP, MF, and CC pathways that were significantly enriched (Figs. S4B– S4C). Pathway analysis showed that those TFs were involved mainly in transcriptional misregulation in cancer, hepatitis B, human T-cell leukemia virus 1 infection, viral carcinogenesis, kaposi sarcoma-associated herpesvirus infection, c-type lectin receptor signaling pathway, lipid and atherosclerosis, acute myeloid leukemia, hepatitis C, Th17 cell differentiation (P < 0.0001; Fig. S4B). Those TFs also enrichment in BPs (intracellular receptor signaling pathway, rhythmic process, pri-miRNA transcription by RNA polymerase II, mononuclear cell differentiation, anterior/posterior pattern specification, gland development, circadian rhythm, regionalization, and lymphocyte differentiation) MFs (DNA-binding transcription activator activity, DNA-binding transcription repressor activity, DNA-binding transcription factor binding, nuclear receptor activity, ligand-activated transcription factor activity, core promoter sequence-specific DNA binding, and chromatin DNA binding), and CCs (tracription regulator complex, RNA polymerase transcrintion regulator complex) (P < 0.0001; Fig. S4C). The results of these enrichment analyses clearly demonstrate the functions and pathways that are enriched by these TFs.

TBX3 promoted the VEGF expression and cell proliferation

We found that TCF7L2, TBX3, FLI1, NFKB2 and EZH2 were most closely related to the corresponding cells by bioinformatics prediction. According to the above study, TCF7L2, TBX3 and EZH2 are closely related to osteoblasts. And among the eight modules and cells, the osteochondral progenitor was closely associated with most of the modules. The TF most closely related to these cells was TBX3. We knocked down the expression of TBX3 in human osteoblast cell lines by treating them with a TBX3 downregulation. We then examined the expression of VEGFA, VEGFB, and VEGFC using qRT-PCR and found that their levels were decreased in the TBX3 knockdown group (Fig. 3A). Furthermore, the downregulation of TBX3 resulted in low expression of VEGFA by western blot (Fig. 3B). Subsequently, we examined the proliferation ability of the osteoblast cell line using CCK-8 and found that TBX3 could promote cell proliferation (Fig. 3C). Thus, we proposed that TBX3 is a transcription factor involved in skeletal microvasculature development by regulating VEGF expression and promoting osteoblast proliferation.