趙銘佳,郭劍,耿女,蘇興,李翠,姚雨宏,晏紅偉,韓寶生,劉美玲,盧文紅,谷翊群*
(1.唐山市婦幼保健院生殖遺傳科,唐山 063000;2.北京協(xié)和醫(yī)學院研究生院,北京 100730;3.國家衛(wèi)生健康委科學技術(shù)研究所 男性生殖健康重點實驗室,北京 100081)
全球范圍內(nèi)大約8%~12%的夫婦受到不孕癥困擾,其中男性因素大約占50%[1]。男性不育以無精子癥最為嚴重,分為梗阻性無精子癥(obstructive azoospermia,OA)和非梗阻性無精子癥(non-obstructive azoospermia,NOA)[2]。NOA病因復雜,大約73%的NOA病因不明[3]。精子發(fā)生是一個極其復雜的細胞分化過程,涉及調(diào)節(jié)生殖細胞發(fā)育和成熟的基因就多達2 300個,其中微小RNA(microRNA,miRNA)在轉(zhuǎn)錄后調(diào)節(jié)基因表達[4]。在以往的研究中,雖然已經(jīng)提出了許多與精子發(fā)生相關(guān)的候選基因,但許多候選基因尚未得到驗證或者在人群中非常罕見。研究表明,miRNA可能通過調(diào)節(jié)靶向基因mRNA的表達介導生精細胞的發(fā)育,參與精子發(fā)生,與男性不育有關(guān)[5]。本研究采用生物信息學的方法進行數(shù)據(jù)挖掘和分析,探討非梗阻性無精子癥相關(guān)的miRNAs和關(guān)鍵基因。
一、數(shù)據(jù)庫檢索
在NCBI的GEO數(shù)據(jù)庫(網(wǎng)址:https://www.ncbi.nlm.nih.gov/geo/)進行基因芯片數(shù)據(jù)檢索,確定檢索關(guān)鍵詞為“非梗阻性無精子癥(non obstructive azoospermia)”和“智人(Homo sapiens)”。在PubMed、Embase和Web of Science中進行檢索,確定檢索詞為“非梗阻性無精子癥(non obstructive azoospermia)”和“miRNA”或“microRNA”或“ncRNA”或“micro ribonucleic acid”。篩選并下載GSE145467和GSE108886表達譜數(shù)據(jù)集。GSE145467數(shù)據(jù)集,由Hodzic等[6]提供,基于Agilent-014850全人類基因組微陣列4x44K G4112F(特征編號版本)的GPL4133平臺,包括10個NOA樣本和10個OA睪丸樣本。GSE108886數(shù)據(jù)集,由Baksi等提供,基于Illumina HumanHT-12 V4.0表達微芯片的GPL10558平臺,包含 8個NOA樣本和4個OA樣本(包括1個睪丸對照樣本)。
本研究中所用數(shù)據(jù)庫及檢索篩選流程見圖1。
圖1 基因數(shù)據(jù)庫篩選及分析流程圖
二、篩選差異表達miRNAs:在GEO數(shù)據(jù)庫搜索已發(fā)表的NOA的miRNA數(shù)據(jù)集合,但是無可用的數(shù)據(jù)。然后,在PubMed、Embase和Web of Science中進行檢索,以文獻中報道非梗阻性無精子癥睪丸組織中存在差異表達的miRNA并且經(jīng)過PCR驗證為文獻篩選標準,共篩選出五項包含miRNA的研究(表1)[7-11]。該五項研究均確定了miRNA的變化趨勢,再以至少兩篇文獻報道且經(jīng)PCR驗證為差異表達的miRNAs作為研究對象的miRNA篩選標準,篩選出差異表達miRNA。為研究關(guān)鍵miRNA在非梗阻性無精子癥中的作用,選擇具有一致表達變化趨勢的miRNAs作為本研究的候選差異表達miRNAs用于下面的分析。
三、差異表達miRNAs的靶基因預測:將上述篩選的miRNAs匯總,通過聯(lián)川生物云平臺的“miRNA靶基因預測云分析”(https://www.omicstudio.cn/analysis/targetGene)預測其靶基因。
四、篩選差異表達基因:利用GEO2R工具(https://www-ncbi-nlm-nih-gov-s.webvpn.cams.cn/geo/geo2r/)在GEO數(shù)據(jù)庫中GSE145467和GSE108886表達譜數(shù)據(jù)集的NOA和匹配的OA中篩選差異表達基因。GSE145467篩選標準是FDR小于0.05和log2FC的絕對值大于2;GSE108886的篩選標準是FDR小于0.05和log2FC的絕對值大于1。對上述兩個數(shù)據(jù)集進行Venn Diagram分析以確定差異表達基因的交集,再應用聯(lián)川生物云平臺(https://www.omicstudio.cn/tool/6)與miRNAs預測靶基因取交集,獲得最終鑒定的差異表達基因并繪制韋恩圖。
五、差異表達基因的GO和KEGG富集分析:利用注釋可視化集成發(fā)現(xiàn)數(shù)據(jù)庫(The Database for Annotation,Visualization and Integrated Discovery DAVID)(https://david.ncifcrf.gov/)進行GO功能和KEGG通路富集分析以確定常見差異表達基因的功能。選擇物種為智人(Homo sapiens)。P<0.05表示差異有統(tǒng)計學意義。
六、差異表達基因的蛋白質(zhì)相互作用網(wǎng)絡和樞紐基因鑒定:利用基因相互作用檢索工具(the Search Tool for the Retrieval of Interacting Genes STRING)數(shù)據(jù)庫(https://cn.string-db.org/)構(gòu)建差異表達基因的蛋白質(zhì)相互作用(protein-protein interaction PPI)網(wǎng)絡,同時要求交互分數(shù)不低于 0.7。然后,應用 Cytoscape 軟件v3.7.1(https://cytoscape.org/)使源自STRING數(shù)據(jù)庫的PPI網(wǎng)絡可視化。應用 Cytoscape 中的 cytoHubba 插件對PPI 網(wǎng)絡中的節(jié)點根據(jù)度數(shù)計算方法進行排序,前20個基因被認為是樞紐基因。P<0.05被認為具有統(tǒng)計學意義的差異。
一、NOA相關(guān)的差異表達miRNAs
在PubMed、Embase和Web of Science中共篩選到5項研究中差異表達并經(jīng)過驗證的miRNAs共56個(表1)。經(jīng)過對文獻結(jié)果中差異表達的miRNAs繪制韋恩圖(圖1),最終篩選出5個關(guān)鍵的候選miRNAs,其中上調(diào)表達1條為miR-10b-5p,下調(diào)表達4條分別為miR-34b-5p,miR-34b-3p,miR-34c-5p,miR-449a(表2)。
表1 篩選出的5項研究基本資料
圖2 五項研究差異表達的56個miRNAs的韋恩圖
表2 篩選出的5個候選miRNAs情況
二、NOA相關(guān)的差異表達基因篩選結(jié)果
通過聯(lián)川生物云平臺的“miRNA靶基因預測云分析”共預測靶基因6 695個。在GSE108886數(shù)據(jù)集中篩選出了2 036個差異表達基因,包括1 674個上調(diào)基因和362個下調(diào)基因。在GSE145467數(shù)據(jù)集中篩選出1 274個差異表達基因,包括1,182個上調(diào)基因和92個下調(diào)基因。對上述兩個數(shù)據(jù)集進行Venn Diagram分析取交集,共鑒定出797個差異,包括776個上調(diào)基因和21個常見的下調(diào)基因(圖3)。再將上述數(shù)據(jù)集的DEGs與預測獲得的miRNAs靶基因取交集,最終鑒定DEGs共180個,其中上調(diào)基因173個,下調(diào)基因7個(圖4)。篩選流程詳見圖1。
A:上調(diào)基因;B:下調(diào)基因。
A:上調(diào)基因;B:下調(diào)基因。
三、差異表達基因的GO和KEGG分析
差異基因GO富集于生物過程(BP),主要包括纖毛組裝,精子軸絲組裝,鞭毛相關(guān)的精子活力,纖毛依賴性細胞運動,頂體組裝,精子發(fā)生,精子細胞發(fā)育,血小板脫顆粒,細胞蛋白質(zhì)代謝過程,蛋白質(zhì)穩(wěn)定化等;細胞組成(CC)主要包括活動纖毛,軸絲,精子鞭毛,中心粒,細胞外空間等;分子功能(MF)主要包括蛋白質(zhì)結(jié)合,纖毛蛋白輕鏈的結(jié)合,整合素結(jié)合等(圖5)。
藍色代表生物過程(BP);桔紅色代表細胞組分(CC);綠色代表分子功能(MF)。
差異表達基因富集的KEGG通路主要包括卵巢類固醇激素生成,多個物種長壽調(diào)控途徑,擴張型心肌病,內(nèi)分泌抵抗,孕酮介導的卵母細胞成熟等(圖6)。
四、差異表達基因的PPI網(wǎng)絡和樞紐基因
通過STRING數(shù)據(jù)庫構(gòu)建出差異表達基因的PPI網(wǎng)絡(圖6)。Cytoscape 中的 cytoHubba 插件篩選出了排名前20的樞紐hub基因,分別為TEKT3、EFHC1、DYNLL2、DNAH2、CETN1、SPATA7、ASRGL1、CCDC146、PLCZ1、SPAG16、DNAL1、EFCAB11、SPA4L、LIN7A、TEKT1、FXR1、RPGRIP1、DPY19L2、DDX25、ZC3H14,及其相互關(guān)系(圖7)。
圖6 差異表達基因的KEGG分析
圓圈代表基因,連線代表基因之間的聯(lián)系。紅線:二者之間相互融合;綠線:二者之間互為接近;藍線:二者之間同時出現(xiàn);紫線:實驗獲得;黃線:文本挖掘獲得;淡藍色線:數(shù)據(jù)庫資料獲得;黑線:共表達。
從紅色到黃色,顏色越深說明相關(guān)程度越高。
現(xiàn)已證明,miRNA在精子發(fā)生和發(fā)育調(diào)節(jié)中發(fā)揮作用[12]。本研究經(jīng)過數(shù)據(jù)庫及文獻的篩選,共確定了5個NOA相關(guān)的差異表達miRNAs,其中上調(diào)表達1個為miR-10b-5p,下調(diào)表達4個分別為miR-34b-5p,miR-34b-3p,miR-34c-5p,miR-449a。
本研究顯示NOA患者睪丸組織中miR-10b-5p表達下調(diào)。Wu等[13]研究發(fā)現(xiàn)miR-10b-5p在拮抗缺氧誘導的心肌細胞凋亡方面發(fā)揮一定的作用,可以推斷miR-10b-5p可能在精子發(fā)生過程中具有一定的作用。miR-34b-5p,miR-34b-3p和miR-34c-5p均屬于miR-34家族,已被證實在精子發(fā)生過程中發(fā)揮重要作用[14-15]。有研究發(fā)現(xiàn)miR-34b和miR-34c的靶基因dazl參與調(diào)控小鼠生殖細胞分化[16],也有研究證明miR-34c通過調(diào)控Nanos2 破壞精原干細胞的穩(wěn)態(tài)[17]。巧合的是,本研究確定的miR-449a所屬的家族與miR-34家族已被證明在結(jié)構(gòu)上相似[12]。同時有研究表明miR-449家族和miR-34b/c 在小鼠睪丸雄性生殖細胞發(fā)育的調(diào)節(jié)中發(fā)揮著相應的作用[18]。Marcet等[19]研究證明,miR-449a的表達受E2F1的正向調(diào)節(jié),從而促進細胞周期進程并誘導受損細胞凋亡,E2F1的過表達導致精母細胞凋亡增加。因此,本研究篩選的5個差異表達miRNAs在精子發(fā)生過程中發(fā)揮重要作用。
本研究發(fā)現(xiàn)差異基因GO富集于生物過程(BP),主要包括纖毛組裝,精子軸絲組裝,細胞蛋白質(zhì)代謝過程,精子發(fā)生等方面。細胞組成(CC)主要包括活動纖毛,軸絲,精子鞭毛,中心粒等。分子功能(MF)主要包括蛋白質(zhì)結(jié)合,纖毛蛋白輕鏈的結(jié)合等。上述結(jié)果表明差異表達基因涉及到精子發(fā)生中的很多方面,上述的某個或多個環(huán)節(jié)出現(xiàn)問題就有可能導致精子發(fā)生障礙,給我們后續(xù)針對5個miRNAs進行功能研究時提供了方向。在NOA相關(guān)的差異表達基因所涉及的生物學功能中,精子發(fā)生嚴重受損是NOA的一個特征[20]。而本研究中差異表達基因的KEGG通路主要包括卵巢類固醇激素生成、多個物種長壽調(diào)控途徑、擴張型心肌病、卵母細胞減數(shù)分裂以及孕酮介導的卵母細胞成熟等,上述大多數(shù)通路參與精子的發(fā)生[21-22],且均有基因ADCY3和IGF1參與。本研究富集到的通路中“卵母細胞減數(shù)分裂”和“孕激素介導的卵母細胞成熟”為女性生殖相關(guān)的途徑,與Hu等[23]進行的類似研究中得到了相同的結(jié)論。因此,“卵巢類固醇激素生成”途徑的某些基因也可能在精子發(fā)生中發(fā)揮作用。
本研究中差異表達基因的PPI網(wǎng)絡前20位的Hub基因中,Tektins(TEKTs)是一個高度保守的微管蛋白質(zhì)家族,已鑒定出五種 TEKT(TEKT1、2、3、4和5)。本研究的hub基因中主要包括TEKT1和TEKT3。有報道發(fā)現(xiàn)TEKT3在精子中主要與線粒體表面和中間部分的外部致密纖維相關(guān),并且TEKT3被發(fā)現(xiàn)存在于精子頭部頂體膜的赤道段區(qū)域[24],說明TEKT3 不僅可以作為精子活力所需的鞭毛成分起作用,還可能參與頂體反應。另外,Oiki等[25]研究發(fā)現(xiàn)TEKT1存在于精子鞭毛和頂體帽的頂端區(qū)域,而且在小鼠精子體外頂體反應后,頂體相關(guān)的TEKT1消失。上述研究均說明了TEKT3和TEKT1在精子發(fā)生過程中的重要角色。本研究中與精子鞭毛相關(guān)的基因還包括EFHC1、DNAH2和FXR1。Li等[26]在研究中證明EFHC1在精子鞭毛組裝相關(guān)蛋白中發(fā)揮重要作用。關(guān)于DNAH2基因即編碼動力蛋白軸絲重鏈2,在最近的一項男性不育畸形精子癥的全外顯子組測序的研究中發(fā)現(xiàn)DNAH2的突變導致精子鞭毛的形態(tài)異常[27]。Huot等[28]在研究中觀察到成年小鼠睪丸中發(fā)現(xiàn)了最高水平表達的FXR1與微管元件相關(guān)。關(guān)于基因DYNLL2,雖然目前還未見其編碼動力蛋白輕鏈LC8型2與男性生殖系統(tǒng)關(guān)系的研究報道,然而,Hu等[23]在類似的研究關(guān)于hub基因的預測中也得出了DYNLL2基因。本研究中與精子纖毛相關(guān)的基因主要包括SPATA7和SPAG16。SPATA7基因編碼一種纖毛蛋白,通常在光感受器和精母細胞中表達[29],SPATA7 基因的突變很可能會導致非梗阻性無精子癥。有研究已經(jīng)證明SPAG16基因功能的喪失會導致雄性不育和嚴重的精子活力缺陷,同時SPAG16基因的雜合突變降低了精子軸突中央裝置的穩(wěn)定性[30]。本研究中發(fā)現(xiàn)中心體相關(guān)的基因主要包括CETN1和CCDC146。Avasthi等[31]在關(guān)于雄性小鼠的中心體蛋白(centrin)的研究中證實CETN1在精子發(fā)生和精子細胞成熟的后期步驟中發(fā)揮重要作用。Firat-Karalar等[32]在通過哺乳動物精子細胞的中體粒蛋白質(zhì)組質(zhì)譜鑒定的研究中發(fā)現(xiàn)CCDC146基因編碼中心體相關(guān)蛋白。PLCZ1已被證明是一種關(guān)鍵的精子卵母細胞激活因子的編碼基因,PLCZ1基因的突變已被證明會導致男性不育[33-34]。DPY19L2是導致精子頭部畸形的主要基因之一[35],Castaneda等[36]研究發(fā)現(xiàn)DPY19L2是小鼠精子頂體發(fā)生和生育力所必需的基因。DDX25是促性腺激素和雄激素作用的靶點,是精子發(fā)生基因的轉(zhuǎn)錄后調(diào)節(jié)關(guān)鍵因子[37]。通過對上述hub基因的分析討論充分說明了這些基因在精子發(fā)生過程中的重要意義,同時為下一步的研究提供了新的方向。
本研究存在一些局限性:(1)對miRNAs的篩選的文獻質(zhì)量有區(qū)別,樣本量和實驗方法也存在一定的區(qū)別。(2)對于篩選差異表達基因的數(shù)據(jù)均是來源于GEO數(shù)據(jù)庫的芯片數(shù)據(jù),樣本量較小,樣本的質(zhì)控等方面均對結(jié)果存在影響。(3)關(guān)鍵的miRNAs和hub基因需要通過實驗來進一步的研究驗證。(4)精子發(fā)生是一個動態(tài)過程,該分析僅在一個時間點提供有關(guān)基因表達和精子發(fā)生狀態(tài)的信息。
本研究通過一系列的生物信息學分析,結(jié)果表明miR-10b-5p、miR-34b-5p、miR-34b-3p、miR-34c-5p及miR-449a5個miRNAs在精子發(fā)生過程中發(fā)揮著重要的作用,為課題組后續(xù)研究miRNAs在精子發(fā)生和男性生育與不育的作用提供研究靶標。NOA相關(guān)的差異表達基因的GO和KEGG分析結(jié)果提示精子發(fā)生過程中的很多環(huán)節(jié)發(fā)生問題都有可能導致非梗阻性無精子癥。此外,ASRGL1、SPA4L、LIN7A、RPGRIP1、ZC3H14等基因未見與精子發(fā)生的相關(guān)研究,可能成為下一步對非梗阻性無精子癥相關(guān)基因進行功能研究的重點。