滕小娟， 林靜霞， 張永軍， 曾鑫權(quán)， 李雪飛， 何仁亮， 周家健

南方醫(yī)科大學皮膚病醫(yī)院，廣東 廣州 510091

皮膚是人體最大的器官，是保護人體免受脫水、紫外線(UV)輻射和傳染性病原體進入的最外層物理屏障[1]。在衰老過程中，人類皮膚組織及其組成細胞會逐漸發(fā)生形態(tài)和結(jié)構(gòu)改變，出現(xiàn)皮膚松弛、皺紋以及彈性等功能降低，由此引起一系列與衰老相關(guān)的皮膚疾病[2],從常見的慢性炎癥性疾病到可能危及生命的疾病(如基底細胞癌、鱗狀細胞癌和黑色素瘤)[3-5]。目前已證實免疫微環(huán)境的異常變化與多種組織和器官衰老有明顯相關(guān)性，其在維持組織器官平衡中發(fā)揮重要作用，而免疫微環(huán)境異常與免疫細胞在組織中的浸潤度有密切聯(lián)系[6-7]。對銀屑病和特應性皮炎患者皮膚組織進行單細胞轉(zhuǎn)錄組數(shù)據(jù)分析發(fā)現(xiàn)，各類免疫細胞類群發(fā)生了急劇變化，免疫微環(huán)境的失衡導致皮膚組織的損害[8-10]。由此推測，隨著皮膚組織衰老，免疫微環(huán)境失去平衡，導致皮膚組織各理化性質(zhì)發(fā)生變化，削弱了皮膚抵御外部環(huán)境變化的能力，但對衰老過程中皮膚免疫微環(huán)境失衡的整體描繪及其分子機制尚待進一步研究。本研究通過整合分析不同年齡段皮膚組織的RNA-seq和scRNA-seq數(shù)據(jù)，探究免疫細胞浸潤水平和免疫炎癥相關(guān)基因隨皮膚衰老的表達動態(tài)變化,為皮膚免疫微環(huán)境的分子機制研究和抗皮膚衰老藥物的開發(fā)提供理論依據(jù)。

1 材料與方法

1.1 材料

從基因型-組織表達研究項目(genotype-tissue expression，GTEx)數(shù)據(jù)庫下載正常皮膚組織的轉(zhuǎn)錄組數(shù)據(jù)(bulk RNA-seq)，包括376個來源于人體暴露皮膚組織和323個來源于非暴露處皮膚組織的轉(zhuǎn)錄組數(shù)據(jù)[11]。隨后，對年輕組和衰老組皮膚組織進行免疫細胞浸潤水平變化分析。從NCBI的基因表達綜合數(shù)據(jù)庫(gene expression omnibus，GEO)中獲取已公開的年輕人、老年人正常皮膚組織的單細胞RNA-seq數(shù)據(jù)，其中年輕皮膚樣本2個，衰老皮膚樣本3個(GSE130973)[12-14]。

1.2 方法

1.2.1 差異表達基因檢測和基因功能富集分析 將來自年齡為20～39歲個體的皮膚樣本定義為年輕組，60～79歲定義為衰老組。在暴露皮膚組織的RNA-seq數(shù)據(jù)集中，年輕組樣本116個，衰老組樣本260個；在非暴露皮膚組織的RNA-seq數(shù)據(jù)集中，年輕組樣本97個，衰老組樣本226個。采用DESeq2算法分析暴露皮膚組織和非暴露處皮膚組織中年輕組與衰老組的轉(zhuǎn)錄組數(shù)據(jù)集，篩選年輕組與衰老組比較后存在差異的表達基因[15]。如果一個基因在年輕組與衰老組的表達量差異∣log2(fold change)∣≥0.58且校正后P<0.05，即檢測為差異表達基因[16]。通過DAVID在線功能注釋工具對衰老組較年輕組差異表達上調(diào)和下調(diào)基因分別進行基因功能富集分析[17]。采用STEM軟件算法檢測基因在皮膚衰老過程中顯著變化的基因簇，并對各基因簇進行基因功能富集分析[18],以驗證差異表達基因是否隨皮膚組織衰老過程動態(tài)變化。

