基于生物信息學技術的非梗阻性無精子癥關鍵microRNAs和基因的研究

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

(1.唐山市婦幼保健院生殖遺傳科，唐山 063000；2.北京協和醫學院研究生院，北京 100730；3.國家衛生健康委科學技術研究所 男性生殖健康重點實驗室，北京 100081)

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

材料和方法

一、數據庫檢索

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

本研究中所用數據庫及檢索篩選流程見圖1。

圖1 基因數據庫篩選及分析流程圖

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

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

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

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

六、差異表達基因的蛋白質相互作用網絡和樞紐基因鑒定：利用基因相互作用檢索工具(the Search Tool for the Retrieval of Interacting Genes STRING)數據庫(https：//cn.string-db.org/)構建差異表達基因的蛋白質相互作用(protein-protein interaction PPI)網絡，同時要求交互分數不低于 0.7。然后，應用 Cytoscape 軟件v3.7.1(https：//cytoscape.org/)使源自STRING數據庫的PPI網絡可視化。應用 Cytoscape 中的 cytoHubba 插件對PPI 網絡中的節點根據度數計算方法進行排序，前20個基因被認為是樞紐基因。P<0.05被認為具有統計學意義的差異。

結 果

一、NOA相關的差異表達miRNAs

在PubMed、Embase和Web of Science中共篩選到5項研究中差異表達并經過驗證的miRNAs共56個(表1)。經過對文獻結果中差異表達的miRNAs繪制韋恩圖(圖1)，最終篩選出5個關鍵的候選miRNAs，其中上調表達1條為miR-10b-5p，下調表達4條分別為miR-34b-5p，miR-34b-3p，miR-34c-5p，miR-449a(表2)。

表1 篩選出的5項研究基本資料

圖2 五項研究差異表達的56個miRNAs的韋恩圖

表2 篩選出的5個候選miRNAs情況

二、NOA相關的差異表達基因篩選結果

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

A：上調基因；B：下調基因。

A：上調基因；B：下調基因。

三、差異表達基因的GO和KEGG分析

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

藍色代表生物過程(BP)；桔紅色代表細胞組分(CC)；綠色代表分子功能(MF)。

差異表達基因富集的KEGG通路主要包括卵巢類固醇激素生成，多個物種長壽調控途徑，擴張型心肌病，內分泌抵抗，孕酮介導的卵母細胞成熟等(圖6)。

四、差異表達基因的PPI網絡和樞紐基因

通過STRING數據庫構建出差異表達基因的PPI網絡(圖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，及其相互關系(圖7)。

圖6 差異表達基因的KEGG分析

圓圈代表基因，連線代表基因之間的聯系。紅線：二者之間相互融合；綠線：二者之間互為接近；藍線：二者之間同時出現；紫線：實驗獲得；黃線：文本挖掘獲得；淡藍色線：數據庫資料獲得；黑線：共表達。

從紅色到黃色，顏色越深說明相關程度越高。

討 論

現已證明，miRNA在精子發生和發育調節中發揮作用[12]。本研究經過數據庫及文獻的篩選，共確定了5個NOA相關的差異表達miRNAs，其中上調表達1個為miR-10b-5p，下調表達4個分別為miR-34b-5p，miR-34b-3p，miR-34c-5p，miR-449a。

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

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

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

本研究存在一些局限性：(1)對miRNAs的篩選的文獻質量有區別，樣本量和實驗方法也存在一定的區別。(2)對于篩選差異表達基因的數據均是來源于GEO數據庫的芯片數據，樣本量較小，樣本的質控等方面均對結果存在影響。(3)關鍵的miRNAs和hub基因需要通過實驗來進一步的研究驗證。(4)精子發生是一個動態過程，該分析僅在一個時間點提供有關基因表達和精子發生狀態的信息。

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