基于轉錄組分析不同溫度下池蝶蚌性別相關基因的表達

曾 起，胡蓓娟，徐紅琴，洪一江*

(南昌大學a.生命科學學院;b.江西省水產動物資源與利用重點實驗室，南昌 330031)

池蝶蚌是中國主要的淡水養殖珍珠蚌之一，因育珠能力強，產出珍珠具備粒徑大，光澤度好和品相佳等特點，養殖規模迅速擴大，成為了我國重要的經濟貝類，應用前景非常巨大[1]。內臟團作為培育珍珠的重要部位，其產出的珍珠光澤好，顆粒大和形狀圓[2]，其中性腺是很重要的因素。

蔡英亞等[3]闡明，產在丹麥的食用牡蠣的性別比例會隨著水溫變化而發生改變，當生存環境的水溫高時，其雌性個體數量處于優勢地位，反之，雄性個體數量占優。金啟增等[4]闡述馬氏珠母貝在培養過程中水溫升高會導致雌貝的數量占比增加；反之，雄貝占比增加。據觀察，水溫與貝類性別的轉變有著很重要的關系，當水平處于13.2 ℃～19.9 ℃時，雄性個體占比高；當水溫在20 ℃～29.3 ℃時，雌性個體和雄性個體占比幾乎相等；但水溫再次降低時，雄性個體的數目又會重新增加。這表明，水溫低，雄性貝數量上是占優勢的。Zapata等人[5]研究發現，歐洲牡蠣配子的發生和性別決定主要受到溫度影響，類固醇不會作為內源性調節因子參與性別發育，而海水溫度的上升和氣候的變暖，可能導致其自然種群中的性別比例傾斜。

最近的研究表明，某些基因會影響牡蠣或淡水貽貝的性別決定和分化。Teaniniuraitemoana等[6]采集了不同發育階段的雌性珠母貝和雄性珠母貝樣本，并進行了轉錄組分析，揭示了其性別決定和性別分化基因，例如影響雄性分化的dmrt和fem1，以及影響雌性發育的foxl2和vit。Zhou等[7]利用了WCGNA(加權基因共表達網絡分析)方法，篩選了扇貝性別分化和決定的關鍵基因，并構建扇貝性別確定和分化的途徑。在這一途徑中，dmrt1具有主導功能。Wang等[8]報道foxl2在三角帆蚌的性別分化中起關鍵作用，而β-catenin[9]在免疫反應和性別決定中起關鍵作用。此外，Hashiyama[10]僅在XX型果蠅性腺的原始生殖細胞中檢測到sxl的表達。

貝類性腺發育過程中，發生了一定程度的性逆轉或者雌性同體的出現。例如，一些珠母貝在受到溫度或者其他因素影響時，在生存過程中表現出連續性的性別逆轉[11]；Christelle[12]通過不同溫度的養殖，發現長牡蠣存在性別逆轉和雌雄同體現象，并伴隨著性別決定相關基因表達的變化；Wu等[13]在研究池蝶蚌發育過程中，發現26～32月齡存在雌性同體現象。因此，為了探究溫度對其性別分化和和性腺發育的影響，本研究采用不同養殖溫度下的池蝶蚌為原材料，并采用Pacbio三代測序技術和Illumina二代測序技術對雌雄蚌進行轉錄組測序，通過差異基因探究池蝶蚌性腺發育和性別分化的原因。

1 材料與方法

1.1 實驗材料

池蝶蚌取自于江西省撫州市池蝶蚌良種場。挑選無插核、健康的29月齡池蝶蚌共50只,放置在充分曝氣的養殖箱進行暫養。1周后，各取雌雄蚌5只作為一組，共分5組，分別在15 ℃，20 ℃，25 ℃，30 ℃和33 ℃下進行為期2周的養殖試驗。養殖結束后，每個溫度組各取3只雌蚌和雄蚌性腺混樣，進行二代RNA測序；另取常溫下(23 ℃)暫養的雌雄蚌性腺作為對照進行三代全長轉錄組測序，性腺樣品液氮速凍，隨后放入-80 ℃冰箱保存備用。

1.2 實驗方法

1.2.1 文庫的構建和測序

