• <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)
    av网站在线播放免费| 久久久久精品国产欧美久久久| 18禁美女被吸乳视频| 国产主播在线观看一区二区| 国产精品亚洲一级av第二区| 美女福利国产在线| 国产成人啪精品午夜网站| 欧美精品一区二区免费开放| 久久人妻熟女aⅴ| 黄片小视频在线播放| 精品亚洲成国产av| 少妇 在线观看| 99久久99久久久精品蜜桃| 欧美精品av麻豆av| 国产国语露脸激情在线看| 人成视频在线观看免费观看| 啦啦啦 在线观看视频| 国产一区有黄有色的免费视频| 麻豆乱淫一区二区| 女人爽到高潮嗷嗷叫在线视频| 高清在线国产一区| 很黄的视频免费| 又大又爽又粗| 黄色视频,在线免费观看| 久久天堂一区二区三区四区| 午夜免费成人在线视频| 成人特级黄色片久久久久久久| 国产成+人综合+亚洲专区| xxx96com| 自拍欧美九色日韩亚洲蝌蚪91| 久久久国产成人免费| 天天操日日干夜夜撸| 精品久久久久久久毛片微露脸| 欧美黑人精品巨大| 丝瓜视频免费看黄片| 亚洲精品成人av观看孕妇| 无限看片的www在线观看| 久久人人爽av亚洲精品天堂| 99久久综合精品五月天人人| 一级作爱视频免费观看| 日韩人妻精品一区2区三区| 国产精品一区二区在线观看99| 他把我摸到了高潮在线观看| 免费黄频网站在线观看国产| 久久久久久免费高清国产稀缺| 久久精品亚洲熟妇少妇任你| 日日爽夜夜爽网站| 国产精品永久免费网站| 两人在一起打扑克的视频| 欧美色视频一区免费| 母亲3免费完整高清在线观看| 人人妻人人澡人人爽人人夜夜| 亚洲精品一二三| 欧美在线黄色| 国产精品影院久久| 国产亚洲一区二区精品| www日本在线高清视频| 啦啦啦视频在线资源免费观看| 国产精品久久久久久精品古装| 免费日韩欧美在线观看| 中文字幕av电影在线播放| 色尼玛亚洲综合影院| 精品欧美一区二区三区在线| 久9热在线精品视频| 国产精品久久视频播放| 久久人妻熟女aⅴ| 亚洲人成77777在线视频| 国产一区二区三区综合在线观看| 欧美黑人精品巨大| 国产成人系列免费观看| 嫩草影视91久久| 美女视频免费永久观看网站| 在线观看免费视频日本深夜| av不卡在线播放| 丝袜美足系列| 亚洲精品av麻豆狂野| 久久热在线av| 多毛熟女@视频| 99久久综合精品五月天人人| 日本vs欧美在线观看视频| 国产极品粉嫩免费观看在线| 男人的好看免费观看在线视频 | 欧美 亚洲 国产 日韩一| 真人做人爱边吃奶动态| 久久人妻av系列| 男女高潮啪啪啪动态图| 成人三级做爰电影| 亚洲免费av在线视频| 人人妻人人澡人人看| 欧美日韩一级在线毛片| 两性夫妻黄色片| 黄片大片在线免费观看| 亚洲精品自拍成人| 交换朋友夫妻互换小说| 久久精品国产亚洲av高清一级| 香蕉丝袜av| 午夜福利在线免费观看网站| 日本黄色日本黄色录像| 99热只有精品国产| 亚洲精品一二三| 日韩欧美一区二区三区在线观看 | 欧美日本中文国产一区发布| 欧美黑人欧美精品刺激| 久久精品国产清高在天天线| av在线播放免费不卡| 搡老岳熟女国产| 日韩成人在线观看一区二区三区| 十八禁高潮呻吟视频| 丝瓜视频免费看黄片| 久久久久久久国产电影| 婷婷成人精品国产| 老汉色∧v一级毛片| 成人特级黄色片久久久久久久| 激情在线观看视频在线高清 | 人人妻人人澡人人爽人人夜夜| 国产精品.久久久| 啦啦啦 在线观看视频| 岛国毛片在线播放| 精品人妻熟女毛片av久久网站| 精品乱码久久久久久99久播| 午夜福利欧美成人| 精品视频人人做人人爽| 国产精品香港三级国产av潘金莲| 美女午夜性视频免费| 免费少妇av软件| 在线观看一区二区三区激情| 欧美精品高潮呻吟av久久| 亚洲美女黄片视频| videos熟女内射| 人人妻,人人澡人人爽秒播| 91字幕亚洲| 亚洲全国av大片| 午夜精品国产一区二区电影| 国产亚洲一区二区精品| 欧美亚洲日本最大视频资源| 国产亚洲一区二区精品| 国产一区二区三区视频了| 日本黄色视频三级网站网址 | 免费看十八禁软件| 视频在线观看一区二区三区| 国产精品久久视频播放| 视频在线观看一区二区三区| a在线观看视频网站| 国产成人影院久久av| 嫁个100分男人电影在线观看| 十八禁人妻一区二区| 两性夫妻黄色片| 亚洲精华国产精华精| 这个男人来自地球电影免费观看| 欧美亚洲 丝袜 人妻 在线| 成人手机av| 在线观看免费午夜福利视频| 久久人人97超碰香蕉20202| 男人舔女人的私密视频| 黑人巨大精品欧美一区二区蜜桃| 国产男靠女视频免费网站| 极品少妇高潮喷水抽搐| 一级a爱片免费观看的视频| 久久亚洲精品不卡| 午夜福利乱码中文字幕| 宅男免费午夜| 99在线人妻在线中文字幕 | 精品少妇久久久久久888优播| 日日夜夜操网爽| 99久久99久久久精品蜜桃| 久久精品亚洲精品国产色婷小说| 亚洲午夜精品一区,二区,三区| 黑人巨大精品欧美一区二区mp4| 热re99久久精品国产66热6| 女人被躁到高潮嗷嗷叫费观| 亚洲精品在线观看二区| 99re在线观看精品视频| 午夜免费观看网址| 在线天堂中文资源库| videosex国产| 国产伦人伦偷精品视频| a级毛片在线看网站| 国产av又大| 国产成人一区二区三区免费视频网站| 精品一区二区三区四区五区乱码| 久久久久视频综合| 亚洲欧洲精品一区二区精品久久久| 脱女人内裤的视频| 欧美乱码精品一区二区三区| 亚洲色图av天堂| 美女视频免费永久观看网站| av天堂在线播放| 午夜福利乱码中文字幕| 亚洲久久久国产精品| 国产精品自产拍在线观看55亚洲 | 一本大道久久a久久精品| 黄色成人免费大全| 午夜视频精品福利| 乱人伦中国视频| 美女午夜性视频免费| 欧美黑人欧美精品刺激| 老司机影院毛片| cao死你这个sao货| 久久精品亚洲av国产电影网| 人人妻人人澡人人看| 高清av免费在线| 久久青草综合色| 又黄又粗又硬又大视频| 黑人巨大精品欧美一区二区mp4| 建设人人有责人人尽责人人享有的| 大型黄色视频在线免费观看| 精品一区二区三卡| 亚洲成av片中文字幕在线观看| 久久久国产精品麻豆| 亚洲黑人精品在线| 99香蕉大伊视频| 欧美成人免费av一区二区三区 | 黄色成人免费大全| 久久人人97超碰香蕉20202| av在线播放免费不卡| 99在线人妻在线中文字幕 | 欧美日韩国产mv在线观看视频| 水蜜桃什么品种好| 国产欧美日韩一区二区三区在线| 久久国产精品男人的天堂亚洲| 黄片小视频在线播放| 99国产综合亚洲精品| 18禁裸乳无遮挡免费网站照片 | 久久99一区二区三区| 欧美成狂野欧美在线观看| 两人在一起打扑克的视频| 国产人伦9x9x在线观看| 熟女少妇亚洲综合色aaa.| xxxhd国产人妻xxx| av不卡在线播放| 日本黄色视频三级网站网址 | 很黄的视频免费| 中文字幕最新亚洲高清| 午夜福利影视在线免费观看| 在线观看www视频免费| 天天躁狠狠躁夜夜躁狠狠躁| 99精品欧美一区二区三区四区| 亚洲欧美一区二区三区黑人| 免费观看a级毛片全部| 久久久国产一区二区| av国产精品久久久久影院| 一进一出抽搐gif免费好疼 | 中文字幕人妻丝袜一区二区| 国产av精品麻豆| 嫁个100分男人电影在线观看| 夜夜夜夜夜久久久久| av天堂久久9| 亚洲一区二区三区不卡视频| 成人特级黄色片久久久久久久| 国产亚洲精品第一综合不卡| 亚洲av美国av| 在线看a的网站| av福利片在线| 黑人猛操日本美女一级片| 黑人操中国人逼视频| 国产欧美日韩一区二区三| 99re6热这里在线精品视频| 日韩欧美免费精品| 国产熟女午夜一区二区三区| 国产激情久久老熟女| 国产精品电影一区二区三区 | 久久99一区二区三区| 国产成人一区二区三区免费视频网站| 欧美精品av麻豆av| 欧美 亚洲 国产 日韩一| 精品电影一区二区在线| 久久精品国产亚洲av香蕉五月 | 美女福利国产在线| 午夜两性在线视频| 日韩 欧美 亚洲 中文字幕| 一个人免费在线观看的高清视频| 大香蕉久久成人网| 亚洲国产中文字幕在线视频| av电影中文网址| 午夜精品国产一区二区电影| a级片在线免费高清观看视频| 亚洲av第一区精品v没综合| 午夜老司机福利片| 黄色怎么调成土黄色| 久久国产精品影院| 亚洲国产毛片av蜜桃av| 久久香蕉国产精品| 欧美日韩成人在线一区二区| 91麻豆av在线| 91九色精品人成在线观看| 欧美av亚洲av综合av国产av| 老鸭窝网址在线观看| 精品少妇久久久久久888优播| tube8黄色片| av超薄肉色丝袜交足视频| 人妻丰满熟妇av一区二区三区 | 少妇被粗大的猛进出69影院| 香蕉丝袜av| 久久这里只有精品19| 50天的宝宝边吃奶边哭怎么回事| 一级黄色大片毛片| 精品一品国产午夜福利视频| 少妇裸体淫交视频免费看高清 | 欧美成人午夜精品| 精品人妻熟女毛片av久久网站| 久久九九热精品免费| svipshipincom国产片| 岛国在线观看网站| 波多野结衣av一区二区av| 国产男女超爽视频在线观看| 人人妻人人澡人人看| 亚洲色图综合在线观看| 丰满饥渴人妻一区二区三| 美女高潮到喷水免费观看| 大片电影免费在线观看免费| 热re99久久精品国产66热6| 精品国产一区二区三区四区第35| 99精国产麻豆久久婷婷| x7x7x7水蜜桃| 丁香六月欧美| 欧美亚洲日本最大视频资源| 午夜福利,免费看| 最新的欧美精品一区二区| 日韩欧美在线二视频 | 中文字幕人妻丝袜制服| 亚洲成a人片在线一区二区| 王馨瑶露胸无遮挡在线观看| 久久人人97超碰香蕉20202| 欧美激情极品国产一区二区三区| 国产av又大| 精品久久久久久久久久免费视频 | 国产在线精品亚洲第一网站| 一边摸一边做爽爽视频免费| 一进一出好大好爽视频| 亚洲一区二区三区欧美精品| 国内久久婷婷六月综合欲色啪| 欧美日韩精品网址| 在线观看免费视频网站a站| 国产成人精品久久二区二区91| 亚洲第一欧美日韩一区二区三区| 高清欧美精品videossex| 在线视频色国产色| 国产伦人伦偷精品视频| 可以免费在线观看a视频的电影网站| 国产在视频线精品| 欧美av亚洲av综合av国产av| 男女之事视频高清在线观看| 久久久国产成人精品二区 | 久久精品亚洲精品国产色婷小说| 国产亚洲精品第一综合不卡| 91av网站免费观看| 精品欧美一区二区三区在线| 久久久久国产精品人妻aⅴ院 | 18禁裸乳无遮挡动漫免费视频| 成人永久免费在线观看视频| 免费黄频网站在线观看国产| 国产亚洲精品久久久久5区| 成年女人毛片免费观看观看9 | 国产日韩欧美亚洲二区| 国内毛片毛片毛片毛片毛片| 首页视频小说图片口味搜索| 久久精品亚洲熟妇少妇任你| 精品一区二区三卡| 亚洲国产欧美一区二区综合| 国产亚洲欧美精品永久| 欧美激情久久久久久爽电影 | 国产成人啪精品午夜网站| 一二三四在线观看免费中文在| 在线永久观看黄色视频| 国产精华一区二区三区| 免费在线观看日本一区| 三上悠亚av全集在线观看| 国产一区二区激情短视频| 免费在线观看亚洲国产| av电影中文网址| 国产真人三级小视频在线观看| 电影成人av| 国产日韩一区二区三区精品不卡| 老汉色∧v一级毛片| 精品亚洲成a人片在线观看| 一个人免费在线观看的高清视频| 国产区一区二久久| 建设人人有责人人尽责人人享有的| 日韩欧美免费精品| 精品少妇久久久久久888优播| 天天躁夜夜躁狠狠躁躁| 大陆偷拍与自拍| 欧美成人免费av一区二区三区 | 国产淫语在线视频| 中国美女看黄片| 99国产精品99久久久久| 精品人妻在线不人妻| 中文欧美无线码| 国产1区2区3区精品| 色在线成人网| 欧美精品人与动牲交sv欧美| av天堂久久9| 飞空精品影院首页| 国产av精品麻豆| 一级片'在线观看视频| 女同久久另类99精品国产91| 久久中文字幕人妻熟女| 成人影院久久| 男人舔女人的私密视频| 黄频高清免费视频| 欧美国产精品一级二级三级| 欧美日韩黄片免| 91成年电影在线观看| 成年女人毛片免费观看观看9 | 久久精品熟女亚洲av麻豆精品| 高清毛片免费观看视频网站 | 国产成人一区二区三区免费视频网站| 日韩免费高清中文字幕av| 两个人免费观看高清视频| 男人操女人黄网站| 国产亚洲精品一区二区www | 婷婷丁香在线五月| 又紧又爽又黄一区二区| 大香蕉久久网| 91麻豆av在线| 妹子高潮喷水视频| 大陆偷拍与自拍| 国产精品 欧美亚洲| 12—13女人毛片做爰片一| 欧美日韩乱码在线| 久久影院123| 十八禁高潮呻吟视频| a级毛片在线看网站| 亚洲精品av麻豆狂野| 国产一区有黄有色的免费视频| 日韩大码丰满熟妇| 亚洲中文日韩欧美视频| 国产精品香港三级国产av潘金莲| 午夜精品在线福利| 欧美乱色亚洲激情| 亚洲视频免费观看视频| 久久这里只有精品19| 九色亚洲精品在线播放| 大香蕉久久网| 亚洲avbb在线观看| 国产一区二区三区视频了| 久热这里只有精品99| 国产黄色免费在线视频| 亚洲国产欧美网| 少妇被粗大的猛进出69影院| 无人区码免费观看不卡| 手机成人av网站| 国产精品.久久久| 法律面前人人平等表现在哪些方面| 悠悠久久av| 精品免费久久久久久久清纯 | 欧美午夜高清在线| 精品久久久久久电影网| 侵犯人妻中文字幕一二三四区| videosex国产| 高清在线国产一区| 99re在线观看精品视频| 国产成人一区二区三区免费视频网站| 777久久人妻少妇嫩草av网站| 男女床上黄色一级片免费看| 久久草成人影院| 涩涩av久久男人的天堂| 免费在线观看亚洲国产| 嫩草影视91久久| 国产精品美女特级片免费视频播放器 | 五月开心婷婷网| 亚洲成国产人片在线观看| 成人特级黄色片久久久久久久| 午夜福利在线免费观看网站| 高清在线国产一区| av中文乱码字幕在线| 欧美日韩成人在线一区二区| 精品久久久精品久久久| 亚洲,欧美精品.| 在线av久久热| 美女视频免费永久观看网站| 亚洲国产精品一区二区三区在线| 亚洲 欧美一区二区三区| 国产精品久久久av美女十八| tocl精华| 不卡av一区二区三区| 国产午夜精品久久久久久| 在线播放国产精品三级| 国产在线一区二区三区精| 麻豆国产av国片精品| 国产精品1区2区在线观看. | 亚洲欧美激情综合另类| 精品国产超薄肉色丝袜足j| 香蕉丝袜av| 久久人人爽av亚洲精品天堂| 欧美av亚洲av综合av国产av| av福利片在线| 午夜激情av网站| 露出奶头的视频| 国产片内射在线| 国产日韩欧美亚洲二区| 国产精华一区二区三区| 国产精品久久久久成人av| 十八禁网站免费在线| 久久久国产一区二区| 18禁观看日本| 制服诱惑二区| 久久香蕉精品热| 老鸭窝网址在线观看| 亚洲av片天天在线观看| 色94色欧美一区二区| 777米奇影视久久| 十分钟在线观看高清视频www| 在线观看免费视频日本深夜| 桃红色精品国产亚洲av| 不卡一级毛片| 99国产综合亚洲精品| 久久午夜综合久久蜜桃| 老司机靠b影院| 男女之事视频高清在线观看| 精品第一国产精品| 国产男靠女视频免费网站| 欧美亚洲 丝袜 人妻 在线| 久久精品国产综合久久久| 少妇被粗大的猛进出69影院| 国产成人一区二区三区免费视频网站| av免费在线观看网站| 露出奶头的视频| 高清在线国产一区| 91成人精品电影| 久久性视频一级片| 两性午夜刺激爽爽歪歪视频在线观看 | 久久狼人影院| 在线观看免费午夜福利视频| 9热在线视频观看99| 人人妻人人澡人人爽人人夜夜| 搡老乐熟女国产| 十八禁高潮呻吟视频| 高清毛片免费观看视频网站 | 国产精品美女特级片免费视频播放器 | 国产精品偷伦视频观看了| av免费在线观看网站| 国产成人精品无人区| 亚洲自偷自拍图片 自拍| 国产精品一区二区在线观看99| 99精品在免费线老司机午夜| 久久久水蜜桃国产精品网| 久久国产亚洲av麻豆专区| 欧美日韩中文字幕国产精品一区二区三区 | 女性被躁到高潮视频| 精品国产乱码久久久久久男人| 欧美激情久久久久久爽电影 | 老司机亚洲免费影院| 老司机影院毛片| 精品国内亚洲2022精品成人 | 国产黄色免费在线视频| 欧美日韩av久久| 69av精品久久久久久| 少妇 在线观看| 黄片大片在线免费观看| 国产高清视频在线播放一区| 午夜老司机福利片| 天天影视国产精品| 高清在线国产一区| 欧美黑人欧美精品刺激| 在线国产一区二区在线| 99re在线观看精品视频| 国产aⅴ精品一区二区三区波| 黄色女人牲交| 欧美国产精品一级二级三级| av欧美777| www.精华液| 一区二区日韩欧美中文字幕| 亚洲国产中文字幕在线视频| 久久狼人影院| 美女扒开内裤让男人捅视频| 成年动漫av网址| www.精华液| 80岁老熟妇乱子伦牲交| 国产人伦9x9x在线观看| 欧美黄色片欧美黄色片| 亚洲 欧美一区二区三区| 国产一区二区三区在线臀色熟女 | av免费在线观看网站| 亚洲熟女精品中文字幕| ponron亚洲| 欧美一级毛片孕妇| 人妻一区二区av| 韩国av一区二区三区四区| 母亲3免费完整高清在线观看| 欧美成人午夜精品| 99久久国产精品久久久| 日韩成人在线观看一区二区三区| 校园春色视频在线观看| 人妻一区二区av| 成人精品一区二区免费| 国产成人精品久久二区二区91| 一二三四在线观看免费中文在| 久久久国产成人免费| 色在线成人网| 久久人妻av系列| 19禁男女啪啪无遮挡网站| 无限看片的www在线观看| a在线观看视频网站| 日韩欧美国产一区二区入口| svipshipincom国产片| 国产欧美亚洲国产| 欧美 亚洲 国产 日韩一| 久久性视频一级片| 免费看a级黄色片| 久久久久精品国产欧美久久久| 国产一区二区激情短视频| 日韩免费高清中文字幕av| 最近最新中文字幕大全免费视频| 欧美日韩精品网址| 嫁个100分男人电影在线观看| 美女视频免费永久观看网站| 18禁黄网站禁片午夜丰满| 午夜影院日韩av| 大片电影免费在线观看免费| 日韩欧美一区二区三区在线观看 | 少妇 在线观看| 香蕉国产在线看| 丝瓜视频免费看黄片| av片东京热男人的天堂|