多组学数据整合分析显示皮肤衰老过程中免疫微环境的改变特征

滕小娟， 林静霞， 张永军， 曾鑫权， 李雪飞， 何仁亮， 周家健

南方医科大学皮肤病医院，广东 广州 510091

皮肤是人体最大的器官，是保护人体免受脱水、紫外线(UV)辐射和传染性病原体进入的最外层物理屏障[1]。在衰老过程中，人类皮肤组织及其组成细胞会逐渐发生形态和结构改变，出现皮肤松弛、皱纹以及弹性等功能降低，由此引起一系列与衰老相关的皮肤疾病[2],从常见的慢性炎症性疾病到可能危及生命的疾病(如基底细胞癌、鳞状细胞癌和黑色素瘤)[3-5]。目前已证实免疫微环境的异常变化与多种组织和器官衰老有明显相关性，其在维持组织器官平衡中发挥重要作用，而免疫微环境异常与免疫细胞在组织中的浸润度有密切联系[6-7]。对银屑病和特应性皮炎患者皮肤组织进行单细胞转录组数据分析发现，各类免疫细胞类群发生了急剧变化，免疫微环境的失衡导致皮肤组织的损害[8-10]。由此推测，随着皮肤组织衰老，免疫微环境失去平衡，导致皮肤组织各理化性质发生变化，削弱了皮肤抵御外部环境变化的能力，但对衰老过程中皮肤免疫微环境失衡的整体描绘及其分子机制尚待进一步研究。本研究通过整合分析不同年龄段皮肤组织的RNA-seq和scRNA-seq数据，探究免疫细胞浸润水平和免疫炎症相关基因随皮肤衰老的表达动态变化,为皮肤免疫微环境的分子机制研究和抗皮肤衰老药物的开发提供理论依据。

1 材料与方法

1.1 材料

从基因型-组织表达研究项目(genotype-tissue expression，GTEx)数据库下载正常皮肤组织的转录组数据(bulk RNA-seq)，包括376个来源于人体暴露皮肤组织和323个来源于非暴露处皮肤组织的转录组数据[11]。随后，对年轻组和衰老组皮肤组织进行免疫细胞浸润水平变化分析。从NCBI的基因表达综合数据库(gene expression omnibus，GEO)中获取已公开的年轻人、老年人正常皮肤组织的单细胞RNA-seq数据，其中年轻皮肤样本2个，衰老皮肤样本3个(GSE130973)[12-14]。

1.2 方法

1.2.1 差异表达基因检测和基因功能富集分析 将来自年龄为20～39岁个体的皮肤样本定义为年轻组，60～79岁定义为衰老组。在暴露皮肤组织的RNA-seq数据集中，年轻组样本116个，衰老组样本260个；在非暴露皮肤组织的RNA-seq数据集中，年轻组样本97个，衰老组样本226个。采用DESeq2算法分析暴露皮肤组织和非暴露处皮肤组织中年轻组与衰老组的转录组数据集，筛选年轻组与衰老组比较后存在差异的表达基因[15]。如果一个基因在年轻组与衰老组的表达量差异∣log2(fold change)∣≥0.58且校正后P<0.05，即检测为差异表达基因[16]。通过DAVID在线功能注释工具对衰老组较年轻组差异表达上调和下调基因分别进行基因功能富集分析[17]。采用STEM软件算法检测基因在皮肤衰老过程中显著变化的基因簇，并对各基因簇进行基因功能富集分析[18],以验证差异表达基因是否随皮肤组织衰老过程动态变化。

1.2.2 免疫细胞浸润度分析 收集GTEx项目暴露皮肤组织和非暴露处皮肤组织的年轻、衰老皮肤组织转录组数据，利用CIBERSORT算法，以perm=1 000, QN=T为参数，对辅助T细胞，CD8阳性的效应T细胞、浆细胞、单核细胞等22种免疫细胞的浸润度水平进行预测[19-20]；随后，对所有样本的22种免疫细胞的免疫细胞浸润水平，采用R程序包进行主成分分析(principal component analysis,PCA)和热图聚类分析(heatmap)，对年轻和年老中免疫细胞浸润度变化进行相关性分析[21]。

