基于生物信息學(xué)方法的小鼠精子發(fā)生關(guān)鍵基因篩選

李 欣 ，田 佳，脫征軍，邵懷峰，張偉新，李委奇，王 琨，徐 寒，溫 萬(wàn)*

（1.中國(guó)農(nóng)業(yè)科學(xué)院北京畜牧獸醫(yī)研究所，北京 100193；2.寧夏畜牧工作站，寧夏銀川 750105）

轉(zhuǎn)錄組是某個(gè)物種或特定細(xì)胞在某一功能狀態(tài)下產(chǎn)生的所有RNA 總和，通過(guò)對(duì)不同樣本的轉(zhuǎn)錄組差異研究，不僅可以了解特定時(shí)間和空間條件下疾病的發(fā)生發(fā)展?fàn)顟B(tài)，也能作為研究疾病診斷和預(yù)后標(biāo)志物的重要研究手段[1]。目前，RNA-seq 已經(jīng)成為轉(zhuǎn)錄組研究中最為常見(jiàn)的手段，來(lái)自不同動(dòng)物的多個(gè)器官的大量RNA-seq數(shù)據(jù)集陸續(xù)出現(xiàn)[1-2]，這些數(shù)據(jù)保存于NCBI 數(shù)據(jù)庫(kù)中。精原細(xì)胞處于雄性生殖細(xì)胞的早期發(fā)育階段，能不斷進(jìn)行有絲分裂，分化為精母細(xì)胞。精母細(xì)胞經(jīng)減數(shù)分裂成為單倍體的精細(xì)胞。單倍體的精細(xì)胞經(jīng)歷形態(tài)變化、細(xì)胞器重排以及代謝途徑和方式的改變等，最終形成終末分化的成熟精子，攜帶著由基因和環(huán)境共同決定的遺傳信息。研究顯示，Sycp3 與精子減數(shù)分裂染色體重組過(guò)程密切相關(guān)[3]，Crem 在精母細(xì)胞到圓形精子的轉(zhuǎn)換過(guò)程起到開(kāi)關(guān)的作用[4]。Zuo 等[5]對(duì)小鼠精子變形過(guò)程中的圓形精細(xì)胞、長(zhǎng)形精細(xì)胞及附睪尾精子進(jìn)行了RNAseq 和生物學(xué)分析，獲得了大量選擇性剪接事件、差異表達(dá)基因和新轉(zhuǎn)錄本信息。雖然對(duì)精子發(fā)生已有大量研究，但是從精原細(xì)胞依次到成熟精子的整體研究很少，且整個(gè)過(guò)程的分子機(jī)制仍不清楚。模式動(dòng)物小鼠容易飼養(yǎng)，繁殖率高，遺傳上有較高的純合度，代謝類型、生理疾病等與人類接近，也具有完整的精子變形過(guò)程，樣本取材方便。因此，本研究通過(guò)整合公共數(shù)據(jù)庫(kù)中的轉(zhuǎn)錄組數(shù)據(jù)，將小鼠精原細(xì)胞和變形期精細(xì)胞和成熟精子進(jìn)行轉(zhuǎn)錄本比較，以期全面揭示小鼠精子發(fā)生過(guò)程中基因表達(dá)的動(dòng)態(tài)性與時(shí)空特異性，為后續(xù)全面研究精子發(fā)生調(diào)控機(jī)制提供理論基礎(chǔ)。

1 材料和方法

1.1 數(shù)據(jù)來(lái)源 測(cè)序數(shù)據(jù)來(lái)自于GEO 數(shù)據(jù)庫(kù)（https://www.ncbi.nlm.nih.gov/geo/），登錄號(hào)為GSE43717。小鼠轉(zhuǎn)錄組原始數(shù)據(jù)信息如表1 所示。

1.2 實(shí)驗(yàn)方法

1.2.1 數(shù)據(jù)預(yù)處理及差異表達(dá)基因篩選 將測(cè)序得到的原始序列（Raw reads）進(jìn)行過(guò)濾去雜，去除測(cè)序接頭序列、重復(fù)冗余序列和低質(zhì)量序列，獲得高質(zhì)量序列數(shù)據(jù)（clean reads）。然后用TopHat 2 將各樣品質(zhì)控后得到的clean reads 與小鼠參考基因組（GRCm38.p6）進(jìn)行比對(duì)，比對(duì)過(guò)程中允許有2 個(gè)堿基的錯(cuò)配。利用FPKM 值反映對(duì)應(yīng)基因的表達(dá)量，使用R 軟件的Ballgown 包進(jìn)行差異分析。差異表達(dá)基因的篩選首先用取對(duì)數(shù)后的倍數(shù)變化（Fold change,FC）和統(tǒng)計(jì)檢驗(yàn)P值初選，然后用FDR（False Discovery Ratio）校驗(yàn)方法對(duì)P值進(jìn)行假陽(yáng)性檢驗(yàn)。當(dāng)FDR≤0.05 且|log2(FC)|＞1 時(shí)，則為差異表達(dá)基因。

