基于高通量測序技術(shù)篩選與山羊精子發(fā)生有關(guān)的lncRNAs

代子諾,何其杰,喬 蕾,何莉娜,黃思藝,趙中權(quán) (西南大學(xué) 動物科技學(xué)院,重慶 北碚400715)

大量的轉(zhuǎn)錄本通常超過200個核苷酸,這些轉(zhuǎn)錄本通常是多聚腺苷酸化的,并且沒有明顯的開放閱讀框(ORF),它們被定義為長鏈非編碼RNA(lncRNAs)[1]。lncRNAs參與了許多重要的生理過程,如X 染色體失活、遺傳印記、染色質(zhì)修飾和細(xì)胞周期調(diào)節(jié)[2-3]。由于lncRNAs的功能眾多且重要,對于lnc RNAs的進一步探索已成為生物信息學(xué)的一個重要方向。

之前的研究已經(jīng)證實,lnc RNAs在睪丸發(fā)育和精子發(fā)生中起著重要的調(diào)節(jié)作用。ANGUERA等[4]的研究結(jié)果表明了lncRNAtsx在粗線期精母細(xì)胞中有特異性表達,但在精原細(xì)胞和圓形精母細(xì)胞中不表達。小鼠tsx基因的敲除促進了更多精子的凋亡,證實了其對精子發(fā)育中的重要作用。

山羊是人類重要的經(jīng)濟家畜,主要提供肉、奶、羊絨、皮等多種產(chǎn)品。與其他山羊相比,大足黑山羊具有較高的繁殖率。大足黑山羊母羊的平均窩產(chǎn)羔數(shù)為2.72只,公山羊具有性成熟早、生長發(fā)育快、抗病性強、遺傳穩(wěn)定等特點。

本研究采用高通量測序方法分析了成熟和不成熟山羊睪丸中差異表達的lncRNAs,并對靶基因進行了預(yù)測。此外,還對這些靶基因進行了分類注釋、GO功能分析和KEGG 通路分析,篩選出了與精子發(fā)生有關(guān)的候選靶基因及l(fā)ncRNA。這些數(shù)據(jù)有助于研究人員深入了解睪丸發(fā)育過程以及生殖過程,為進一步了解lncRNAs的功能提供了許多有價值的信息。

1 材料與方法

1.1 試驗樣品及初情期鑒定方法本研究選取9只大足黑山羊按照年齡平均分為3組。ET 組包括3只出生時的山羊(ET1、ET2和ET3),PT 組包括3只初情期山羊(PT1、PT2和PT3),C 組包括3只成年期山羊(AT1、AT2和AT3)。山羊麻醉后,立即收集睪丸,一部分送入深圳華大基因公司進行測序,一部分置入液氮中進行保存。按初情期的定義來確定公畜的初情期在實踐中較為困難。目前常采用體質(zhì)量、睪丸大小的評定來做估測。初情期的鑒定方法:初情期與體質(zhì)量的關(guān)系通常比年齡更為密切,當(dāng)公羊初次表現(xiàn)出求偶、勃起、爬跨、交配等行為時,對這些公羊進行稱重發(fā)現(xiàn)其體質(zhì)量為成年體質(zhì)量的40%。此外,對這些公羊的睪丸大小進行測量發(fā)現(xiàn),睪丸直徑約為6 cm。綜合以上特征,判斷其為處于初情期的公羊。

1.2 蘇木精伊紅染色采用4%甲醛溶液和常規(guī)組織學(xué)方法對3個不同發(fā)育階段的睪丸組織進行固定處理。之后移入二甲苯與酒精混合液處理5 min,依次浸入100%,95%,85%,70%酒精,分別處理2～5 min,蒸餾水轉(zhuǎn)入蘇木精染液進行染色5～15 min,蒸餾水洗滌后,轉(zhuǎn)入0.1%～0.5%伊紅染液染色1～5 min后經(jīng)脫水、脫色、封片后加蓋玻片封固,在顯微鏡下進行觀察。

