趙銘佳，郭劍，耿女，蘇興，李翠，姚雨宏，晏紅偉，韓寶生，劉美玲，盧文紅，谷翊群*

(1.唐山市婦幼保健院生殖遺傳科，唐山 063000；2.北京協(xié)和醫(yī)學(xué)院研究生院，北京 100730；3.國(guó)家衛(wèi)生健康委科學(xué)技術(shù)研究所 男性生殖健康重點(diǎn)實(shí)驗(yàn)室，北京 100081)

全球范圍內(nèi)大約8%～12%的夫婦受到不孕癥困擾，其中男性因素大約占50%[1]。男性不育以無(wú)精子癥最為嚴(yán)重，分為梗阻性無(wú)精子癥(obstructive azoospermia，OA)和非梗阻性無(wú)精子癥(non-obstructive azoospermia，NOA)[2]。NOA病因復(fù)雜，大約73%的NOA病因不明[3]。精子發(fā)生是一個(gè)極其復(fù)雜的細(xì)胞分化過(guò)程，涉及調(diào)節(jié)生殖細(xì)胞發(fā)育和成熟的基因就多達(dá)2 300個(gè)，其中微小RNA(microRNA，miRNA)在轉(zhuǎn)錄后調(diào)節(jié)基因表達(dá)[4]。在以往的研究中，雖然已經(jīng)提出了許多與精子發(fā)生相關(guān)的候選基因，但許多候選基因尚未得到驗(yàn)證或者在人群中非常罕見(jiàn)。研究表明，miRNA可能通過(guò)調(diào)節(jié)靶向基因mRNA的表達(dá)介導(dǎo)生精細(xì)胞的發(fā)育，參與精子發(fā)生，與男性不育有關(guān)[5]。本研究采用生物信息學(xué)的方法進(jìn)行數(shù)據(jù)挖掘和分析，探討非梗阻性無(wú)精子癥相關(guān)的miRNAs和關(guān)鍵基因。

材料和方法

一、數(shù)據(jù)庫(kù)檢索