1.2.2 免疫細胞浸潤度分析 收集GTEx項目暴露皮膚組織和非暴露處皮膚組織的年輕、衰老皮膚組織轉(zhuǎn)錄組數(shù)據(jù)，利用CIBERSORT算法，以perm=1 000, QN=T為參數(shù)，對輔助T細胞，CD8陽性的效應T細胞、漿細胞、單核細胞等22種免疫細胞的浸潤度水平進行預測[19-20]；隨后，對所有樣本的22種免疫細胞的免疫細胞浸潤水平，采用R程序包進行主成分分析(principal component analysis,PCA)和熱圖聚類分析(heatmap)，對年輕和年老中免疫細胞浸潤度變化進行相關(guān)性分析[21]。

1.2.3 單細胞轉(zhuǎn)錄組數(shù)據(jù)分析 采用10×genomics的CellRanger軟件對從GEO數(shù)據(jù)庫下載的每個皮膚樣本的scRNA-seq數(shù)據(jù)進行初步處理[22]。先利用STAR2比對軟件對二代測序短序列片段進行比對，用每一個短序列片段的標簽序列，獲得每一細胞的轉(zhuǎn)錄組，并獲得樣本中每一個基因在每一個細胞的表達矩陣[23]。隨后，利用Seurat程序包中的NormalizeData、ScaleData將高通量測序比對所得的基因表達矩陣進行標準化；用FindVariableFeatures找出整個單細胞數(shù)據(jù)中差異變化最大的2 000個基因，并用findCluster對所有細胞進行聚類分析；最后使用RunUMAP將所有細胞映射到二維平面，進而對細胞分群及聚類，識別年輕皮膚和衰老皮膚的免疫細胞群的差異。通過Seurat包提取皮膚中特定的免疫細胞群，并對特定免疫細胞類群進行亞群分析[24-25]。

1.2.4 統(tǒng)計學處理 年輕組與衰老組的RNA-seq數(shù)據(jù)分析采用kruskal wallis檢驗進行基因差異表達比較；兩組scRNA-seq數(shù)據(jù)分析采用Seurat軟件包的t檢驗進行基因差異表達比較。P<0.05為差異具有統(tǒng)計學意義。

2 結(jié)果

2.1 皮膚組織中炎癥相關(guān)基因在衰老皮膚中轉(zhuǎn)錄表達水平上升

通過整合分析衰老皮膚非暴露皮膚組織的轉(zhuǎn)錄組數(shù)據(jù)與年輕非暴露皮膚組織的轉(zhuǎn)錄組數(shù)據(jù)，篩選得到402個顯著表達上調(diào)的基因和619個顯著表達下調(diào)的基因(圖1A)，表明衰老皮膚相對年輕皮膚的基因轉(zhuǎn)錄水平有明顯差異及其在衰老過程中扮演重要角色。隨后對非暴露衰老皮膚表達上調(diào)的基因進行基因本體(GO)功能富集分析，發(fā)現(xiàn)其主要富集在免疫反應和炎癥相關(guān)分子信號通路,例如ERK1/2信號通路、IL-1信號通路等(圖1B)，提示人體皮膚組織的炎癥水平隨衰老逐漸增強；而對表達下調(diào)的基因進行功能富集分析，發(fā)現(xiàn)其主要富集在B細胞受體信號通路活化和補體系統(tǒng)激活(圖1C)，表明衰老皮膚抵御外來微生物侵擾的能力下降。通過比較人體暴露處的年輕皮膚和衰老皮膚轉(zhuǎn)錄組數(shù)據(jù)，篩選出282個表達上調(diào)的基因和272個表達下調(diào)的基因(圖1D)。通過基因本體功能富集分析，發(fā)現(xiàn)上調(diào)基因主要富集在炎癥反應，免疫反應細胞趨化性等(圖1E)，而下調(diào)基因富集在與衰老相關(guān)的通路，表明暴露皮膚組織炎癥水平在人體衰老中逐漸增強(圖1F)。用STEM算法檢測基因隨暴露皮膚衰老變化的基因簇時，發(fā)現(xiàn)皮膚差異表達基因隨衰老而表達量逐步上升，顯示49個差異表達基因隨年齡增長表達水平升高(圖1G)，進一步基因本體功能富集分析發(fā)現(xiàn)其富集在趨化因子信號通路及免疫反應等通路(圖1H)?？傊狙芯客ㄟ^整合分析轉(zhuǎn)錄組學數(shù)據(jù)，發(fā)現(xiàn)暴露皮膚組織與非暴露皮膚組織存在一定的差異，從轉(zhuǎn)錄表達水平展示了炎癥相關(guān)基因隨人體衰老而逐漸增強。