1.3 RNA提取以及文庫的構(gòu)建通過TRIzol Plus RNA 純化試劑盒(Invitrogen,美國)提取總RNA。用NanoDrop-1000分光光度計測定總RNA 的濃度和純度。利用生物素標(biāo)記的特異性探針去除核糖體r RNA。純化之后,將RNA 在一定的溫度和離子環(huán)境下片段化。隨后使用TruSeq?標(biāo)準(zhǔn)試劑盒中隨機引物和反轉(zhuǎn)錄酶合成cDNA 第1鏈,之后使用DNA聚合酶I和核糖核酸酶H 來合成雙鏈的cDNA。雙鏈的cDNA 產(chǎn)物隨后進行了加“A”堿基和接頭連接。連接產(chǎn)物將被擴增,純化后即得到了最終的cDNA 文庫。最后將構(gòu)建好的測序文庫進行測序。

1.4 數(shù)據(jù)的過濾組裝以及編碼能力預(yù)測過濾掉含有接頭,大量低質(zhì)量堿基,以及含N 過多的讀數(shù),得到清潔讀數(shù)。之后使用比對軟件HISAT[5]將清潔讀數(shù)比對到山羊參考基因組上并用String Tie[6]進行組裝。在獲得每個樣品所有轉(zhuǎn)錄本的序列后用Cuffcompare[7]將這些轉(zhuǎn)錄本與已知的m RNA 及l(fā)ncRNA 進行比較,獲得它們相互位置關(guān)系的信息。并使用蛋白質(zhì)數(shù)據(jù)庫pfam[8]以及預(yù)測軟件CPC[9],txCdsPredict,CNCI[10]對轉(zhuǎn)錄本的編碼進行打分以區(qū)分m RNA 和lncRNA。匹配上pfam 的轉(zhuǎn)錄本則認(rèn)為是mRNA,否則為lncRNA。4種方法至少有3種方法結(jié)果一致,才會確定該轉(zhuǎn)錄本的類型。3款預(yù)測軟件的打分閾值分別為CPC threshold＝0,CNCI threshold＝0,txCdsPredict threshold＝500。

1.5 基因定量分析及組間差異分析使用RSEM[11]計算轉(zhuǎn)錄本的表達量并對基因的表達量進行FPKM 標(biāo)準(zhǔn)化處理。利用差異分析軟件DEGseq[12]進行組件差異分析。顯著性差異基因的過濾條件為:Fold Change≥2.000和FDR≤0.001。

1.6 IncRNA的靶基因以及GO 功能和KEGG 富集分析計算lnc RNA 與m RNA 的2個相關(guān)系數(shù),斯皮爾曼相關(guān)性系數(shù)與皮爾森相關(guān)性系數(shù),要求斯皮爾曼相關(guān)性系數(shù)≥0.6且皮爾森相關(guān)性系數(shù)≥0.6。lncRNA 在mRNA 上游10 000內(nèi)或在mRNA 下游20 000內(nèi)則判定為順式調(diào)控作用。超出這范圍的,將用RNAplex[13]分析lncRNA 與m RNA 的結(jié)合能,若結(jié)合能＜－30,則判定為反式調(diào)控。利用BLAST2GO[14]和Diamond[15]對lnc RNA 顯著性差異的靶基因進行GO 和KEGG 功能注釋,顯著富集的篩選條件:FDR≤0.01。

1.7 逆轉(zhuǎn)錄與定量聚合酶鏈反應(yīng)根據(jù)lnc RNAs序列設(shè)計引物避開與靶基因有重疊的區(qū)域,引物由上海生工生物公司合成。GAPDH 被用作內(nèi)參基因。我們使用TaKaRa PrimescriptTMRT 試劑盒從基因組中去除gDNA,并進行RNA 的逆轉(zhuǎn)錄。在Bio-Rad CFX96 實時PCR 檢測系統(tǒng)中使用TBGrenTMPremix ExTaqTMⅡ(Ta KaRa)進行反應(yīng),如下所示:95℃持續(xù)30 s,隨后40 個95℃循環(huán)持續(xù)5 s,60℃持續(xù)30 s。用2－ΔΔCt法測定lncRNAs的相對表達水平。每個選擇鑒定的lncRNAs均進行3次生物學(xué)重復(fù),試驗所用引物如表1所示。