總RNA的提取，cDNA文庫的構建，轉錄組測序工作由北京諾禾致源科技股份有限公司完成。二代轉錄組測序平臺為Illunima Hiseq 2500，三代全長轉錄組測序平臺為Pacbio Sequel。

1.2.3 基因表達水平分析

利用CD-HIT軟件將雌雄蚌轉錄本合并去冗余后，得到的轉錄本作為參考序列，使用bowtie2軟件將clean data比對到全長轉錄組，并通過RSEM軟件獲得單個樣品中每個轉錄本的reads數量。考慮到序列深度和基因長度對片段的影響，并將所有reads計數標準化為FPKM，以將FPKM大于0.3的相對水平定義為顯著表達基因，并用其做基因表達水平統計。FPKM計算公式如下：

注：Xi:單個轉錄本的讀取數;Li:轉錄本長度；N:轉錄本的讀取總數。

1.2.4 基因差異表達分析和聚類分析

實驗采用了DEseq2 R包分析有重復組的數據，對各個溫度下雌雄之間進行差異表達分析。其中將差異倍數(FoldChange)大于2及調整后的p值(padj)小于0.05作為標準，篩選符合條件的基因作為差異表達基因。FoldChange是兩個樣品間某一轉錄本表達量的比值，padj是通過對p值進行校正后得到的以不同實驗條件下的差異基因的FPKM值為表達水平，做層次聚類分析。

1.2.5 差異表達基因功能注釋

使用NT，NR，Swiss-Prot，KEGG和GO數據庫對差異轉錄本進行注釋。

1.2.6 差異基因功能富集分析

差異表達基因的功能富集分析主要包括GO富集分析和KEGG通路富集分析，其中GO富集分析和KEGG富集分析是通過clusterProfiler R包完成。

1.2.7 qRT-PCR定量驗證

為了驗證轉錄組數據，選擇了6個差異表達基因的mRNA進行qRT-PCR。qRT-PCR引物設計使用Oligo7.0工具完成，引物序列由北京擎科生物公司合成，具體引物見表1。β-actin作為內參基因，采用2-ΔΔCt相對定量法定量，每個樣本每個基因進行3次重復實驗。

表1 qRT-PCR選用基因及引物序列

2 結果

2.1 測序結果與序列分析

PacBio單分子熒光測序為環狀測序，測序過程的單分子產出的有效片段為Subreads，經過自我矯正和去冗余后，獲得雌雄全長轉錄本數目分別為96 867和111 520條，平均長度分別為4 095和3 876 bp，主要分布在2.5～3.5 kbp之間(圖1)。

2.2 轉錄本的表達定量

將clean data的reads比對到全長轉錄本上發現，雌性性腺轉錄本比對上的序列在88%左右，雄性性腺轉錄本比對上的序列在75%左右，雌性轉錄本比對上數目多于雄性轉錄本。用FPKM統計不同表達水平下基因的數量以及單個基因的表達水平，發現在所有轉錄組中，FPKM分布最多的在0～1之間，占了50%左右，其中FPKM>5的大約占25%圖2(A)。為了直觀的比較各樣本整體的表達水平，以及對每個處理溫度組表達水平的離散程度，采用箱線圖展示FPKM的分布顯示10個組樣本表達水平整體上一致，其中雌性25 ℃處理組樣本表達水平略低，雌性30 ℃處理組表達水平略高圖2(B)。

轉錄本長度

轉錄本長度圖1 轉錄本長度分布圖

2.3 差異表達基因統計

采用DEseq2分析鑒定5個溫度下不同性別處理組間的差異表達基因，在15 ℃雌雄比較組至33 ℃雌雄比較組中，分別檢測到6 311，3 013，6 673，6 798，5 685個差異表達基因(表2)，結果發現，雌雄之間的差異基因非常豐富。在20 ℃和33 ℃時差異基因數量比其它組更少，15 ℃時雌雄對比組含3 013個差異表達基因，其中1 045個基因上調，1 968個基因下調；33 ℃對比組含5 703個差異表達基因，2426個上調，3277下調。在5個處理組中，共有17696條轉錄本在所有組間表達。用層次聚類分析5個溫度脅迫處理組間差異表達基因的分布情況發現，雌蚌組與雄蚌組具有類似的表達情況，其中15 ℃和20 ℃聚為一支，25 ℃，30 ℃和33 ℃聚為一支(圖3)。