圖1 皮膚組織中炎癥相關(guān)基因在衰老皮膚中轉(zhuǎn)錄表達水平上升 1A:火山圖展示在年輕和衰老非暴露皮膚的顯著差異表達基因；1B:在衰老皮膚較年輕非暴露皮膚的表達上調(diào)基因的GO富集分析結(jié)果；1C:在衰老皮膚較年輕非暴露皮膚的表達下調(diào)基因的GO富集分析結(jié)果；1D:火山圖展示年輕和衰老暴露皮膚的顯著差異表達基因；1E:在衰老皮膚較年輕暴露皮膚的表達上調(diào)基因的GO富集分析結(jié)果；1F:在衰老皮膚較年輕暴露皮膚的表達下調(diào)基因的GO富集分析結(jié)果；1G:箱線圖展示49個基因在暴露皮膚中隨年齡增加而基因表達量逐漸上調(diào)；1H:隨衰老逐漸上調(diào)的49個基因的GO富集分析結(jié)果

2.2 衰老皮膚組織中免疫細胞浸潤水平較年輕皮膚組織的異質(zhì)性大

利用CIBERSORT算法以暴露(圖2A)和非暴露(圖2B)皮膚組織的標準化基因表達矩陣為輸入預測皮膚組織中免疫細胞亞群比例。比較年輕組與衰老組的免疫細胞浸潤水平時，PCA顯示各免疫細胞浸潤水平數(shù)據(jù)并不能區(qū)分年輕組和衰老組里的個體，表明皮膚組織中免疫細胞的浸潤水平隨著衰老沒有發(fā)生明顯的變化,但衰老組可觀察到更多的離群值，表明免疫細胞浸潤水平在衰老皮膚里異質(zhì)性更強(圖2A、2B)。熱圖聚類分析顯示暴露和非暴露皮膚組織中不同免疫細胞亞群的浸潤程度與年齡段及性別分布均無顯著差異(圖2C、2D)。以上結(jié)果表明皮膚免疫細胞亞群的浸潤程度在不同年齡段改變程度不明顯，衰老皮膚組織的免疫浸潤異質(zhì)性較年輕皮膚組織要大。

圖2 衰老皮膚組織中免疫細胞浸潤水平較年輕皮膚組織的異質(zhì)性大 2A:主成分分析顯示暴露皮膚的年輕和衰老皮膚標本的免疫浸潤度差異不明顯，但是衰老皮膚標本的異質(zhì)性較年輕樣本大；2B:主成分分析顯示非暴露皮膚的年輕和衰老皮膚標本不能分開，但是衰老皮膚的標本異質(zhì)性較年輕樣本大；2C:熱圖分析顯示22種不同類型免疫細胞浸潤水平在暴露皮膚標本的分布情況；2D:熱圖分析顯示22種不同類型免疫細胞浸潤水平在非暴露皮膚標本的分布情況

2.3 衰老皮膚組織中T細胞類群數(shù)量較年輕皮膚組織顯著減少

對皮膚組織的單細胞轉(zhuǎn)錄組數(shù)據(jù)進行分析，根據(jù)年輕和衰老皮膚組織單細胞轉(zhuǎn)錄組數(shù)據(jù)進行細胞聚類，利用UMAP算法映射細胞類群到二維平面；在17個不同的皮膚細胞類群中發(fā)現(xiàn)0、13、16細胞群為巨噬細胞/樹突狀免疫細胞，第6細胞群為T細胞類群(圖3A)。進一步分析T細胞和巨噬細胞/樹突狀免疫細胞在年輕和衰老皮膚組織中的群體變化時，發(fā)現(xiàn)T細胞的相對比例在衰老皮膚組織中比年輕皮膚組織顯著下降；而巨噬細胞/樹突狀免疫細胞在衰老皮膚組織中的相對比例較年輕皮膚組織中高，說明T細胞的群體變化有可能使年老皮膚組織的炎癥水平顯著提高(圖3B)。為了進一步探究皮膚細胞亞群的組成，提取0、6、13、16細胞群進行重新細胞聚類分析，得到14個細胞類群并對細胞亞群信息注釋，發(fā)現(xiàn)多個免疫T細胞和巨噬細胞亞群(圖3C、3D)，且大多數(shù)的免疫細胞亞群發(fā)生了變化。特別引人注意的是，天然免疫細胞在衰老皮膚組織里急劇下降；調(diào)控T細胞(regulatory T cells)、輔助T細胞和CD8陽性的效應T細胞在衰老皮膚中也急劇下降，表明隨著皮膚衰老，效應免疫細胞類群數(shù)量不斷下降，降低了衰老皮膚應對外部刺激和侵染的能力(圖3E、3F)?？傊Y(jié)果表明人體皮膚衰老過程中，T免疫細胞類群和天然免疫細胞類群會不斷下降。