1.8 統(tǒng)計分析所有qPCR 數(shù)據(jù)均表示為±s。使用單因素方差(one-way ANOVA)對所有數(shù)據(jù)進行分析,通過Levene檢驗數(shù)據(jù)方差的均勻性,然后進行最小顯著性差異(s x) 檢驗。

2 結(jié)果

2.1 不同發(fā)育階段睪丸組織形態(tài)我們將收集到的一部分睪丸組織用橫切和縱切的方式制作切片,然后用蘇木精染色,放置顯微鏡下進行形態(tài)學(xué)觀察(圖1)。我們觀察到這3組睪丸的形態(tài)學(xué)差異很大。在放大200倍的情況下,出生期睪丸的曲細(xì)精管直徑明顯小于初情期和成年期;曲細(xì)精管排列較近,空間內(nèi)結(jié)締組織較小;精原細(xì)胞數(shù)量稀少,未有精子細(xì)胞的產(chǎn)生。初情期內(nèi)精原細(xì)胞數(shù)量明顯增多,已有精母細(xì)胞和精子細(xì)胞的產(chǎn)生。在成年期可見大量精原細(xì)胞已經(jīng)發(fā)育為初級精母細(xì)胞和次級精母細(xì)胞,次級精母細(xì)胞又發(fā)育為精子細(xì)胞,部分精子細(xì)胞已逐漸分化為了精子。結(jié)果表明,3 組處在不同發(fā)育階段的睪丸在形態(tài)學(xué)上有顯著性差異,收集的3組睪丸與我們預(yù)期的發(fā)育時期相吻合。

圖1 不同發(fā)育階段睪丸的顯微形態(tài)觀察 A.出生期山羊睪丸組織橫切切片形態(tài)(200×);B.出生期山羊睪丸組織縱切切片形態(tài)(200×);C.初情期山羊睪丸組織橫切切片形態(tài)(200×);D.初情期山羊睪丸組織縱切切片在形態(tài)(200×);E.成年山羊睪丸組織橫切切片形態(tài)(200×);F.成年山羊睪丸組織縱切切片形態(tài)(200×)

2.2 測序數(shù)據(jù)概述ET、PT 和AT 組共計分別產(chǎn)生了654 200 616,658 526 698和720 412 344個原始讀數(shù),過濾掉含接頭和大量低質(zhì)量堿基以及含N過多的讀數(shù)后,分別產(chǎn)生了619 375 206(94.68%),620 322 606(94.20%)和688 471 996(95.57%)的清潔讀數(shù),這些清潔讀數(shù)將用于后續(xù)分析?？偣矊⒔?0%的清潔讀數(shù)被定位到山羊的參考基因組中,表明文庫質(zhì)量良好。

2.3 山羊睪丸中l(wèi)ncRNAs的識別及結(jié)構(gòu)特征我們從CPC,txCdsPredict,CNCI和pfam 分析結(jié)果的交叉點中鑒定產(chǎn)生了8 183個新發(fā)現(xiàn)的lncRNA(圖2)。同時,我們計算了這些lncRNAs轉(zhuǎn)錄本的平均長度為8 676.50 bp。然而,計算m RNA 的平均長度為10 376.40 bp。此外,lncRNAs外顯子的平均數(shù)為2.76,其中62.65%的轉(zhuǎn)錄本僅含有1個或2個外顯子。然而,47 193個蛋白質(zhì)編碼轉(zhuǎn)錄本中的外顯子平均數(shù)為11.99,遠(yuǎn)遠(yuǎn)大于lncRNAs的轉(zhuǎn)錄本。表明lnc RNA 編碼能力較弱,主要是從表觀修飾的層面上對機體進行調(diào)控(圖3)。

圖2 8 183個新的lncRNAs從CPC、txCdsPredict、CNCI和pfam 分析結(jié)果的交集中獲得

圖3 外顯子分布圖