1.2.3 单细胞转录组数据分析 采用10×genomics的CellRanger软件对从GEO数据库下载的每个皮肤样本的scRNA-seq数据进行初步处理[22]。先利用STAR2比对软件对二代测序短序列片段进行比对，用每一个短序列片段的标签序列，获得每一细胞的转录组，并获得样本中每一个基因在每一个细胞的表达矩阵[23]。随后，利用Seurat程序包中的NormalizeData、ScaleData将高通量测序比对所得的基因表达矩阵进行标准化；用FindVariableFeatures找出整个单细胞数据中差异变化最大的2 000个基因，并用findCluster对所有细胞进行聚类分析；最后使用RunUMAP将所有细胞映射到二维平面，进而对细胞分群及聚类，识别年轻皮肤和衰老皮肤的免疫细胞群的差异。通过Seurat包提取皮肤中特定的免疫细胞群，并对特定免疫细胞类群进行亚群分析[24-25]。

1.2.4 统计学处理 年轻组与衰老组的RNA-seq数据分析采用kruskal wallis检验进行基因差异表达比较；两组scRNA-seq数据分析采用Seurat软件包的t检验进行基因差异表达比较。P<0.05为差异具有统计学意义。

2 结果

2.1 皮肤组织中炎症相关基因在衰老皮肤中转录表达水平上升

通过整合分析衰老皮肤非暴露皮肤组织的转录组数据与年轻非暴露皮肤组织的转录组数据，筛选得到402个显著表达上调的基因和619个显著表达下调的基因(图1A)，表明衰老皮肤相对年轻皮肤的基因转录水平有明显差异及其在衰老过程中扮演重要角色。随后对非暴露衰老皮肤表达上调的基因进行基因本体(GO)功能富集分析，发现其主要富集在免疫反应和炎症相关分子信号通路,例如ERK1/2信号通路、IL-1信号通路等(图1B)，提示人体皮肤组织的炎症水平随衰老逐渐增强；而对表达下调的基因进行功能富集分析，发现其主要富集在B细胞受体信号通路活化和补体系统激活(图1C)，表明衰老皮肤抵御外来微生物侵扰的能力下降。通过比较人体暴露处的年轻皮肤和衰老皮肤转录组数据，筛选出282个表达上调的基因和272个表达下调的基因(图1D)。通过基因本体功能富集分析，发现上调基因主要富集在炎症反应，免疫反应细胞趋化性等(图1E)，而下调基因富集在与衰老相关的通路，表明暴露皮肤组织炎症水平在人体衰老中逐渐增强(图1F)。用STEM算法检测基因随暴露皮肤衰老变化的基因簇时，发现皮肤差异表达基因随衰老而表达量逐步上升，显示49个差异表达基因随年龄增长表达水平升高(图1G)，进一步基因本体功能富集分析发现其富集在趋化因子信号通路及免疫反应等通路(图1H)。总之，本研究通过整合分析转录组学数据，发现暴露皮肤组织与非暴露皮肤组织存在一定的差异，从转录表达水平展示了炎症相关基因随人体衰老而逐渐增强。

图1 皮肤组织中炎症相关基因在衰老皮肤中转录表达水平上升 1A:火山图展示在年轻和衰老非暴露皮肤的显著差异表达基因；1B:在衰老皮肤较年轻非暴露皮肤的表达上调基因的GO富集分析结果；1C:在衰老皮肤较年轻非暴露皮肤的表达下调基因的GO富集分析结果；1D:火山图展示年轻和衰老暴露皮肤的显著差异表达基因；1E:在衰老皮肤较年轻暴露皮肤的表达上调基因的GO富集分析结果；1F:在衰老皮肤较年轻暴露皮肤的表达下调基因的GO富集分析结果；1G:箱线图展示49个基因在暴露皮肤中随年龄增加而基因表达量逐渐上调；1H:随衰老逐渐上调的49个基因的GO富集分析结果

2.2 衰老皮肤组织中免疫细胞浸润水平较年轻皮肤组织的异质性大

利用CIBERSORT算法以暴露(图2A)和非暴露(图2B)皮肤组织的标准化基因表达矩阵为输入预测皮肤组织中免疫细胞亚群比例。比较年轻组与衰老组的免疫细胞浸润水平时，PCA显示各免疫细胞浸润水平数据并不能区分年轻组和衰老组里的个体，表明皮肤组织中免疫细胞的浸润水平随着衰老没有发生明显的变化,但衰老组可观察到更多的离群值，表明免疫细胞浸润水平在衰老皮肤里异质性更强(图2A、2B)。热图聚类分析显示暴露和非暴露皮肤组织中不同免疫细胞亚群的浸润程度与年龄段及性别分布均无显著差异(图2C、2D)。以上结果表明皮肤免疫细胞亚群的浸润程度在不同年龄段改变程度不明显，衰老皮肤组织的免疫浸润异质性较年轻皮肤组织要大。