表1 小鼠轉(zhuǎn)錄組原始數(shù)據(jù)信息

1.2.2 差異基因的GO 分析和KEGG 富集分析 通過(guò)WebGestalt（WEB-based GEne SeT AnaLysis Toolkit，http://www.webgestalt.org/）分析工具對(duì)篩選出的差異表達(dá)基因進(jìn)行GO 和KEGG 富集分析。

1.2.3 關(guān)鍵基因篩選 篩選出關(guān)鍵通路，將所包含的基因進(jìn)行FPKM 值排序，取前18 名的基因進(jìn)行詳細(xì)的功能分析。

2 結(jié)果

2.1 差異表達(dá)基因篩選結(jié)果 對(duì)雄性小鼠不同階段生殖細(xì)胞的轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行預(yù)處理，并通過(guò)閾值篩選得到2個(gè)不同階段的差異基因，從精原細(xì)胞到圓形精細(xì)胞共獲得10 785 個(gè)差異表達(dá)基因，其中上調(diào)基因6 176 個(gè)，下調(diào)基因4 609 個(gè)（圖1-A）；從圓形精細(xì)胞到成熟精子共獲得8 545 個(gè)基因，其中上調(diào)基因3 770 個(gè)，下調(diào)基因4 775 個(gè)（圖1-B）。

2.2 差異基因通路富集分析 由圖2 可知，精原細(xì)胞到圓形精細(xì)胞過(guò)程的差異基因主要參與了生物合成過(guò)程（生物過(guò)程）、質(zhì)膜和囊泡部分（細(xì)胞組分）及陰離子配位（細(xì)胞功能）等過(guò)程。由圖3 可知，圓形精細(xì)胞到成熟精子的差異基因主要參與了細(xì)胞器組織和芳香族化合物生物合成（生物過(guò)程）、水解酶活性（細(xì)胞功能）等過(guò)程，同時(shí)，可以明顯看出圓形精細(xì)胞變?yōu)槌墒炀拥幕蚨嗉杏谏镞^(guò)程。

圖1 差異基因火山圖

精原細(xì)胞到圓形精細(xì)胞的差異基因，富集的KEGG信號(hào)通路主要是肌動(dòng)蛋白細(xì)胞骨架的調(diào)節(jié)和絲裂原活化蛋白激酶（Mitogen-activated Protein Kinase,MAPK）信號(hào)通路，而圓形精細(xì)胞到成熟精子上調(diào)差異基因主要是緊密連接和MAPK 信號(hào)通路。

2.3 關(guān)鍵基因的篩選 鑒于GO 和KEGG 富集分析獲得的條目較多，本研究根據(jù)精子發(fā)生過(guò)程的關(guān)鍵點(diǎn)篩選出相關(guān)性較強(qiáng)的8 個(gè)通路，并對(duì)這些基因進(jìn)行了查重去冗余，獲得628 個(gè)差異基因；再將638 個(gè)基因在3 個(gè)階段的FPKM 值進(jìn)行排序，取FPKM 值最大的前18 個(gè)基因進(jìn)行研究，該18 個(gè)基因所屬的GO 條目詳細(xì)信息見(jiàn)表2 和圖4。表2 所示的7 個(gè)通路是精子發(fā)生過(guò)程的關(guān)鍵條目。該條目中，Reproduction（GO:0000003）和Reproductive process（GO:0022414）條目中的基因可能與繁殖相關(guān)；Motile cilium（GO:0031514）、Sperm flagellum（GO:0036126）和Cilium（GO:0005929）條目中的基因可能與精子尾部形成相關(guān)；Golgi apparatus（GO:0005794）條目中的基因可能與精子頂體的形成相關(guān)；Spermatogenesis（GO:0007283）條目中的基因可能與精子發(fā)生相關(guān)。

圖2 精原細(xì)胞到圓形精細(xì)胞過(guò)程的上下調(diào)基因GO/KEGG 富集分析

圖3 圓形精細(xì)胞到成熟精子過(guò)程的上下調(diào)基因GO/KEGG 富集分析

表2 精子發(fā)生過(guò)程的關(guān)鍵條目

3 討 論