2.4 各組間lncRNAs的差異分析通過定量分析,我們發(fā)現(xiàn)在預(yù)測的8 183個新的lncRNAs中,7 808個lncRNAs在樣本中表達。3組間lncRNAs的統(tǒng)計分析表明,AT 和ET 組相比共有7 508個lnc RNAs表達差異,其中5 465 個表達差異顯著。AT和PT 組相比共有7 513個不同表達的lncRNAs,其中5 392個表達差異顯著。PT 和ET 組相比共有5 593個lnc RNAs差異表達,其中只有1 167個表達差異顯著(圖4)。

2.5 qPCR驗證差異表達的lncRNAs我們從差異表達的數(shù)據(jù)中隨機選擇9個lnc RNAs,用qPCR 驗證其相對表達量。如圖5所示,大多數(shù)lncRNAs在體內(nèi)3個階段表達,表達趨勢以及差異顯著與我們的測序分析數(shù)據(jù)高度一致。

2.6 IncRNAs的靶基因預(yù)測IncRNA 的功能主要通過順式或反式作用于靶基因來實現(xiàn)。順式作用靶基因預(yù)測的基本原理認(rèn)為lncRNAs的功能與其臨近的編碼蛋白基因有關(guān)[16],對lnc RNAs 以及mRNA 的位置關(guān)系進行了詳細(xì)的分類并進行統(tǒng)計(圖6),計算lncRNAs與m RNA 的相關(guān)系數(shù)并判定在m RNA 上游10 000 內(nèi)或在mRNA 下游20 000 內(nèi)的lncRNAs。共有4 111個lnc RNAs參與到了順式調(diào)控作用,5 477個m RNA 被lnc RNAs順式調(diào)控。在此范圍之外,用RNAplex分析lncRNAs和mRNA 的結(jié)合能,最終發(fā)現(xiàn)只有445個lnc RNAs參與了反式作用。

2.7 差異表達lncRNAs靶基因的富集分析及與精子發(fā)生相關(guān)lncRNA的篩選與所有可能被調(diào)控的靶基因相比,我們更關(guān)注的是組件顯著性差異的靶基因功能。為了探索這些靶基因的功能,我們進行了GO 功能和KEGG 通路富集分析。在GO 功能分析中,精子發(fā)生、精細(xì)胞發(fā)育、精細(xì)胞分化等功能模塊顯著富集(圖7),提示這些差異表達的靶基因在雄性生殖中可能扮演著重要角色。此外對于KEGG 通路的富集分析發(fā)現(xiàn)靶基因主要富集于MAPK 信號通路,HTLV-I信號通路,醛固酮的合成分泌、酮體信號通路,其中MAPK 信號通路,醛固酮的合成分泌、酮體信號通路均與精子發(fā)生密切相關(guān)。57個候選m RNA 富集在精子發(fā)生功能中共靶向到了50個新發(fā)現(xiàn)的lnc RNAs(表2),表明這50個lncRNAs可能是調(diào)控精子發(fā)生的關(guān)鍵基因。

2.8 與精子發(fā)生相關(guān)lncRNAs的初步鑒定從篩選出的50個與精子發(fā)生相關(guān)的lncRNAs選擇出3個進行qPCR定量分析(圖8),結(jié)合mRNA 測序與已知的研究結(jié)果,對這些lncRNAs及其靶基因進行關(guān)聯(lián)分析。選出的3 個lncRNAs LTCONS_00002707、LTCONS_00052201、LTCONS_00072416在成年期的表達量比其他時期極顯著增加,其表達趨勢與其靶基因和測序數(shù)據(jù)基本一致,表明在出生期和初情期精子發(fā)生相關(guān)的基因表達量比較低,在成年期性成熟時則顯著表達。推測這3個lncRNAs有可能是通過促進精子發(fā)生相關(guān)的mRNA的表達而發(fā)揮作用。

圖4 不同組中差異表達的lncRNAs和m RNA 數(shù)量 A.過濾前差異表達基因數(shù)量;B.過濾后差異表達基因數(shù)量;1.AT-VSET;2.AT-VS-PT;3.PT-VS-ET

