人胸腺基质发育及退化的单细胞转录组分析

李雨晨，白丽莲，黄荷凤，刘欣梅

上海交通大学医学院附属国际和平妇幼保健院妇产科，上海市胚胎源性疾病重点实验室，上海 200030

胸腺是建立适应性免疫及自身免疫耐受的关键器官。T细胞在胸腺内的分化有赖于以胸腺上皮细胞（thymic epithelial cell，TEC）、成纤维细胞（fibroblast，Fb）和树突状细胞（dendritic cell，DC）等为代表的胸腺基质细胞形成功能健全的微环境［1］。在胸腺内，分化中的T细胞在TEC诱导下完成T细胞受体（T cell receptor，TCR）基因重排形成TCR组库，并识别由TEC和DC等抗原提呈细胞（antigen-presenting cells，APC）通过主要组织相容性复合物（major histocompatibility complex，MHC）提呈的自身抗原肽，作出适当强度的免疫反应，完成阳性选择和阴性选择，获得成熟的免疫反应功能，并清除自身反应性T细胞，建立免疫系统稳态［2］。

人类胸腺皮质和髓质TEC分化大约始于胚胎第8周，10～12周时皮、髓质分界逐渐清晰；12～15周时，胸腺淋巴细胞数量达到高峰，T细胞在趋化因子作用下向外周淋巴组织迁移［3］。胸腺的体积于青春期达到高峰，并在成年后逐渐退化［4］。了解胸腺的发育和退化，对于研究机体免疫稳态的维持、自身免疫疾病的发生规律和免疫功能重建等方面均有重要意义［5］。

由于正常人类胸腺样本不易获取，过去对胸腺发育和退化的研究主要来自小鼠等模式生物，并且通常局限于整个组织水平或分选的特定细胞类型［6-7］。近年来，单细胞RNA测序（single-cell RNA-sequencing，scRNA-seq）技术的发展使得我们能从少量组织样本中获得单细胞分辨率的转录组图谱，以揭示生物学过程中的细胞异质性［8-10］。尽管目前已有针对人类胸腺的scRNA-seq研究发表，且覆盖了从胚胎至成年的各个生命时期，但这些研究倾向于将来自不同发育时间点的样本进行整合分析，而忽略了由发育阶段变化造成的生物学差异［11-12］。

为了探索人类胸腺基质细胞在不同发育阶段的转录动态，我们利用ArrayExpress数据库中的人类胸腺单细胞转录组数据（E-MATB-8581）［12］进行生物信息学分析。根据样本来源的孕周（post-conception weeks，PCW）和出生后年龄将其划分为胚胎期、婴儿期和成人期，分析不同类型的基质细胞在各个发育阶段的基因表达差异，并结合富集分析推测潜在的功能变化。

1 资料与方法

1.1 单细胞转录组数据准备

从zenodo数据库下载DOI编号为10.5281/zenodo.3572422的E-MATB-8581［12］的“细胞-基因”表达矩阵，根据原数据提供的细胞类型注释，选取了DC细胞、TEC细胞、内皮细胞（endothelial cell，Endo）和Fb细胞等亚群进行后续分析。该公共数据代表的胸腺样本来自人类胚胎期（孕7～17周）、婴儿期（3、6、10、15、30月龄）及成人期（13～35岁），包含的生命阶段较完整，符合本研究的需要。

1.2 细胞分群

主要使用R语言Seurat对软件包表达矩阵进行分析［13］。保留符合以下条件的细胞数据：①检测到基因数>500。②检测到转录本数>2 000。通过Log-Normalize方法，设置10 000为放缩因子，对转录本计数进行归一化。使用R语言RSCORE软件包［14］，参考蛋白质-蛋白质相互作用（protein-protein interaction，PPI）数据库BioGRID（4.2.192）［15］，对10 000个高可变基因构建加权的相关性网络，并根据相关性将多个基因模块化，最终每个细胞的基因表达信息被转换为基因模块的评分值。基于该“细胞-模块”矩阵进行主成分分析（principal components analysis，PCA）降维。根据碎石图选取前25个主成分（principal components，PC），使用统一流形逼近与投影（uniform manifold approximation and projection，UMAP）进行非线性降维，在二维平面反映细胞间的相似程度［16］。

1.3 差异表达分析