圖3 衰老皮膚組織中T細胞類群數(shù)量較年輕皮膚組織顯著減少 3A:UMAP算法展示人體年輕和衰老個體皮膚的不同細胞類群的分布情況，其中紅圈為免疫細胞類群；3B:柱形圖展示巨噬細胞和T細胞在衰老皮膚組織中細胞群體數(shù)量變化；3C:UMAP圖展示人體年輕和衰老個體皮膚的免疫細胞群體的分布圖譜；3D:各免疫細胞亞群特異性表達基因在各免疫細胞中轉(zhuǎn)錄表達水平；3E:UMAP圖展示免疫細胞亞群在年輕和衰老皮膚中的分布情況；3F:皮膚組織中T細胞亞群數(shù)量在衰老皮膚組織中顯著減少

2.4 衰老皮膚組織中免疫相關(guān)基因在非免疫細胞顯現(xiàn)表達上調(diào)趨勢

為了從單細胞水平上探索皮膚衰老過程中免疫相關(guān)基因的表達變化，本研究探討這些因子在年輕皮膚和衰老皮膚組織表達情況。小提琴圖分析結(jié)果顯示：IRF8、IL-1RN、IL-33、CCL21、CXCL3和CXCL1的表達量水平在衰老皮膚組織較年輕皮膚組織中升高，表明這些基因的轉(zhuǎn)錄表達差異可能對衰老皮膚的免疫水平造成影響(圖4A～4C)。進一步觀察IL-1RN、IL-33、CXCL1的細胞來源，發(fā)現(xiàn)這些基因的表達來源主要是皮膚組織中的非免疫細胞，如角質(zhì)形成細胞、成纖維細胞等。IL-1RN的表達主要來源于皮膚組織的角質(zhì)形成細胞和巨噬細胞/樹突狀細胞(圖4D)；IL-33的表達主要來源于皮膚組織的血管上皮細胞和成纖維細胞(圖4E)；CXCL1的表達主要來源于皮膚組織的成纖維細胞和巨噬細胞/樹突狀細胞(圖4F),表明皮膚組織中非免疫細胞表達多種免疫相關(guān)因子，可能參與皮膚組織免疫微環(huán)境平衡調(diào)節(jié)。炎癥相關(guān)因子表達水平分析顯示，皮膚組織中非免疫細胞的免疫相關(guān)基因表達水平在衰老皮膚組織中顯著升高。

圖4 衰老皮膚組織中非免疫細胞的炎癥相關(guān)基因顯現(xiàn)表達上調(diào)趨勢 4A～4C:小提琴圖展示IRF8、IL-1RN、IL-33、CCL21、CXCL3和CXCL1在年輕和衰老皮膚各細胞類群中的轉(zhuǎn)錄表達水平；4D～4F:UMAP算法展示IL-1RN、IL-33和CXCL1在年輕和衰老皮膚組織各組織類群的轉(zhuǎn)錄表達水平

3 討論

皮膚由三個主要層組成：表皮、真皮和皮下組織。最外層的表皮是抵御外部環(huán)境的物理屏障[26]。表皮主要由角質(zhì)形成細胞組成，包括常駐基底細胞(BC)、棘細胞(SC)和顆粒細胞(GC)。此外，表皮含有抗紫外線的黑素細胞和朗格漢斯細胞(LC)，它們介導外來病原體的免疫識別[27-28]。而真皮組織充滿復雜的細胞外基質(zhì)和真皮成纖維細胞，為皮膚提供強度和彈性，也充滿了免疫細胞和內(nèi)皮細胞[29]。連接表皮和真皮是一個復雜的基底膜區(qū)域，對細胞間通訊和連接至關(guān)重要[30]。因此，皮膚組織中各細胞類群和細胞外基質(zhì)共同組成皮膚組織微環(huán)境，特別是由免疫細胞類群介導的免疫微環(huán)境對皮膚組織穩(wěn)態(tài)維持發(fā)揮重要作用。隨著皮膚的衰老，表皮變薄，黑素細胞和LCs數(shù)量逐步減少，皮膚組織細胞微環(huán)境失衡，免疫微環(huán)境遭到破壞，但是隨皮膚衰老細胞類群的變化圖譜尚待進一步研究[10,31-32]。