(A) 不同表達水平區間的基因數量統計圖

(B) 基因表達水平對比圖圖2

2.4 不同溫度下差異基因的KEGG富集

差異基因的富集分析是通過超幾何分布，其富集分析以通路為單位，將轉錄本中的基因作為背景，在差異表達基因中顯著性富集的通路。在所有富集中，均展現富集程度最高的前20個通路。對5個不同溫度下雌雄間所鑒定到的6311，3013，6673，6798和5685個差異表達基因進行KEGG注釋，分別完成4585，2203，4896，4932和4069個。其中在15 ℃組，腎素分泌(ko04924)，甘氨酸、絲氨酸和蘇氨酸代謝(ko00260)，淀粉和蔗糖代謝(ko00500)顯著富集等；在20 ℃組中，顯著富集的有卵母細胞減數分裂(ko04114)，MAPK信號通路(ko04016)，醛固酮合成與分泌(ko04925)和GnRH信號傳導途徑(ko04912)等，腎素分泌(ko04924)及磷脂酰肌醇信號系統(ko04070)也被顯著富集；在25 ℃情況下，主要富集到卵母細胞減數分裂(ko04114)，細胞周期(ko04110)和Notch信號通路和腎素分泌(ko04924)等；在30 ℃的時，富集到細胞周期(ko04110)，細胞周期-酵母(ko04111)和孕激素介導的卵母細胞成熟(ko04914)。33 ℃時，可以看出顯著富集的包括丙酮酸代謝(ko00620)，腎素分泌(ko04924)，胰高血糖素信號通路(ko04922)。KEGG富集顯示在各溫度的處理組均含有卵母細胞分裂，但是在20 ℃和25 ℃的時呈現出顯著富集(qvalue<0.05)；在30 ℃出現了孕激素介導的卵母細胞成熟；在25 ℃的時候開始具有細胞周期。

表2 相同溫度下雌雄性腺之間的差異基因結果

注：聚類分析圖，將log10(FPKM+1)值進行歸一化轉換并進行聚類，深色表示高表達，淺色表示低表達。顏色從深到淺，表示log10(FPKM+1)從大到小。

2.5 差異基因中性別相關基因的篩選和分析

對5個溫度處理組之間的差異表達基因的進行組間比較，維恩圖有802個差異表達基因在5個對比組中重疊(見圖5)。其中共上調的基因為319個數，共下調的基因個數為483個。將5個組共表達差異的差異基因進行GO富集分析發現，下調基因的GO富集顯示在幾丁質的結合，幾丁質代謝過程和蛋白結合，鈣離子的結合，細胞外區域，膜，胞外空間上，這些均與與珍珠質的分泌和形成有關；上調基因的GO富集顯示在蛋白質結合，ATP的結合和酶的活動等，這些與精子的生成和活動有關(見圖6)。

注：A：15 ℃；B：20 ℃；C：25 ℃；D：30 ℃；E：33 ℃。圖4 各個溫度下差異基因KEGG富集結果

2.6 差異基因中性別相關基因的篩選和分析

在差異表達基因中篩選出大量卵巢發生相關基因，其中包括中蛋黃鐵蛋白(FRIY)，卵黃囊蛋白(SIAE)，核黃素結合蛋白(RBP)，卵黃囊衍生的胚胎大分子(ELYS)等，還有卵細胞發生相關基因，包括卵母細胞表達核仁蛋白(ANO39)，卵母細胞成熟因子(MOS)，卵母細胞特異性母體影響因子(ZAR1)和卵黃原蛋白(VIT)等。

在差異表達基因中鑒別到大量精子發生相關基因，精子發生相關蛋白4，7，17，22(spat4，spat7，spt17，spt22)，精子發生和中心粒相關蛋白1(sper1)和精子特異性接頭組蛋白H1樣蛋白(HILS1)；精子運動相關蛋白精子鞭毛蛋白1、2(spef1，spef2)，精子表面蛋白17(sp17)等；還發現大量精巢相關基因，例如，睪丸特異性絲氨酸/蘇氨酸蛋白激酶2(TSSK2)，T復合物睪丸表達蛋白1(Tcte-1)和精巢表達蛋白10，11，30(Tex10，Tex11，Tex30)等。

