{'fraction_non_alphanumeric': 0.03135732417500373, 'fraction_numerical': 0.017370962122343337, 'mean_word_length': 5.106990881458967, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Single-Cell RNA Sequencing of Bone Marrow Mesenchymal Stem Cells from the Elderly People

Background and Objectives Bone marrow mesenchymal stem cells (BMSCs) show considerable promise in regenerative medicine. Many studies demonstrated that BMSCs cultured in vitro were highly heterogeneous and composed of diverse cell subpopulations, which may be the basis of their multiple biological characteristics. However, the exact cell subpopulations that make up BMSCs are still unknown. Methods and Results In this study, we used single-cell RNA sequencing (scRNA-Seq) to divide 6,514 BMSCs into three clusters. 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%). The gene expression profile and function of the cells in the same cluster were similar. The vast majority of cells expressed the markers defining BMSCs by flow cytometry and gene expression analysis. Each cluster had at least 20 differentially expressed genes (DEGs). We conducted Gene Ontology enrichment analysis on the top 20 DEGs of each cluster and found that the three clusters had different functions, which were related to self-renewal, multilineage differentiation and cytokine secretion, respectively. In addition, the function of the top 20 DEGs of each cluster was checked by the National Center for Biotechnology Information gene database to further verify our hypothesis. Conclusions This study indicated that scRNA-Seq can be used to divide BMSCs into different subpopulations, demonstrating the heterogeneity of BMSCs.


Introduction
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.

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 microfluidics. 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).
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/).

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).

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.

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 ubiquiti- nation, innate immune response and response to cytokine; cellular components such as extracellular organelle and vesicle; molecular functions such as protein binding and enzyme binding.

Discussion
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 immunomodulatory 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 pa-tients' 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 pro-liferation 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 immunemodulating 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 dif-ferentiating 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.