本研究通過(guò)分析GSE43717 轉(zhuǎn)錄組測(cè)序數(shù)據(jù)，從精原細(xì)胞到圓形精細(xì)胞篩選到的上調(diào)基因比下調(diào)基因多1 576 個(gè)，而從圓形精細(xì)胞到成熟精子篩選到的下調(diào)基因比上調(diào)基因多1 005 個(gè)，這說(shuō)明在早期單倍體細(xì)胞、圓形精子細(xì)胞中合成了許多新的mRNA，但在晚期單倍體細(xì)胞、長(zhǎng)形精細(xì)胞中，染色質(zhì)固縮進(jìn)而導(dǎo)致結(jié)構(gòu)的改變使轉(zhuǎn)錄失活[7]。通過(guò)GO 分析發(fā)現(xiàn)差異基因主要參與精子質(zhì)膜和囊泡形成、生物合成、水解酶活性過(guò)程。差異表達(dá)基因主要富集的KEGG 信號(hào)通路主要是肌動(dòng)蛋白細(xì)胞骨架的調(diào)節(jié)、緊密連接和MAPK 信號(hào)傳導(dǎo)通路。其中MAPK 信號(hào)傳導(dǎo)通路是貫穿精原細(xì)胞到精子整個(gè)過(guò)程的共有部分，該通路不但參與了精子生成與凋亡、運(yùn)動(dòng)、獲能和頂體反應(yīng)等過(guò)程，還是附睪中精子成熟的重要通路。這些重要生物學(xué)過(guò)程以及生物通路的發(fā)現(xiàn)，是根據(jù)精子發(fā)生過(guò)程關(guān)鍵點(diǎn)的轉(zhuǎn)錄組測(cè)序、差異基因和FPKM 值較大基因篩選而來(lái)的，可見(jiàn)在研究生物發(fā)育過(guò)程中，選擇合適的時(shí)間點(diǎn)、差異基因篩選標(biāo)準(zhǔn)以及基因表達(dá)量對(duì)機(jī)理的研究至關(guān)重要。

圖4 關(guān)鍵基因的熱圖分析

3.1 參與精子發(fā)生相關(guān)基因Prm2、Selenof、Tnp1、Tnp2、Ropn1l、Pafah1b2、Spink2、Prm1參與精子發(fā)生過(guò)程。在精子發(fā)生過(guò)程中，mRNA 的定位和翻譯被認(rèn)為是以一種階段性的方式調(diào)節(jié)的。有報(bào)道顯示，Prm2mRNA 位于圓形精細(xì)胞的擬染色體，在長(zhǎng)形精細(xì)胞的胞漿進(jìn)行翻譯，說(shuō)明Prm2在整個(gè)精子發(fā)生過(guò)程中具有時(shí)空調(diào)節(jié)的特性[8]。而且Prm1、Prm2、Tnp1 和Tnp2均是在染色質(zhì)重塑過(guò)程中取代組蛋白的基本染色體蛋白質(zhì)，它們?cè)诰幼冃吻捌诒晦D(zhuǎn)錄并以惰性狀態(tài)存在，并在長(zhǎng)形精子細(xì)胞中被激活翻譯[9-10]。絲氨酸蛋白酶抑制劑Kazal 2 型（Spink2）基因?qū)儆赟PINK 家族。研究表明，Spink2在小鼠睪丸中表達(dá)量較高，粗線期精母細(xì)胞階段首次顯著表達(dá)，而在附睪中表達(dá)量顯著降低[11]。Ropn1l是構(gòu)成精子鞭毛纖維鞘的關(guān)鍵基因之一，直接影響到精子活力和生育力[12]，這也是精子發(fā)生過(guò)程的重要一環(huán)。另外，Selenof可能參與二硫鍵形成相關(guān)的氧化還原反應(yīng)，質(zhì)控經(jīng)內(nèi)質(zhì)網(wǎng)折疊所產(chǎn)生的蛋白質(zhì)[13-14]，進(jìn)而影響精子發(fā)生過(guò)程。