目前已證實免疫微環(huán)境的異常變化與多種組織和器官衰老有明顯相關(guān)性，其在維持組織器官平衡中發(fā)揮重要作用，而免疫微環(huán)境異常與免疫細胞在組織中的浸潤度有密切聯(lián)系[6-7]。多項研究利用GTEx數(shù)據(jù)集發(fā)現(xiàn)多個與衰老相關(guān)的基因，例如在肌肉組織、皮膚組織、脂肪組織和腦組織等[33-35]。有研究通過單細胞轉(zhuǎn)錄組測序技術(shù)發(fā)現(xiàn)衰老皮膚相對于年輕皮膚組織中4個類群成纖維細胞的分泌特性和炎癥因子的表達發(fā)生了變化，顯示衰老過程中皮膚組織的細胞微環(huán)境發(fā)生了異常[14，36]。通過對眼瞼皮膚組織進行單細胞轉(zhuǎn)錄組分析發(fā)現(xiàn)，隨皮膚衰老，光老化損傷逐漸積累、皮膚炎癥水平不斷提高，而轉(zhuǎn)錄因子KLF6可能導致角質(zhì)形成細胞表現(xiàn)出炎癥因子表達升高和衰老的表型[37]。

本研究通過整合分析皮膚組織的RNA-seq和scRNA-seq數(shù)據(jù)，探究免疫細胞浸潤水平和免疫炎癥相關(guān)基因表達隨皮膚衰老過程的動態(tài)變化。利用GTEx項目中的皮膚轉(zhuǎn)錄組數(shù)據(jù)進行差異基因分析，發(fā)現(xiàn)暴露皮膚組織與非暴露皮膚組織存在一定的基因表達差異，其中暴露皮膚中隨衰老表達上調(diào)的基因與免疫炎癥水平相關(guān)；STEM分析也顯示隨衰老上調(diào)的基因與免疫炎癥水平相關(guān)，表明衰老過程中皮膚微環(huán)境發(fā)生了急劇變化，呈現(xiàn)出持續(xù)免疫炎癥狀態(tài)，可能與皮膚功能隨年齡的增長而引起的功能衰退和老年性皮膚疾病的發(fā)生發(fā)展有關(guān)。但基于轉(zhuǎn)錄組數(shù)據(jù)預測年輕和衰老的皮膚組織的免疫浸潤水平方法依賴于各細胞類群的特異表達基因表達程度，預測結(jié)果可能未反映免疫細胞亞群隨衰老過程的變化圖譜；通過scRNA-seq數(shù)據(jù)分析發(fā)現(xiàn)，衰老皮膚相對年輕皮膚巨噬細胞/樹突細胞，T細胞的類群發(fā)生了明顯的變化；通過對免疫細胞的分群，發(fā)現(xiàn)輔助T細胞、調(diào)節(jié)T細胞、自然殺傷細胞和CD8+T細胞在衰老皮膚里的群體顯著下降。本研究結(jié)果表明，人體皮膚衰老過程中T免疫細胞類群和天然免疫細胞類群會不斷下降，可能與衰老皮膚抗外部刺激、微生物侵染以及創(chuàng)傷后修復能力下降有關(guān)；炎癥相關(guān)因子表達水平分析表明，部分免疫相關(guān)因子在皮膚衰老過程中發(fā)生了表達變化，而這些基因的表達變化既來源于皮膚組織的免疫細胞，同時也來源于非免疫細胞，因此在皮膚組織衰老過程中，非免疫細胞的炎癥相關(guān)基因表達變化可能在皮膚組織免疫微環(huán)境平衡調(diào)節(jié)發(fā)揮重要的作用。

綜上所述，本研究表明在衰老過程中，輔助T細胞、調(diào)節(jié)T細胞、自然殺傷細胞和CD8+T細胞急劇下降，而巨噬細胞在衰老過程中上升，表明皮膚組織內(nèi)免疫細胞失衡在皮膚衰老過程中發(fā)揮重要作用。本研究描繪了皮膚組織衰老過程中免疫細胞轉(zhuǎn)錄組的精細景觀圖譜，反映出相關(guān)基因隨皮膚衰老的表達動態(tài)變化，有助于闡明復雜的人體皮膚老化生物學過程，為抗皮膚老化及相關(guān)疾病的新治療靶點或化合物的開發(fā)提供理論基礎(chǔ)。

本研究中詳細數(shù)據(jù)及資料請見https://doi.org/10.5281/zenodo.6638573。