Bone marrow mesenchymal stem cells (BMSCs) show considerable promise for regenerative medicine, and tissue engineering because they have the competence of self-renewal and differentiation into various types of tissue-specific cells, including osteoblasts, chondrocytes and adipocytes (1, 2). Many studies indicated that BMSCs was a highly heterogeneous cell population composed of different cell subpopulations, which may be the basis of their multiple biological characteristics (3, 4). However, it is not clear which cell subpopulations with different properties or functions they consist of. In addition, most of the BMSCs in current studies came from younger donors, not elderly people.

By collectively sequencing tens of thousands of cell samples, bulk RNA sequencing only reveals the average effect of the cell population or information about the larger number of cells, masking the heterogeneity among cells. Therefore, the ideal method to analyze which cells the tissue consists of or to explore the heterogeneity among cells with the same type is to analyze individual cells in the sample. In 2009, Tang et al. (5) combined single-cell RNA amplification technology with second-generation sequencing technology to conduct single-cell RNA sequencing (scRNA-seq) experiments on mature mouse egg cells for the first time. With the development of single-cell isolation, gene amplification, sequencing, and other technologies, scRNA-seq has been gradually applied to study tumor cells, immune cells, etc., but few studies have used it to analyze human BMSCs (6, 7). It is thus of great interest for us to dissect the cellular heterogeneity of BMSCs by scRNA-seq profiling.

Isolation and culture of BMSCs

The experimental sample was taken from an 85-year-old female patient with intertrochanteric fracture admitted to our hospital. Informed consent was obtained from all participants. This study was approved by the ethics committee of the institutional review board at the Seventh Medical Center of Chinese PLA General Hospital (Number. 2017-87).

This patient was treated with intramedullary nailing. During the operation, after opening the apex of the trochiter with a drill, the aspirator (about 20 cm in length) was inserted into the medullary cavity. An injector was used to connect the tail of the aspirator, and about 10 ml of bone marrow was extracted from the fractured femur (Supplementary Fig. S1). Pay attention to avoid mixing with peripheral blood, and avoid violent oscillations during transportation (8).

This experiment adopted density gradient centrifuging to separate bone marrow mononuclear cells of the donor. The cells were cultured in dulbecco’s modified eagle medium (DMEM) (StemRD, USA) containing 10% fetal bovine serum (FBs) (StemRD, USA) and 1% penicillin and streptomycin (Corning, USA). The cells were maintained in a 5% CO2 humidified incubator at a constant temperature of 37℃. The medium was changed every 3 days. The cells were observed by inverted phase contrast microscope. When the adherent cells reached 80% confluence, they were digested with 0.25% trypsin (Corning, USA) and passaged. The BMSCs at passage 2 were harvested for the further experiments.

Flow cytometry analysis

The BMSCs at passage 2 were digested with trypsin, washed with phosphate-buffered saline (PBS), and stained with CD29-APC (BD Biosciences, USA), CD44-FITC (BD Biosciences, USA), CD90-FITC (BD Biosciences, USA), CD105-APC (BD Biosciences, USA), CD34-PE (BD Biosciences, USA), CD45-PE (BD Biosciences, USA), and HLA-DR-PE (BD Biosciences, USA). The cells were washed and resuspended with PBS. Flow cytometry analysis was performed on the flow cytometer (BD Biosciences, USA).

ScRNA-seq

The BMSCs at passage 2 were used to prepare a single-cell suspension. 10x Genomics Chromium was used to construct scRNA-seq response system, and all experiments were carried out according to the official operating manual (9). The single cells were separated by droplet micro-fluidics. A Gel Bead with Read1, cell label (Barcode), Unique Molecular Identifiers (UMI), and Poly (dT) was mixed with a single cell through oil phase to form the Gel Bead in Emulsions (GEMs). Barcode could distinguish different cells, UMI could distinguish different cDNA in the same cell, and Poly (dT) was used to enrich mRNA (10, 11). After lysing the cells, free mRNA underwent reverse transcriptional reaction and cDNA with Barcode and UMI information was generated. The Emulsions were then broken and the appropriate length of cDNA was selected as the template for PCR amplification. The Single Cell 30 Reagent Kit (v2 chemistry) was used to construct library. Agilent 2100 bioanalyzer was used to detect the peak type and fragment size of library. The resulting libraries were sequenced on the Illumina System.

Bioinformatics analysis