在NCBI的GEO數(shù)據(jù)庫(kù)(網(wǎng)址：https：//www.ncbi.nlm.nih.gov/geo/)進(jìn)行基因芯片數(shù)據(jù)檢索，確定檢索關(guān)鍵詞為“非梗阻性無(wú)精子癥(non obstructive azoospermia)”和“智人(Homo sapiens)”。在PubMed、Embase和Web of Science中進(jìn)行檢索，確定檢索詞為“非梗阻性無(wú)精子癥(non obstructive azoospermia)”和“miRNA”或“microRNA”或“ncRNA”或“micro ribonucleic acid”。篩選并下載GSE145467和GSE108886表達(dá)譜數(shù)據(jù)集。GSE145467數(shù)據(jù)集，由Hodzic等[6]提供，基于A(yíng)gilent-014850全人類(lèi)基因組微陣列4x44K G4112F(特征編號(hào)版本)的GPL4133平臺(tái)，包括10個(gè)NOA樣本和10個(gè)OA睪丸樣本。GSE108886數(shù)據(jù)集，由Baksi等提供，基于Illumina HumanHT-12 V4.0表達(dá)微芯片的GPL10558平臺(tái)，包含 8個(gè)NOA樣本和4個(gè)OA樣本(包括1個(gè)睪丸對(duì)照樣本)。

本研究中所用數(shù)據(jù)庫(kù)及檢索篩選流程見(jiàn)圖1。

圖1 基因數(shù)據(jù)庫(kù)篩選及分析流程圖

二、篩選差異表達(dá)miRNAs：在GEO數(shù)據(jù)庫(kù)搜索已發(fā)表的NOA的miRNA數(shù)據(jù)集合，但是無(wú)可用的數(shù)據(jù)。然后，在PubMed、Embase和Web of Science中進(jìn)行檢索，以文獻(xiàn)中報(bào)道非梗阻性無(wú)精子癥睪丸組織中存在差異表達(dá)的miRNA并且經(jīng)過(guò)PCR驗(yàn)證為文獻(xiàn)篩選標(biāo)準(zhǔn)，共篩選出五項(xiàng)包含miRNA的研究(表1)[7-11]。該五項(xiàng)研究均確定了miRNA的變化趨勢(shì)，再以至少兩篇文獻(xiàn)報(bào)道且經(jīng)PCR驗(yàn)證為差異表達(dá)的miRNAs作為研究對(duì)象的miRNA篩選標(biāo)準(zhǔn)，篩選出差異表達(dá)miRNA。為研究關(guān)鍵miRNA在非梗阻性無(wú)精子癥中的作用，選擇具有一致表達(dá)變化趨勢(shì)的miRNAs作為本研究的候選差異表達(dá)miRNAs用于下面的分析。

三、差異表達(dá)miRNAs的靶基因預(yù)測(cè)：將上述篩選的miRNAs匯總，通過(guò)聯(lián)川生物云平臺(tái)的“miRNA靶基因預(yù)測(cè)云分析”(https：//www.omicstudio.cn/analysis/targetGene)預(yù)測(cè)其靶基因。

四、篩選差異表達(dá)基因：利用GEO2R工具(https：//www-ncbi-nlm-nih-gov-s.webvpn.cams.cn/geo/geo2r/)在GEO數(shù)據(jù)庫(kù)中GSE145467和GSE108886表達(dá)譜數(shù)據(jù)集的NOA和匹配的OA中篩選差異表達(dá)基因。GSE145467篩選標(biāo)準(zhǔn)是FDR小于0.05和log2FC的絕對(duì)值大于2；GSE108886的篩選標(biāo)準(zhǔn)是FDR小于0.05和log2FC的絕對(duì)值大于1。對(duì)上述兩個(gè)數(shù)據(jù)集進(jìn)行Venn Diagram分析以確定差異表達(dá)基因的交集，再應(yīng)用聯(lián)川生物云平臺(tái)(https：//www.omicstudio.cn/tool/6)與miRNAs預(yù)測(cè)靶基因取交集，獲得最終鑒定的差異表達(dá)基因并繪制韋恩圖。

五、差異表達(dá)基因的GO和KEGG富集分析：利用注釋可視化集成發(fā)現(xiàn)數(shù)據(jù)庫(kù)(The Database for Annotation，Visualization and Integrated Discovery DAVID)(https：//david.ncifcrf.gov/)進(jìn)行GO功能和KEGG通路富集分析以確定常見(jiàn)差異表達(dá)基因的功能。選擇物種為智人(Homo sapiens)。P<0.05表示差異有統(tǒng)計(jì)學(xué)意義。

六、差異表達(dá)基因的蛋白質(zhì)相互作用網(wǎng)絡(luò)和樞紐基因鑒定：利用基因相互作用檢索工具(the Search Tool for the Retrieval of Interacting Genes STRING)數(shù)據(jù)庫(kù)(https：//cn.string-db.org/)構(gòu)建差異表達(dá)基因的蛋白質(zhì)相互作用(protein-protein interaction PPI)網(wǎng)絡(luò)，同時(shí)要求交互分?jǐn)?shù)不低于 0.7。然后，應(yīng)用 Cytoscape 軟件v3.7.1(https：//cytoscape.org/)使源自STRING數(shù)據(jù)庫(kù)的PPI網(wǎng)絡(luò)可視化。應(yīng)用 Cytoscape 中的 cytoHubba 插件對(duì)PPI 網(wǎng)絡(luò)中的節(jié)點(diǎn)根據(jù)度數(shù)計(jì)算方法進(jìn)行排序，前20個(gè)基因被認(rèn)為是樞紐基因。P<0.05被認(rèn)為具有統(tǒng)計(jì)學(xué)意義的差異。

結(jié) 果

一、NOA相關(guān)的差異表達(dá)miRNAs

在PubMed、Embase和Web of Science中共篩選到5項(xiàng)研究中差異表達(dá)并經(jīng)過(guò)驗(yàn)證的miRNAs共56個(gè)(表1)。經(jīng)過(guò)對(duì)文獻(xiàn)結(jié)果中差異表達(dá)的miRNAs繪制韋恩圖(圖1)，最終篩選出5個(gè)關(guān)鍵的候選miRNAs，其中上調(diào)表達(dá)1條為miR-10b-5p，下調(diào)表達(dá)4條分別為miR-34b-5p，miR-34b-3p，miR-34c-5p，miR-449a(表2)。

表1 篩選出的5項(xiàng)研究基本資料

圖2 五項(xiàng)研究差異表達(dá)的56個(gè)miRNAs的韋恩圖

表2 篩選出的5個(gè)候選miRNAs情況

二、NOA相關(guān)的差異表達(dá)基因篩選結(jié)果

通過(guò)聯(lián)川生物云平臺(tái)的“miRNA靶基因預(yù)測(cè)云分析”共預(yù)測(cè)靶基因6 695個(gè)。在GSE108886數(shù)據(jù)集中篩選出了2 036個(gè)差異表達(dá)基因，包括1 674個(gè)上調(diào)基因和362個(gè)下調(diào)基因。在GSE145467數(shù)據(jù)集中篩選出1 274個(gè)差異表達(dá)基因，包括1，182個(gè)上調(diào)基因和92個(gè)下調(diào)基因。對(duì)上述兩個(gè)數(shù)據(jù)集進(jìn)行Venn Diagram分析取交集，共鑒定出797個(gè)差異，包括776個(gè)上調(diào)基因和21個(gè)常見(jiàn)的下調(diào)基因(圖3)。再將上述數(shù)據(jù)集的DEGs與預(yù)測(cè)獲得的miRNAs靶基因取交集，最終鑒定DEGs共180個(gè)，其中上調(diào)基因173個(gè)，下調(diào)基因7個(gè)(圖4)。篩選流程詳見(jiàn)圖1。

A：上調(diào)基因；B：下調(diào)基因。

A：上調(diào)基因；B：下調(diào)基因。

三、差異表達(dá)基因的GO和KEGG分析

差異基因GO富集于生物過(guò)程(BP)，主要包括纖毛組裝，精子軸絲組裝，鞭毛相關(guān)的精子活力，纖毛依賴(lài)性細(xì)胞運(yùn)動(dòng)，頂體組裝，精子發(fā)生，精子細(xì)胞發(fā)育，血小板脫顆粒，細(xì)胞蛋白質(zhì)代謝過(guò)程，蛋白質(zhì)穩(wěn)定化等；細(xì)胞組成(CC)主要包括活動(dòng)纖毛，軸絲，精子鞭毛，中心粒，細(xì)胞外空間等；分子功能(MF)主要包括蛋白質(zhì)結(jié)合，纖毛蛋白輕鏈的結(jié)合，整合素結(jié)合等(圖5)。

藍(lán)色代表生物過(guò)程(BP)；桔紅色代表細(xì)胞組分(CC)；綠色代表分子功能(MF)。

差異表達(dá)基因富集的KEGG通路主要包括卵巢類(lèi)固醇激素生成，多個(gè)物種長(zhǎng)壽調(diào)控途徑，擴(kuò)張型心肌病，內(nèi)分泌抵抗，孕酮介導(dǎo)的卵母細(xì)胞成熟等(圖6)。

四、差異表達(dá)基因的PPI網(wǎng)絡(luò)和樞紐基因

通過(guò)STRING數(shù)據(jù)庫(kù)構(gòu)建出差異表達(dá)基因的PPI網(wǎng)絡(luò)(圖6)。Cytoscape 中的 cytoHubba 插件篩選出了排名前20的樞紐hub基因，分別為T(mén)EKT3、EFHC1、DYNLL2、DNAH2、CETN1、SPATA7、ASRGL1、CCDC146、PLCZ1、SPAG16、DNAL1、EFCAB11、SPA4L、LIN7A、TEKT1、FXR1、RPGRIP1、DPY19L2、DDX25、ZC3H14，及其相互關(guān)系(圖7)。

圖6 差異表達(dá)基因的KEGG分析

圓圈代表基因，連線(xiàn)代表基因之間的聯(lián)系。紅線(xiàn)：二者之間相互融合；綠線(xiàn)：二者之間互為接近；藍(lán)線(xiàn)：二者之間同時(shí)出現(xiàn)；紫線(xiàn)：實(shí)驗(yàn)獲得；黃線(xiàn)：文本挖掘獲得；淡藍(lán)色線(xiàn)：數(shù)據(jù)庫(kù)資料獲得；黑線(xiàn)：共表達(dá)。

從紅色到黃色，顏色越深說(shuō)明相關(guān)程度越高。

討 論

現(xiàn)已證明，miRNA在精子發(fā)生和發(fā)育調(diào)節(jié)中發(fā)揮作用[12]。本研究經(jīng)過(guò)數(shù)據(jù)庫(kù)及文獻(xiàn)的篩選，共確定了5個(gè)NOA相關(guān)的差異表達(dá)miRNAs，其中上調(diào)表達(dá)1個(gè)為miR-10b-5p，下調(diào)表達(dá)4個(gè)分別為miR-34b-5p，miR-34b-3p，miR-34c-5p，miR-449a。

本研究顯示NOA患者睪丸組織中miR-10b-5p表達(dá)下調(diào)。Wu等[13]研究發(fā)現(xiàn)miR-10b-5p在拮抗缺氧誘導(dǎo)的心肌細(xì)胞凋亡方面發(fā)揮一定的作用，可以推斷miR-10b-5p可能在精子發(fā)生過(guò)程中具有一定的作用。miR-34b-5p，miR-34b-3p和miR-34c-5p均屬于miR-34家族，已被證實(shí)在精子發(fā)生過(guò)程中發(fā)揮重要作用[14-15]。有研究發(fā)現(xiàn)miR-34b和miR-34c的靶基因dazl參與調(diào)控小鼠生殖細(xì)胞分化[16]，也有研究證明miR-34c通過(guò)調(diào)控Nanos2 破壞精原干細(xì)胞的穩(wěn)態(tài)[17]。巧合的是，本研究確定的miR-449a所屬的家族與miR-34家族已被證明在結(jié)構(gòu)上相似[12]。同時(shí)有研究表明miR-449家族和miR-34b/c 在小鼠睪丸雄性生殖細(xì)胞發(fā)育的調(diào)節(jié)中發(fā)揮著相應(yīng)的作用[18]。Marcet等[19]研究證明，miR-449a的表達(dá)受E2F1的正向調(diào)節(jié)，從而促進(jìn)細(xì)胞周期進(jìn)程并誘導(dǎo)受損細(xì)胞凋亡，E2F1的過(guò)表達(dá)導(dǎo)致精母細(xì)胞凋亡增加。因此，本研究篩選的5個(gè)差異表達(dá)miRNAs在精子發(fā)生過(guò)程中發(fā)揮重要作用。

本研究發(fā)現(xiàn)差異基因GO富集于生物過(guò)程(BP)，主要包括纖毛組裝，精子軸絲組裝，細(xì)胞蛋白質(zhì)代謝過(guò)程，精子發(fā)生等方面。細(xì)胞組成(CC)主要包括活動(dòng)纖毛，軸絲，精子鞭毛，中心粒等。分子功能(MF)主要包括蛋白質(zhì)結(jié)合，纖毛蛋白輕鏈的結(jié)合等。上述結(jié)果表明差異表達(dá)基因涉及到精子發(fā)生中的很多方面，上述的某個(gè)或多個(gè)環(huán)節(jié)出現(xiàn)問(wèn)題就有可能導(dǎo)致精子發(fā)生障礙，給我們后續(xù)針對(duì)5個(gè)miRNAs進(jìn)行功能研究時(shí)提供了方向。在NOA相關(guān)的差異表達(dá)基因所涉及的生物學(xué)功能中，精子發(fā)生嚴(yán)重受損是NOA的一個(gè)特征[20]。而本研究中差異表達(dá)基因的KEGG通路主要包括卵巢類(lèi)固醇激素生成、多個(gè)物種長(zhǎng)壽調(diào)控途徑、擴(kuò)張型心肌病、卵母細(xì)胞減數(shù)分裂以及孕酮介導(dǎo)的卵母細(xì)胞成熟等，上述大多數(shù)通路參與精子的發(fā)生[21-22]，且均有基因ADCY3和IGF1參與。本研究富集到的通路中“卵母細(xì)胞減數(shù)分裂”和“孕激素介導(dǎo)的卵母細(xì)胞成熟”為女性生殖相關(guān)的途徑，與Hu等[23]進(jìn)行的類(lèi)似研究中得到了相同的結(jié)論。因此，“卵巢類(lèi)固醇激素生成”途徑的某些基因也可能在精子發(fā)生中發(fā)揮作用。

本研究中差異表達(dá)基因的PPI網(wǎng)絡(luò)前20位的Hub基因中，Tektins(TEKTs)是一個(gè)高度保守的微管蛋白質(zhì)家族，已鑒定出五種 TEKT(TEKT1、2、3、4和5)。本研究的hub基因中主要包括TEKT1和TEKT3。有報(bào)道發(fā)現(xiàn)TEKT3在精子中主要與線(xiàn)粒體表面和中間部分的外部致密纖維相關(guān)，并且TEKT3被發(fā)現(xiàn)存在于精子頭部頂體膜的赤道段區(qū)域[24]，說(shuō)明TEKT3 不僅可以作為精子活力所需的鞭毛成分起作用，還可能參與頂體反應(yīng)。另外，Oiki等[25]研究發(fā)現(xiàn)TEKT1存在于精子鞭毛和頂體帽的頂端區(qū)域，而且在小鼠精子體外頂體反應(yīng)后，頂體相關(guān)的TEKT1消失。上述研究均說(shuō)明了TEKT3和TEKT1在精子發(fā)生過(guò)程中的重要角色。本研究中與精子鞭毛相關(guān)的基因還包括EFHC1、DNAH2和FXR1。Li等[26]在研究中證明EFHC1在精子鞭毛組裝相關(guān)蛋白中發(fā)揮重要作用。關(guān)于DNAH2基因即編碼動(dòng)力蛋白軸絲重鏈2，在最近的一項(xiàng)男性不育畸形精子癥的全外顯子組測(cè)序的研究中發(fā)現(xiàn)DNAH2的突變導(dǎo)致精子鞭毛的形態(tài)異常[27]。Huot等[28]在研究中觀(guān)察到成年小鼠睪丸中發(fā)現(xiàn)了最高水平表達(dá)的FXR1與微管元件相關(guān)。關(guān)于基因DYNLL2，雖然目前還未見(jiàn)其編碼動(dòng)力蛋白輕鏈LC8型2與男性生殖系統(tǒng)關(guān)系的研究報(bào)道，然而，Hu等[23]在類(lèi)似的研究關(guān)于hub基因的預(yù)測(cè)中也得出了DYNLL2基因。本研究中與精子纖毛相關(guān)的基因主要包括SPATA7和SPAG16。SPATA7基因編碼一種纖毛蛋白，通常在光感受器和精母細(xì)胞中表達(dá)[29]，SPATA7 基因的突變很可能會(huì)導(dǎo)致非梗阻性無(wú)精子癥。有研究已經(jīng)證明SPAG16基因功能的喪失會(huì)導(dǎo)致雄性不育和嚴(yán)重的精子活力缺陷，同時(shí)SPAG16基因的雜合突變降低了精子軸突中央裝置的穩(wěn)定性[30]。本研究中發(fā)現(xiàn)中心體相關(guān)的基因主要包括CETN1和CCDC146。Avasthi等[31]在關(guān)于雄性小鼠的中心體蛋白(centrin)的研究中證實(shí)CETN1在精子發(fā)生和精子細(xì)胞成熟的后期步驟中發(fā)揮重要作用。Firat-Karalar等[32]在通過(guò)哺乳動(dòng)物精子細(xì)胞的中體粒蛋白質(zhì)組質(zhì)譜鑒定的研究中發(fā)現(xiàn)CCDC146基因編碼中心體相關(guān)蛋白。PLCZ1已被證明是一種關(guān)鍵的精子卵母細(xì)胞激活因子的編碼基因，PLCZ1基因的突變已被證明會(huì)導(dǎo)致男性不育[33-34]。DPY19L2是導(dǎo)致精子頭部畸形的主要基因之一[35]，Castaneda等[36]研究發(fā)現(xiàn)DPY19L2是小鼠精子頂體發(fā)生和生育力所必需的基因。DDX25是促性腺激素和雄激素作用的靶點(diǎn)，是精子發(fā)生基因的轉(zhuǎn)錄后調(diào)節(jié)關(guān)鍵因子[37]。通過(guò)對(duì)上述hub基因的分析討論充分說(shuō)明了這些基因在精子發(fā)生過(guò)程中的重要意義，同時(shí)為下一步的研究提供了新的方向。

本研究存在一些局限性：(1)對(duì)miRNAs的篩選的文獻(xiàn)質(zhì)量有區(qū)別，樣本量和實(shí)驗(yàn)方法也存在一定的區(qū)別。(2)對(duì)于篩選差異表達(dá)基因的數(shù)據(jù)均是來(lái)源于GEO數(shù)據(jù)庫(kù)的芯片數(shù)據(jù)，樣本量較小，樣本的質(zhì)控等方面均對(duì)結(jié)果存在影響。(3)關(guān)鍵的miRNAs和hub基因需要通過(guò)實(shí)驗(yàn)來(lái)進(jìn)一步的研究驗(yàn)證。(4)精子發(fā)生是一個(gè)動(dòng)態(tài)過(guò)程，該分析僅在一個(gè)時(shí)間點(diǎn)提供有關(guān)基因表達(dá)和精子發(fā)生狀態(tài)的信息。

本研究通過(guò)一系列的生物信息學(xué)分析，結(jié)果表明miR-10b-5p、miR-34b-5p、miR-34b-3p、miR-34c-5p及miR-449a5個(gè)miRNAs在精子發(fā)生過(guò)程中發(fā)揮著重要的作用，為課題組后續(xù)研究miRNAs在精子發(fā)生和男性生育與不育的作用提供研究靶標(biāo)。NOA相關(guān)的差異表達(dá)基因的GO和KEGG分析結(jié)果提示精子發(fā)生過(guò)程中的很多環(huán)節(jié)發(fā)生問(wèn)題都有可能導(dǎo)致非梗阻性無(wú)精子癥。此外，ASRGL1、SPA4L、LIN7A、RPGRIP1、ZC3H14等基因未見(jiàn)與精子發(fā)生的相關(guān)研究，可能成為下一步對(duì)非梗阻性無(wú)精子癥相關(guān)基因進(jìn)行功能研究的重點(diǎn)。