对同一细胞亚型，以不同发育阶段作为分组变量，使用Wilcoxon秩和检验鉴定各个时间点的差异表达基因。为了减少样本间批次效应的影响，将发育轨迹划分为7～8孕周组（07-08PCW）、9～11孕周组（09-11PCW）、12～15孕周组（12-15PCW）、16～17孕周组（16-17PCW）、婴儿期组（Infant）和成人期组（Adult）。使用Bonferroni校正法进行多重检验校正，修正后P值<0.05认为差异具有统计学意义。选取各个阶段中差异表达满足|log2FC|>0.5的基因用于基因本体论（Gene Ontology，GO）功能分析和京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genome，KEGG）通路富集分析。

1.4 GO功能及KEGG通路富集分析

使用R语言clusterProfiler软件包对各个阶段富集表达的差异基因进行GO功能注释和KEGG通路富集分析［17］。使用Benjamini-Hochberg校正法校正多重检验P值，修正后P值<0.05时认为富集分析结果有统计学意义。

2 结果

2.1 人胸腺基质细胞的分群

基于UMAP结果，几种主要基质细胞在二维平面上能相互区分（图1A）。我们根据经典的细胞亚型标记基因验证了分群的可靠性（图1B）。TEC高表达泛上皮细胞标记基因EPCAM和TEC标记基因FOXN1［3］，其中皮质TEC（cortical TEC，cTEC）表达PSMB11。DC的亚群分别高表达DC1标记基因CLEC9A、DC2标记基因CLEC10A和活化DC（activated DC，aDC）标记基因LAMP3［18］。Endo高表达标记基因CDH5［19］。Fb则高表达PDGFRA和COLEC11，以及2型Fb标记基因FBN1［20］。

图1 人类胸腺基质细胞的UMAP分群图Fig 1 UMAPplots of human thymic stromacells

Fb主要来自胚胎阶段和成人阶段，并依据不同发育阶段分成离散的2群细胞。然而，来自同一阶段内不同时间点、不同建库策略的Fb细胞几乎完全混合，提示这种离散分群可能来源于胚胎期和成人期的生物学差异，而非批次效应。TEC主要以cTEC和Ⅰ型髓质TEC（medullary TEC，mTEC）细胞为主，其次是高表达自身免疫调控基因AIRE（autoimmune regulator）的Ⅱ型mTEC。尽管相同TEC亚型的细胞形成了相对集中且连续的分群，但来自不同时间点的细胞并未完全混合。类似地，来自胚胎期和成年期的Endo也形成了连续但不混合的分群。DC1、DC2和aDC亚群之间能清晰地区分，但来自不同发育阶段的同一亚群细胞则被均匀混合。以上结果提示，基于PPI构建基因相关性网络的分析方式，既有效地区分了不同细胞类型，同时揭示了可能由发育阶段造成的差异，且相比于DC，这种差异在TEC、Endo和Fb中更为显著。

2.2 APC随发育进程的转录动态

对差异基因的GO和KEGG分析显示，cTEC在胚胎早期主要富集包括核糖体参与的蛋白合成和加工等功能，而自9周龄开始直到成年则主要富集抗原提呈、自身免疫疾病相关的信号通路（图2A）。此外，在出生后，cTEC下调表达与黏着斑（focal adhesion）、PI3K-Akt信号通路相关的基因（图2B）。与cTEC略有不同，Ⅰ型mTEC在胚胎和婴儿期主要富集抗原提呈相关的信号通路，在成年期则富集与核糖体、mRNA代谢和细胞凋亡等有关的信号通路，并相应下调多种自身免疫疾病相关通路的基因（图3A）。Ⅱ型mTEC在胚胎期和婴儿期的功能富集与Ⅰ型mTEC相似，但在成年期的差异表达基因未发现显著上调的信号通路，仅发现少数下调的基因可能与糖代谢相关（图3B）。

图2 人类cTECs差异基因的KEGG通路分析Fig 2 The KEGG pathway analysis of differentially expressed genesof human cTECs

图3人类mTECs差异基因的KEGG通路分析Fig 3 KEGG pathway analysis of differentially expressed genes of human mTECs

2类DC细胞自胚胎期至成人期均富集包括抗原提呈或多种自身免疫疾病相关的通路。其中，与其他时期相比，孕7～8周时的DC1细胞下调表达1型糖尿病等疾病相关的基因，而婴儿期DC1细胞和胚胎期相比，还富集细胞凋亡相关的通路（图4A）；出生后的DC2细胞和胚胎期相比，还富集肿瘤坏死因子（tumor necrosis factor，TNF）信号相关的通路，并且下调表达和糖类或蛋白质代谢相关的基因（图4B）。

图4 人类胸腺DCs差异基因的KEGG通路分析Fig 4 KEGG pathway analysis of differentially expressed genes of human thymic DCs

