亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        類風(fēng)濕關(guān)節(jié)炎的潛在生物標(biāo)志物及其免疫調(diào)控機(jī)制:基于GEO數(shù)據(jù)庫(kù)

        2024-11-03 00:00:00陳莉莉吳天宇張銘丁子夏張妍楊依清鄭佳倩張小楠
        關(guān)鍵詞:生物標(biāo)志物類風(fēng)濕關(guān)節(jié)炎機(jī)器學(xué)習(xí)

        摘要:目的 探討類風(fēng)濕關(guān)節(jié)炎(RA)的潛在發(fā)病機(jī)制,挖掘可以作為RA早期診斷的生物標(biāo)志物,并探究其可能的免疫調(diào)控機(jī)制。方法 依托GEO數(shù)據(jù)庫(kù)多個(gè)基因組數(shù)據(jù),利用limma包、RRA方法、批次矯正方法與clusterProfiler包篩選差異表達(dá)基因并預(yù)測(cè)其作用功能;通過STRING在線數(shù)據(jù)庫(kù)得到蛋白質(zhì)互作網(wǎng)絡(luò),利用Cytoscape3.8.0軟件與GeneMANIA數(shù)據(jù)庫(kù)篩選差異表達(dá)關(guān)鍵基因并預(yù)測(cè)其互作機(jī)制;構(gòu)建ROC工作曲線驗(yàn)證關(guān)鍵基因診斷模型的信效度,并通過機(jī)器學(xué)習(xí)方法篩選疾病特征免疫細(xì)胞;使用corrplot包分析特征免疫細(xì)胞與疾病關(guān)鍵基因的相關(guān)性;使用GSEA富集方法探究關(guān)鍵基因生物學(xué)功能和相關(guān)性;通過構(gòu)建膠原誘導(dǎo)性關(guān)節(jié)炎(CIA)大鼠模型探究具有較高診斷價(jià)值的差異基因在RA滑膜組織中的表達(dá)情況。結(jié)果 得到與RA密切相關(guān)核心關(guān)鍵基因9個(gè),分別為:CD3G、CD8A、SYK、LCK、IL2RG、STAT1、CCR5、ITGB2、ITGAL,富集分析結(jié)果表明差異基因可能主要通過細(xì)胞因子相關(guān)通路調(diào)控滑膜炎癥;ROC工作曲線分析表明9個(gè)核心基因的預(yù)測(cè)模型信效度較高,其中STAT1預(yù)測(cè)模型的AUC值最高(0.909);相關(guān)性分析結(jié)果顯示,CD3G、ITGAL、LCK、CD8A和STAT1等5個(gè)核心基因與疾病特征免疫細(xì)胞相關(guān)性較強(qiáng),其中STAT1與M1型巨噬細(xì)胞的相關(guān)性聯(lián)系最強(qiáng)(R=0.68,Plt;0.001);CIA大鼠膝關(guān)節(jié)的HE、番紅-O固綠染色結(jié)果表明關(guān)節(jié)炎大鼠造模成功且藥物治療有效,免疫熒光結(jié)果表明STAT1 和磷酸化的STAT1 在CIA大鼠踝關(guān)節(jié)滑膜組織中表達(dá)高于正常組和治療組;免疫印跡實(shí)驗(yàn)結(jié)果表明,STAT1 蛋白表達(dá)量在滑膜成纖維細(xì)胞的細(xì)胞核(P=0.0002)與細(xì)胞質(zhì)(Plt;0.0001)中存在差異,且經(jīng)過治療,細(xì)胞核內(nèi)p-STAT1 與STAT1 的蛋白表達(dá)量顯著降低(Plt;0.0001)。結(jié)論 CD3G、CD8A、SYK、LCK、IL2RG、STAT1、CCR5、ITGB2 與ITGAL等9 個(gè)基因可作為RA的早期診斷生物標(biāo)志物,CD3G/CD8A/LCK-γδ T細(xì)胞、ITGAL-Tfh細(xì)胞、STAT1-M1型巨噬細(xì)胞等基因-免疫細(xì)胞途徑可能與RA的發(fā)生發(fā)展密切相關(guān)。

        關(guān)鍵詞:類風(fēng)濕關(guān)節(jié)炎;生物標(biāo)志物;免疫調(diào)控;生物信息學(xué);機(jī)器學(xué)習(xí)

        類風(fēng)濕關(guān)節(jié)炎(RA)是一種發(fā)病機(jī)制復(fù)雜的進(jìn)行性肌肉骨骼損傷疾病,遺傳因素、環(huán)境因素和微生物群組的交互作用誘導(dǎo)機(jī)體固有免疫的異?;罨钦T發(fā)RA的主要原因[1-3]。RA多發(fā)病于中老年人群體,以女性患者居多,近年來其全球患病率和流行率呈現(xiàn)上升趨勢(shì),且其流行趨勢(shì)出現(xiàn)年輕化[4, 5]。盡管提倡疾病的早發(fā)現(xiàn)、早診斷、早治療有助于患者預(yù)后,但由于RA的高度異質(zhì)性,部分患者發(fā)病前期的生化指標(biāo)和體征并不符合RA的診斷金標(biāo)準(zhǔn),在確診時(shí)疾病往往已經(jīng)進(jìn)展為中后期,嚴(yán)重降低了患者的生命質(zhì)量[6, 7]。因此探究疾病新型生物標(biāo)志物對(duì)于RA的臨床前期診斷有重要意義。此外,針對(duì)RA免疫系統(tǒng)異?;罨奶攸c(diǎn),目前臨床常用治療方案多以具有免疫抑制效果的藥物為主,但長(zhǎng)期服藥引起的非靶組織的毒副作用和藥物耐受限制了傳統(tǒng)抗風(fēng)濕藥物的應(yīng)用[8, 9]。因此,為了進(jìn)一步推進(jìn)RA的精準(zhǔn)診療,探究潛在生物標(biāo)志物在RA中的免疫調(diào)控機(jī)制是有必要的。

        基因表達(dá)綜合(GEO)數(shù)據(jù)庫(kù)是國(guó)內(nèi)外常用的疾病基因分析公共數(shù)據(jù)庫(kù),其涵蓋了多種疾病的不同類型的基因組數(shù)據(jù),可以用于探索病因未明疾病的表觀遺傳改變并篩選新型疾病診斷和預(yù)后生物標(biāo)志物[10]。本研究基于GEO數(shù)據(jù)庫(kù)多個(gè)數(shù)據(jù)集的基因組數(shù)據(jù),篩選得到了與RA發(fā)生發(fā)展密切相關(guān)的9個(gè)核心基因,通過ROC工作曲線預(yù)測(cè)了核心基因作為生物標(biāo)志物的信效度,并初步探究了核心基因與疾病特征免疫細(xì)胞之間的關(guān)系。隨后通過構(gòu)建膠原誘導(dǎo)性關(guān)節(jié)炎(CIA)大鼠模型驗(yàn)證了相關(guān)性較強(qiáng)的核心基因在RA滑膜組織中的表達(dá)情況。本研究通過整合分析多個(gè)數(shù)據(jù)集的基因組數(shù)據(jù),從多方面探究了RA發(fā)生發(fā)展過程中具有重要生物學(xué)意義的潛在生物標(biāo)志物的表達(dá)和診斷價(jià)值,并結(jié)合疾病特征免疫細(xì)胞進(jìn)一步探討了相關(guān)的生物標(biāo)志物在RA發(fā)病過程中的免疫調(diào)控過程,對(duì)挖掘RA精準(zhǔn)診療新靶點(diǎn)具有一定意義。

        1 材料和方法

        1.1 數(shù)據(jù)資料來源

        從GEO數(shù)據(jù)庫(kù)中篩選含有滑膜組織健康對(duì)照組和類風(fēng)濕關(guān)節(jié)炎患者組的訓(xùn)練集和驗(yàn)證集數(shù)據(jù),共納入5個(gè)訓(xùn)練集:GSE1919,GSE12021,GSE55235,GSE55457與GSE77298;一個(gè)驗(yàn)證集:GSE89408。

        1.2 實(shí)驗(yàn)動(dòng)物

        5~7周齡雄性SD大鼠25只,體質(zhì)量180±20 g,飼養(yǎng)于蚌埠醫(yī)科大學(xué)科研中心SPF級(jí)動(dòng)物房,給予充足飲食和飲水,適應(yīng)性飼養(yǎng)1周后建立CIA大鼠模型。所有動(dòng)物實(shí)驗(yàn)經(jīng)蚌埠醫(yī)科大學(xué)實(shí)驗(yàn)動(dòng)物管理和倫理委員會(huì)批準(zhǔn)(倫動(dòng)科批號(hào)[2021]第192號(hào))。

        1.3 主要試劑

        牛Ⅱ 型膠原(Chondrex)、弗氏不完全佐劑(Chondrex);甲氨蝶呤(國(guó)藥準(zhǔn)字H31020644,上海上藥信誼藥廠有限公司);4%多聚甲醛通用型組織固定液(Biosharp);鼠抗信號(hào)轉(zhuǎn)導(dǎo)和轉(zhuǎn)錄激活因子1(STAT1,Proteintech);TRIzol(Ambion);細(xì)胞核蛋白與細(xì)胞漿蛋白抽提試劑盒(碧云天)。

        1.4 差異表達(dá)基因篩選

        1.4.1 差異表達(dá)基因獲得

        對(duì)單個(gè)數(shù)據(jù)集的基因表達(dá)矩陣進(jìn)行差異分析,利用Robust Rank Aggreg(RRA)方法篩選在多個(gè)數(shù)據(jù)集均表現(xiàn)差異的基因。此外,為剔除數(shù)據(jù)集間的批次效應(yīng),將多個(gè)數(shù)據(jù)集合并后使用BatchNormalization 方法進(jìn)行批次矯正,對(duì)批次矯正后的數(shù)據(jù)集進(jìn)行差異分析。最后將兩種方法得到的差異表達(dá)基因取交集,得到共差異表達(dá)基因。以上兩種差異分析篩選標(biāo)準(zhǔn)均為:|logFC|gt;1,adj.P.Vallt;0.05。

        1.4.2 GO和KEGG功能富集分析

        對(duì)共差異表達(dá)基因進(jìn)行GO功能富集分析和KEGG通路富集分析。并對(duì)GO富集分析和KEGG富集分析結(jié)果中各條目相關(guān)性統(tǒng)計(jì)學(xué)意義排列的前5個(gè)條目進(jìn)行可視化展示。

        1.4.3 免疫細(xì)胞浸潤(rùn)分析

        為量化浸潤(rùn)性免疫細(xì)胞在RA中的相對(duì)比例,利用CIBERSORT方法分析矯正合并后的共差異表達(dá)基因在每個(gè)樣品中的表達(dá)矩陣,以Plt;0.05過濾免疫浸潤(rùn)結(jié)果。計(jì)算樣品中每種免疫細(xì)胞類型的百分比,并以柱形圖顯示。使用pheatmap 包繪制22種免疫細(xì)胞的聚類熱圖,使用corrplot軟件包繪制22種免疫細(xì)胞的相關(guān)性熱圖。

        1.4.4 Hub 基因篩選

        從STRING在線數(shù)據(jù)庫(kù)(https://string-db.org/)檢索關(guān)鍵基因編碼蛋白間可能的潛在相互作用,同時(shí)構(gòu)建蛋白質(zhì)互作網(wǎng)絡(luò)。隨后將蛋白互作網(wǎng)絡(luò)導(dǎo)入Cytoscape3.8.0軟件,利用cytohubba插件中的10種計(jì)算程序(MCC, DMNC, MNC, Degree, BottleNeck,EcCentricity, Closeness, Radiality, Betweenness, Stress)對(duì)基因的表達(dá)情況進(jìn)行打分,將每種打分分?jǐn)?shù)排名前50的基因視為有效打分基因,同時(shí)滿足10 種有效打分的基因認(rèn)為是蛋白互作網(wǎng)絡(luò)的關(guān)鍵基因。隨后從GEO數(shù)據(jù)庫(kù)中獲取含有正常組與RA患病早期組和晚期組滑膜組織的基因表達(dá)矩陣(GSE89408),用于驗(yàn)證關(guān)鍵基因在病程早期和晚期在滑膜組織中的表達(dá)情況。

        1.4.5 疾病特征免疫細(xì)胞篩選

        對(duì)每種免疫細(xì)胞的含量在正常組和RA組間循環(huán)進(jìn)行差異分析,將Plt;0.05的免疫細(xì)胞視為差異免疫細(xì)胞。同時(shí)利用glmnet包對(duì)免疫細(xì)胞浸潤(rùn)結(jié)果構(gòu)建LASSO回歸模型,利用交叉驗(yàn)證和LASSO回歸模型篩選特征免疫細(xì)胞。將上述實(shí)驗(yàn)得到的差異免疫細(xì)胞和特征免疫細(xì)胞取交集,得到疾病核心免疫細(xì)胞。

        1.4.6 核心基因與疾病核心免疫細(xì)胞相關(guān)性分析

        將核心基因的表達(dá)量與核心免疫細(xì)胞的含量進(jìn)行相關(guān)性檢驗(yàn),設(shè)置過濾條件為相關(guān)系數(shù)(R)的絕對(duì)值大于0.4,相關(guān)性檢驗(yàn)的Plt;0.001,滿足條件說明核心基因的表達(dá)與免疫細(xì)胞的含量是具有較高相關(guān)性的。然后利用corrplot包繪制相關(guān)性熱圖和單個(gè)基因與核心免疫細(xì)胞散點(diǎn)圖,若得到的散點(diǎn)圖相關(guān)系數(shù)大于0,則說明兩者間存在正調(diào)控關(guān)系,反之為負(fù)調(diào)控關(guān)系。

        1.4.7 單基因ssGSEA分析

        將關(guān)鍵基因表達(dá)分為高表達(dá)和低表達(dá)兩組,將基因按照在RA組和對(duì)照組中的差異表達(dá)程度排序,引用org.Hs.eg.db人類數(shù)據(jù)庫(kù)轉(zhuǎn)化ID和logFoldChange方法進(jìn)行差異表達(dá)分析。對(duì)關(guān)鍵基因的KEGG通路進(jìn)行GSEA富集分析,富集通路過濾條件為:GSSize≤200amp;GSSize≥10,Plt;0.05。為了得到單基因在疾病中的富集通路功能差異,從GSEA分子特征數(shù)據(jù)庫(kù)(www.gsea-msigdb.org/gsea/msigdb)中下載標(biāo)志基因集文件,計(jì)算ssGSEA功能富集評(píng)分并對(duì)得分結(jié)果進(jìn)行差異分析,繪制各功能通路在RA組與對(duì)照組間差異。

        1.5 膠原誘導(dǎo)性關(guān)節(jié)炎(CIA)大鼠模型的構(gòu)建及干預(yù)

        根據(jù)參考文獻(xiàn)和本實(shí)驗(yàn)室沿用的方法構(gòu)建CIA大鼠模型[11, 12]。將等量的弗氏不完全佐劑和牛II 型膠原(2 mg/mL)在冰上充分混勻,直至完全乳化得到造模所需的乳化劑。為保證動(dòng)物成模,采取兩次加強(qiáng)免疫方法造模。初次免疫時(shí)在SD大鼠尾根部皮內(nèi)注射200 μL乳化劑,1 周后加強(qiáng)免疫,在SD大鼠尾根部皮內(nèi)注射100 μL乳化劑。加強(qiáng)免疫5~8 d 后,SD大鼠均出現(xiàn)關(guān)節(jié)紅腫和活動(dòng)受限,通過關(guān)節(jié)炎評(píng)分(AI)、足爪紅腫狀態(tài)和足跖厚度證明成模。AIgt;4 分時(shí)認(rèn)為模型構(gòu)建成功,開始給藥。

        采用隨機(jī)數(shù)法將造模成功的大鼠分為3組,對(duì)照組(vehicle)、模型組(CIA組)和甲氨蝶呤組(MTX組),MTX組每?jī)芍芙?jīng)口灌胃給藥1次,給藥濃度為2 mg/kg,對(duì)照組和模型組大鼠經(jīng)口灌胃2 mL的PBS溶液,連續(xù)給藥4周。

        1.6 組織病理學(xué)染色與免疫熒光染色

        大鼠脫頸處死,取下肢骨組織后用4%多聚甲醛固定,慢脫鈣處理后脫水透明、浸蠟包埋切片,厚度為5 μm。使用蘇木精伊紅(HE)染色和番紅O-固綠染色觀察大鼠膝關(guān)節(jié)組織炎性細(xì)胞浸潤(rùn)、血管翳形成、滑膜襯里層變化和關(guān)節(jié)軟骨破環(huán)情況;使用免疫熒光染色觀察大鼠踝關(guān)節(jié)滑膜組織中目標(biāo)蛋白表達(dá)水平。

        1.7 滑膜組織細(xì)胞胞質(zhì)與胞核蛋白提取

        收集對(duì)照組、模型組和MTX治療組大鼠滑膜組織的滑膜成纖維細(xì)胞至1.5 mL EP管中,使用預(yù)冷的PBS洗2次,加入200 L預(yù)冷的細(xì)胞漿蛋白抽提試劑A (pH 7.910 mmol/L HEPES、10 mmol/L KCl、0.1 mmol/LEDTA、1 mmol/L dithiothreitol、0.4% Igepal CA-630、5 mol/L leupeptin、2 mol/L pepstatin A、1 mol/L aprotinin和1 mmol/L phenylmethyl sulfonyl fluoride)懸浮細(xì)胞,冰上放置15 min,4℃ 12000 r/min離心5 min,上清液為細(xì)胞質(zhì)含胞漿蛋白-70 ℃保存?zhèn)溆?,沉淀物為?xì)胞核。隨后將沉淀物中殘余的上清液吸凈,加入50 L 預(yù)冷的細(xì)胞核蛋白抽提試劑 B(pH 7.9 20 mmol/L HEPES、0.4 mol/L NaCl、1 mmol/L EDTA、1 mmol/L二硫蘇糖醇和1 mmol/L苯甲基磺酰氯),4℃渦旋30 min,12 000 r/min離心10 min,隨后立即吸取上清液至預(yù)冷的EP管中,即得到細(xì)胞核蛋白。

        1.8 Western blotting

        使用上述提取的細(xì)胞核蛋白和細(xì)胞質(zhì)蛋白進(jìn)行免疫印跡實(shí)驗(yàn)。首先使用BCA蛋白定量法測(cè)定蛋白質(zhì)濃度。調(diào)整蛋白上樣量使其符合實(shí)驗(yàn)要求。120 V電泳1.5 h,200 mA轉(zhuǎn)膜2 h。使用5%脫脂牛奶封閉2 h,隨后將條帶置于STAT1 抗體(Abcam)與p-STAT1 抗體(Abcam)中4 ℃搖床孵育過夜。使用TBST沖洗后,二抗室溫封閉2 h,隨后用超敏ECL化學(xué)發(fā)光底物在凝膠成像系統(tǒng)曝光,收集可用圖像并用EXCEL做統(tǒng)計(jì)學(xué)分析。

        1.9 統(tǒng)計(jì)學(xué)分析

        使用R軟件處理實(shí)驗(yàn)數(shù)據(jù),多組間數(shù)據(jù)使用單因素方差分析,兩組間數(shù)據(jù)采用Wilcoxon秩和檢驗(yàn),Plt;0.05認(rèn)為差異有統(tǒng)計(jì)學(xué)意義,所有實(shí)驗(yàn)均獨(dú)立重復(fù)3次。

        2 結(jié)果

        2.1 差異表達(dá)基因篩選

        利用R軟件篩選得到不同數(shù)據(jù)集的差異表達(dá)基因(圖1)。基于RRA方法獲得差異表達(dá)基因共454個(gè),其中上調(diào)基因292個(gè),下調(diào)基因162個(gè)。利用批次標(biāo)準(zhǔn)化方法得到差異基因349個(gè)。對(duì)通過兩種方法得到的差異表達(dá)基因取交集,共得到338 個(gè)共差異表達(dá)基因,以供后續(xù)分析。

        2.2 交集基因的GO和KEGG富集分析結(jié)果

        GO富集分析結(jié)果表明,共差異表達(dá)基因在BP中主要富集在白細(xì)胞的細(xì)胞間粘附、單核細(xì)胞分化、淋巴細(xì)胞分化、白細(xì)胞遷移、調(diào)節(jié)白細(xì)胞的細(xì)胞間粘附性;在CC層面主要集中于細(xì)胞質(zhì)膜外側(cè)、凝集素包裹的內(nèi)泡膜、凝集素包裹的內(nèi)吞囊泡、凝集素包裹的囊泡膜、膜筏;在MF層面主要富集于趨化因子受體結(jié)合、趨化因子活性、免疫受體活性、G蛋白偶聯(lián)受體結(jié)合、細(xì)胞因子活性(圖2A)。KEGG富集分析結(jié)果表明,關(guān)鍵基因主要富集于趨化因子信號(hào)通路、細(xì)胞因子-細(xì)胞因子受體信號(hào)通路、類風(fēng)濕關(guān)節(jié)炎信號(hào)通路、NF-κB信號(hào)通路、Th1和Th2免疫細(xì)胞分化(圖2B)。

        2.3 核心基因篩選和表達(dá)水平分析

        通過Cytoscape軟件共篩選得到9個(gè)且均在RA中存在差異表達(dá)核心基因(圖3A、B),分別為CD3G、CD8A、脾酪氨酸激酶(SYK)、淋巴細(xì)胞特異性蛋白酪氨酸激酶(LCK)、白細(xì)胞介素2受體γ(IL2RG)、信號(hào)傳導(dǎo)及轉(zhuǎn)錄激活蛋白1(STAT1)、趨化因子受體5(CCR5)、整合素亞基β2(ITGB2)、整合素亞基αL(ITGAL)。

        選擇含有正常對(duì)照組、RA早期患者和晚期患者滑膜組織樣本的數(shù)據(jù)集分析核心基因的表達(dá)情況。相較于正常組的低表達(dá),9個(gè)核心基因在RA發(fā)病早期的滑膜組織中表達(dá)均升高(Plt;0.001);與RA 早期相比,ITGB2和ITGBAL兩個(gè)核心基因的中位值降低,差異具有統(tǒng)計(jì)學(xué)意義(Plt;0.05),但仍高于正常組的表達(dá)水平,其余7個(gè)基因的表達(dá)水平在RA早晚期的滑膜組織中沒有統(tǒng)計(jì)學(xué)差異(Pgt;0.05,圖3C~K)。

        2.4 核心基因功能預(yù)測(cè)和ROC工作曲線結(jié)果

        將通過GeneMANIA數(shù)據(jù)庫(kù)預(yù)測(cè)得到的基因功能結(jié)果根據(jù)FDR值進(jìn)行排序,發(fā)現(xiàn)核心基因可能通過細(xì)胞因子或其表面受體通路調(diào)控細(xì)胞異常激活和炎癥反應(yīng)發(fā)揮作用(圖4A)。

        ROC曲線分析表明,9 個(gè)核心基因的AUC值均大于0.65,其中,ITGB2 預(yù)測(cè)模型的AUC 值為0.667;ITGAL、LCK、CD8A、CD3G、SYK、CCR5和IL2RG等7個(gè)預(yù)測(cè)模型的AUC值介于0.70~0.85(圖4B~I(xiàn));STAT1預(yù)測(cè)模型的AUC值為0.909(圖4J)。

        2.5 免疫細(xì)胞浸潤(rùn)分析與疾病特征免疫細(xì)胞篩選

        CIBERSORT分析結(jié)果顯示,多種免疫細(xì)胞在RA中存在差異表達(dá)(圖5A)。其中RA組的記憶B細(xì)胞、CD8 T細(xì)胞、活化CD4 T細(xì)胞、濾泡性輔助性T細(xì)胞、γδT細(xì)胞、M0型巨噬細(xì)胞和M1型巨噬細(xì)胞的表達(dá)數(shù)量高于正常組(Plt;0.01),而正常組的靜息CD4 T細(xì)胞、漿細(xì)胞、調(diào)節(jié)性T細(xì)胞、活化的NK細(xì)胞、單核細(xì)胞、靜息樹突狀細(xì)胞、靜息肥大細(xì)胞和嗜酸性粒細(xì)胞的表達(dá)數(shù)量高于RA組(Plt;0.01,圖5B)。此外,免疫細(xì)胞相關(guān)性熱圖結(jié)果表明,多對(duì)免疫細(xì)胞在RA的發(fā)生發(fā)展中具有協(xié)同或拮抗作用(圖5C)。LASSO回歸分析共得到16個(gè)與RA發(fā)生發(fā)展密切相關(guān)的免疫細(xì)胞。與上文中篩選的差異免疫細(xì)胞取交集后共得到11 個(gè)疾病特征免疫細(xì)胞(圖5D)。分別為:漿細(xì)胞、靜息CD4 T細(xì)胞、濾泡性輔助性T細(xì)胞、γδ T細(xì)胞、單核細(xì)胞、M0型巨噬細(xì)胞、M1型巨噬細(xì)胞、靜息樹突狀細(xì)胞、靜息肥大細(xì)胞、活化的肥大細(xì)胞和嗜酸性粒細(xì)胞。

        2.6 核心基因與疾病特征免疫細(xì)胞的相關(guān)性分析

        9 個(gè)核心基因均與疾病特征免疫細(xì)胞存在不同程度的相關(guān)性(圖6A)。其中CD3G、ITGAL、LCK、CD8A和STAT1 等5 個(gè)核心基因與相關(guān)的疾病特征免疫細(xì)胞相關(guān)性較強(qiáng)(圖6B~F),CD3G(R=0.54,Plt;0.001)、CD8A(R=0.65,Plt;0.001)與LCK(R=0.59,Plt;0.001)和γδ T細(xì)胞存在正相關(guān),ITGAL與Tfh細(xì)胞存在正相關(guān)(R=0.56,Plt;0.001),STAT1 與M1 型巨噬細(xì)胞存在正相關(guān),且其相關(guān)性系數(shù)最高(R=0.68,Plt;0.001)。

        2.7 STAT1在各組大鼠中的表達(dá)情況

        GSEA富集分析結(jié)果展示了STAT1表達(dá)升高情況下可能介導(dǎo)的5條生物學(xué)過程,分別為同種異體移植排斥,移植物抗宿主病,產(chǎn)生IgA的腸道免疫網(wǎng)絡(luò),原發(fā)性免疫缺陷和病毒蛋白與細(xì)胞因子及細(xì)胞因子受體的相互作用(圖7A)。

        膝關(guān)節(jié)HE結(jié)果表明,CIA組大鼠與正常組大鼠相比滑膜組織增生明顯,骨質(zhì)侵蝕程度較重,可見大量炎性細(xì)胞浸潤(rùn)和血管翳產(chǎn)生,經(jīng)MTX治療后,鏡下炎細(xì)胞浸潤(rùn)程度明顯降低且關(guān)節(jié)表面光滑無糜爛(圖7B);膝關(guān)節(jié)番紅-O固綠結(jié)果顯示,CIA組大鼠軟骨面和軟骨下骨組織界限不清晰,軟骨細(xì)胞減少,而經(jīng)MTX處理的大鼠的軟骨表面光滑且軟骨細(xì)胞數(shù)量明顯增加,軟骨潮線清晰、骨組織無明顯破壞(圖7C)。

        免疫熒光結(jié)果表明,CIA組大鼠的踝關(guān)節(jié)滑膜組織中STAT1表達(dá)水平明顯高于正常組,而經(jīng)過MTX治療后,STAT1的表達(dá)水平明顯降低(圖7D)。

        2.8 STAT1磷酸化和細(xì)胞核遷移

        CIA組大鼠踝關(guān)節(jié)滑膜組織中p-STAT1的表達(dá)水平顯著增高,而MTX 組大鼠踝關(guān)節(jié)滑膜組織中p-STAT1的表達(dá)水平下調(diào)(圖8A)。Western blotting實(shí)驗(yàn)結(jié)果表明CIA組大鼠滑膜組織細(xì)胞核內(nèi)的STAT1與p-STAT1的蛋白表達(dá)量高于細(xì)胞質(zhì),且顯著高于對(duì)照組大鼠細(xì)胞質(zhì)與細(xì)胞核中STAT1(Cytosolo:P=0.0002;Nuclear:Plt;0.0001)與p-STAT1(Cytosolo:P=0.0023;Nuclear:Plt;0.0001)蛋白的表達(dá)水平,而經(jīng)MTX治療后發(fā)現(xiàn),細(xì)胞核內(nèi)p-STAT1 的蛋白表達(dá)量明顯減少(Plt;0.0001,圖8B~D)。

        3 討論

        借助計(jì)算機(jī)生物信息學(xué)的方法可以探究不同基因間的互作過程和其調(diào)控疾病的可能發(fā)生機(jī)制[13]。RA作為一種高度異質(zhì)性疾病,具體發(fā)病機(jī)制尚不明晰[1, 14]。臨床常通過患者遠(yuǎn)端小關(guān)節(jié)的病理變化和血清類風(fēng)濕因子、C-反應(yīng)蛋白和紅細(xì)胞沉降率等指標(biāo)判斷是否發(fā)病,但傳統(tǒng)生化指標(biāo)不適用于部分血清學(xué)檢驗(yàn)陰性的RA患者[15, 16]。因此本研究依托生信分析探究RA的潛在生物標(biāo)志物,初步探究有益于RA臨床早期診斷的潛在標(biāo)志物。

        本研究通過整合GEO數(shù)據(jù)庫(kù)中5 個(gè)RA患者滑膜組織數(shù)據(jù)集的基因組數(shù)據(jù),初步篩選得到338個(gè)與RA發(fā)病密切相關(guān)的共差異表達(dá)基因。KEGG和GO富集分析表明差異基因可能通過激活細(xì)胞內(nèi)外質(zhì)膜表面受體發(fā)揮相關(guān)功能,主要調(diào)控炎性細(xì)胞因子及其受體通路以及免疫細(xì)胞的增殖分化。隨后進(jìn)一步篩選與RA發(fā)生發(fā)展密切相關(guān)的核心基因,發(fā)現(xiàn)CD3G、CD8A、SYK、LCK、IL2RG、STAT1、CCR5、ITGB2 和ITGAL等9 個(gè)核心基因在疾病差異基因拓?fù)渚W(wǎng)絡(luò)作為樞紐節(jié)點(diǎn)與其他基因關(guān)聯(lián)緊密且節(jié)點(diǎn)較多。此外,我們發(fā)現(xiàn)9個(gè)核心基因在RA患者的滑膜組織中表達(dá)水平均升高,且CD3G、CD8A、SYK、LCK、IL2RG、STAT1 和CCR5 這7 個(gè)核心基因在RA早期和晚期的滑膜組織中表達(dá)水平差異不具有統(tǒng)計(jì)學(xué)意義(Pgt;0.05),這提示以上7 個(gè)核心基因在RA早期可能作為疾病的診斷標(biāo)志物。ROC曲線結(jié)果結(jié)果提示9個(gè)核心基因用于RA早期診斷具有一定的意義,其中STAT1 作為RA預(yù)測(cè)模型的AUC值最高(0.909),提示其作為RA發(fā)病的生物診斷標(biāo)志物具有較高信效度。

        RA作為一種自身免疫性疾病,免疫細(xì)胞的異?;罨⑶忠u與互作在疾病發(fā)展中發(fā)揮重要作用[17, 18]。在既往的一些研究中,CD3G、CD8A、ITGAL、LCK和STAT1已被證明與RA的發(fā)生發(fā)展關(guān)系密切[19-21],但大多研究?jī)H通過觀察其蛋白質(zhì)或mRNA的表達(dá)水平說明這些基因在RA的病程進(jìn)展中有一定的價(jià)值,而并未與疾病相關(guān)免疫細(xì)胞相聯(lián)系進(jìn)而說明其在RA中的可能調(diào)控機(jī)制。本研究篩選得到16 種疾病特征免疫細(xì)胞,初步探究了差異表達(dá)核心基因與疾病特征免疫細(xì)胞的潛在相關(guān)性。結(jié)果表明CD3G、CD8A、ITGAL、LCK和STAT1這5個(gè)疾病差異表達(dá)核心基因主要與不同亞型的T淋巴細(xì)胞和M1型巨噬細(xì)胞存在密切關(guān)聯(lián)。

        γδ T細(xì)胞是具有γ鏈和δ鏈組成的異二聚體T細(xì)胞受體復(fù)合物(TCR)的一類T淋巴細(xì)胞亞群[22, 23]。研究表明,γδ T細(xì)胞在RA患者滑液中高表達(dá),與滑膜炎癥的嚴(yán)重程度和疾病活動(dòng)度呈正相關(guān)[24, 25]。本研究表明,CD3G、CD8A 與LCK 基因與γδ T 細(xì)胞密切相關(guān),CD3G、CD8A與LCK表達(dá)水平的增加可能與γδ T細(xì)胞的異常生物學(xué)行為有關(guān)。盡管在過去的研究中這3 個(gè)基因作為RA的生物標(biāo)志均被提及[20, 21, 26],但尚未有研究明確表明其與γδ T細(xì)胞之間的潛在互作關(guān)系,因此在后續(xù)的研究中將CD3G、CD8A與LCK作為調(diào)控γδ T細(xì)胞的潛在靶點(diǎn)在RA的靶向治療和明晰其發(fā)病機(jī)制方面可能具有一定意義。Tfh細(xì)胞主要通過分泌炎性介質(zhì)或介導(dǎo)B 淋巴細(xì)胞的增殖分化促進(jìn)RA 的疾病進(jìn)展[27]。本研究發(fā)現(xiàn)ITGAL與Tfh細(xì)胞的水平呈現(xiàn)正相關(guān),提示ITGAL的表達(dá)水平與Tfh細(xì)胞的增殖分化可能存在一定關(guān)系。研究表明ITGAL(CD11A)與自身免疫性疾病中CD4+T細(xì)胞的過度活化密切相關(guān)[26, 28],但關(guān)于ITGAL調(diào)控Tfh 細(xì)胞在RA中的具體作用機(jī)制仍需進(jìn)一步探究。M1型促炎巨噬細(xì)胞可以通過分泌多種炎性因子或激活滑膜細(xì)胞異常生物學(xué)行為加重滑膜炎癥[29, 30]。研究表明STAT1的表達(dá)增加能夠促進(jìn)巨噬細(xì)胞向M1 型促炎表型極化[31],本研究發(fā)現(xiàn)STAT1 與M1 型巨噬細(xì)胞的水平存在正相關(guān),同時(shí)STAT1可能參與Th1與Th2免疫細(xì)胞的增殖分化。GSEA富集分析結(jié)果提示STAT1可能通過腸道免疫網(wǎng)絡(luò)和細(xì)胞因子及其受體的激活這兩種生物學(xué)過程發(fā)揮免疫調(diào)控作用,而這兩種信號(hào)通路已被證明在RA的發(fā)生發(fā)展中具有重要作用[32-34]。

        為進(jìn)一步驗(yàn)證生信分析結(jié)果的可信性,本研究選擇AUC值最大且與疾病特征免疫細(xì)胞聯(lián)系最為緊密的基因(STAT1)進(jìn)行動(dòng)物實(shí)驗(yàn)。本研究結(jié)果表明,在CIA大鼠中,踝關(guān)節(jié)滑膜組織的STAT1 與p-STAT1 的表達(dá)明顯高于對(duì)照組;p-STAT1在CIA大鼠踝關(guān)節(jié)滑膜成纖維細(xì)胞細(xì)胞核中的表達(dá)量高于細(xì)胞質(zhì),而經(jīng)過MTX治療后,細(xì)胞核內(nèi)p-STAT1 的表達(dá)量顯著降低,這提示STAT1 在RA中通過磷酸化以及核遷移介導(dǎo)信號(hào)通路的激活,發(fā)揮轉(zhuǎn)錄因子的作用。以上結(jié)果說明STAT1作為RA的潛在生物標(biāo)志物具有重要意義,其可能通過介導(dǎo)M1 型巨噬細(xì)胞炎性活化進(jìn)而調(diào)節(jié)細(xì)胞因子及其受體通路或腸道免疫網(wǎng)絡(luò)途徑在RA的病程進(jìn)展中發(fā)揮重要作用,但具體調(diào)控機(jī)制還需進(jìn)一步探討。

        JAK-STAT信號(hào)通路的激活與RA微環(huán)境中炎癥介質(zhì)的釋放、破骨細(xì)胞的分化以及骨質(zhì)的破壞密切相關(guān),針對(duì)這一通路所研發(fā)并應(yīng)用于臨床的靶點(diǎn)阻滯劑也較成熟[35, 36]。目前臨床應(yīng)用的阻滯劑主要以靶向抑制JAK1/2/3為主,目前報(bào)道的靶向STAT3的一些藥物,如Stattic、OPB-111077和HL-237等也已處于臨床試驗(yàn)階段,但是關(guān)于靶向STAT1的藥物研究較少,針對(duì)于臨床STAT1特異性靶向藥也被認(rèn)為有較大的探索空間[37]。研究表明,STAT1能夠驅(qū)動(dòng)趨化因子誘導(dǎo)白細(xì)胞募集到關(guān)節(jié)炎,加重疾病活動(dòng)度[38];抑制STAT1能夠促進(jìn)Th1細(xì)胞分化,緩解關(guān)節(jié)炎癥狀[39]。因此,結(jié)合本研究結(jié)果和文獻(xiàn)報(bào)道,我們推測(cè)STAT1可以作為RA臨床潛在治療靶點(diǎn)。

        綜上所述,本研究結(jié)果表明CD3G、CD8A、SYK、LCK、IL2RG、STAT1、CCR5、ITGB2與ITGAL這9個(gè)基因可能作為RA的潛在生物標(biāo)志物用于疾病早期診斷,且CD3G/CD8A/LCK?γδ T 細(xì)胞、ITGAL-Tfh 細(xì)胞、STAT1-M1 型巨噬細(xì)胞等基因-免疫細(xì)胞途徑在RA的發(fā)生發(fā)展中具有重要作用。本研究不足之處在于本研究缺乏對(duì)調(diào)控機(jī)制的深入探究。在未來的實(shí)驗(yàn)中,本團(tuán)隊(duì)將依托動(dòng)物實(shí)驗(yàn)、細(xì)胞實(shí)驗(yàn)以及大型公共數(shù)據(jù)庫(kù)進(jìn)一步探究本研究中獲得的多個(gè)具有重要意義的潛在生物標(biāo)志物在RA中的具體調(diào)控機(jī)制及其與免疫細(xì)胞相互作用的生物學(xué)過程。

        參考文獻(xiàn):

        [1] Alivernini S, Firestein GS, McInnes IB. The pathogenesis ofrheumatoid arthritis[J]. Immunity, 2022, 55(12): 2255-70.

        [2] Xu Q, Ni JJ, Han BX, et al. Causal relationship between gutmicrobiota and autoimmune diseases: a two-sample Mendelianrandomization study[J]. Front Immunol, 2021, 12: 746998.

        [3] Kronzer VL, Davis JM. Etiologies of rheumatoid arthritis: update onmucosal, genetic, and cellular pathogenesis[J]. Curr Rheumatol Rep, 2021, 23(4): 21.

        [4] Chiang PH, Ju PC, Chiang YC, et al. Correspondence on‘ Incidencetrend of five common musculoskeletal disorders from 1990 to 2017at the global, regional and national level: results from the globalburden of disease study 2017'[J]. Ann Rheum Dis, 2023, 82(2): e46.

        [5] Guan SY, Zheng JX, Sam NB, et al. Global burden and risk factors ofmusculoskeletal disorders among adolescents and young adults in204 countries and territories, 1990-2019[J]. Autoimmun Rev, 2023,22(8): 103361.

        [6] Lenti MV, Rossi CM, Melazzini F, et al. Seronegative autoimmunediseases: a challenging diagnosis[J]. Autoimmun Rev, 2022, 21(9):103143.

        [7] Tidblad L, Westerlind H, Delcoigne B, et al. Comorbidities andtreatment patterns in early rheumatoid arthritis: a nationwideSwedish study[J]. RMD Open, 2022, 8(2): e002700.

        [8] Ben Mrid R, Bouchmaa N, Ainani H, et al. Anti-rheumatoid drugsadvancements: new insights into the molecular treatment ofrheumatoid arthritis[J]. Biomed Pharmacother, 2022, 151: 113126.

        [9] Tan Y, Buch MH. 'Difficult to treat' rheumatoid arthritis: currentposition and considerations for next steps[J]. RMD Open, 2022, 8(2): e002387.

        [10]Clough E, Barrett T. The gene expression omnibus database[J].Methods Mol Biol, 2016, 1418: 93-110.

        [11] Zhao T, Xie Z, Xi Y, et al. How to model rheumatoid arthritis inanimals: from rodents to non-human Primates[J]. Front Immunol,2022, 13: 887460.

        [12]Zhang X, Zhang X, Wang X, et al. Efficient delivery of triptolideplus a miR-30-5p inhibitor through the use of near infrared laserresponsive or CADY modified MSNs for efficacy in rheumatoidarthritis therapeutics[J]. Front Bioeng Biotechnol, 2020, 8: 170.

        [13]Rincón-Riveros A, Morales D, Rodríguez JA, et al. Bioinformatictools for the analysis and prediction of ncRNA interactions[J]. Int JMol Sci, 2021, 22(21): 11397.

        [14] Jang S, Kwon EJ, Lee JJ. Rheumatoid arthritis: pathogenic roles ofdiverse immune cells[J]. Int J Mol Sci, 2022, 23(2): 905.

        [15]Yu R, Zhang J, Zhuo Y, et al. Identification of diagnostic signaturesand immune cell infiltration characteristics in rheumatoid arthritisby integrating bioinformatic analysis and machine-learningstrategies[J]. Front Immunol, 2021, 12: 724934.

        [16]Cheng Q, Chen X, Wu H, et al. Three hematologic/immune systemspecificexpressed genes are considered as the potential biomarkersfor the diagnosis of early rheumatoid arthritis throughbioinformatics analysis[J]. J Transl Med, 2021, 19(1): 18.

        [17]Firestein GS, McInnes IB. Immunopathogenesis of rheumatoidarthritis[J]. Immunity, 2017, 46(2): 183-96.

        [18]Wu D, Luo Y, Li T, et al. Systemic complications of rheumatoidarthritis: focus on pathogenesis and treatment[J]. Front Immunol,2022, 13: 1051082.

        [19]Matsumoto H, Fujita Y, Onizawa M, et al. Increased CEACAM1expression on peripheral blood neutrophils in patients withrheumatoid arthritis[J]. Front Immunol, 2022, 13: 978435.

        [20]Li Z, Xu M, Li R, et al. Identification of biomarkers associated withsynovitis in rheumatoid arthritis by bioinformatics analyses[J].Biosci Rep, 2020, 40(9): BSR20201713.

        [21]Zhou S, Lu H, Xiong M. Identifying immune cell infiltration andeffective diagnostic biomarkers in rheumatoid arthritis bybioinformatics analysis[J]. Front Immunol, 2021, 12: 726747.

        [22]Vantourout P, Hayday A. Six-of-the-best: unique contributions of γδT cells to immunology[J]. Nat Rev Immunol, 2013, 13: 88-100.

        [23]Hayday AC. γδ T cell update: adaptate orchestrators of immunesurveillance[J]. J Immunol, 2019, 203(2): 311-20.

        [24]Mo WX, Yin SS, Chen H, et al. Chemotaxis of Vδ2 T cells to thejoints contributes to the pathogenesis of rheumatoid arthritis[J].Ann Rheum Dis, 2017, 76(12): 2075-84.

        [25] Jacobs MR, Haynes BF. Increase in TCR gamma delta Tlymphocytes in synovia from rheumatoid arthritis patients withactive synovitis[J]. J Clin Immunol, 1992, 12(2): 130-8.

        [26]Pascual-García S, Martínez-Peinado P, López-Jaén AB, et al.Analysis of novel immunological biomarkers related to rheumatoidarthritis disease severity[J]. Int J Mol Sci, 2023, 24(15): 12351.

        [27]Lu J, Wu J, Xia XL, et al. Follicular helper T cells: potentialtherapeutic targets in rheumatoid arthritis[J]. Cell Mol Life Sci,2021, 78(12): 5095-106.

        [28]Zhao S, Wang Y, Liang Y, et al. MicroRNA-126 regulates DNAmethylation in CD4+ T cells and contributes to systemic lupuserythematosus by targeting DNA methyltransferase 1[J]. ArthritisRheum, 2011, 63(5): 1376-86.

        [29]Cutolo M, Campitiello R, Gotelli E, et al. The role of M1/M2macrophage polarization in rheumatoid arthritis synovitis[J]. FrontImmunol, 2022, 13: 867260.

        [30]Tardito S, Martinelli G, Soldano S, et al. Macrophage M1/M2polarization and rheumatoid arthritis: a systematic review[J].Autoimmun Rev, 2019, 18(11): 102397.

        [31]Ye Q, Luo F, Yan T. Transcription factor KLF4 regulated STAT1 topromote M1 polarization of macrophages in rheumatoid arthritis[J].Aging: Albany NY, 2022, 14(14): 5669-80.

        [32]Guo DG, Lv JH, Chen X, et al. Study of miRNA interactome inactive rheumatoid arthritis patients reveals key pathogenic roles ofdysbiosis in the infection-immune network[J]. Rheumatology:Oxford, 2021, 60(3): 1512-22.

        [33]Zaiss MM, Wu HJ J, Mauro D, et al. The gut-joint axis in rheumatoidarthritis[J]. Nat Rev Rheumatol, 2021, 17: 224-37.

        [34]Kondo N, Kuroda T, Kobayashi D. Cytokine networks in thepathogenesis of rheumatoid arthritis[J]. Int J Mol Sci, 2021, 22(20):10922.

        [35]Simon LS, Taylor PC, Choy EH, et al. The Jak/STAT pathway: afocus on pain in rheumatoid arthritis[J]. Semin Arthritis Rheum,2021, 51(1): 278-84.

        [36]Banerjee S, Biehl A, Gadina M, et al. JAK-STAT signaling as atarget for inflammatory and autoimmune diseases: current and futureprospects[J]. Drugs, 2017, 77(5): 521-46.

        [37]Hu L, Liu RJ, Zhang LL. Advance in bone destruction participatedby JAK/STAT in rheumatoid arthritis and therapeutic effect of JAK/STAT inhibitors[J]. Int Immunopharmacol, 2022, 111: 109095.

        [38]Karonitsch T, Saferding V, Kieler M, et al. Amino acids fuelingfibroblast-like synoviocyte activation and arthritis by regulatingchemokine expression and leukocyte migration[J]. ArthritisRheumatol, 2024, 76(4): 531-40.

        [39]Mai YP, Yu XT, Gao T, et al. Autoantigenic peptide andimmunomodulator codelivery system for rheumatoid arthritistreatment by reestablishing immune tolerance[J]. ACS Appl MaterInterfaces, 2024: Online ahead of print.

        (編輯:余詩(shī)詩(shī))

        基金項(xiàng)目:國(guó)家級(jí)大學(xué)生創(chuàng)新訓(xùn)練項(xiàng)目(202310367024,202310367027)

        猜你喜歡
        生物標(biāo)志物類風(fēng)濕關(guān)節(jié)炎機(jī)器學(xué)習(xí)
        基于UPLC—Q—TOF—MS技術(shù)的牛血清白蛋白誘導(dǎo)過敏反應(yīng)的代謝組學(xué)研究
        基于UPLC—Q—TOF—MS技術(shù)的牛血清白蛋白誘導(dǎo)過敏反應(yīng)的代謝組學(xué)研究
        類風(fēng)濕關(guān)節(jié)炎蒙醫(yī)藥治療現(xiàn)狀與展望
        類風(fēng)濕關(guān)節(jié)炎的診斷及其藥物治療研究進(jìn)展
        中西醫(yī)結(jié)合治療類風(fēng)濕關(guān)節(jié)炎的臨床療效觀察
        基于網(wǎng)絡(luò)搜索數(shù)據(jù)的平遙旅游客流量預(yù)測(cè)分析
        康復(fù)護(hù)理在類風(fēng)濕關(guān)節(jié)炎治療中的應(yīng)用分析
        今日健康(2016年12期)2016-11-17 14:30:20
        前綴字母為特征在維吾爾語(yǔ)文本情感分類中的研究
        基于支持向量機(jī)的金融數(shù)據(jù)分析研究
        海洋環(huán)境監(jiān)測(cè)中生物標(biāo)志物的研究進(jìn)展
        日本免费一区二区三区影院| 中文字幕乱偷乱码亚洲| 黑丝美女被内射在线观看| 国产精品视频白浆免费视频| 色一情一乱一伦一视频免费看| 久久人人爽人人爽人人片av麻烦 | 一本久久精品久久综合桃色| 亚洲av在线观看播放| 亚洲精品无码专区| 亚洲综合一区无码精品| 亚洲av中文无码乱人伦在线咪咕| 亚州av高清不卡一区二区| 亚洲一区二区三区影院| 国产亚洲视频在线观看网址| av无码特黄一级| 日本九州不卡久久精品一区| 亚洲国产成人久久三区| 色拍拍在线精品视频| 亚洲精品白浆高清久久| 亚洲中文字幕久久精品色老板 | 中文字幕中文字幕在线中二区 | 国产91色综合久久免费| 国产好大好硬好爽免费不卡| 91麻豆精品激情在线观看最新| 日韩av在线不卡一二三区| 国产婷婷色一区二区三区深爱网| 97久久人人超碰超碰窝窝| 美女一级毛片免费观看97| 亚洲av毛片在线播放| 一区二区三区亚洲视频| 娜娜麻豆国产电影| 97se在线| 一本色道久久88综合| 久久综合伊人77777麻豆| 中国a级毛片免费观看| 激情 一区二区| 老女人下面毛茸茸的视频| 国产aⅴ无码专区亚洲av麻豆 | 精品 无码 国产观看| 亚洲发给我的在线视频| 亚洲熟妇色自偷自拍另类|