The official 10x Genomics pipeline Cell Ranger was used to deal with raw data. The BCL format file generated by Illumina sequencers was converted to Fastq format file by using the Cell Ranger mkfastq pipeline. Next, the reads library was processed by using the Cell Ranger count pipeline to generate a gene-barcode matrix. Based on the STAR algorithm (12), the reads were aligned to a human reference genome. Cell barcodes and UMIs were calibrated and filtered. The expected number of recovered cells was set to 5,000, as a parameter related to cell barcode filtering.

The dimension of the data on these genes was reduced through principal component analysis (PCA). The cells were clustered into diverse clusters through the k-means clustering approach. Next, the clustered cells were visualized in two-dimensional projection of t-distributed stochastic neighbour embedding (t-SNE).

The Loupe Browser was used to analyze the expression of the International Society for Cellular Therapy (ICST)-proposed positive marker genes CD29 (ITGB1), CD44, CD90 (THY1), CD105 (ENG), and negative marker genes CD34, CD45 (PTPRC), HLA-DRA (1). Besides, the expression of markers for early passage, FGFR2, and late passage BMSCs, PLAT, were researched (13).

Based on the results of cell clustering, the expression difference of each gene between that cluster and the average of the rest of clusters were computed to identify the differentially expressed genes (DEGs) of each cluster. Meanwhile, the corresponding fold change (FC) and p-value of each gene were obtained (14, 15). DEGs must meet the following conditions: (1) mean UMI counts>1; (2) log2FC>1; (3) p<0.05. The top 20 DEGs were selected according to the log2FC sequence from large to small. Kobas 3.0 was used to conduct Gene Ontology (GO) enrichment analysis on the top 20 DEGs of each cluster, so as to determine their main functions. The function of the top 20 DEGs of each cluster was checked by the National Center for Biotechnology Information (NCBI) gene database (https://www.ncbi.nlm.nih.gov/gene/).

Expression of BMSCs markers

Almost all the cells expressed CD29 (99.95%), CD44 (96.18%), CD90 (99.32%), and CD105 (99.50%). Inversely, very few cells expressed CD34 (1.69%), CD45 (0.11%), or HLA-DR (0.63%) (Fig. 1a).

Figure 1. Expression of BMSCs markers by flow cytometry analysis and scRNA-seq. (a) Expression of surface marker antigens of the BMSCs. (b) Expression patterns of characteristic marker genes were projected on the t-SNE plot. (c) Expression patterns of marker gene for early passage BMSCs, FGFR2, and late passage BMSCs, PLAT, were projected on the t-SNE plot. Each cell is coloured based on the expression level of the indicated marker gene.

The expression of a variety of characteristic marker genes in BMSCs was examined. We found that the expression of the ISCT-proposed positive marker genes, such as CD29 (ITGB1), CD44, CD90 (THY1), and CD105 (ENG) was very high. On the contrary, the cells lacked expression of the negative marker genes, such as CD34, CD45 (PTPRC), and HLA-DRA (Fig. 1b). The proportion of cells expressing the marker genes is as follows: CD29 (ITGB1): 99.94%, CD44: 99.74%, CD90 (THY1): 99.75%, CD105 (ENG): 87.27%, CD34: 0.12%, CD45 (PTPRC): 0.23%, HLA-DRA: 0.03%. In addition, more than half of the cells expressed PLAT (52.26%), the marker for late passage BMSCs, but very few cells expressed FGFR2 (0.09%), the marker for early passage BMSCs (Fig. 1c).

Sequencing data statistics

High quality scRNA-seq data were got (Table 1). The number of effective cells detected was assessed as 6,514 (Fig. 2).

Table 1 . Summary statistics for the sequencing data of the BMSCs

Number of cells	6,514
Number of reads	889,767,619
Mean reads per cell	136,593
Median genes per cell	5,969
Fraction reads in cells	81.7%
Valid barcodes	97.1%
Q30 bases in barcode	97.9%
Q30 bases in RNA read	93.2%
Q30 bases in UMI	97.9%

Figure 2. The barcode on the left of the green marker line is the valid cell, and the Barcode on the right is the background. The figure on the right of the marker line has a steep downward trend, indicating that the effective cells are well differentiated from the background. The intersection point between the marker line and the abscissa is the effective cell number.

Different cell subpopulations

The cells that have similar gene expression profile were clustered. The 6,514 cells were divided into three clusters through the k-means clustering approach. The axes correspond to the 2-dimensional embedding produced by the t-SNE algorithm (Fig. 3). The number and corresponding proportion of cells in clusters 1 to 3 were 3,766 (57.81%), 1,720 (26.40%), and 1,028 (15.78%), respectively. Meanwhile, the top 20 DEGs of each cluster were obtained (Table 2), and the clustering effect of them in each cluster was ideal (Fig. 4). The expression of the top 20 DEGs of three clusters in all cells was revealed in Supplementary Fig. S2∼4.

Table 2 . Log2FC and p-value corresponding to the top 20 DEGs of the three clusters

Cluster 1				Cluster 2				Cluster 3
Gene name	log2FC	p		Gene name	log2FC	p		Gene name	log2FC	p
CDCA5	1.48	2e-08		RP11-400N13.3	2.19	1e-12		ACTA2	2.16	8e-12
FAM64A	1.43	7e-08		TGM2	1.83	4e-11		CLU	2.11	8e-12
MYBL2	1.31	9e-07		NEAT1	1.72	2e-10		PTGDS	1.95	2e-08
CENPA	1.27	3e-06		RP11-47I22.3	1.62	2e-09		AKR1C3	1.92	2e-10
PAQR4	1.27	3e-06		MXD4	1.53	1e-08		TUBA1A	1.91	3e-11
HIST1H4C	1.26	4e-06		TNS1	1.50	6e-08		CYBA	1.84	5e-09
TCF19	1.22	7e-06		BAALC	1.43	1e-07		C1GALT1C1	1.83	5e-10
ASF1B	1.22	7e-06		FTL	1.43	1e-07		UBC	1.82	3e-10
H2AFX	1.21	7e-06		PNRC1	1.33	2e-06		LUM	1.81	2e-09
NR2C2AP	1.20	1e-05		IFITM10	1.28	5e-06		ERP44	1.76	2e-09
HMGB2	1.19	1e-05		TMEM158	1.25	1e-05		MVP	1.72	7e-09
PKMYT1	1.18	2e-05		A1BG	1.25	7e-06		DPT	1.70	5e-08
CHAF1A	1.18	2e-05		COL11A1	1.23	1e-05		TIMP1	1.64	3e-08
MKI67	1.17	2e-05		FTH1	1.22	9e-06		HSPA5	1.60	1e-07
CEP55	1.16	2e-05		MYEOV	1.22	2e-05		UTP3	1.58	2e-07
TMPO	1.16	2e-05		RHOQ	1.19	3e-05		PTX3	1.55	1e-06
MND1	1.15	3e-05		THBS2	1.16	7e-05		BLVRA	1.53	7e-07
AURKA	1.15	3e-05		ADAMTS2	1.16	5e-05		ANXA1	1.52	6e-07
PSMC3IP	1.14	3e-05		HIST1H2BK	1.15	6e-05		KPNA2	1.51	1e-06
GINS2	1.13	3e-05		COL5A1	1.13	6e-05		PSMD6	1.50	1e-06

Figure 3. T-SNE projection of cells colored by clustering. Each dot in the figure represents a cell. The legend on the right is the name of the cluster. Different colors represent different clusters.

Figure 4. Clustering heat map of top 20 DEGs corresponding to the three clusters. Each column in the figure represents a gene and each row represents a cluster. Each small area color in the heat map represents the average expression level of the corresponding gene in the cluster. As shown in the right legend, the color from red to blue represents the gene expression from high to low. Different clusters are distinguished by different colors in the first line.

GO analysis

The top 20 DEGs corresponding to the three clusters were analyzed by GO analysis, and we found that the top 20 DEGs of different clusters were significantly enriched with different molecular functions, biological processes and cellular components (Supplementary Fig. S5∼7).

The top 20 DEGs of cluster 1 were significantly enriched in biological processes such as DNA conformation change and DNA packaging; cellular components such as chromosome and histone binding; molecular functions such as histone binding.

The top 20 DEGs of cluster 2 were significantly enriched in biological processes such as collagen fibril organization and endodermal cell differentiation; cellular components such as ferritin complex and fibrillar collagen trimer; molecular functions such as extracellular matrix structural constituent.

The top 20 DEGs of cluster 3 were significantly enriched in biological processes such as prostaglandin biosynthetic process, positive regulation of protein ubiquitination, innate immune response and response to cytokine; cellular components such as extracellular organelle and vesicle; molecular functions such as protein binding and enzyme binding.

BMSCs show considerable potential in clinical applications in regenerative medicine, especially because they are easily accessible and hold no ethical concerns. They possess the capability for tissue repair and immunomo-dulatory through cell-to-cell contact concomitant with paracrine secretion of large amounts of bioactive macromolecules (1). Various studies suggest that autologous BMSCs transplantation can be used to treat diseases such as myocardiopathy, osteoarthritis, and nerve injury (2, 16). In addition, the patients requiring this therapy are usually elderly. In this study, we successfully isolated and cultured BMSCs from an 85-year-old patient with hip fracture, which were consistent with the characteristics of BMSCs in terms of cell morphology, growth characteristics and phenotype.

Many studies demonstrated that invasive methods, such as bone marrow puncture, were used to collect bone marrow from the iliac crest, which not only increased the patients’ pain, but also increased the risk of complications such as infection (17). The volunteer in this study was an elderly patient with intertrochanteric fracture who was treated with intramedullary nailing. To avoid additional damage, we extracted bone marrow from the fractured femoral bone marrow cavity and successfully separated BMSCs. BMSCs from the elderly often show reduced proliferation and differentiation (18). In our study, when the cells were passed to the third generation, there was a clear decrease in proliferation. Therefore, we finally chose the BMSCs at passage 2 for the experiment.

By using flow cytometry, we found that nearly all the BMSCs expressed CD29, CD44, CD90, and CD105, but very few cells expressed CD34, CD45 or HLA-DR. We detected the expression of these markers by Loupe Browser and obtained similar results. On the one hand, these BMSCs fulfill the criteria suggested by the ICST (1); on the other hand, it shows that scRNA-Seq quantification of BMSCs populations corresponds with immunophenotyping by flow cytometry. At the same time, it demonstrates the reliability of the sequencing technology. Similarly, Oetjen et al. (19) used scRNA-Seq to evaluate over 76,000 cells from 20 healthy adult human donors, all the major subpopulations of bone marrow mononuclear cells were identified, and they found overall subpopulation frequencies were resemble to flow cytometry of the matched samples.

The BMSCs cultured in vitro may contain cell types at different stages of differentiation (20). Therefore, they are highly heterogeneous, with different morphology, phenotype, and differentiation ability (3, 4). In addition, they may be composed of different subpopulations of cells with different biological functions. However, studies on the biological characteristics of BMSCs subpopulations are very limited, possibly because there is no appropriate analysis method, or the cells cannot be classified into reasonable subpopulations. ScRNA-seq is a new method to study the heterogeneity among cells. It can amplify and sequence the transcriptome at the level of a single cell, and then divide the cell population into different clusters by gene expression difference among cells, so that cell type characteristics can be more fully defined. However, there are very few studies reporting this technology to analyze human BMSCs. Liu et al. (13) analyzed different passages of BMSCs from a healthy donor, and they described the changes that occurred in different passages of BMSCs. Their study demonstrated that FGFR2 was likely to be an important marker of early passage BMSCs. PLAT is a key gene associated with a subpopulation more prevalent in late passage BMSCs and that it is a potential marker for BMSC senescence. However, most of the early passage (passage 2) BMSCs probed in our study, expressed PLAT rather than FGFR2, which might mean that, unlike the BMSCs from young people, the BMSCs from elderly people show signs of aging at early passage. Passage 2 BMSCs may be considered early passage cells for the young people, but late passage cells for the elderly people.

Currently, due to the high cost of scRNA-seq, few large sample studies have been conducted. To overcome this shortcoming, we examined 6,514 cells and obtained a set of high-quality sequencing data (Table 1). Similarly, Zhang et al. (21) obtained valuable experimental results by sequencing 5,330 umbilical cord mesenchymal stem cells. The 6,514 BMSCs were divided into three clusters by scRNA-seq and bioinformatics analysis, and the gene expression profile and function of the cells in the same cluster were similar. By performing GO analysis to the top 20 DEGs in each cluster, we found that the biological functions of the three clusters were different, and preliminarily believed that cluster 1 was related to cell proliferation, cluster 2 was related to cell differentiation, and cluster 3 was related to various factors secreted by cells, which also corresponded to the characteristics of BMSCs. Actually, there were more than 20 DEGs in each cluster. In order to facilitate the study, we selected the top 20 DEGs for analysis, which was also the method commonly used in the previous studies (13). The function of the top 20 DEGs of each cluster was checked by NCBI gene database to further verify our hypothesis. And many genes’ function has been experimentally verified by scholars. Some functions of the top 20 DEGs in each cluster are as following:

Functional analysis of some genes in cluster 1: Cell division cycle associated 5 (CDCA5) is a vital regulator of sister chromatid cohesion and separation, and it is involved in the generation and maintenance of sister chromatid cohesion (22). FAM64A is a mitotic regulator which can foster cell metaphase-anaphase transition (23). MYBL2 is a physiological regulator of cell cycle progression, cell survival, and cell differentiation (24). The integrity of the human centromere DNA repeats is protected by CENP-A (25). PAQR4 contributes to modulation of cell proliferation and knockdown of PAQR4 leads to reduction of cell proliferation (26). During DNA replication, the histone chaperone Asf1b is important for Importin-histone recognition (27). HMGB2 as transcriptional activators through the modulation of chromatin structure (28). Chromatin assembly factor-1 (CAF-1) complex can adjust chromatin assembly in DNA replication and repair (29). It can be seen that many DEGs in cluster 1 are related to cell self-renewal.

Functional analysis of some genes in cluster 2: TGM2 is involved in the regulation of chondrogenesis (30). NEAT1 is thought to promote osteogenic differentiation in BMSCs by regulating miR-29b-3p/BMP1 axis (31). COL11A1 plays an essential role during the cartilage development process (32). Type V collagen might be considered as a stromal component that impairs osteogenesis (33). It can be seen that many top 20 DEGs in cluster 2 are related to multidirectional differentiation of BMSCs.

Functional analysis of some genes in cluster 3: Cyba gene encoding p22^phox are essential for microbial killing and innate immunity (34). TIMP-1 is a key effector molecule responsible for the antiangiogenic effects of MSC, which is a mechanism to tune down immunity and inflammation (35). Annexin A1 (ANXA1) is an immune-modulating protein with various functions in the immune system (36). Lumican (LUM) affects cell migration and promotes wound repair (37). Dermatopontin (DPT) possesses significant pro-angiogenic properties and plays vital roles in cutaneous wound healing (38). In addition, some genes are involved in regulating cytokine secretion. ERp44 is a pH-regulated chaperone of the secretory pathway and an essential regulator of protein secretion at the ER-Golgi interface (39). The HSPA5 gene encodes the binding immunoglobulin protein is a necessary component of the translocation machinery for protein import into endoplasmic reticulum (40). It can be seen that proteins encoded by some top 20 DEGs in cluster 3 can be used as secretory factors to participate in immune regulation, damage repair, and other processes.

Liu et al. (13) analyzed 36 BMSCs from passage 3 by scRNA-seq and identified three BMSCs subpopulations, too. In their study, the BMSCs were separated into subpopulation A, B, and C. Subpopulation A was characterized by the strong expression of FGFR2. Since FGFR2 was involved in osteogenesis and more likely to be expressed by early passage cells, this subpopulation might include skeletal stem cells. Subpopulation B expressed higher levels of FGF5. FGF5 is a secreted heparin-binding growth factor, which could increase osteogenic differentiation of MSCs. This subpopulation may represent differentiating skeletal stem cell progenitors. Subpopulation C was characterized by strong expression of PLAT and VCAM1. Because VCAM1 promoted angiogenesis, this subpopulation may be involved with angiogenesis. Although both studies divided cells into three subpopulations, each subpopulation corresponds to different characteristic genes. The reasons are multifaceted. Firstly, we studied the elderly, while they studied relatively young individuals. Secondly, it was passage 2 BMSCs that we studied，which was different from the passage 3 BMSCs they studied. And their study proved that BMSC subpopulations change with passage. In addition, they grouped only 36 cells, perhaps because of different sequencing techniques.

This study was limited by the fact that the BMSCs from a single subject. Whether the top 20 DEGs of each subpopulation can be used as its characteristic gene is still up for discussion. With the further work, such as expanding sample size and functional verification, the characteristic genes of each subpopulation will be identified. The verification of DEGs and the sorting of cell subpopulations by these genes will be the next work. In the future, we may analyze each subpopulation in more detail, and further explore the relationship between DEGs and cell phenotypes.

In conclusion, this study demonstrated that BMSCs with high quality could be isolated and cultured from the bone marrow of the affected side femur of elderly hip fracture patients. The BMSCs could be divided into three subpopulations by using scRNA-seq, and the functions of subpopulations were related to self-renewal, multilineage differentiation, and cytokine secretion, respectively. The heterogeneity of BMSCs has been demonstrated. More studies including BMSCs from a wider range of subjects are needed.

Supplementary data including seven figures can be found with this article online at https://doi.org/10.15283/ijsc21042.

This work was supported by the special project for innovation of military medical research plan (16CXZ002) and the Beijing Nova Program (xx2018021).

The authors have no conflicting financial interest.