图2 衰老皮肤组织中免疫细胞浸润水平较年轻皮肤组织的异质性大 2A:主成分分析显示暴露皮肤的年轻和衰老皮肤标本的免疫浸润度差异不明显，但是衰老皮肤标本的异质性较年轻样本大；2B:主成分分析显示非暴露皮肤的年轻和衰老皮肤标本不能分开，但是衰老皮肤的标本异质性较年轻样本大；2C:热图分析显示22种不同类型免疫细胞浸润水平在暴露皮肤标本的分布情况；2D:热图分析显示22种不同类型免疫细胞浸润水平在非暴露皮肤标本的分布情况

2.3 衰老皮肤组织中T细胞类群数量较年轻皮肤组织显著减少

对皮肤组织的单细胞转录组数据进行分析，根据年轻和衰老皮肤组织单细胞转录组数据进行细胞聚类，利用UMAP算法映射细胞类群到二维平面；在17个不同的皮肤细胞类群中发现0、13、16细胞群为巨噬细胞/树突状免疫细胞，第6细胞群为T细胞类群(图3A)。进一步分析T细胞和巨噬细胞/树突状免疫细胞在年轻和衰老皮肤组织中的群体变化时，发现T细胞的相对比例在衰老皮肤组织中比年轻皮肤组织显著下降；而巨噬细胞/树突状免疫细胞在衰老皮肤组织中的相对比例较年轻皮肤组织中高，说明T细胞的群体变化有可能使年老皮肤组织的炎症水平显著提高(图3B)。为了进一步探究皮肤细胞亚群的组成，提取0、6、13、16细胞群进行重新细胞聚类分析，得到14个细胞类群并对细胞亚群信息注释，发现多个免疫T细胞和巨噬细胞亚群(图3C、3D)，且大多数的免疫细胞亚群发生了变化。特别引人注意的是，天然免疫细胞在衰老皮肤组织里急剧下降；调控T细胞(regulatory T cells)、辅助T细胞和CD8阳性的效应T细胞在衰老皮肤中也急剧下降，表明随着皮肤衰老，效应免疫细胞类群数量不断下降，降低了衰老皮肤应对外部刺激和侵染的能力(图3E、3F)。总之，结果表明人体皮肤衰老过程中，T免疫细胞类群和天然免疫细胞类群会不断下降。

图3 衰老皮肤组织中T细胞类群数量较年轻皮肤组织显著减少 3A:UMAP算法展示人体年轻和衰老个体皮肤的不同细胞类群的分布情况，其中红圈为免疫细胞类群；3B:柱形图展示巨噬细胞和T细胞在衰老皮肤组织中细胞群体数量变化；3C:UMAP图展示人体年轻和衰老个体皮肤的免疫细胞群体的分布图谱；3D:各免疫细胞亚群特异性表达基因在各免疫细胞中转录表达水平；3E:UMAP图展示免疫细胞亚群在年轻和衰老皮肤中的分布情况；3F:皮肤组织中T细胞亚群数量在衰老皮肤组织中显著减少

2.4 衰老皮肤组织中免疫相关基因在非免疫细胞显现表达上调趋势

为了从单细胞水平上探索皮肤衰老过程中免疫相关基因的表达变化，本研究探讨这些因子在年轻皮肤和衰老皮肤组织表达情况。小提琴图分析结果显示：IRF8、IL-1RN、IL-33、CCL21、CXCL3和CXCL1的表达量水平在衰老皮肤组织较年轻皮肤组织中升高，表明这些基因的转录表达差异可能对衰老皮肤的免疫水平造成影响(图4A～4C)。进一步观察IL-1RN、IL-33、CXCL1的细胞来源，发现这些基因的表达来源主要是皮肤组织中的非免疫细胞，如角质形成细胞、成纤维细胞等。IL-1RN的表达主要来源于皮肤组织的角质形成细胞和巨噬细胞/树突状细胞(图4D)；IL-33的表达主要来源于皮肤组织的血管上皮细胞和成纤维细胞(图4E)；CXCL1的表达主要来源于皮肤组织的成纤维细胞和巨噬细胞/树突状细胞(图4F),表明皮肤组织中非免疫细胞表达多种免疫相关因子，可能参与皮肤组织免疫微环境平衡调节。炎症相关因子表达水平分析显示，皮肤组织中非免疫细胞的免疫相关基因表达水平在衰老皮肤组织中显著升高。

图4 衰老皮肤组织中非免疫细胞的炎症相关基因显现表达上调趋势 4A～4C:小提琴图展示IRF8、IL-1RN、IL-33、CCL21、CXCL3和CXCL1在年轻和衰老皮肤各细胞类群中的转录表达水平；4D～4F:UMAP算法展示IL-1RN、IL-33和CXCL1在年轻和衰老皮肤组织各组织类群的转录表达水平

