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

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

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

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

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

1 材料與方法

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

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

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

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

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

1.6 IncRNA的靶基因以及GO 功能和KEGG 富集分析計(jì)算lnc RNA 與m RNA 的2個(gè)相關(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]對(duì)lnc RNA 顯著性差異的靶基因進(jìn)行GO 和KEGG 功能注釋,顯著富集的篩選條件:FDR≤0.01。

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

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

2 結(jié)果

2.1 不同發(fā)育階段睪丸組織形態(tài)我們將收集到的一部分睪丸組織用橫切和縱切的方式制作切片,然后用蘇木精染色,放置顯微鏡下進(jìn)行形態(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ā)育為初級(jí)精母細(xì)胞和次級(jí)精母細(xì)胞,次級(jí)精母細(xì)胞又發(fā)育為精子細(xì)胞,部分精子細(xì)胞已逐漸分化為了精子。結(jié)果表明,3 組處在不同發(fā)育階段的睪丸在形態(tài)學(xué)上有顯著性差異,收集的3組睪丸與我們預(yù)期的發(fā)育時(shí)期相吻合。

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

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

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

圖3 外顯子分布圖

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

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

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

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

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

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

表1 試驗(yàn)所用引物

圖5 qPCR 驗(yàn)證9個(gè)差異表達(dá)的lncRNAs ?.P＜0.05;??.P＜0.01。下同

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

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

3 討論

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

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

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

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

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

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