我们进一步关注MHC分子，即人白细胞抗原（human leukocyte antigen，HLA）基因在各个时间点的表达差异。cTEC中的HLA基因大多在孕12周后表达上调，此后的发育阶段中无显著差异（图5A）。而早在孕9～11周时，Ⅰ型mTEC中的HLA基因即呈现表达上调，并持续到婴儿期；进入成年期后，Ⅰ型mTEC中的多种MHC-Ⅱ类分子基因，包括HLA-DQB1、HLA-DQA1、HLADRB5、HLA-DPA1、HLA-DPB1和HLA-DRB1等均显著下调，MHC-Ⅰ类分子基因（HLA-A、HLA-B、HLA-C等）的下调则较轻微（图5A）。HLA基因的表达模式在成年Ⅱ型mTEC细胞中未发现显著改变（图5A）。2类DC细胞的HLA分子均自胚胎早期就开始表达，在出生后达到峰值，并持续至成年（图5B）。

图5 人类胸腺APCs中HLA基因的表达差异Fig 5 Differential expression of HLA genesin human thymic APCs

2.3 Fb细胞随发育进程的转录动态

Fb 1型细胞在胚胎期主要上调胰岛素样生长因子2（insulin-like growth factor 2，IGF2）、转化生长因子β-1（transforming growth factor-β1，TGF-β1），在胚胎期和婴儿期上调COL1A1和COL3A1等胶原α链基因（图6A）。Fb 2型细胞的胶原基因表达谱和Fb 1型细胞类似。然而，2种Fb自婴儿期开始即上调MHC-Ⅰ类分子基因（HLAB、HLA-C）（图6A），并在成人期下调表达多种胶原基因（图6B）。此外，2种Fb在出生后阶段均上调多种趋化因子配体基因，包括CXCL14、CCL19和CCL21等。作为比较，我们发现Ⅰ型mTEC的CCL19和CCL21在成年期较胚胎和婴儿期下调（图6C）。

图6 人类胸腺Fb1和Fb2在不同发育阶段的差异表达基因Fig 6 Differentially expressed genes of human thymic Fb1 and Fb2 across development stages

GO和KEGG分析提示，2种Fb细胞在胚胎期的功能以核糖体相关的蛋白质合成或加工为主。Fb 1型细胞在婴儿期主要富集细胞外基质（extracellular matrix，ECM）相关的通路，而在成人期则富集抗原提呈相关的通路（图7A）。与Fb 1型细胞略有不同，Fb 2型细胞在婴儿期和成人期均富集抗原提呈和免疫反应相关的通路（图7B）。而2种Fb在成人期的下调基因均富集了以蛋白质消化吸收和Prion病等为主的信号通路。以上结果提示，在出生后的胸腺发育中，Fb从以一般基质细胞的支持功能为主，逐渐产生参与T细胞抗原识别、分化和激活的功能。

图7 人类胸腺Fb 1型和Fb 2型差异基因的KEGG通路分析Fig 7 KEGG pathway analysis of differentially expressed genes of human thymic Fb1 and Fb2

2.4 Endo随发育进程的转录动态

在胚胎期和婴儿期，Endo主要富集ECM和黏着斑等相关的通路。而自出生后开始，Endo还富集抗原提呈相关通路，特别是在成年期同时富集包括1型糖尿病、自身免疫性甲状腺疾病相关的通路，并且下调细胞黏附、ECM或PI3K-Akt信号通路相关的基因。HLA基因的表达谱显示Endo在出生后显著上调MHC-Ⅰ类和多种MHC-Ⅱ类分子基因（HLA-DRA、HLA-DRB1、HLA-DRB5、HLADPA1和HLA-DPB1等）（图8）。

图8 人类胸腺Endo在不同发育阶段的差异表达基因Fig 8 Differentially expressed genes of human thymic Endo across development stages

3 讨论

在2类MHC分子中，通常认为MHC-Ⅰ类分子主要与CD8分子结合，而MHC-Ⅱ主要结合CD4［21］。在胸腺内，cTEC主要通过提呈2类自身抗原肽-MHC复合物，刺激T细胞前体作出适当强度的免疫反应，使CD4/CD8双阳性（double positive，DP）细胞获得MHC限制性免疫反应，激活生存信号并分化为CD4或CD8单阳性（single positive，SP）细胞［2］。我们的分析发现，自胚胎期分化成熟以来，人类cTEC的抗原提呈功能一直持续到成年，同时2类MHC分子的表达均未见明显下调。这些结果提示，胸腺的生理性退化可能对主要由cTEC完成的阳性选择功能影响较小。