3 讨论

皮肤由三个主要层组成：表皮、真皮和皮下组织。最外层的表皮是抵御外部环境的物理屏障[26]。表皮主要由角质形成细胞组成，包括常驻基底细胞(BC)、棘细胞(SC)和颗粒细胞(GC)。此外，表皮含有抗紫外线的黑素细胞和朗格汉斯细胞(LC)，它们介导外来病原体的免疫识别[27-28]。而真皮组织充满复杂的细胞外基质和真皮成纤维细胞，为皮肤提供强度和弹性，也充满了免疫细胞和内皮细胞[29]。连接表皮和真皮是一个复杂的基底膜区域，对细胞间通讯和连接至关重要[30]。因此，皮肤组织中各细胞类群和细胞外基质共同组成皮肤组织微环境，特别是由免疫细胞类群介导的免疫微环境对皮肤组织稳态维持发挥重要作用。随着皮肤的衰老，表皮变薄，黑素细胞和LCs数量逐步减少，皮肤组织细胞微环境失衡，免疫微环境遭到破坏，但是随皮肤衰老细胞类群的变化图谱尚待进一步研究[10,31-32]。

目前已证实免疫微环境的异常变化与多种组织和器官衰老有明显相关性，其在维持组织器官平衡中发挥重要作用，而免疫微环境异常与免疫细胞在组织中的浸润度有密切联系[6-7]。多项研究利用GTEx数据集发现多个与衰老相关的基因，例如在肌肉组织、皮肤组织、脂肪组织和脑组织等[33-35]。有研究通过单细胞转录组测序技术发现衰老皮肤相对于年轻皮肤组织中4个类群成纤维细胞的分泌特性和炎症因子的表达发生了变化，显示衰老过程中皮肤组织的细胞微环境发生了异常[14，36]。通过对眼睑皮肤组织进行单细胞转录组分析发现，随皮肤衰老，光老化损伤逐渐积累、皮肤炎症水平不断提高，而转录因子KLF6可能导致角质形成细胞表现出炎症因子表达升高和衰老的表型[37]。

本研究通过整合分析皮肤组织的RNA-seq和scRNA-seq数据，探究免疫细胞浸润水平和免疫炎症相关基因表达随皮肤衰老过程的动态变化。利用GTEx项目中的皮肤转录组数据进行差异基因分析，发现暴露皮肤组织与非暴露皮肤组织存在一定的基因表达差异，其中暴露皮肤中随衰老表达上调的基因与免疫炎症水平相关；STEM分析也显示随衰老上调的基因与免疫炎症水平相关，表明衰老过程中皮肤微环境发生了急剧变化，呈现出持续免疫炎症状态，可能与皮肤功能随年龄的增长而引起的功能衰退和老年性皮肤疾病的发生发展有关。但基于转录组数据预测年轻和衰老的皮肤组织的免疫浸润水平方法依赖于各细胞类群的特异表达基因表达程度，预测结果可能未反映免疫细胞亚群随衰老过程的变化图谱；通过scRNA-seq数据分析发现，衰老皮肤相对年轻皮肤巨噬细胞/树突细胞，T细胞的类群发生了明显的变化；通过对免疫细胞的分群，发现辅助T细胞、调节T细胞、自然杀伤细胞和CD8+T细胞在衰老皮肤里的群体显著下降。本研究结果表明，人体皮肤衰老过程中T免疫细胞类群和天然免疫细胞类群会不断下降，可能与衰老皮肤抗外部刺激、微生物侵染以及创伤后修复能力下降有关；炎症相关因子表达水平分析表明，部分免疫相关因子在皮肤衰老过程中发生了表达变化，而这些基因的表达变化既来源于皮肤组织的免疫细胞，同时也来源于非免疫细胞，因此在皮肤组织衰老过程中，非免疫细胞的炎症相关基因表达变化可能在皮肤组织免疫微环境平衡调节发挥重要的作用。

综上所述，本研究表明在衰老过程中，辅助T细胞、调节T细胞、自然杀伤细胞和CD8+T细胞急剧下降，而巨噬细胞在衰老过程中上升，表明皮肤组织内免疫细胞失衡在皮肤衰老过程中发挥重要作用。本研究描绘了皮肤组织衰老过程中免疫细胞转录组的精细景观图谱，反映出相关基因随皮肤衰老的表达动态变化，有助于阐明复杂的人体皮肤老化生物学过程，为抗皮肤老化及相关疾病的新治疗靶点或化合物的开发提供理论基础。

本研究中详细数据及资料请见https://doi.org/10.5281/zenodo.6638573。