3.2 與精子尾部形成相關(guān)基因 精子尾部結(jié)構(gòu)與精子運(yùn)動(dòng)功能密切相關(guān)，具有正常尾部結(jié)構(gòu)的精子才能游動(dòng)，尾部每一個(gè)部分的發(fā)育均受到相關(guān)基因的精確調(diào)控，如果發(fā)生基因缺失或突變會(huì)導(dǎo)致相應(yīng)編碼蛋白的改變，從而導(dǎo)致精子尾部結(jié)構(gòu)異?；虍a(chǎn)生功能障礙[15]。Ldhc、Hsp90aa1、Gstm5、Csnk2b、Spa17均與精子尾部形成相關(guān)。Hsp90aa1 屬于分子伴侶，參與ATP 酶活性相關(guān)的功能循環(huán)，從而誘導(dǎo)目的蛋白的構(gòu)象變化，導(dǎo)致其活化[16]。Hsp90aa1 位于人精子頸部、尾部中段等區(qū)域，在精子獲能時(shí)表達(dá)量增加，在胞內(nèi)鈣穩(wěn)態(tài)以及蛋白酪氨酸磷酸化方面發(fā)揮重要作用；在氧化應(yīng)激條件下，可通過(guò)激活一氧化氮合酶來(lái)保持精子活力[17]。Gstm5 屬于谷胱甘肽（GST）家族，在纖維鞘中表達(dá)，促進(jìn)生殖細(xì)胞的增殖和分化，同時(shí)，在氧化應(yīng)激過(guò)程中可能發(fā)生改變，以對(duì)抗自由基對(duì)睪丸和生精細(xì)胞的攻擊[18]。Spa17是一種在睪丸組織中特異性表達(dá)的蛋白，定位于初級(jí)精母細(xì)胞、次級(jí)精母細(xì)胞和成熟精子中，而精原細(xì)胞、支持細(xì)胞和間質(zhì)細(xì)胞中無(wú)表達(dá)[19]。成熟精子中Spa17 主要表達(dá)于精子尾部[20]，推測(cè)其與精子運(yùn)動(dòng)有關(guān)。另外，Ldhc和Csnk2b未見(jiàn)在精子中的相關(guān)研究。

3.3 與繁殖相關(guān)的基因 與繁殖相關(guān)的基因包括Txnrd3、Tcp1、Gpx4、Dkkl1和Dbil5。動(dòng)物硫氧還蛋白還原酶Txnrd 是還原型輔酶Ⅱ（NADPH）依賴性蛋白質(zhì)，屬于吡啶核苷酸二硫化物氧化還原酶家族，含有黃素腺嘌呤二核苷酸（FAD），已知的有Txnrd1、Txnrd2 和Txnrd3。其中，由Txnrd3基因編碼的同工酶主要存在于睪丸和發(fā)育早期的精子細(xì)胞中，在結(jié)構(gòu)上有一個(gè)谷胱甘肽還原型結(jié)構(gòu)域，推測(cè)它參與精子成熟過(guò)程中的氧化還原反應(yīng)[21]。Tcp1 即t 復(fù)合多肽1，在精子發(fā)生過(guò)程中表達(dá)量較高，定位于精子頂體囊泡中，影響精細(xì)胞形成過(guò)程中頂體的形成[22]。Gpx4 即谷胱甘肽過(guò)氧化物酶4，是一種膜相關(guān)的谷胱甘肽過(guò)氧化物酶，在哺乳動(dòng)物睪丸中表達(dá)量較高，其亞型均參與精子發(fā)育[23]，線粒體Gpx4（mGpx4）敲除小鼠精子質(zhì)量下降，進(jìn)而導(dǎo)致不育，該基因的所有亞型全部敲除后引起胚胎致死[24]。Dkkl1的mRNA 主要產(chǎn)生于精母細(xì)胞，其蛋白定位于頂體，可促進(jìn)精子穿透透明帶[25]。而且，Dkkl1在弱精子癥患者的精子中表達(dá)明顯降低[26]，推測(cè)其在精子運(yùn)動(dòng)中起重要作用。另外，研究顯示，Dbil5 是一種睪丸特異性蛋白，在精子發(fā)生的后期表達(dá)[27]，大鼠Dbil5 的mRNAs 翻譯在精細(xì)胞早期受到強(qiáng)烈抑制，在精細(xì)胞后期被激活[28]，推測(cè)其可能通過(guò)影響精子后期的變形，進(jìn)而影響繁殖。

4 結(jié) 論

本研究通過(guò)對(duì)小鼠精子形成期3 個(gè)階段的差異表達(dá)基因進(jìn)行分析，發(fā)現(xiàn)FPKM 值靠前的差異基因主要集中在精子發(fā)生、精子尾部形成以及繁殖等條目。尤其與精子發(fā)生相關(guān)，包含魚精蛋白在內(nèi)的8 個(gè)基因是學(xué)者研究比較全面的，體現(xiàn)了這些基因的重要性。此外，對(duì)繁殖相關(guān)的基因研究也較多，主要關(guān)系到精子受精能力引起的不育現(xiàn)象。總之，精子形成期間發(fā)生了大量基因的變動(dòng)，關(guān)注期間關(guān)鍵基因可為精子變形機(jī)制的闡明提供理論基礎(chǔ)。