由于mTEC的抗原提呈主要用于筛选出非自身反应性T细胞，其抗原提呈功能下降可能使某些自身反应性T细胞从阴性选择中逃逸，进入外周造成自身免疫疾病［22］。过去对自身免疫疾病模型NZB小鼠的研究［23］发现，其胸腺中MHC-Ⅱ阳性的mTEC比例在生理性退化时的数量衰减较对照组更显著，提示了这类mTEC在自身免疫疾病发生中的潜在作用。我们的分析在单细胞分辨率水平更精确地推断了人类胸腺生理性退化时最易影响的潜在mTEC亚型（Ⅰ型mTEC），其在成年期显著下调MHC-Ⅱ类分子，抗原提呈功能较胚胎期和婴儿期有所下降。考虑不同MHC分子的亲和性，这一改变可能更容易影响CD4单阳性细胞，即诸如辅助性T细胞或调节性T细胞等的前体，进而影响青春期后的机体免疫稳态［24］。值得注意的是，过去普遍认为表达AIRE基因的Ⅱ型mTEC对诱导中心耐受、调控mTEC表达自身抗原肽等过程尤为关键［25-26］，然而，根据我们的功能富集分析，在成年期，该类细胞未见显著的转录改变，且AIRE基因本身和HLA基因的表达也无明显下调。因此，胸腺退化对Ⅱ型mTEC的影响有待设计实验进一步验证。

对于胸腺中的非APC细胞，我们主要分析了Fb和Endo。单细胞测序的细胞分群揭示了两者在出生前后的转录差异。胚胎期的Fb除了表达多种胶原基因参与构建ECM之外，还同时上调IGF2基因，这可能与其通过胰岛素样生长因子受体（insulin-like growth factor receptor，IGF-R）诱导胚胎期T细胞前体的发育有关［27］。IGF2在出生后的Fb中下调，可能提示这种诱导作用的下降。在成年期，Fb中多种胶原基因的下调，则可能使胸腺难以维持有序排列的小叶结构。由于MHC-Ⅰ类分子通常在人类大多数细胞类型中均有表达，其表达上调不足以说明Fb在出生后分化出抗原提呈的功能。然而，CCL19和CCL21是mTEC通过CCR7诱导SP细胞进入胸腺髓质的主要趋化因子［28-29］。与Ⅰ型mTEC逐渐下调CCL19和CCL21相对应，出生后的Fb上调这2个基因提示其具有类似的诱导T细胞在胸腺内迁移的功能，也可能与出生后诱导T细胞归巢或退化时的炎性反应有关。

与Fb不同，胸腺中的Endo在出生后显著上调2类MHC分子的表达，并可能由此参与T细胞的抗原识别。与MHC-Ⅰ类分子不同，MHC-Ⅱ类分子通常由专职的APC（professional APC）表达［30］。MHC-Ⅱ类分子在Endo中上调可由γ干扰素（γ-interferon，IFN-γ）诱导，并参与调控效应-记忆T细胞（effector memory T cells）的抗原识别［31］，可能是与出生后适应性免疫系统的构建相协调的改变。

综上所述，我们基于单细胞转录组数据，利用生物信息学方法，揭示了多种人类胸腺基质细胞在不同发育阶段中潜在的功能改变，推断了胸腺在出生后发育和退化过程中易受影响的细胞类型。胸腺生理性退化可能主要影响T细胞分化的阴性选择。这些结果为后续研究人体免疫系统稳态在生命早、中期的动态改变和胸腺微环境的体外重建提供了新的依据。

参·考·文·献

[1] Han J,Zúñiga-Pflücker JC.A 2020 view of thymus stromal cells in T cell development[J].JImmunol,2021,206(2):249-256.

[2] Starr TK,Jameson SC,Hogquist KA.Positive and negative selection of T cells[J].Annu Rev Immunol,2003,21:139-176.

[3] Farley AM,Morris LX,Vroegindeweij E,et al.Dynamics of thymus organogenesis and colonization in early human development[J].Development,2013,140(9):2015-2026.

[4] Nishino M,Ashiku SK,Kocher ON,et al.The thymus:a comprehensive review[J].Radiographics,2006,26(2):335-348.

[5] Palmer S,Albergante L,Blackburn CC,et al.Thymic involution and rising disease incidence with age[J].Proc Natl Acad Sci U S A,2018,115(8):1883-1888.

[6] Dik WA,Pike-Overzet K,Weerkamp F,et al.New insights on human T cell development by quantitative T cell receptor gene rearrangement studies and geneexpression profiling[J].JExp Med,2005,201(11):1715-1723.

