肋软骨是主要的自体软骨来源[1]，肋软骨细胞作为种子细胞[2]，在细胞治疗[3,4,5,6]、组织工程等领域都有广泛应用[7,8,9]。因此，对肋软骨细胞进行深度测序，表征细胞异质性图谱，甚至发现其中全新的细胞类型，对其临床应用具有重要指导意义。

复杂的生物系统由单个细胞的协同功能共同决定，传统方法虽然提供了大量的基因组或转录组数据，但却无法揭示造成这种复杂性的细胞异质性[10]。新一代的高通量测序方法——单细胞转录组测序(single-cellRNA-sequencing,scRNA-seq)可以对单个细胞的DNA或RNA分子进行深度测序，是研究单个目标细胞中基因组、转录组、非编码RNA的测序方法[11]。单细胞转录组测序是使用整个转录组的数千个细胞去鉴定细胞类型，它的分析有助于深入理解干细胞多能性与发育分化[12]、细胞重编程[13]、疾病发生与治疗等过程的基因调节网络[14]。骨关节置换术后的软骨样本的单细胞测序结果表明，关节软骨中有7个不同状态的软骨细胞群：增殖性软骨细胞群、肥大软骨细胞群、肥大前软骨细胞群、纤维软骨细胞、效应软骨细胞、调控软骨细胞及自我平衡软骨细胞，拟时间序列分析发现增殖软骨细胞比肥大软骨细胞发育更早期，而纤维软骨在这些细胞中处于发育的最晚期[15]。对健康的半月板组织的单细胞测序，聚类得到7个细胞亚群，分别为软骨细胞祖细胞、纤维软骨细胞、预肥大软骨细胞及纤维软骨细胞祖细胞、增殖纤维软骨细胞、调节性软骨细胞、内皮细胞[16]。这些文章不仅明确了每个亚群的标记基因以及主要的生物功能，深度挖掘了软骨发育轨迹，还筛选出了相关疾病的候选基因。

该研究对肋软骨组织来源细胞的表达数据进行单细胞转录组测序分析，将不同的细胞聚类为不同的亚群，并对其进行注释，得到各亚群细胞的标记基因。随后，对于每个细胞亚群的标记基因进行GO分析和KEGG信号通路富集分析，完成肋软骨细胞的亚群鉴定和功能分析。

1、材料和方法

1.1 设计

单细胞转录组测序分析。

1.2 时间及地点

实验于2019年11月至2021年1月在中国医学科学院整形外科医院研究中心完成。

1.3 材料

样本来源于1例31岁女性患者经小耳重建手术后废弃肋软骨组织，该研究获得了中国医学科医学院整形外科医院伦理委员会机构审查委员会的批准，患者签署知情同意书。

实验试剂和仪器：PBS、1%青霉素-链霉素、胰酶、4%锥虫蓝(Gibco，美国)；Ⅱ型胶原酶(Sigma，美国)；恒温摇床(上海智诚分析仪器有限公司)；S1000TMTouchThermalCycler(BioRad，美国)；Agilent4200(CapitalBioTechnology，北京)；Countess®ⅡAutomatedCellCounter(Invitrogen，美国)。

1.4 方法

1.4.1 人软骨细胞的分离和扩增

首先，从手术室获取新鲜肋软骨，并将其置于冰上的PBS和1%青霉素-链霉素的溶液中，在20min内运送到实验室。将肋软骨切成1mm3的小块，在含100mg/L青霉素和100mg/L链霉素的磷酸盐缓冲盐水溶液中洗涤，添加0.25%胰酶，在37℃，80r/min下摇动30min，弃去上清液，添加0.2%Ⅱ型胶原酶，在37℃，80r/min下摇动8h。使用锥虫蓝染色计数活细胞，如果原代细胞活力高于80%，将用10XGenomics平台对样品进行上机处理。

1.4.2 细胞捕获和单细胞文库的构建