注：m15 vs f15:在15 ℃下雄雌間差異基因的數量，其他類似。

注：A:上調差異表達基因的GO富集；B:下調差異表達基因的GO富集。

為了進一步探究性別決定和性別分化，從各個組間中挑選出差異基因，列出出部分性別決定和性別分化相關基因如表4.7，其中sox2，wnt4，fem1，rspo1，vit在5個溫度對比組均顯著表達差異，但是sox5，sox9和folx2在20 ℃時沒有表達差異。其中tra1,wnt4，fem1，rspo1，sox9和folx2等，文獻記載它們與性別決定和性別分化有關。

表3 差異基因中性別決定基因的篩選

2.7 qRT-PCR驗證分析

隨機挑選了6個差異表達基因pif，tex43，wnt4，rsop1，tra1和fem1，進行qRT-PCR驗證。其中wnt4，rspo1，pif，calm在雌性中顯著表達(P<0.05，log2FC<-1)，tra1和fem1在雄性中顯著表達(P<0.05，log2FC>1)。通過qRT-PCT進行驗證，篩選的6個基因雌雄間表達差異極顯著(P<0.01)，且差異表達基因與轉錄組測序結果中的表達量變化倍數趨勢基本一致(圖7)。結果表明，轉錄組結果可靠，利用RNA-Seq技術對池蝶蚌雌雄株進行轉錄組測序以及篩選與性別相關的差異基因是可行的。根據熒光定量結果，fem1,tra1,rspo1和wnt4的mRNA表達水平如圖8，在雌蚌中fem1在15 ℃和20 ℃時高于其他水平；在雌蚌和雄蚌中rspo1和tra1在15 ℃和20 ℃時均低于其他溫度水平；在雌蚌中wnt4在25 ℃時表達量最高。

3 討論

SMRT測序技術可以提升無參轉錄組的組裝效果[14]。在這項研究中，我們首次采用RNA-seq和SMRT測序結合的方式，生成了池蝶蚌全長轉錄組。經分析后獲得了96 867和111 520個高質量的全長序列。在這項研究中，我們首次采用RNA-seq和SMRT測序結合的方式，生成了池蝶蚌全長轉錄組。經分析后獲得了96867和111520個高質量的全長序列，為無需PCR的基因結構探索和遺傳功能探究提供了寶貴數據。與之前使用Illumina平臺的池蝶蚌轉錄組研究相比[15]，二代的雌雄轉錄本的平均長度分別為484和472 bp，而本實驗的轉錄本的長達4 095和3 876 bp。研究調查了5個不同溫度下池蝶蚌雌雄性腺的基因表達。對差異表達基因研究發現，在20 ℃和33 ℃的時候，差異基因數量最少，推測這兩個溫度下雌雄表達差異比其他組小。各雌雄對比組只有802個基因在各個溫度下以不同方式共同表達。KEGG富集分析發現，25 ℃時出現細胞周期通路，說明細胞活動相對開始活躍；同時雌性個體生殖相關的活動變得活躍起來，在25 ℃時出現卵母細胞減速分裂通路和notch信號通路等。KEGG表明25 ℃時，雌性相關活動更加活躍，有利于雌性生長發育。

注：Log2FC:雄雌之間的差異倍數。圖7 qRT-PCR定量結果

注：fem1，tra1基因在15 ℃下雌蚌性腺中的表達量分別定義1;rspo1，wnt4基因在15 ℃下雄蚌性腺中的表達量分別定義1;**表示P<0.01。