表1 試驗所用引物

圖5 qPCR 驗證9個差異表達的lncRNAs ?.P＜0.05;??.P＜0.01。下同

圖6 lncRNAs和靶基因的位置分布

圖7 靶基因的GO 功能富集圖譜

3 討論

在山羊睪丸的研究中,更多的研究集中在蛋白質(zhì)編碼RNA 和miRNA 上,而lncRNAs作為一種重要的調(diào)節(jié)因子鮮有研究進行測序發(fā)掘。在本研究中,我們利用深度測序方法分析了不同發(fā)育階段山羊睪丸中l(wèi)ncRNAs的表達,通過靶基因預(yù)測以及組件差異分析尋找出組件差異顯著的lncRNAs及其靶基因,將這些顯著性差異的靶基因通過GO 功能與KEGG 富集分析以篩選出與精子發(fā)生相關(guān)的lncRNAs。

表2 篩選出的部分與精子發(fā)生有關(guān)的lncRNAs及靶基因

圖8 與精子發(fā)生相關(guān)lncRNAs的qPCR 驗證

本研究共獲得了12.68 Gb的數(shù)據(jù),90%的基因映射到了山羊參考基因組上。發(fā)掘出新的lncRNAs 8 183個,并篩選獲得了50個新發(fā)現(xiàn)與精子發(fā)生有關(guān)的lncRNAs及對應(yīng)的57個候選靶基因,大多數(shù)候選靶基因都位于lnc RNAs附近,通過順式作用被lncRNA 所調(diào)控,提示順式調(diào)控可能是lncRNAs對靶基因的主要調(diào)控方式,這與之前報道的一致[16]。這些候選靶基因在精子的形成過程中扮演著重要角色。例如tdtp(睪丸發(fā)育相關(guān)蛋白)之前發(fā)現(xiàn)在睪丸精原細(xì)胞中表達量很高,與生殖細(xì)胞的減數(shù)分裂有關(guān)。當(dāng)tdtp缺失時將會導(dǎo)致精子發(fā)生功能停止,此外臨床研究也表明,在患有精子生成障礙癥的男性中,tdtp基因幾乎不表達[17]。zpbp2(透明帶結(jié)合蛋白)是在精子頂體中特異表達并高度保守的蛋白。之前的研究認(rèn)為zpbp2 是透明帶在精子上的重要受體,主要在精子與卵子的相互作用中行使功能,最近的研究表明zpbp2在精子發(fā)生過程中對精子頂體的形成起重要作用[18]。敲除zpbp2會導(dǎo)致生育力下降和精子頭部細(xì)微畸形,這與zpbp2沿精子吻脊的離散定位相對應(yīng)。odf2(外周致密纖維2)是精子尾部的促成成分,是維持精子結(jié)構(gòu)的骨架蛋白,在精子運動過程中可以保護鞭毛不受剪切力的作用從而維持精子的運動[19]。此外,tssk家族也是我們篩選出的候選靶基因,它可以動態(tài)磷酸化odf2來維持精子的完整性與odf2協(xié)同調(diào)節(jié)精子結(jié)構(gòu)和運動[20]。

從篩選出的lnc RNAs中經(jīng)qPCR 驗證發(fā)現(xiàn)它們的表達趨勢與靶基因一致,即在出生期和初情期表達量較低,在成年性成熟期表達量顯著增高,推測這些lnc RNAs可能促進它們靶基因的表達,即通過促進精子發(fā)生相關(guān)的m RNA 的表達發(fā)揮作用。本研究的睪丸組織染色結(jié)果也表示在出生期和初情期幾乎未發(fā)現(xiàn)明顯的精子形態(tài),而在成年期觀察到了大量精子的產(chǎn)生,表明這些lncRNAs可能在精子發(fā)生中扮演著重要角色。

總之,本研究的結(jié)果極大地縮小了精子發(fā)生相關(guān)lnc RNAs的研究范圍,為進一步了解lnc RNAs在睪丸發(fā)育和精子發(fā)生中的作用提供了重要的數(shù)據(jù)支持。