• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Spatially defined single-cell transcriptional profiling characterizes diverse chondrocyte subtypes and nucleus pulposus progenitors in human intervertebral discs

    2021-10-27 04:06:00YiboGanJianHeJunZhuZhengyangXuZhongWangJingYanOuHuZhijieBaiLinChenYangliXieMinJinShuoHuangBingLiuandPengLiu
    Bone Research 2021年3期
    關(guān)鍵詞:時需敵手會話

    Yibo Gan,Jian He ,Jun Zhu,Zhengyang Xu,Zhong Wang,Jing Yan,Ou Hu,Zhijie Bai,Lin Chen,Yangli Xie,Min Jin,Shuo Huang,Bing Liu,5,6?and Peng Liu ?

    INTRODUCTION

    Degenerative disc disease(DDD)is regarded as the primary cause of low back pain,resulting in a global healthcare burden and significant socioeconomic costs.1It may lead to a severe impact on the quality of life of patients.2The current treatment of DDD,mainly including bed rest,rehabilitation,medication,interventional therapy,and surgery,3provides only symptomatic relief but fails to reestablish the homeostasis of the intervertebral disc(IVD).4Furthermore,the deterioration of the health of the compromised spine cannot be prevented.5Thus,the unrelenting threat posed by DDD to human health has motivated the search for an increased understanding of human IVD physiology and pathology.

    The IVD has a well-confined structure,including three components:the central hydrated nucleus pulposus(NP),the surrounding lamellar annulus fibrosus(AF),and the cartilage endplate(CEP)that is adjoining to the vertebra.6The confined structure of the IVD plays a part in the mechanical function.7Unfortunately,alterations in the cellular composition and microenvironment cause the IVD to undergo a slow but relentless program that causes the confined structure to be compromised during the degenerative process.8–10The origin of the IVD is heterologous,where the NP is believed to be derived from the notochord,11–12and the AF and CEP are derived from the sclerotome.13–14Consequently,the cells in the IVD are also heterogeneous,composed of NP cells,and notochord cells in the NP,AF cells in the AF,and chondrocytes in the CEPs.15However,classification based on spatial location cannot uncover the highly heterogeneous cell populations in regard to phenotype and function.Although previous studies have revealed phenotypes of IVD cells by bulk RNA sequencing,16–18the search for molecular mechanisms underlying degeneration has been complicated by the large amount of heterogeneity in cellular compositions and the subsequently highly complex cellular microenvironment of the IVD.To further examine the cellular heterogeneity,some efforts were made to distinguish the critical cell types in IVD.The hypothesis regarding cellular heterogeneity in the IVD was initially supported by Hunter CJ et al.,as evidenced by the existence of large vacuolar notochordal cells in the NP and small rounded chondrocytes.12,19Notochordal cells are thought to disappear starting in adolescence in the human IVD,20–21which has been questioned becausebrachyury(TBXT),a notochord lineage marker,continued to be expressed in the IVD.22Thus,notochord cells are thought to be the precursors of all NP cells regardless of variations in morphology and size at different stages.23In addition,mesenchymal stem cells(MSCs)are thought to exist in the IVD due to the expression of the MSC markersENG(CD105),CD44,THY1(CD90),NT5E(CD73),andNGFR(CD271).24–25NP progenitor cells are characterized by clonogenicity,pluripotency,and NP reorganization properties.26However,the different lineages remain largely unknown due to the lack of high-precision and unbiased resolution for distinguishing cell populations in the human IVD,although its importance is widely acknowledged.

    Single-cell RNA sequencing(scRNA-seq)is considered as a powerful tool for resolving cellular heterogeneity and hierarchical factors forming a complicated cell niche.27–28Here,we performed scRNA-seq to obtain an unbiased picture of IVD cell populations.Our findings provide a better understanding of the inherent heterogeneity and reshape the existing classifications of chondrocytes in the IVD.Notably,we also confirmed the existence of progenitor cells in the IVD marked byPDGFRAandPROCR.Thus,our study reveals the cellular landscape of the human IVD and provides insights that could help to identify therapeutic targets for human DDD.

    RESULTS

    Comprehensive scRNA-seq analyses resolve the major cell types in the human IVD

    To determine the cellular composition of the human IVD,we employed droplet-based single-cell transcriptomic profiling(10X Genomics Chromium System)of cells from the NP,AF,and CEP from five healthy human IVDs(Pfirrmann I)(Fig.1a and Supplementary Table 1),as evaluated by magnetic resonance imaging(MRI)according to the Pfirrmann grading system29(Supplementary Fig.1a).The integrity of the IVD was confirmed because the sagittal cross-section showed that it met the criteria of high hydration and ordered organization with increased deposition of chondroitin sulfate based on hematoxylin & eosin and safranin O/fast green staining(Supplementary Fig.1a).Because it was difficult to distinguish the boundary between the NP and inner AF,we harvested gelatinous tissue from the central region as the NP.Thus,the tissue origins of harvested cells were identified clearly due to the strict criteria of sampling.A total of 128 833 individual human IVD cells were profiled,and 108 108 cells were retained for subsequent analysis after rigorous quality control and doublet exclusion(Supplementary Table 1).The resulting cells were sequenced to a median depth of 5 367 unique molecular identifiers(UMIs)per cell,with a median of 1 569 genes detected per cell(Supplementary Fig.1b and Supplementary Table 1).Similarities between samples determined by Pearson’s correlations and the sequencing depth suggested that all samples were comparable(Supplementary Fig.1c).

    We performed fastMNN30to correct batch effects among different data sets.Unbiased clustering based on t-distributed stochastic neighbor embedding(tSNE)identified nine putative root clusters in the healthy human IVD(Fig.1b and Supplementary Fig.1d),including(1–3)three clusters ofSOX9+chondrocytes(Chond1,Chond2,and Chond3);(4)notochord cells;(5)stromal cells;(6)pericytes;(7)endothelial cells(ECs);(8)nucleus pulposus progenitor cells(NPPCs);and(9)blood cells.The chondrogenic marker geneSOX9+and chondrocyte-specific ECM genes(COL2A1andACAN)were ubiquitously expressed in the three chondrocyte clusters(Fig.1c,d).The notochord origin marker gene,TBXT,31was dominantly expressed in the notochord cell cluster,along with notochord-derived cytokeratin genes,such asKRT832(Fig.1c,d).FOXC2,GJA1,andHES4,which are essential for stromal cell differentiation,33–35that were mainly expressed in the stromal cell cluster.Pericyte and EC clusters were identified by feature gene expression(ACTA2,TAGLN,andMCAMfor pericytes36–38andPECAM1,CD34,CDH5,ERG,andVWFfor EC39–42)(Fig.1c,d and Supplementary Table 2).We found thatPDGFRA(the mesenchymal progenitor marker),43PRRX1(which is restricted to the mesodermal origin and regulator of mesenchymal precursors)44andIGF1(a growth factor that effectively differentiates MSCs into NP-like cells)45were specifically expressed in the NPPC clusters(Fig.1c,d).Thus,we speculated thatPDGFRA+NPPCs could be a mesoderm-derived progenitor cell cluster in the IVD.A total of 2 651 differentially expressed genes(DEGs)were identified that distinguished human IVD cell populations(Fig.1d and Supplementary Table 2).Spatially,chondrocytes and stromal cells were abundant in the NP,CEP,and AF,while notochord cells were mainly found in the NP(Fig.1e and Supplementary Fig.1e).The expression of some widely reported marker genes of the IVD was also detected in these cell populations(Supplementary Fig.1f).We then performed an immunohistochemistry assay to validate the spatial distribution of major cell types(Fig.1f).We found that most SOX9+chondrocytes were detected in the NP,AF,and CEP,as expected.PDGFRA+NPPCs were mainly distributed in the NP and rarely found in the AF and CEP.ACTA2+pericytes and PECAM1+ECs were sporadically distributed in the NP and were present in the tube-like CEP,in line with previous findings on capillaries in the CEP.46Moreover,immunofluorescence staining of the human IVD(Pfirrmann I and II)validated the presence of scattered PECAM1+CD34+cells and ACTA2+cells in the IVD(Supplementary Fig.1g).

    Pairwise correlation analysis clearly distinguished the chondrocyte and nonchondrocyte subsets(Supplementary Fig.1h).Gene ontology(GO)analysis revealed distinct functional enrichment in these cell types(Supplementary Fig.1i).For example,Chond1 was enriched for signaling regulation and stimulus-response,while Chond2 was enriched for ECM synthesis and organization.As expected,pericyte and ECs were enriched for genes involved in regulating vasculature development,cell adhesion,and junctions.Interestingly,the NPPC cluster was enriched for terms that regulated skeletal development and ossification.

    Fig.1 Single-cell transcriptomic landscape of human intervertebral disc(IVD)cells.a Schematic workflow of the experimental strategy.Cells isolated from the NP,AF,and CEP of the human IVD were subjected to droplet-based scRNA-seq.NP nucleus pulposus,CEP cartilage endplate,AF annulus fibrosus,IVD intervertebral disc,scRNA-seq single-cell RNA sequencing.b Distribution of 108 108 cells from human intervertebral disks.Eight cell clusters were visualized by a tSNE plot.Cell numbers for each cluster are indicated in brackets.NPPC nucleus pulposus progenitor cell,EC endothelial cell,Chond chondrocyte,tSNE t-distributed stochastic neighbor embedding.c The average expression of curated feature genes for cell clusters defined in b on the tSNE map.d Heatmap revealing the scaled expression of DEGs for each cell cluster.DEGs differentially expressed genes.e Fraction of cell clusters in the NP,CEP,and AF.f Representative immunohistochemistry staining of signature markers of the indicated cell clusters in the AF,NP,and CEP of healthy human IVD tissues and quantification of positive cells displayed with a box plot(n=3).Scale bar,100 μm

    To validate the conserved cell heterogeneity of the IVD across species,we compared the transcriptome of IVD cells between humans and rats by reanalyzing the scRNA-seq data from a recent rat study.47As expected,most of the cell clusters identified in the human IVD were also found in the rat IVD and showed gene expression conservation across cell types,including NPPCs,ECs,and pericytes(Supplementary Fig.2a,b).In particular,NPPCs in rats also highly expressedPDGFRA,PRRX1,andIGF1and shared distinct gene expression patterns with their counterparts in humans(Supplementary Fig.2c,d).

    Overall,these results revealed the cellular diversity in the human IVD,and we identified a set of markers that can potentially be used to recognize the cell clusters in the human IVD.

    The functional definition of chondrocyte subpopulations in the IVD

    As chondrocytes are known to play a pivotal role in ECM homeostasis and the degeneration of the IVD,48we sought to determine their composition.Each of the three chondrocyte clusters was divided into two subclusters(Fig.2a).The distribution of subclusters exhibited apparent distinctions in the three compartments of the IVD(Fig.2b).The subclusters of C1 and C2 were mostly located in the AF and CEP,while C5 was mainly located in the NP.Subclusters of C3,C4,and C6 were relatively evenly distributed in the NP,AF,and CEP.A total of 912 DEGs were found among the six chondrocyte subclusters(Fig.2c and Supplementary Table 3).We found that C1 preferentially expressed growth factor(GF)genes such asBMP2,TGFB1,andFGF2.Subclusters C3 and C4 preferentially expressed the genes of the main ECM components of the IVD,such asACANandCOL2A1.Subclusters C5 and C6 preferentially expressedPRG4andCNMD,suggesting that they may play a protective role and stabilize the chondrocyte phenotype.49–50

    Fig.2 Characterization of chondrocytes in the human IVD.a tSNE plot of the six subclusters of 93 495 chondrocytes defined in the IVD.b Fraction of each chondrocyte subcluster in the NP,CEP,and AF.c Heatmap revealing the scaled expression of DEGs for each chondrocyte subcluster.d Heatmap showing pairwise Pearson correlations in the global transcriptome between IVD chondrocytes and articular chondrocytes(Ji et al.,51).FC fibrocartilage chondrocyte,HomC homeostatic chondrocyte,HTC hypertrophic chondrocyte,preHTC prehypertrophic chondrocyte,ProC proliferative chondrocyte,RegC regulatory chondrocyte,EC effector chondrocyte.e Dot plot showing the mean expression of selected chondrocyte function-associated genes among the six chondrocyte subclusters.Dot size indicates the percentage of cells in subclusters with detected expression.f The fraction of each chondrocyte subcluster arrested in the different cell-cycle phases.g Radar map showing the performance of six gene sets associated with the indicated function and metabolic pathway among each chondrocyte subcluster.h Heatmap showing pairwise Pearson correlations of expressed matrisome genes in chondrocytes.Two signature patterns(matrisome-associated and core matrisome)were identified by hierarchical clustering.i The number of expressed genes associated with six matrisome patterns in each chondrocyte subcluster.ECM extracellular matrix.j Violin plots showing the expression levels of representative genes associated with six matrisome patterns in each chondrocyte subcluster

    To better understand the specific characteristics of IVD chondrocytes,we compared the transcriptomic differences between these chondrocyte subclusters and articular cartilage chondrocytes at different stages of osteoarthritis(stages 0-4)(Fig.2d).51There was no characteristic correspondence among subclusters C1,C2,and articular chondrocytes.We found that subclusters C1 and C2 shared some DEGs with articular regulatory chondrocytes(RegCs),includingCKS2andHMOX1(Fig.2e),and highly expressedIBSP and CYTL1,which are NP-negative biomarkers(Supplementary Table 3).17,52Moreover,subclusters C1 and C2 showed a higher percentage of cells arrested in the G2/M phase than that in others,which was indicative of relatively higher proliferative activity(Fig.2f).We also evaluated these subclusters using a gene set related to chondrocyte function(Fig.2g).The results showed that oxidative phosphorylation played a role in the metabolic pattern of C1 and C2,which could be explained by the fact that C1 and C2 were mainly located in the vascularized AF and CEP.Notably,subcluster C1 was enriched for genes related to chondrogenic differentiation.Therefore,we hypothesized that subclusters C1 and C2 represent regulatory chondrocytes that stimulated surrounding cells by secreting GFs.

    Pairwise correlation analysis revealed close relationships among C3,C4,and homeostatic chondrocytes(HomCs,Fig.2d)with a similar pattern of gene expression as articular HomCs,such asCCNL1andWSB1(Fig.2e).51In contrast with regulatory chondrocytes C1 and C2,fewer cells in C3 and C4 were arrested in the G2/M phase(Fig.2f).In particular,both C3 and C4 exhibited strong enrichment of circadian regulation genes and moderate enrichment of chondrogenic differentiation(Fig.2g).C3 was also enriched for cellular adhesion genes,which are critical for forming chondrocyte clonal columns within an ordered,three-dimensional cell array.53Considering that they preferentially expressed ECMrelated genes,we chose to classify C3 and C4 as homeostatic chondrocytes,which function in maintaining ECM homeostasis and circadian rhythm.

    The C6 subcluster was relatively similar to hypertrophic chondrocytes(HTCs)and prehypertrophic chondrocytes(preHTCs,Fig.2d).COL5A1andEPYCwere expressed in the subclusters of C5 and C6(Fig.2e),similar to articular HTCs and preHTCs.51We also found that C5 and C6 highly expressed genes reflecting protective characteristics(KLF2andCHI3L1)(Fig.2e)54–55and existed in the dormant stage of proliferation(Fig.2f).Interestingly,subclusters C5 and C6 preferentially performed metabolic processes,including oxidative phosphorylation and glycolysis,showing the traits of a high metabolism(Fig.2g)and the characteristics of articular effector chondrocytes.51Unlike the resident quiescent chondrocytes with low metabolism,56these subclusters were possibly adapted to anaerobic metabolism because C5 was mainly located in the NP,which has an avascular and hypoxic microenvironment,consistent with previous study showing that the NP that is predominantly glycolytic due to vigorousHIF1activity.57Collectively,we inferred that C5 and C6 were effector chondrocytes with high metabolic rates and protective/repair functions.

    To reveal the core function of chondrocytes in modulating ECM homeostasis,we detected the expression of matrisomerelated genes.Matrisome genes were categorized into the core matrisome(collagens,proteoglycans,and ECM glycoproteins)and matrisome-associated(ECM regulators,ECM affiliation,and secreted factors)according to a matrisome classification database(matrisomeproject.mit.edu).58We first evaluated the average expression of six modules in eight clusters(Supplementary Fig.3a)and compared the expression of matrisome genes distinctly expressed in the NP,CEP,and AF(Supplementary Fig.3b and Supplementary Table 4).Correlation analysis of matrisome-related genes in chondrocytes revealed two patterns:the core matrisome and matrisome-associated(Fig.2h).To clarify the primary function of matrisome-related gene subsets in six chondrocyte subclusters,we compared the expression abundance of these genes(Fig.2i,j).We found that secreted factors were predominantly expressed in the regulatory C1 subset,while the homeostatic C3 subset preferentially expressed genes of the core matrisome(Fig.2i,j).In contrast,the effector C5 subset exhibited high expression of ECM regulators,reflecting its regulatory role in ECM homeostasis(Fig.2i,j).

    Taken together,these data add to the knowledge on the functions of chondrocyte subclusters in human IVD.

    Delineating nucleus pulposus progenitor cells and their signature genes

    NP progenitor/stem cells are critical in the physiological and pathological processes of the IVD.59–60We identified NPPCenriched genes related to bone development,bone morphogenesis,connective tissue development,and endochondral bone growth(Supplementary Fig.1i).To better understand the role of the NPPC cluster,we sought to determine their composition in the human IVD.We partitioned NPPCs into four subclusters(Fig.3a).The localization of the discogenic markerPAX1confirmed the physical presence of the NPPC-1 subcluster.PAX1is expressed in the sclerotome,which is critical for the formation of vertebrae and IVDs,61indicating the potential role of discogenic differentiation in NPPCs.Subcluster NPPC-2 specifically expressedANGPT1,which is critical for the survival of nucleus pulposus cells.26PRG4,the signature gene of NPPC-3,was also highly expressed in articular cartilage progenitor cells.62SOX9expression indicated the chondrogenic priming of NPPC-4(Fig.3b).These NPPC subclusters were also distinguished by the indicated DEGs(Supplementary Fig.4a and Supplementary Table 5).GO analysis of these DEGs showed that NPPC-1 and NPPC-3 were enriched for genes regulating ECM organization,while NPPC-4 was enriched for genes involved in mRNA catabolic metabolism(Supplementary Fig.4b).Gene set enrichment analysis(GSEA)showed that NPPC-1 was enriched for the calcium signaling pathway,which played a vital role in modulating NP homeostasis by regulatingAQP2.63NPPC-2 was enriched for the MAPK signaling pathway,potentially playing a protective role in cell survival in the NP.64NPPC-3 preferentially expressed theSMAD2/3pathway,and NPPC-4 was enriched forNOTCHsignaling,which plays a role in cell growth(Supplementary Fig.4c).65

    To explore the regulatory networks that determine cell fate specification in the NPPC subclusters,we utilized single-cell regulatory network inference and clustering(SCENIC)to infer the regulatory activity(regulon)from the coexpression of transcription factors(TFs)and their downstream target genes.66We filtered 21 core regulons out of 227 regulons that were used to discriminate the four NPPC clusters(Fig.3c,Supplementary Fig.4d,and Supplementary Table 6).The highly enriched regulons in NPPC-1 includedHOXA10andHOXA7.TheSOX4,RARA,andMEIS1regulons were specific to NPPC-2.NPPC-3 exhibited strong enrichment ofZFP14andSMAD3.NPPC-4 was enriched for regulons such asGLI1,EGR2,andNR2F1(Fig.3c and Supplementary Table 6).Some important regulons,includingHOXA10,SOX4,SMAD3,andGLI1,together with their downstream target genes,such as the abovementionedPAX1,PRG4,andANGPT,had the potential to regulate the function of NPPCs(Fig.3d and Supplementary Table 5).Specifically,HOXA10is a critical regulator of osteogenesis.67SOX4is highly expressed in osteoblast progenitors,and its expression is increased during osteoblast differentiation.68GLI1marks mesenchymal progenitors responsible for bone formation and fracture repair and regulates chondrocyte differentiation.69SMAD3,the downstream target ofTGF-β,plays a dominant role in chondrogenesis and maintaining the phenotype of chondrocytes.70

    Fig.3 Characterization of NPPC in human IVD.a Four subclusters of 2 157 NPPCs were visualized by a tSNE plot.b tSNE plot of signature gene expression in NPPCs.c Heatmap revealing binary regulon activities analyzed with SCENIC in each subcluster of NPPCs.“ON”indicates active regulons;“OFF”indicates inactive regulons.SCENIC Single-Cell Regulatory Network Inference and Clustering.d The HOXA10,SOX4,SMAD3,and GLI1 regulon networks in NPPC subclusters.The TFs are in dark red,and the corresponding target genes are in light green.The line thickness indicates the level of GENIE3 weights.The dot size indicates the number of enriched TF motifs.TFs transcription factors,GENIE3 GEne Network Inference with Ensemble of trees.e Dot plot showing differentially expressed genes encoding surface markers in each subcluster of NPPCs.f Violin plots showing the expression levels of the signature genes of NPPC-3 in the IVD.g Immunofluorescence staining showing the coexpression of PDGFRA,PROCR,and PRRX1 in human IVD cells in situ(n=3).Scale bar,5 μm.h Flow cytometry gating strategies for sorting PDGFRA+PROCR+in the human IVD.i Representative crystal violet staining of CFU-F colonies generated by sorted primary PROCR+cells of the human IVD(left,n=3).Scale bar,5 mm.Quantification of the number of CFU-F colonies(right).The statistical significance of differences was determined using one-way ANOVA with multiple comparison tests(LSD).**P<0.000 1.Error bars indicate the SEM.CFU-F colony-forming unit-fibroblast,ANOVA Analysis of Variance,LSD least significant difference,SEM Standard Error of the Mean.j Immunofluorescence staining of SMAD3 and p-SMAD3 in the PROCR+and PROCR?cells of the human IVD(n=3).Scale bar,40 μm

    To immunophenotype these NPPC subclusters,we screened for cell surface marker genes that were differentially expressed among the four NPPC subclusters.Among them,PDGFRAshowed higher expression in NPPC-1,NPPC-2,and NPPC-3 than in NPPC-4.Interestingly,we found that NPPC-3 preferentially expressedPROCR(Fig.3e),a widely reported signature gene for progenitor cells in multiple organs,including the hematopoietic and vascular systems,71–74pancreas,75ovaries,76etc.Thus,the specific expression ofPROCRsuggested the potential stemness capacity of NPPC-3.We then combined the expression of the membranous marker genesPDGFRAandPROCRand the transcription factorPRRX1as a signature for the identification of NPPC-3(Fig.3f and Supplementary Fig.4e)and performed immunofluorescence staining to examine their coexpression(Fig.3g).Immunostaining of a healthy IVD(Pfirrmann I,Supplementary Table 1)showed that PDGFRA+PROCR+NPPCs were mainly located in the NP zone(Supplementary Fig.4f).To assess the proportion of PDGFRA+PROCR+NPPCs in the NP,primary PDGFRA+PROCR+cells were flow cytometrically sorted from the human IVD(Pfirrmann II,Supplementary Table 1).The results showed that the frequency of PDGFRA+PROCR+cells was 0.36%(Fig.3h),and PDGFRA was enriched in almost all PROCR+cells in the IVD.To test the clonogenicity of NPPC-3,primary PROCR+cells were sorted by flow cytometry for a colonyforming unit-fibroblast(CFU-F)colony formation assay.The counts of typical colonies derived from primary PROCR+cells were 25.9±3.3 per 1 000 cells,which was comparable to that for PDGFRA+MSCs77and significantly higher than that for PROCR?cells(2.9±1.4 per 1 000 cells)(Fig.3i),indicating that NPPCs exhibited enhanced colony-formation ability.To verify the in silico finding of the enriched regulatory activity of SMAD3 in PROCR+NPPCs,we detected the expression level of p-SMAD3 in P2 PROCR+and PROCR?cells from the human IVD(Pfirrmann II,Supplementary Table 1).As expected,the expression of p-SMAD3 in the nucleus was higher in PROCR+cells than in PROCR?cells(Fig.3j).These results indicated that SMAD3 was highly activated in PROCR+cells,suggesting the potential role of PROCR+cells in the chondrogenesis of IVD.

    Taken together,these data elucidated the cellular heterogeneity in NPPCs,which was highly regulated and comprised the population with clonogenicity that could be enriched by PDGFRA and PROCR.

    Reconstruction of the bilineage trajectory of PDGFRA+PROCR+NPPCs

    Connective tissue comprised stromal cells with phenotypic and functional complexity,78which provided support during NP development and repair.79We collected 1 372 stromal cells from the NP that were divided into six subclusters(Supplementary Fig.5a and Supplementary Table 7),including three subclusters of fibroblasts(Fib1,Fib2,and Fib3)that expressed high levels of fibroblast signature genes,such asCEMIP,AKR1C1,MGP,COMP,DNER,andMELTF,80–81two subclusters of neurogenic cells(Neu1 and Neu2)with high expression of the neurogenic markersSOX2,NGFR,NCMAP,andCLDN19,82–85and osteogenic cells that expressed high levels of the osteogenic regulatorsRUNX2andDLX5(Supplementary Fig.5b).86–88To further verify the existence of the cell clusters in the IVD,immunofluorescence staining of human IVDs(Pfirrmann I and II)showed a few RUNX+SP7+cells and SOX2+cells in the IVD(Supplementary Fig.5c),consistent with the findings from the scRNA-seq analysis.

    We next sought to investigate the differentiation trajectories that determined the cellular hierarchy in NP cells.All NP cells,including four subclusters of NPPCs,three subclusters of fibroblasts,three subclusters of chondrocytes,and osteogenic cells,were involved in reconstructing the differentiation trajectories using Monocle 3(Fig.4a),an algorithm for the reconstruction of lineage programs based on similarity at the transcriptional level.89We set NPPC-3 as the starting point of the differentiation trajectories due to its high expression of pluripotent genes and progenitor potential identified above,and then computed pseudotime for cells along the inferred developmental axis(Fig.4a,b).More specifically,NPPC-3 was predicted to differentiate into two distinct cell lineages,including the chondrogenic branch,which includes NPPC-3,NPPC-1,Chond2,and Chond3,and the osteogenic branch,which includes NPPC-3,NPPC-2,NPPC-4,osteogenic cells,and Fib3(Fig.4a,b).

    To explore gene expression dynamics along the trajectories,we grouped genes that varied between cell clusters into 16 modules using Louvain community analysis(Supplementary Table 8).A heatmap showed the aggregated expression in each module across cell clusters(Fig.4c).We found that the expression of genes in module 4 was downregulated along both trajectories,such as the differential regulator genesTWIST1andFOXP1,90–91which were enriched for genes related to ECM organization(Fig.4d–f).In contrast,the expression of chondrogenic genes was gradually elevated along the chondrogenic trajectory,such asCOL2A1andACANin module 7,and remained at high expression levels until terminal differentiation(Fig.4d,e),which was also evidenced by the expression of module 9(e.g.,the chondrogenicCNMD and FGFBP2)(Supplementary Fig.5d,e).The expression of the osteogenic gene set was elevated along the osteogenic trajectory,such asRUNX2andALPLin module 11(Fig.4d,e)and the expression ofSP7,BGLAP,andMMP11in module 16(Supplementary Fig.5d,e).

    Fig.4 Reconstruction of bilineage trajectories in NP cells.a UMAP visualization of NPPC,chondrocyte,fibroblast,and osteogenic subclusters.The principal graph of trajectories reported by Monocle 3 was rooted in NPPC-3,as indicated by a cycle.UMAP,Uniform manifold approximation and projection.b Developmental pseudotime for cells present along the trajectory inferred by Monocle 3,with osteogenic and chondrogenic branches coming from NPPC subclusters.c Heatmap showing the scaled mean expression of modules of coregulated genes grouped by Louvain community analysis across the subclusters.d UMAP plots showing the relative expression level of representative gene modules in NP cell subclusters.e Pseudotime kinetics of the indicated genes in the modules in d from the root along the trajectories to chondrogenic and osteogenic differentiation.f Histogram showing the pathways enriched by ReactomePA for each module indicated in d.g Representative alizarin red(top left),oil red O(bottom left),alcian blue(top right),and safranine O/fast green(bottom right)staining after the osteogenic,adipogenic,and chondrogenic differentiation of PROCR+cells(n=3).Magnified images of the boxed areas are shown on the right.White scale bars,400 μm.Black scale bars,200 μm

    Taken together,these data depicted the trajectories of NP cells,in which PROCR+cells were enriched for multipotent NPPCs that generate three lineages,consequently revealing the successive activation of transcriptional programs in NP homeostasis.

    Putative signaling network for the intercellular crosstalk regulating the homeostasis of the NP

    To seek further insights into the critical factors involved in the NP cell niche of the human IVD,we investigated the signaling network among the main cell types in the NP.CellChat analysis of these 14 subclusters in the NP identified the signaling network for intercellular crosstalk.Relative active bidirectional signaling interactions among these cell subclusters revealed highly regulated cellular communications(Fig.5a and Supplementary Table 9).ECs,pericytes,fibroblasts,and neurogenic cells identified as niche components in the NP played distinct roles in signaling interactions to regulate the differential process.To determine the important factors,we further analyzed the intercellular signaling networks of VEGF,TGFB,PDGF,and FGF(Fig.5b–e).

    Interestingly,Fib3 was involved in VEGF signaling,both autocrine and paracrine(Fig.5b,f).ECs were the leading receiver of VEGF signals,as expected,and NPPC subclusters functioned as regulators of the communication(Fig.5b).Moreover,the TGF-β pathway was involved in many signaling interactions among chondrocyte subclusters and NPPC clusters viaTGFB3-TGFBRorTGFB3-ACVR1(Fig.5c,f).As shown above,NPPC-3 was enriched forSMAD3,the key downstream target of TGF-β,which prompted us to further investigate the role of TGF-β3 in chondrogenesis in NPPCs.The results showed that 10 ng·mL?1TGF-β3 effectively induced chondrogenesis and the formation of dense cartilage extracellular matrix(ECM)compared with that in the negative control group after 28 days of differentiation.However,supplementation with 10 μmol·L?1SB505124,a TGF-β receptor inhibitor,blocked the chondrogenesis of PROCR+cells both with and without TGF-β3(Fig.5g).The results demonstrated that the TGF-β family plays an important role in the chondrogenic regulation of PROCR+cells.

    Fig.5 Overview of the crosstalk networks among the clusters in the NP.a Overview of the cellular network regulating the homeostasis of the NP.Dots indicate cell clusters.The dot size indicates the relative quantity of each cluster.The thickness of the directed line indicates the relative quantity of significant ligand-receptor pairs between any two pairs of cell clusters.b-e Circle plots showing the inferred VEGF(b),TGFβ(c),PDGF(d),and FGF(e)signaling networks.f Dot plot showing the communication probability of the indicated ligand-receptor pairs between EC,Pericyte,Neu,and Fib3 subclusters(sending signals)and four NPPC subclusters(accepting signals).g Representative alcian blue and toluidine blue staining for the chondrogenic effect of TGF-β3 supplementation(10 ng·mL?1)for 28 days on PROCR+cells from the human IVD.Scale bars,400 μm.h Histogram showing the proliferation of PDGF-AA(20 ng·mL?1)on PROCR+cells from the human IVD detected by a CCK-8 kit(n=3).The statistical significance of differences was determined using one-way ANOVA with multiple comparison tests(LSD).**P<0.01 compared with the control group at day 10.##P<0.01 compared with the PDGF-AA group at day 10.Error bars indicate the SEM

    In the PDGF signaling network,the NPPC clusters acted as critical contributors by secreting PDGFA ligand,leading to the paracrine activity of NPPCs to osteogenic cells,pericytes,and Fib1 and the autocrine activity of NPPCs to themselves.Specifically,NPPC-3 was the key population that dominated the PDGF signaling network(Fig.5d,f).Previous studies have reported that PDGF-AA is involved in the regulation of cell proliferation.92–93Therefore,we explored the effect of PDGF-AA on the proliferation of PROCR+cells from the human IVD(Fig.5h).The results showed that 20 ng·mL?1PDGF-AA significantly promoted the proliferation of PROCR+cells on the 10th day of expansion,and 100 nmol·L?1crenolanib,a PDGFR α/β inhibitor,significantly inhibited the proliferation of PROCR+cells in the presence or absence of PDGFAA after treatment for 10 days.The FGF signaling network exhibited intensive exchanges among almost all the cell types withFGFligands that were mainly secreted by Fib3(Fig.5e,f).

    By comprehensively predicting signaling networks for intercellular crosstalk,large numbers of ligand-receptor pairs participated in ligand-receptor pairs ofVEGF,TGFB,SEMA3,PDGF,NGF,LAMC,FGF,BMP,andANGPTbetween NPPCs and other cell types(Fig.5f).Interestingly,Fib3 was involved in almost all the above pathways,suggesting its significance in NP homeostasis.We further revealed that NPPC-1,NPPC-2,NPPC-3,and pericytes sent communications to other cells viaIGFandPDGF(Supplementary Fig.6a).As expected,Fib3 was exclusively dominated byFN1in regard to outgoing communication.In addition,Fib3 and Chond1 received incoming communication byBMP,GDF,andANGPT,which reportedly played a prominent essential role in the IVD(Supplementary Fig.6b).26,94–95EGF,used by NPPC-2,NPPC-4,and Neu1 for incoming signaling,could be a protective factor in IVD regeneration(Supplementary Fig.6b).96

    CellChat analysis of NP,AF,and CEP cells also revealed a large number of signaling networks among cell subclusters from the three substructures of the IVD(Supplementary Fig.6c).For example,NPPCs interacted with Chond1 from the AF and CEP.In particular,the GAS signaling pathways were intensively regulated between NPPCs and Chond3 from the AF(Supplementary Fig.6d),possibly protecting the IVD from inflammatory factors.97SPP1,an osteogenesis-related factor,was highly involved in the interaction among NPPCs and stromal cells in the CEP(Supplementary Fig.6d).This was in line with our above hypothesis that osteogenic cells might play a role in IVD homeostasis and/or degenerative processes.However,these proposed signaling pathways should be considered as multiple biological cascades rather than a sole event because the three substructures always work as a whole.

    Taken together,these results indicated that there is a complicated relationship among the distinct cell types and described a cellular crosstalk network with a hierarchical signaling pathway that regulates NP homeostasis in a coordinated manner.

    DISCUSSION

    The severe threat of DDD to human health prompted us to seek an innovative treatment that reestablishes IVD homeostasis.Inadequate knowledge of IVD physiology,and pathology poses a challenge to the development of novel treatment strategies.Due to the cellular heterogeneity and resulting complex microenvironment in the human IVD,an in-depth understanding of specific markers and their roles in IVD homeostasis is urgently needed.Here,we resolved the cellular diversity at a single-cell level using transcriptomic profiling and identified the cell types with a set of specific markers in the human IVD.We classified IVD chondrocytes into three subtypes based on their potential roles in ECM homeostasis.Notably,we identified new subtypes of progenitor cells with signature genes,spatial distribution in situ,and progenitor potential.Moreover,we analyzed the intercellular crosstalk based on the signaling network and uncovered key factors,such as the PDGF and TGF-β cascades,as important cues for regulating the NP microenvironment.Together with previous studies,12,98–99a better understanding of the cellular heterogeneity of the human IVD is developing,with the aim of contributing to new therapeutic strategies for DDD.

    同2.1類似,本協(xié)議中節(jié)點(diǎn)可直接利用私鑰計算與交易相關(guān)的隨機(jī)數(shù)。不同的是:此協(xié)議是交互性的,雙方進(jìn)行密鑰協(xié)商時需選取隨機(jī)數(shù)wA、wB作為與會話相關(guān)的秘密值并做ZKP,對方通過ZKP驗(yàn)證隨機(jī)數(shù)的真實(shí)性后方可計算κ。假設(shè)IoT節(jié)點(diǎn)私鑰的安全性被攻陷,敵手截獲了交易后的某次會話。由于每次會話都重新選取隨機(jī)數(shù),則實(shí)現(xiàn)了后續(xù)會話的不可追蹤,提供了前向安全性。

    The cellular heterogeneity of IVD cells has been a long-debated controversy due to the complexity of the IVD ontogeny,a tricomponent organization with distinct origins.100Multiple developmental origins lead to the inhomogeneity of the cell composition.Although some scholars have attempted to examine the IVD at the single-cell level,a highly precise and unbiased description of cell populations in the human IVD remains to be elucidated.52,98Previously,notochord cells and chondrocytes were recognized in the NP,which was regarded as the notochordal lineage,evidenced by the constant expression ofTBXT.101In line with previous findings,we found a minor cluster that expressed high levels of the markersTBXTandKRT8,which could be a rare but distinct notochord cell cluster.As expected,we found three major clusters of chondrocytes,which are always regarded as core players in ECM homeostasis in the human IVD.Although the expression ofTBXTwas not detected,another notochord marker,NOG,was expressed in the majority of chondrocytes(Supplementary Fig.1f).This interesting finding coincides with a previous theory that distinctive cellular morphology in the NP is due to the various phases along the notochord lineage during aging and degeneration.5,102–105

    Apart from the leading role of notochord lineage cells,the supporting role of minor cell clusters is more notable because of their unclear function,which has been infrequently reported.First,SOX2+NGFR+neurogenic cells,one of the stromal subclusters,were also found in the NP(Supplementary Fig.5a–c).Although the healthy disc was regarded as an aneural tissue,106the pattern of nerve endings has been previously confirmed in healthy and degenerative IVDs,107–110which were small in diameter and relatively sparse.111Thus,sporadic SOX2+neurogenic cells were probably related to neural ingrowth.Furthermore,RUNX2played a part in postnatal IVD development and regulated the notochordal transition into chondrocyte-like cells.112UpregulatedRUNX2expression was also found in the degenerated IVD,which led to IVD calcification.113–114In addition,the stem cells in the IVD exhibited osteogenic potential during ex vivo culture.25These studies may have indicated that the homeostasis of bone formation is important for the physiological and pathological processes of IVD.Our scRNA-seq analysis and immunofluorescence staining revealed the existence of a rare cell cluster that differentially expressed the osteogenic genesRUNX2,DLX5,andSP7,86–88which were defined as osteogenic cells(Supplementary Fig.5a–c).This finding suggested that osteogenic cells exist in healthy IVDs.We hypothesized that osteogenic cells likely contribute to the homeostasis of the IVD or are involved in the pathological process of early degeneration,which began as early as during the teenage years.115–116Finally,the dynamics of vascularization,represented by ECs and pericytes,play a role in disc homeostasis.Previous studies showed that blood vessels penetrated the AF and CEP during the early postnatal years but regressed later,leaving an avascular microenvironment,which accounted for the poor ability for remodeling and repair in IVDs.117–119However,blood vessels are present in the human IVD until even the third decade of life.120During the slow process of vascular regression,it is reasonable that some remnants are left behind,such as ECs.A recent study reported that cross bridges after vascular regression are indeed present in both healthy and degenerated human disks.The cross-bridges of the IVD stained positively for PECAM1 in adult sheep,although the PECAM1+cross-bridges declined with aging.121In line with scRNA-seq analysis,ACTA2+MCAM+pericytes and PECAM1+CD34+ECs were scattered in the IVD(Fig.1b–f and Supplementary Fig.1g).Our data showed that ECs and pericytes communicated with NPPCs via the VEGF,PDGF,and TGF-β signaling pathways,suggesting that they played a role in NP homeostasis(Fig.5).Notably,MCAM is regarded as a classical surface marker of pericytes/MSCs.122Previously,periosteal and meniscal MCAM+cells were shown to exhibit canonical features of skeletogenesis,123–124and MCAM+or ACTA2+cells were also detected in the disc.47,125–127Interestingly,MCAM was specifically expressed in the cell population with migration and repopulating potential in degenerative IVDs.125The functional characteristics of these cell types should be investigated in future studies.The highly conserved cellular heterogeneity across cell clusters between human and rat IVDs(Supplementary Fig.2)suggested that the rat is an ideal animal model to study the role of the above cell clusters in IVD homeostasis.

    Cells in the IVD are generally referred to as“chondrocyte-like”cells or“IVD chondrocytes”.Traditionally,chondrocytes in the IVD are classified into NP,AF,and CEP chondrocytes based on their spatial distribution.However,the spatial-based classification of the cell population was insufficient because of the cellular heterogeneity and possible cell migration among the three sites of the IVD.128Thus,the precise roles of IVD chondrocytes in ECM homeostasis are still largely unknown.16–18,129Therefore,a deeper understanding of the roles of IVD chondrocytes in ECM homeostasis is necessary.Taking advantage of the high throughput nature of analysis at the single-cell level with scRNA-seq,we were able to identify six subclusters of IVD chondrocytes with three functional patterns(Fig.2).First,we identified a new population of regulatory chondrocytes with active GF expression and chondrogenic pathway regulators,implying its regulatory role in chondroid ECM homeostasis.In contrast,homeostatic chondrocytes showed high similarity to classical chondrocytes,which were quiescent,fully differentiated,and responsible for ECM deposition.130Interestingly,homeostatic chondrocytes were enriched in circadian regulation genes,which involved key pathways regulating the homeostasis of IVDs.131This finding suggests that homeostatic chondrocytes could be a potential therapeutic target for circadian rhythm in the human IVD.It is noteworthy that the effector chondrocytes were metabolically active,which is important in maintaining the ECM biogenesis of the IVD.132In addition,the high expression ofPRG4(lubricin)also implies that they play a protective role in reducing shear stress and inflammation and keeping the joint healthy.133In contrast,effector chondrocytes were characterized by ossification and shared expression patterns with articular HTCs.51Thus,the definitive function of effector chondrocytes is certainly worth future investigation.Overall,the six transcriptomically defined populations of chondrocytes exhibited distinct roles in ECM homeostasis,providing new perspectives for exploring the mechanism of IVD chondrocytes.

    The IVD possesses the capability of spontaneous regeneration,as evidenced by self-healing after disc degeneration,134probably due to the presence of in situ progenitor cells.Progenitor cells expressing different marker gene sets existed in three compartments of the IVD.59–60The progenitor cells exhibited certain plasticity and the ability to slow down disc degeneration.135–136Thus,it is a promising strategy to activate endogenous progenitor cells or transplant exogenous progenitor cells for DDD therapy.However,a comprehensive understanding of their in vivo characteristics,including discriminable identity,lineage,spatial distribution,and functional role,is still lacking.We sought to help to increase the understanding of progenitor cells at a single-cell resolution.Surprisingly,we found a cluster of cells that exclusively expressedPDGFRA,a signature of MSCs,77,137–138and was mainly distributed in the NP(Supplementary Fig.4).Notably,the PDGFRA+PROCR+NPPC subcluster was enriched for genes in the SMAD3 signaling pathway and exhibited higher activation of p-SMAD3(Fig.3),which determines the TGF-β-induced chondrogenesis139and cell fate decisions of stem cells by participating in the cell-cycle process and binding of m6A methyltransferase.140–141Moreover,PROCRwas used to sort rare progenitor/stem cells with high efficacy.For example,PROCR(encoding CD201)was used as a sorting marker to harvest isolated 1%of islet cells,which robustly formed islet-like organoids.75Applications in the hematopoietic system showed thatPROCRenriched T1 prehematopoietic stem cells at a resolution of 68 parts per million and functional HSCs in the human fetal liver.71,73In this study,we identified an NPPC cluster that highly expressedPROCR,which exhibited pluripotency with colony-formation capacity and osteochondrogenic potentials(Figs.3 and 4),similar to the characteristics of multipotent mesenchymal stromal cells.142Thus,we characterized these cells as resident progenitor cells in the human IVD.It is possible that the alternative cell fate in NPPCs determines the outcome of the IVD when a degenerative program is initiated.On the one hand,the chondrogenic fate could help rebalance IVD homeostasis via cell replenishment.143On the other hand,the osteogenic fate could lead to DDD by inducing heterotopic ossification.144Accordingly,these results have new implications for innovative therapeutic strategies targeting NPPCs.

    The two branches of the cell fate of NPPCs motivated us to explore the key regulatory factors.Resident progenitor cells are exhausted or altered during degeneration,26,145indicating that the microenvironment has a significant influence on cell fate.To identify the key factors regulating the fate of NPPCs,CellChat analysis was used to dissect the intercellular crosstalk based on the signaling network in the human IVD(Fig.5).We found that GFrelated signaling pathways were involved in the crosstalk network,mainly including the previously reported FGF family,143,146TGF-β family,147–148BMP family,149–150and PDGF family.151Among them,TGF-β was important due to the high activation of SMAD3 in NPPCs.TGF-β directs embryonic matrix development within the notochord and promotes the differentiation of the sclerotome into the AF,61,152–153suggesting that it is an inherent regulator of the human IVD.Previous studies have shown that the TGF-β family plays an important role in the development and protection of the IVD,especially in maintaining the phenotype of chondrocytes.154Moreover,the loss of TGF-β signaling in growth plate chondrocytes and inner AF cells led to the loss of matrix tissue and endplate cartilage cells and abnormal growth plate cartilage morphology in Tgfbr2 conditional knockout mice.155The critical role of TGF-β was also evidenced by the observation that the knockout of SMAD3,the key downstream target of TGF-β,led to the spontaneous development of IVD degeneration in 30-day-old mice.156In addition,TGF-β has been shown to have a beneficial effect on chondrogenic anabolism in MSCs.157In this study,TGFBwas involved in regulating NPPCs,as evidenced by TGF-β3 promoting the chondrogenesis of PROCR+cells(Fig.5).Meanwhile,the secretory role of chondrocyte clusters on TGF-β should not be neglected in the human IVD(Fig.5).Furthermore,PDGFwas found to engage in regulating NPPCs,probably due to the exclusive expression of its receptor genePDGFRAin NPPCs.Previous studies have reported that PDGF-AA is involved in the regulation of cell proliferation.92–93In line with the CellChat analysis,we found that PDGF-AA significantly promoted the proliferation of PROCR+cells(Fig.5).Interestingly,all the minor clusters in the NP are involved in interacting with NPPCs,suggesting their potential role in regulating NPPCs and subsequently maintaining IVD homeostasis.Moreover,further investigations need to elucidate their roles and establish an innovative strategy to optimize the microenvironment and benefit IVD stem/progenitor cells.

    Although we validated the existence of identified cell populations by flow cytometry,immunofluorescence staining,and scRNA-seq evidence from the rat IVD,we surprisingly found that Sample 1 was from a 16-year-old donor who suffered from vertebral fracture exhibited obvious variability in the proportion of cell clusters(Supplementary Fig.1e).Acute trauma has been shown to stimulate resident cells to regenerate in previous studies.158–160Interestingly,a recent study reported that NP cells derived from trauma patients showed higher adipogenic and chondrogenic potential than those derived from degenerated IVDs.161Thus,we are more inclined to hypothesize that ECs,pericytes,and NPPCs are rare in the IVD,and acute trauma may induce local regeneration,which accounts for the unwanted distribution variability across donors.Due to the scarcity of desirable samples of healthy disks from young patients with vertebral fractures,this needs to be explored in future studies.

    In summary,our study described the cell atlas of the human IVD,providing a valuable resource for further investigation of IVD homeostasis at the mechanistic level.The cellular heterogeneity and signaling network we uncovered help to increase the understanding of the human IVD at a single-cell level and provide crucial clues for establishing new therapeutic strategies for DDD treatment.

    MATERIALS AND METHODS

    For full methods,see the Supplementary Methods.

    Human IVD tissue specimens

    This study was approved by the Institutional Ethics Review Board of Daping Hospital[Ethics Committee(2019-127)]and the Chinese Clinical Trial registry(ChiCTR1900028201).All procedures were performed in accordance with the ethical standards of the committee responsible for human experimentation and with the Declaration of Helsinki of 1975,as revised in 2000.Informed consent was obtained from all patients for inclusion in the study.Eleven human IVDs were carefully dissected from nine donors in this study(Supplementary Table 1).The gelatinous tissue from the central region was harvested as the NP.The peripheral lamellar structure of the outer IVD was harvested as AF.The superior and inferior homogeneous cartilage tissue was harvested as CEP.The sampling areas of NP,AF,and CEP are indicated(Supplementary Fig.1a).

    Single-cell RNA sequencing

    The cells were washed with PBS three times and concentrated to 700–1 200 cells per μL.The suspension was then loaded on a Chromium Controller(10X Genomics).For scRNA-seq library construction,a Chromium Single Cell 3′Library and Gel Bead Kit V2(10X Genomics,PN120237)was utilized to generate single-cell gel beads in emulsion(GEM)within barcoded,full-length cDNA from polyadenylated mRNA.The captured cells were lysed in GEM,and the released RNA was reverse-transcribed with primers containing poly-T,a barcode,UMIs,and the read 1 primer sequence,in that order.Barcoded,full-length cDNA was PCR amplified for library construction.After enzymatic fragmentation,an adapter ligation reaction was performed to add a sample index and read 2 primer sequences to the cDNA fragment.After quality control,the libraries were sequenced on an Illumina NovaSeq 6000 platform to generate 150-bp paired-end reads,according to the manufacturer’s instructions(Berry Genomics).

    DATA AVAILABILITY

    All data from the study are available in online supplementary files.The scRNA-seq data have been deposited in GEO(GSE160756).All other relevant data from this study are available from the corresponding authors upon reasonable request.

    ACKNOWLEDGEMENTS

    The authors thank Hongxia Hu and Haitao Hu from Berry Genomics for technical support with scRNA-seq.The authors acknowledge Zhigang Zhou and Feng Wei from Peking University Third Hospital and Fei Luo,Bo Yu,Bo Huang,Qizhao Huang,Qinghua Ma,and Ruili Cai from Army Medical University for their excellent technical support with sample collection and single-cell preparation.The authors are grateful to Wenxia Zheng and Hengsheng Tao from Olympus and Xue Yang and Qing Zhou from Army Medical University for help with section staining and imaging.This study was supported by grants from the National Natural Science Foundation of China(81802165 and 31930054),National Key Research and Development Program of China(2017YFA0103401 and 2019YFA0110201),Training Plan of Talents’Innovation of Army Medical Center of PLA(2019CXJSB013),Postdoctoral Innovative Talent Support Program in Chongqing(2019-298),and Fund for Excellent Young Scholars of the State Key Laboratory of Trauma,Burns and Combined Injury(SKLYQ201902).

    AUTHOR CONTRIBUTIONS

    P.L.and B.L.designed and supervised the study;Y.G.,J.Z.,and O.H.performed the sample collection and preparation with help from P.L.,L.C.,and Y.X.;Y.G.performed single-cell RNA sequencing with the help of Z.X.,Z.W.,and J.H.;Y.G.performed the histological,immunohistochemistry,and immunofluorescence staining with the help of M.J.and S.H.;Y.G.,J.H.,J.Y.,Z.B.,and O.H.performed the FAC sorting,CFU-F,and trilineage differentiation experiments with the help of Z.X.;J.H.performed the bioinformatic analysis with help from B.L.,P.L.,and Y.G.;P.L.,B.L.,Y.G.,and J.H.wrote the manuscript.All authors read and approved the manuscript.

    ADDITIONAL INFORMATION

    Supplementary informationThe online version contains supplementary material available at https://doi.org/10.1038/s41413-021-00163-z.

    Competing interests:The authors declare no competing interests.

    猜你喜歡
    時需敵手會話
    辨別數(shù)列的特征,選用合適的方法求數(shù)列的和
    湖北
    生活哲理
    不帶著怒氣做任何事
    有意冒犯性言語的會話含義分析
    漢語教材中的會話結(jié)構(gòu)特征及其語用功能呈現(xiàn)——基于85個會話片段的個案研究
    我畫我心
    沖突語的會話分析研究
    對外漢語課堂英語通用語的會話調(diào)整功能
    不帶著怒氣作戰(zhàn)
    老司机福利观看| 看片在线看免费视频| 成人18禁在线播放| 人妻久久中文字幕网| 成人亚洲精品av一区二区| 黄色丝袜av网址大全| 人人妻人人看人人澡| 国产精品99久久久久久久久| 成人国产综合亚洲| 亚洲专区国产一区二区| 国产高清视频在线播放一区| 99国产极品粉嫩在线观看| 黄色女人牲交| 欧美性猛交╳xxx乱大交人| a级毛片在线看网站| 给我免费播放毛片高清在线观看| x7x7x7水蜜桃| 国产av一区在线观看免费| 亚洲 欧美一区二区三区| 亚洲中文日韩欧美视频| 久久热在线av| 亚洲一区二区三区不卡视频| 不卡av一区二区三区| 999精品在线视频| 免费看a级黄色片| 亚洲国产精品sss在线观看| 中文字幕人成人乱码亚洲影| 丁香欧美五月| 精品国产三级普通话版| 黄频高清免费视频| 亚洲欧美精品综合久久99| 国产午夜精品论理片| 波多野结衣高清无吗| 床上黄色一级片| 香蕉国产在线看| 国产又黄又爽又无遮挡在线| 欧美日韩黄片免| 91麻豆精品激情在线观看国产| 欧美zozozo另类| 久久久久国产一级毛片高清牌| 成人永久免费在线观看视频| 色噜噜av男人的天堂激情| 日韩欧美免费精品| 亚洲avbb在线观看| 日本熟妇午夜| 国产精品久久久人人做人人爽| 搡老妇女老女人老熟妇| 蜜桃久久精品国产亚洲av| 亚洲va日本ⅴa欧美va伊人久久| 久久久久亚洲av毛片大全| 久久久久久国产a免费观看| 亚洲国产色片| 黑人操中国人逼视频| 久久欧美精品欧美久久欧美| 少妇裸体淫交视频免费看高清| 老司机午夜十八禁免费视频| 一个人看视频在线观看www免费 | 天天躁狠狠躁夜夜躁狠狠躁| 精品免费久久久久久久清纯| 亚洲精品乱码久久久v下载方式 | 99精品欧美一区二区三区四区| 99riav亚洲国产免费| 亚洲第一电影网av| av在线蜜桃| www.999成人在线观看| 久久久久九九精品影院| 国产欧美日韩精品一区二区| 精品久久久久久,| 一个人免费在线观看电影 | 免费电影在线观看免费观看| 一级作爱视频免费观看| 9191精品国产免费久久| 18禁黄网站禁片免费观看直播| 观看美女的网站| 首页视频小说图片口味搜索| 网址你懂的国产日韩在线| 国产一区二区在线观看日韩 | 久久精品综合一区二区三区| 成年女人毛片免费观看观看9| 成人永久免费在线观看视频| 丰满的人妻完整版| 亚洲国产色片| 搞女人的毛片| 国产精品久久久人人做人人爽| 日韩高清综合在线| 可以在线观看毛片的网站| 日本免费一区二区三区高清不卡| 国产精品精品国产色婷婷| 欧美黄色淫秽网站| 网址你懂的国产日韩在线| 国产蜜桃级精品一区二区三区| 国产精品 国内视频| svipshipincom国产片| av福利片在线观看| 亚洲精品美女久久av网站| 在线观看日韩欧美| 波多野结衣巨乳人妻| 此物有八面人人有两片| 国产成人精品无人区| 欧洲精品卡2卡3卡4卡5卡区| 99久久精品一区二区三区| 天天一区二区日本电影三级| 宅男免费午夜| 欧美在线一区亚洲| 欧美乱妇无乱码| 午夜精品久久久久久毛片777| 99久久国产精品久久久| 欧美另类亚洲清纯唯美| 脱女人内裤的视频| 欧美精品啪啪一区二区三区| av在线蜜桃| 丁香欧美五月| 99精品欧美一区二区三区四区| 亚洲av第一区精品v没综合| 嫁个100分男人电影在线观看| 欧美黄色淫秽网站| 日日夜夜操网爽| 在线a可以看的网站| 免费看a级黄色片| 麻豆国产97在线/欧美| 国产日本99.免费观看| 欧美乱色亚洲激情| 国产伦精品一区二区三区视频9 | 成人亚洲精品av一区二区| 精品不卡国产一区二区三区| 女生性感内裤真人,穿戴方法视频| xxxwww97欧美| 中文亚洲av片在线观看爽| 女生性感内裤真人,穿戴方法视频| xxxwww97欧美| 亚洲中文字幕日韩| 偷拍熟女少妇极品色| 狂野欧美激情性xxxx| 国产精品1区2区在线观看.| 五月玫瑰六月丁香| 亚洲第一欧美日韩一区二区三区| 亚洲专区国产一区二区| 国产一区二区三区在线臀色熟女| 亚洲天堂国产精品一区在线| 女人高潮潮喷娇喘18禁视频| 亚洲五月婷婷丁香| 欧美绝顶高潮抽搐喷水| 欧美黄色片欧美黄色片| 欧美3d第一页| 欧美三级亚洲精品| 黄频高清免费视频| 一级毛片女人18水好多| 五月伊人婷婷丁香| 男人的好看免费观看在线视频| 免费一级毛片在线播放高清视频| 国产精品久久久人人做人人爽| 久久久国产成人免费| 亚洲无线观看免费| 九九久久精品国产亚洲av麻豆 | 成人特级黄色片久久久久久久| 久久久水蜜桃国产精品网| 亚洲成a人片在线一区二区| 色老头精品视频在线观看| 黄片小视频在线播放| 国产成人aa在线观看| 国产免费男女视频| www.999成人在线观看| e午夜精品久久久久久久| 曰老女人黄片| a级毛片a级免费在线| 欧美日韩亚洲国产一区二区在线观看| 欧美+亚洲+日韩+国产| 国产黄色小视频在线观看| 天堂影院成人在线观看| 久久国产乱子伦精品免费另类| 国产精品 国内视频| 免费在线观看日本一区| 18禁黄网站禁片免费观看直播| 日本一二三区视频观看| 亚洲最大成人中文| 国产三级黄色录像| 99久久成人亚洲精品观看| 国产乱人视频| 十八禁人妻一区二区| 99国产精品99久久久久| 免费电影在线观看免费观看| 欧美成人免费av一区二区三区| 床上黄色一级片| 首页视频小说图片口味搜索| ponron亚洲| 最近最新中文字幕大全免费视频| 国产激情欧美一区二区| 在线观看一区二区三区| 亚洲人成伊人成综合网2020| 变态另类成人亚洲欧美熟女| 国内毛片毛片毛片毛片毛片| 国内精品久久久久精免费| 日韩成人在线观看一区二区三区| 一级作爱视频免费观看| 在线观看免费视频日本深夜| 两个人视频免费观看高清| 天堂网av新在线| 日韩成人在线观看一区二区三区| 国产成人aa在线观看| 国语自产精品视频在线第100页| 亚洲熟妇熟女久久| 特大巨黑吊av在线直播| 欧美+亚洲+日韩+国产| 国产精品 国内视频| 成人精品一区二区免费| 国产高清视频在线观看网站| 午夜成年电影在线免费观看| 亚洲专区中文字幕在线| 禁无遮挡网站| 婷婷精品国产亚洲av| 国产亚洲欧美在线一区二区| 国产高清三级在线| 最好的美女福利视频网| 99精品久久久久人妻精品| 99国产精品一区二区三区| 97超视频在线观看视频| 成人三级做爰电影| 婷婷精品国产亚洲av在线| a级毛片在线看网站| 亚洲午夜理论影院| 欧美黑人欧美精品刺激| 丰满人妻熟妇乱又伦精品不卡| 嫩草影视91久久| 男女视频在线观看网站免费| 小蜜桃在线观看免费完整版高清| 88av欧美| 亚洲欧美精品综合久久99| 免费在线观看影片大全网站| 日韩高清综合在线| 一个人看的www免费观看视频| 亚洲av日韩精品久久久久久密| 色尼玛亚洲综合影院| 国产精品久久久久久亚洲av鲁大| 亚洲精品久久国产高清桃花| 老熟妇乱子伦视频在线观看| 搡老妇女老女人老熟妇| 久久精品国产综合久久久| 18禁国产床啪视频网站| 成年女人毛片免费观看观看9| 国产一区二区三区视频了| 性色av乱码一区二区三区2| 中文字幕久久专区| 国产91精品成人一区二区三区| 美女高潮的动态| 我的老师免费观看完整版| 亚洲国产精品成人综合色| 免费一级毛片在线播放高清视频| 在线免费观看的www视频| 欧美丝袜亚洲另类 | 成年免费大片在线观看| 悠悠久久av| 成人午夜高清在线视频| 一本综合久久免费| 日韩欧美一区二区三区在线观看| 后天国语完整版免费观看| 又黄又爽又免费观看的视频| 丰满人妻一区二区三区视频av | 很黄的视频免费| 亚洲电影在线观看av| 欧美在线一区亚洲| www日本在线高清视频| 亚洲中文av在线| 亚洲天堂国产精品一区在线| 国产伦精品一区二区三区视频9 | 窝窝影院91人妻| 亚洲激情在线av| www国产在线视频色| 日本免费一区二区三区高清不卡| 最近在线观看免费完整版| 亚洲专区国产一区二区| 欧美黑人巨大hd| 国产精品99久久久久久久久| 色综合站精品国产| 天堂网av新在线| 国产蜜桃级精品一区二区三区| 性欧美人与动物交配| 日本免费一区二区三区高清不卡| 亚洲国产色片| 午夜福利在线观看免费完整高清在 | 午夜两性在线视频| 国产精品久久视频播放| 亚洲成av人片在线播放无| 精华霜和精华液先用哪个| 色综合亚洲欧美另类图片| 精品久久久久久,| 制服丝袜大香蕉在线| 国产亚洲欧美98| 中文资源天堂在线| 无遮挡黄片免费观看| 亚洲av熟女| 国产aⅴ精品一区二区三区波| 久久精品91蜜桃| 亚洲,欧美精品.| 午夜福利视频1000在线观看| 天堂影院成人在线观看| 一级黄色大片毛片| 丰满的人妻完整版| 给我免费播放毛片高清在线观看| 久久热在线av| 成人av一区二区三区在线看| 成人特级av手机在线观看| 午夜激情欧美在线| 中文字幕av在线有码专区| 日韩精品青青久久久久久| 综合色av麻豆| 女人被狂操c到高潮| 久久久久免费精品人妻一区二区| 久久99热这里只有精品18| 国产美女午夜福利| 亚洲最大成人中文| 99精品欧美一区二区三区四区| xxxwww97欧美| 国产伦在线观看视频一区| 丁香六月欧美| 午夜a级毛片| 欧美黑人巨大hd| 少妇的逼水好多| 男人舔奶头视频| 成人特级av手机在线观看| 国产精品永久免费网站| 日韩欧美国产一区二区入口| 中亚洲国语对白在线视频| 免费观看精品视频网站| 美女大奶头视频| 麻豆国产av国片精品| 99国产精品一区二区三区| 88av欧美| 欧美黄色淫秽网站| 九九在线视频观看精品| 露出奶头的视频| 级片在线观看| 亚洲精品一区av在线观看| 18禁裸乳无遮挡免费网站照片| 无限看片的www在线观看| 亚洲国产欧洲综合997久久,| 欧美极品一区二区三区四区| 亚洲一区二区三区色噜噜| 黄色视频,在线免费观看| 级片在线观看| 国产男靠女视频免费网站| 久久久久国内视频| 亚洲国产日韩欧美精品在线观看 | 精品国产美女av久久久久小说| 在线观看日韩欧美| 757午夜福利合集在线观看| 免费高清视频大片| 中出人妻视频一区二区| 身体一侧抽搐| 两个人视频免费观看高清| 国产又黄又爽又无遮挡在线| 国产野战对白在线观看| 国产成人影院久久av| 欧美高清成人免费视频www| 久久这里只有精品中国| 很黄的视频免费| 999精品在线视频| 欧美在线一区亚洲| 黄片大片在线免费观看| 午夜激情福利司机影院| 麻豆成人av在线观看| 精品久久久久久久毛片微露脸| 黄色丝袜av网址大全| 一区二区三区高清视频在线| 日韩成人在线观看一区二区三区| 美女大奶头视频| cao死你这个sao货| 啦啦啦韩国在线观看视频| 无遮挡黄片免费观看| 十八禁人妻一区二区| 亚洲av成人一区二区三| 久久精品国产亚洲av香蕉五月| 国产麻豆成人av免费视频| 国产伦一二天堂av在线观看| 人人妻人人澡欧美一区二区| 亚洲国产欧美一区二区综合| 小蜜桃在线观看免费完整版高清| 神马国产精品三级电影在线观看| 日日干狠狠操夜夜爽| 黑人欧美特级aaaaaa片| 黑人巨大精品欧美一区二区mp4| 国产欧美日韩精品亚洲av| 久久草成人影院| 动漫黄色视频在线观看| 欧美性猛交╳xxx乱大交人| 熟女少妇亚洲综合色aaa.| 精品国产超薄肉色丝袜足j| 亚洲无线在线观看| 校园春色视频在线观看| 午夜视频精品福利| 一边摸一边抽搐一进一小说| 天天一区二区日本电影三级| 女生性感内裤真人,穿戴方法视频| 欧美国产日韩亚洲一区| 热99re8久久精品国产| 99久久99久久久精品蜜桃| av国产免费在线观看| 可以在线观看的亚洲视频| 精品不卡国产一区二区三区| 亚洲av第一区精品v没综合| 国产1区2区3区精品| 丰满的人妻完整版| 中出人妻视频一区二区| 999久久久国产精品视频| www.www免费av| 无限看片的www在线观看| 长腿黑丝高跟| 亚洲国产中文字幕在线视频| 日韩人妻高清精品专区| 午夜福利在线观看吧| 亚洲国产欧美网| 国产精品一及| 91麻豆av在线| 亚洲av成人不卡在线观看播放网| 久久久久九九精品影院| 一本综合久久免费| 99在线人妻在线中文字幕| 天堂影院成人在线观看| 中文字幕人成人乱码亚洲影| 精品一区二区三区视频在线 | 国产成人av激情在线播放| 日韩 欧美 亚洲 中文字幕| 国产成人aa在线观看| 夜夜看夜夜爽夜夜摸| 亚洲黑人精品在线| 搡老熟女国产l中国老女人| 男女床上黄色一级片免费看| 美女扒开内裤让男人捅视频| 日本五十路高清| 99精品久久久久人妻精品| 在线a可以看的网站| 日日干狠狠操夜夜爽| 午夜日韩欧美国产| 五月玫瑰六月丁香| 欧美成人免费av一区二区三区| 国产午夜福利久久久久久| 色精品久久人妻99蜜桃| 欧美乱码精品一区二区三区| 一二三四社区在线视频社区8| 成人18禁在线播放| 日韩成人在线观看一区二区三区| 久99久视频精品免费| 日本黄色视频三级网站网址| 69av精品久久久久久| 香蕉久久夜色| 国产一区二区在线av高清观看| 一级作爱视频免费观看| 欧美极品一区二区三区四区| 2021天堂中文幕一二区在线观| 最近在线观看免费完整版| 麻豆成人av在线观看| 日韩av在线大香蕉| 成人永久免费在线观看视频| 亚洲av电影在线进入| 可以在线观看毛片的网站| 免费电影在线观看免费观看| 99久国产av精品| 久久久国产欧美日韩av| 日本撒尿小便嘘嘘汇集6| 久久久久久人人人人人| 香蕉国产在线看| 动漫黄色视频在线观看| 久久精品91无色码中文字幕| 免费在线观看视频国产中文字幕亚洲| 久久久久久久午夜电影| 亚洲欧美日韩高清在线视频| 国产精品久久久久久精品电影| 成熟少妇高潮喷水视频| 精品久久久久久,| 在线看三级毛片| 午夜视频精品福利| 毛片女人毛片| 亚洲 欧美 日韩 在线 免费| 亚洲精华国产精华精| 白带黄色成豆腐渣| 国产高清有码在线观看视频| 亚洲av成人一区二区三| 久久精品影院6| 真实男女啪啪啪动态图| 香蕉av资源在线| 伦理电影免费视频| 男女做爰动态图高潮gif福利片| 欧美日韩亚洲国产一区二区在线观看| 久久久国产欧美日韩av| 国产69精品久久久久777片 | 国模一区二区三区四区视频 | a级毛片在线看网站| 亚洲av五月六月丁香网| 狠狠狠狠99中文字幕| 亚洲欧美日韩东京热| 我要搜黄色片| 亚洲自偷自拍图片 自拍| 99国产综合亚洲精品| netflix在线观看网站| 视频区欧美日本亚洲| 淫妇啪啪啪对白视频| 日韩欧美三级三区| 两个人看的免费小视频| 午夜免费观看网址| 级片在线观看| 他把我摸到了高潮在线观看| 级片在线观看| 99热只有精品国产| 日本一本二区三区精品| 97超视频在线观看视频| 九色国产91popny在线| 久久中文看片网| 国产aⅴ精品一区二区三区波| 精品乱码久久久久久99久播| 久久热在线av| 色综合亚洲欧美另类图片| 国产av一区在线观看免费| 他把我摸到了高潮在线观看| 伦理电影免费视频| 国产精品免费一区二区三区在线| 特大巨黑吊av在线直播| 国内揄拍国产精品人妻在线| 熟女少妇亚洲综合色aaa.| 男女下面进入的视频免费午夜| 精品久久久久久成人av| 在线观看一区二区三区| 日本在线视频免费播放| 九九久久精品国产亚洲av麻豆 | 亚洲国产欧美网| 国产av不卡久久| 国产激情偷乱视频一区二区| 久久久国产精品麻豆| 欧美高清成人免费视频www| 偷拍熟女少妇极品色| 99国产精品一区二区三区| 亚洲av熟女| 欧美不卡视频在线免费观看| 最近最新免费中文字幕在线| 在线观看66精品国产| 高清毛片免费观看视频网站| 国产视频一区二区在线看| 青草久久国产| 国产高清视频在线播放一区| 热99re8久久精品国产| 久久人人精品亚洲av| 成人三级做爰电影| 男女做爰动态图高潮gif福利片| 久久久久精品国产欧美久久久| 国产精品,欧美在线| 欧美一级a爱片免费观看看| 宅男免费午夜| 国内精品久久久久精免费| 国产精品日韩av在线免费观看| 亚洲成av人片在线播放无| 日本免费a在线| avwww免费| 在线播放国产精品三级| 亚洲激情在线av| 精品一区二区三区av网在线观看| 免费高清视频大片| 国产伦精品一区二区三区四那| 90打野战视频偷拍视频| 中文字幕高清在线视频| 日韩欧美免费精品| 免费看光身美女| 一个人观看的视频www高清免费观看 | 成人欧美大片| 中文字幕熟女人妻在线| 两性午夜刺激爽爽歪歪视频在线观看| 成人午夜高清在线视频| av福利片在线观看| 又大又爽又粗| 身体一侧抽搐| 午夜激情欧美在线| 日韩成人在线观看一区二区三区| 欧美另类亚洲清纯唯美| 黄色女人牲交| 精品国内亚洲2022精品成人| 欧美黑人欧美精品刺激| 黄色丝袜av网址大全| 亚洲片人在线观看| 久久人妻av系列| 午夜免费观看网址| 精品99又大又爽又粗少妇毛片 | 一二三四社区在线视频社区8| 中文字幕av在线有码专区| 少妇裸体淫交视频免费看高清| 欧美不卡视频在线免费观看| 一本久久中文字幕| 亚洲 国产 在线| 亚洲第一欧美日韩一区二区三区| 国产伦精品一区二区三区视频9 | 女人被狂操c到高潮| 最新中文字幕久久久久 | 不卡一级毛片| 老司机午夜十八禁免费视频| 啦啦啦免费观看视频1| 色尼玛亚洲综合影院| 99热只有精品国产| 波多野结衣巨乳人妻| 日韩欧美免费精品| 亚洲欧美精品综合一区二区三区| 真人一进一出gif抽搐免费| 亚洲精品在线观看二区| 99久久综合精品五月天人人| 欧美午夜高清在线| 可以在线观看毛片的网站| 香蕉丝袜av| 欧美日韩综合久久久久久 | 亚洲欧美一区二区三区黑人| 婷婷亚洲欧美| 国产亚洲av嫩草精品影院| 日韩精品青青久久久久久| 成人三级做爰电影| 国产成年人精品一区二区| 欧美高清成人免费视频www| 岛国在线免费视频观看| 手机成人av网站| 两性夫妻黄色片| 国产男靠女视频免费网站| 色尼玛亚洲综合影院|