[7] Mingueneau M,Kreslavsky T,Gray D,et al.The transcriptional landscape ofαβT cell differentiation[J].Nat Immunol,2013,14(6):619-632.

[8] Tang F,Barbacioru C,Wang Y,et al.mRNA-Seq whole-transcriptome analysis of asinglecell[J].Nat Methods,2009,6(5):377-382.

[9] Macosko EZ,Basu A,Satija R,et al.Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets[J].Cell,2015,161(5):1202-1214.

[10] Schier AF.Single-cell biology:beyond the sum of its parts[J].Nat Methods,2020,17(1):17-20.

[11] Zeng Y,Liu C,Gong Y,et al.Single-cell RNA sequencing resolves spatiotemporal development of pre-thymic lymphoid progenitors and thymus organogenesis in human embryos[J].Immunity,2019,51(5):930-948.e6.

[12] Park JE,Botting RA,Domínguez Conde C,et al.A cell atlas of human thymic development defines T cell repertoire formation[J].Science,2020,367(6480):eaay3224.

[13] Stuart T,Butler A,Hoffman P,et al.Comprehensive integration of singlecell data[J].Cell,2019,177(7):1888-1902.e21.

[14] Dong J,Zhou P,Wu Y,et al.Enhancing single-cell cellular state inference by incorporating molecular network features[EB/OL].(2019-10-15)[2021-03-20].https://www.biorxiv.org/content/10.1101/699959v2.

[15] Oughtred R,Rust J,Chang C,et al. The BioGRID database:a comprehensive biomedical resource of curated protein,genetic,and chemical interactions[J].Protein Sci,2021,30(1):187-200.

[16] Becht E,McInnes L,Healy J,et al.Dimensionality reduction for visualizing single-cell data using UMAP[J].Nat Biotechnol,2018:38-44.

[17] Yu G,Wang LG,Han Y,et al.clusterProfiler:an R package for comparing biological themes among gene clusters[J].OMICS,2012,16(5):284-287.

[18] Jardine L,Haniffa M. Reconstructing human DC,monocyte and macrophage development in utero using single cell technologies[J].Mol Immunol,2020,123:1-6.

[19] De Val S,Black BL.Transcriptional control of endothelial cell development[J].Dev Cell,2009,16(2):180-195.

[20] Thiele BJ,Doller A,Kähne T,et al.RNA-binding proteins heterogeneous nuclear ribonucleoprotein A1,E1,and K are involved in post-transcriptional control of collagenⅠandⅢsynthesis[J].Circ Res,2004,95(11):1058-1066.

[21] van der Merwe PA,Davis SJ.Molecular interactions mediating T cell antigen recognition[J].Annu Rev Immunol,2003,21:659-684.

[22] Thomas R,Wang W,Su DM.Contributions of age-related thymic involution to immunosenescenceand inflammaging[J].Immun Ageing,2020,17:2.

[23] Fletcher AL,Seach N,Reiseger JJ,et al.Reduced thymic aire expression and abnormal NF-κB2 signaling in a model of systemic autoimmunity[J].JImmunol,2009,182(5):2690-2699.

[24] Oh J,Wang W,Thomas R,et al.Capacity of tTreg generation is not impaired in the atrophied thymus[J].PLoSBiol,2017,15(11):e2003352.

[25] Cavadini P,Vermi W,Facchetti F,et al.AIRE deficiency in thymus of 2 patientswith omenn syndrome[J].JClin Invest,2005,115(3):728-732.

[26] Passos GA,Speck-Hernandez CA,Assis AF,et al.Update on aire and thymic negative selection[J].Immunology,2018,153(1):10-20.

[27] Kecha O,Brilot F,Martens H,et al.Involvement of insulin-like growth factors in early T cell development:a study using fetal thymic organ cultures[J].Endocrinology,2000,141(3):1209-1217.

[28] Kozai M,Kubo Y,Katakai T,et al.Essential role of CCL21 in establishment of central self-tolerancein T cells[J].JExp Med,2017,214(7):1925-1935.

[29] Halkias J,Melichar HJ,Taylor KT,et al.Tracking migration during human T cell development[J].Cell Mol Life Sci,2014,71(16):3101-3117.

[30] Wieczorek M,Abualrous ET,Sticht J,et al.Major histocompatibility complex(MHC)classⅠand MHC classⅡproteins:conformational plasticity in antigen presentation[J].Front Immunol,2017,8:292.

[31] Pober JS,Merola J,Liu R,et al.Antigen presentation by vascular cells[J].Front Immunol,2017,8:1907.