将样本细胞悬液送至博奥晶典(北京)公司进行测序，使用微流体芯片将单个细胞与单个凝珠(beads)形成“油包水”的液滴，在液滴中裂解捕获的细胞，通过Oligo(dT)磁珠富集带有polyA尾的mRNA，随后在FragmentationBuffer中用二价阳离子将得到的mRNA随机打断。以片段化的mRNA为模版，寡核苷酸为引物，在M-MuLV反转录酶体系中合成cDNA第一条链，随后用RNaseH降解RNA链，并在DNApolymeraseⅠ体系下，以dNTPs为原料合成cDNA第二条链。纯化后的双链cDNA经过末端修复、加A尾并连接测序接头，用AMPureXPbeads筛选370-420bp的cDNA进行PCR扩增并再次使用AMPureXPbeads纯化PCR产物，最终获得文库。在测序的flowcell中加入4种荧光标记的dNTP、DNA聚合酶以及接头引物进行扩增，在每一个测序簇延伸互补链时，每加入一个被荧光标记的dNTP就能释放出相对应的荧光，测序仪通过捕获荧光信号，并用计算机软件将光信号转化为测序峰，从而获得待测片段的序列信息。

1.4.3 CellRanger测序

CellRanger软件是从10XGenomics网站(https://support.10Xgenomics.com/single-cell-geneexpression/software/downloads/)获得的。使用cellranger计数模块进行比对、过滤、条形码计数和单分子标签(uniquemolecularidentifier,UMI)计数，以生成特征条形码矩阵并确定聚类。使用主成分分析进行降维，T分布随机邻域嵌入可视化后通过K-means算法和Graph-based算法生成聚类。

1.4.4 主成分分析

主成分分析是一种减少数据集的维数，同时保持数据集之间最大差异的一种算法。这是通过保留低阶主成分，忽略高阶主成分做到的，这样低阶成分往往能够保留住数据的最重要方面。主成分分析限于线性维数，它不能解释特征之间的复杂多项式关系，但它降维相对较快，当与稀疏的矩阵一起使用时，它可以扩展到非常大的数据集。

1.4.5 T分布随机邻域嵌入

T分布随机邻域嵌入非线性降维算法是通过识别具有多个相似特征的数据点，投射到不同数据簇中的一种算法，它也是一种用于大型高维数据可视化的统计方法，它使用概率分布来估计嵌入的情况，将数据投射到各个孤立的簇中，实现细胞群的可视化。

1.4.6 K-meas算法

K-means聚类是一种划分聚类算法，计算原理是先随机选择k个对象作为凝聚点，按照就近原则将其余对象向凝聚点凝集，并计算其余对象与凝集点的距离，计算距离的方法一般选择欧氏距离。凝集结束后，计算出各个初始分类的中心位置，用计算出的中心位置重新进行聚类，K-means聚类通常适用于样本量较大的连续型变量，并且要求已知类别数。

1.4.7 簇细胞类型注释和标记基因鉴定

使用SingleR突出显示已知的标记基因可识别簇，分别对每个亚群与其他所有亚群的基因表达平均值进行比较，获得每个亚群较高表达的基因集合，即标记基因，使用VlnPlot函数生成标记基因的小提琴图。

1.5 主要观察指标

(1)碱基含量的分布；(2)基因数和单分子标签数比例；(3)特异性基因在亚群中表达量分布。

1.6 统计学分析

采用SPSS25.0统计软件分析。计量资料以表示，行独立样本t检验，P<0.05为差异有显著性意义。

2、结果

2.1 10XGenomics样品处理和cDNA文库制备

人肋软骨组织制备成单细胞悬液，使用Countess®ⅡAutomatedCellCounter进行细胞计数和细胞活率测定，细胞活率为85%，将细胞浓度调整为1×109L-1。利用IlluminaHiSeq测序平台进行测序和初级分析，进行数据可视化，在单细胞水平上生成了肋软骨组织中基因表达的深层转录图，获得6634个细胞进行分析，检测到单个细胞的单分子标签和基因中位数分别为16405和2956。81.8%的测序饱和率表明对可用转录本进行了全面取样。以上数据见图1A。

对数据分析之前进行了严格的数据质控，首先检测了碱基含量的分布，见图1B，其用于检测有无AT、GC分离现象。理论上普通文库的G和C碱基及A和T碱基含量每个测序循环上应分别相等，且整个测序过程稳定不变，呈水平线。由于随机引物扩增偏差等原因，常常会导致在测序得到的每个read前几个碱基有较大的波动，这种波动属于正常情况。由图1B中GC含量分布看出，未检测到AT、GC分离现象。在文库构建过程中可能会掺入死细胞，或者多个细胞被捕获在同一个液滴中，产生潜在的多重细胞，因此利用基因数和单分子标签数进行质控，见图1C，基因数大于5000、大部分细胞的单分子标签数大于10000，是一组高质量的测序数据。随着测序数据量的增加，单细胞检测到的基因数量和单分子标签的数量也随之增加，两者之间呈明显正相关，见图1D。细胞的单分子标签数与基因数的相关性可以确定是否有大量死亡或濒死细胞的线粒体污染，此次数据得分为0.89，符合要求，可以进行后续分析。

2.2 根据单细胞转录组测序数据进行细胞分群

为了深入剖析肋软骨细胞的异质性，首先对表达数据单分子标签归一化，然后主成分分析降维，见图2A，选取K-means聚类算法进行细胞分群，同时利用T分布随机邻域嵌入进行可视化的聚类分析，获得了4个聚类(Cluster)的细胞群，不同类的细胞亚群标记为不同的颜色，见图2B。每个亚群都有特异的标记基因，图中显示了各亚群中前20名的差异基因，见图2C。结合软骨细胞特征性基因在各群中的表达情况，发现Cluster1是以COL10A1、S100A2为标记基因的细胞群，这群细胞在该样本中占比36%;Cluster2是以BMP2、COL2A1基因为标记基因的细胞群，细胞数占总样本的36%;Cluster3是以FOS、JUN基因为标记基因的细胞群，细胞数占总样本的25%;Cluster4是以CD90、CD146基因为标记基因的细胞群，细胞数占总样本的数量最少，仅为3%，见图2D。

为了进一步区分每个聚类，基于基因平均表达量和亚群之间的差异表达基因进行协方差统计分析，研究发现标记基因均具有较高的特异性，可以识别细胞类型，见图3。Cluster1中肥大软骨基因COL10A1、钙调蛋白S100A2、抑制基质降解蛋白SERPINE2、软骨标记基因SOX9表达比较高，成软骨基因COL2A1、COL9A1、增殖调控基因FOS、JUN、BCL2，热调蛋白基因HSPA1A表达相对较低，因此这群细胞高表达肥大软骨细胞标记基因，将Cluster1注释为肥大软骨细胞群。Cluster2中成软骨基因COL2A1、COL9A1、SOX9、BMP2，钙调蛋白S100A2的基因表达相对较高，增殖性基因FOS、JUN，热调蛋白HSPA1A，干细胞基因CD90、CD146的表达均比较低，因此Cluster2是一群向软骨分化完全，但增殖能力不强的软骨细胞群。Cluster3中增殖性基因FOS、JUN，抑制凋亡的基因BCL2，调控细胞生长的基因HSPA1A、HAPA1B特异性高表达，软骨分化基因SOX9、COL2A1，干细胞基因CD146在Cluster3的表达相对较少，因此Cluster3被定义为一类增殖性细胞群。Cluster4中肥大软骨基因COL10A1，成软骨基因SOX9、COL2A1表达均较少，但调控细胞发育相关的基因MYLK、干细胞标记基因CD146、CD90等在Cluster4中特异性表达，将Cluster4注释为干细胞群。

2.3 细胞亚群标记基因参与的生物过程

为了进一步研究肋软骨组织中各细胞亚群相关的功能状态和潜在的调控因子，对每个细胞群中富集的差异表达基因进行GO和KEGG分析，见图4，结果显示：在被注释为肥大软骨细胞群的Cluster1中，生物学功能方面主要富集细胞间的黏附和对细胞迁移的调控；在细胞组分方面主要富集在细胞外基质、细胞外泌体之中；在分子功能方面主要富集蛋白质结合和细胞黏附等功能。在被注释为软骨细胞群的Cluster2中，生物学功能方面主要富集软骨细胞分化、Wnt信号通路的调节、细胞增殖和凋亡的调节；在细胞组分方面主要富集在细胞外泌体，细胞外基质中；在分子功能方面主要富集胶原结合、蛋白质结合、细胞外基质结合等功能。在被注释为增殖性软骨细胞群的Cluster3中，生物学功能主要富集信号传导、细胞周期调控等方面；细胞组分方面主要富集在细胞质、细胞核、细胞外基质等部位；在分子功能方面主要富集蛋白结合、ATP结合等方面。在被注释为干细胞的Cluster4中，生物学功能主要富集在细胞增殖、细胞分化、细胞黏附、胶原分解代谢等方面；细胞组分主要在细胞外基质、细胞外泌体、细胞膜等部位；分子功能主要富集在激活金属酶等方面。综上所述，在细胞组分方面4个亚群都主要富集在细胞外基质，在分子功能方面都主要与蛋白质结合有关。

3、讨论

肋软骨是人体组织中一类重要的软骨来源，与耳软骨、关节软骨相比，其来源更为丰富，取出一部分后对正常生理功能影响较小，常作为软骨自体移植和组织工程构建的种子细胞来源[17]。针对肋软骨细胞的深入研究对其作为种子细胞应用及疾病认识都有积极意义。

该研究通过单细胞转录组测序分析，初步将细胞分为4个亚群，包括肥大软骨细胞群、软骨细胞群、增殖性细胞群和干细胞群。与JI等[15]报道的有关骨关节炎单细胞测序的分群结果基本一致，其在研究中发现关节软骨中也含有肥大软骨细胞群、增殖性细胞群、前体细胞群。但在骨关节炎的细胞分群中有一群标记为稳态软骨细胞的亚群，在肋软骨中未发现。稳态软骨细胞群中特异性高表达调节生物节律的PER1基因和调节细胞周期的SIRT1基因[18,19]，但这2个基因在肋软骨各亚群中分布比较少且无特异性。骨关节炎测序结果还将软骨细胞细分出一群调节软骨细胞群，其基因集富集分析结果显示调节软骨细胞群具有大量与抗原处理和抗原呈递相关的基因表达，表明这些细胞可能具有免疫细胞功能。作者分析，肋软骨与关节软骨同样作为透明软骨，因此在细胞类型上有一致性，但文献中患骨关节炎的关节软骨是病损组织，因此在细胞分群上会有区别。

目前关于小耳畸形患者的肋软骨是否有发育异常的相关研究证实：小耳畸形患者除耳郭发育异常外，常伴有其他器官系统畸形，其中肋软骨畸形达32.6%[20]。该研究样本来源于31岁小耳畸形患者的肋软骨组织，在术前对其进行胸部CT和X射线片检查，排除患者肋软骨发育不全、肋软骨畸形、肋骨或肋软骨骨折、肋软骨高度钙化等情况，此样本肋软骨有轻微钙化，结合其年龄属于正常情况。从软骨细胞向肥大细胞分化进而骨化是透明软骨随年龄变化的一个特征[21]。COL2A1和SOX9是已知的软骨细胞标记基因，COL10A1是肥大软骨特异性基因，可促进细胞外基质的钙化[22]。根据这些特征基因，看到在31岁女性肋软骨中的成熟软骨细胞和肥大细胞各占到36%左右。软骨细胞肥大化是从成软骨分化到软骨基质矿化的中间过程，软骨细胞进入增殖、肥大阶段，表达COL10A1使软骨降解，会进一步钙化[23]。随着年龄增加，肥大软骨细胞亚群比例是否会进一步增加，还需要后续不同年龄样本的测序比较研究。该研究同时也为这些细胞发现新的标志物，比如在肥大软骨细胞中特异性高表达SERPINE2。该基因是丝氨酸蛋白酶抑制剂家族的一员，可以多种方式影响细胞外基质蛋白的代谢，调节细胞外基质微环境[24]。结合肥大软骨细胞的特性，SERPINE2可以作为肥大软骨细胞的候选标记基因。COMP在软骨细胞亚群中特异性高表达，该基因通过与其他细胞外基质蛋白(例如胶原蛋白和纤连蛋白)相互作用，在软骨的结构完整性中发挥作用[25],COMP可作为透明软骨细胞的标记基因。

在细胞亚群中，看到有一群软骨表型表达不明显，但是其增殖能力较强，在功能上与细胞周期相关的增殖性软骨细胞群。这群细胞不仅高表达调节细胞增殖的基因如BCL2、FOS，热休克蛋白家族基因如HSPA1A、HSPA1B等表达也具有特异性。该家族不仅是机体应激反应迅速产生的特异性蛋白，同时也可与端粒酶协同作用，具有重要的抗凋亡作用[26]，可作为增殖性细胞亚群的候选标记基因。增殖性软骨细胞占比达1/4左右，其中富集大量调控RNA代谢过程的基因，该群细胞在肋软骨细胞作为种子细胞应用时的作用应受到重视。

细胞分群中，还有极少比例的干细胞。作者发现注释为干细胞的亚群中既有CD146标记的内皮来源干细胞，也有CXCR4标记的神经来源干细胞，还有COL1A1标记的纤维母细胞，CD90标记的间质来源干细胞，提示在肋软骨发育中可能有多种干细胞参与。

胶原蛋白家族在标记细胞亚群中发挥了重要作用。一直以来，Ⅱ型胶原被认为是软骨细胞的特征性标志，肋软骨中软骨基质以Ⅱ型胶原蛋白为主，也含有少量的Ⅰ、Ⅸ型胶原蛋白，肋软骨胶原蛋白特别是Ⅱ型胶原蛋白的表达量及分布决定了其生物学特性[27]，Ⅱ型胶原向Ⅰ型胶原的转换，是被公认为去分化过程的标志[18]，此过程使软骨细胞转变为类似成纤维样细胞，进而失去软骨形成能力[28]。该研究中，COL2A1在Cluster2中表达量最多，在Cluster4中表达最少，COL1A1基因的表达情况正好相反，在Cluster2中表达最少，在Cluster4中表达最多，这个去分化过程是如何进行的，Cluster1和Cluster3在此过程中发挥了怎样的作用，可以后续通过Monocle轨迹分析进行发育轨迹方面的探索。

GO分析和Pathway通路分析显示，Cluster1-4的细胞组分大部分都在细胞外基质中，且生物功能也都与蛋白聚糖的合成有关，在组织工程化软骨的构建研究中，细胞外基质是细胞附着的支架[29]，为细胞提供生长和代谢的场所[30]。蛋白多糖聚集体是软骨两大特异性基质成分之一[31]，细胞外基质中蛋白多糖聚集体的含量可作为判断体外大量培养扩增软骨细胞分泌基质功能的可靠依据[32]，因此，可以从单细胞角度上判断肋软骨细胞是种植于组织工程支架的良好种子细胞。

该研究利用单细胞转录组测序技术对1例31岁女性肋软骨细胞进行分群，初步探讨了肋软骨细胞的异质性。既往研究表明，肋软骨随年龄的增长会出现不同程度的钙化[33,34]，不同年龄和性别会导致软骨细胞群体生物学特性有差异性[35,36]，吴迎等[37]对70例女性肋软骨进行双源CT三维重建发现肋软骨的钙化率随年龄增加，但并不呈正相关。GUO等[38]用320例连续小耳畸形患者的CT影像学数据进行分析，发现女性患者的软骨钙化率高于男性(P<0.05)。在6-15岁组中，男性和女性的钙化率相似，而其他3组中女性的钙化率均高于男性。此外，女性的中至重度钙化率高于男性，因此下一步将增加不同年龄、不同性别来源软骨细胞分群研究，进一步明确细胞亚群间的转归和影响，深入了解各细胞亚群的特性和潜能，对后续绘制肋软骨细胞的发育轨迹及作为种子细胞构建组织工程软骨提供研究基础，也为明确软骨相关疾病发生机制提供线索。

参考文献:

[8]靳小兵.自体组织工程生长板修复兔长骨生长板损伤的实验研究[D].西安:第四军医大学,2004.

[9]差疆.在体构建组织工程软骨的实验研究[D].西安:第四军医大学,2012.

[17]胡坤,陈竹,罗栩伟,等.组织工程软骨的研究新进展[J].西部医学,2020,32(6):927-932.

[20]杨美蓉.先天性小耳畸形伴发胸廓畸形的临床研究及初步遗传学分析[DJ].北京:协和医学院,2018.

[22]薛松,姜亚飞,桑伟林,等.肥大软骨细胞在骨关节炎发病中的作用[J].中国矫形外科杂志,2020,28(6).522-526.

[23]陈海涛,倪曲波,陈廖斌.关节软骨的发育及调控[J].武汉大学学报(医学版),2021,42(2):333-337.

[27]杜俊杰.罗卓荆,胡蕴玉,等.rhBMP-2在体内诱导成骨中I,II型胶原及碱性磷酸酶的表达[J].第四军医大学学报,2001,22(11)981-983.

[28]张艳,柴岗,刘伟,等.人肋软骨细胞体外培养中生长代谢及功能的变化[J].中华整形外科杂志,2004,20(5);372-376.

[33]赵琴琴,蔡震,游晓波,等.不同性别不同层次肋软骨的生物力学性能研究[J].中华整形外科杂志,2020,36(9):1042-1046.

[36]周志文,鞠黎,楼跃.不同年龄软骨组织多肽的定量分析[J].中国组织工程研究,2018,22(8):1178-1183.

[37]吴迎,王先成,熊祥,等.70例女性第6-8肋软骨组织量及钙化特点分析[J].中华整形外科杂志,2019,35(8):764-771.

文章来源：师航,李佳,刘霞,蒋海越.基于单细胞转录组测序分析人肋来源软骨细胞的异质性[J].中国组织工程研究,2022,26(19):3155-3161.