轉錄組分析結果表明，wnt4，rspo1，fem1和tra1等4個性別決定或性別分化基因在不同溫度下的池蝶蚌性腺中始終出現差異。wnt4和rsop1的表達量雌性高于雄性，fem1和tra1的表達量雄性高于雌性。wnt4是性別決定基因，且在雌雄哺乳動物的形態發育中起到關鍵的作用[16]。小鼠性別決定前的性腺中wnt4存在表達，性別分化后在卵巢中仍持續表達，但精巢中的表達顯著下降[17]。李海龍等人[18]研究發現wnt4在櫛孔扇貝的性腺發育的各個時期均存在表達，尤其成熟性腺表達最高，推測其在性腺生殖細胞成熟的整個過程中起著不可忽視的作用。Rspo1在雌性性別決定通路上游發揮作用，研究發現能啟動wnt/β-catenin通路，對雌性卵巢的分化是不可或缺的，另外，rspo1能促進卵巢發育并抑制精巢發育[19]。在青鳉中研究發現rspo1基因激活的通路，對雌性卵巢的分化是不可或缺的[20]，rspo1在青鳉雄性胚胎中過表達能夠引起完全的性逆轉(雄轉雌)，證明rspo1通路在青鳉雌性卵巢發育中發揮著關鍵的作用。在線蟲中，fem1是雄性個體發育，以及雄性和雌雄同體精子發生所必須的[21]，它可以通過阻止tra1的作用，從促進雄性發育產生。Doniach等[22]發現在櫛孔扇貝中，fem1只在成熟的雄性個體內表達，而不是雌性或者雌雄同體內，表明fem1在維持雄性化可能發揮作用。周祖陽等[23]在長牡蠣中成功克隆出fem1，并對長牡蠣各個發育時期的fem1進行表達分析，推測其參與了性別決定和性別分化的過程。fem1，rspo1和wnt4的qRT-PCR與轉錄測序的表達量變化倍數一致，且rspo1和wnt4在雌性中的表達水平高于雄性，而fem1在雄性中的表達水平高于雌性。推測fem1，rspo1和wnt4在性別調控中發揮作用，rspo1和wnt4在維持雌性化發生作用，而fem1在維持雄性化發揮作用。

Berkseth等[24]表明tra1在發育過程中通過參與阻止與雄性發育的基因來促進雌性發育方面有普遍作用。Ronald[25]總結出tra1可以被視為是雌性動物中的一種阻遏物，它可以關閉雄性相關基因然后發育成雌性。在此次研究中，通過差異基因分析發現，tra1在雄性中顯著表達，而且熒光定量結果分析也顯示tra1表達量在雄性中高于雌性。Wu等[13]發現池蝶蚌在26～32月齡時，出現雌雄同體的現象，且大部分的雌性同體蚌從雄蚌中發育而來。此次研究的池蝶蚌處于30月齡，推測tra1的出現，可能關閉池蝶蚌雄性相關基因，進而產生雌雄同體的現象。熒光定量結果顯示，25 ℃以上時，tra1的表達量高于10 ℃和15 ℃水平，溫度的升高，可能進一步促進池蝶蚌雌雄同體的形成。

Holleley等[26]發現當孵化溫度在22 ℃～32 ℃之間時，鬃獅蜥由性染色體決定性別；當溫度超過32 ℃時，后代中雌性的比例越來越高。Weber C和Zhou等[27]在巴西紅耳龜中發現，早期的原始性腺在31 ℃時分化為卵巢，在26 ℃時分化為精巢。在具有溫度依賴性性別決定的紅耳龜中，rspo1的上調發生在溫度敏感期，即性腺的發育受到某個特定溫度的影響。因此，rspo1表達是對溫度敏感的，并且在處于胚胎從雌性轉變為雄性的孵化溫度時，rspo1基因就會下調[28]。熒光定量結果表明，rspo1在15 ℃和20 ℃時，雌蚌和雄蚌中表達量較25 ℃，30 ℃和33 ℃下顯著(P<0.01)降低，而推測在15 ℃和20 ℃溫度下有利于雄性的發育；雌蚌中的wnt4在25 ℃時表達量最高，而在雄蚌中的wnt4隨著溫度升高而升高；雌蚌中的fem1在25 ℃之后表達量出現顯著(P<0.01)降低，而雄蚌中的fem1受溫度影響不大。

綜上，根據KEGG富集分析和rspo1，wnt4，tra1和fem1的表達情況判斷，推測25 ℃以上的溫度，適合雌性的生殖發育。本研究為探索池蝶蚌的性別分化和性腺發育做出一定的參考，關于溫度對池蝶蚌性腺發育和性別分化的影響，還涉及諸多基因與其產物之間的相互作用，依然有待進一步的深入發掘和探索。