中華鱉性別差異表達基因比較研究

周先文，王曉清，王 佩，江 輝，曾 丹，彭英海

(1.湖南農業大學，湖南 長沙 410128；2.湘西州畜牧水產事務中心，湖南 吉首 416000)

中華鱉(Trionyxsinensis)是龜鱉目龜鱉科的代表動物，為中國的特色養殖品種，養殖產量位居世界第一，近年來我國學者對中華鱉的生長[1]、營養[2-4]、病害[5-6]和育種[7]等方面開展了大量的研究。中華鱉雄性生長普遍快于雌性，其雄性裙邊厚，市場價格高于雌性，實行中華鱉雄性單性養殖能獲得更好的經濟效益[8]。根據水產動物的特性有針對性地進行育種研究，有利于提高相關產業的經濟效益。在水產動物的性別研究中，硬骨魚類的性成熟受腦-垂體-性腺軸的調節，下丘腦通過分泌激素到垂體促性腺細胞上，促進卵泡激素等的合成與釋放，對性腺發育和配子形成等產生作用[9-10]。近來的研究表明，Dmrt1和Amh是中華鱉雄性的重要調控因子[8，11]，Cyp19a1和Foxl2為中華鱉卵巢的調節因子。然而，中華鱉的性別決定機制尚未闡明，中華鱉在自然生產過程中3齡達到性成熟，此階段為性別相關基因高度表達階段。通過差異表達基因的研究，可進一步驗證中華鱉性別差異候選基因，探索其對性別決定的影響，對明晰中華鱉的性別決定機制有著十分重要的意義。

1 材料與方法

1.1 試驗材料

2017年9月25日從常德市鼎城區河州龜鱉養殖專業合作社購買3齡中華鱉雌雄各6只，體重994.05～1 286.48 g。分別取中華鱉雌性的卵巢和雄性的精巢，每個個體取2份性腺組織樣品放置于凍存管中，液氮速凍后置-80℃冰箱保存。

1.2 RNA提取及cDNA制備

解剖獲得中華鱉精巢和卵巢組織，每個樣品取2份，其中一份用于轉錄組分析，另一份進行RNA提取用于熒光定量PCR驗證試驗，RNA提取采用RNA提取試劑盒(OMGA，廣州)，購自廣州飛揚生物工程有限公司。cDNA制備采用賽默飛的cDNA試劑盒，按試劑盒說明合成cDNA。

1.3 文庫構建及測序

文庫構建及測序委托武漢菲沙基因信息有限公司完成，樣品檢測合格后，通過試劑盒去除rRNA富集mRNA。加入fragementation buffer將mRNA打斷，以mRNA為模板，用隨機引物合成一鏈cDNA，然后加入dNTPs等合成cDNA，再用AMpure XP beads進行片段大小選擇、擴增及純化PCR產物，最后獲得文庫；對文庫采取Qubit 2.0初步定量，然后用Agilent 2100對插入片段長度進行檢測。庫建合格后，利用Illumina Hiseq X平臺進行HiSeq雙端測序。

1.4 生物信息學分析

數據下機得到raw data，經過數據處理，去接頭信息，低質量堿基等得到clean data，將clean data與中華鱉參考基因組比對，進行基因注釋，注釋后的基因作表達分析、樣品相關性分析、差異表達分析、GO富集、KEGG富集、轉錄因子分析、蛋白質互作網絡分析等。轉錄組數據與參考基因組序列比對利用TopHat 2[12]，調用Bowtie 2[13]將每個樣本質控后的二代序列比對到中華鱉的參考基因組序列，中華鱉的參考基因組下載地址為https：//www.ncbi.nlm.nih.gov/genome/14578?genome_assembly_id=44853。

表達水平分析利用RSEM[14]，調用Bowtie 2進行比對，將樣品比對到每個轉錄本上的Reads數目，進行FPKM(Fragments per kilobase per Million bases)轉換，最后，獲得基因和轉錄本的表達水平。基因差異表達分析的輸入數據為基因表達水平分析中得到的expected_count數據，用R語言包DEseq 2[15]進行差異分析，采用FDR(False discovery rate)和log2FC(Fold change)進行篩選，其中上調表達的，滿足條件FDR<0.05、log2FC>1；下調的差異表達的基因數目，滿足條件FDR<0.05、log2FC<-1。差異表達基因的Gene Ontology[16](簡稱GO，http：//www.geneontology.org/)是基因功能國際標準分類體系。分析獲得的差異基因在Gene Ontology中的分布狀況及其在中華鱉性腺發育中的作用，GO采用超幾何分析，選取FDR≤0.05的GO term作為顯著富集的GO條目。同時對差異表達基因進行了KEGG分析，差異表達基因的KEGG[17-18](Kyoto Encyclopedia of Genes and Genomes)是有關Pathway的公共數據庫，Pathway顯著性富集分析以KEGG Pathway為單位，以獲得差異基因相對于所有注釋基因顯著富集的通路。為分析差異基因蛋白互作情況，研究通過STRING蛋白質互作數據庫(http：//string-db.org/)進行蛋白互作網絡的構建，對差異表達蛋白以Cytoscape 3.6.1軟件進行可視化編輯。

1.5 RT-qPCR驗證

選取12個差異表達基因進行驗證，用Primer 6.0設計引物，以GAPDH為內參基因，采用TB GreenTMPre mix EXTapTM熒光定量表達試劑盒，采用2-ΔΔCt計算相對表達量，每個樣品設置3表重復。

2 結果與分析

2.1 測序數據結果與分析

對3齡中華鱉精巢和卵巢分別進行轉錄組測序，得到精巢和卵巢clean data各為21 494 396 400 bp和24 605 142 600 bp。通常以過濾后得到的clean data的百分比來說明數據的可靠性和精確性；以Q30作為堿基質量檢測的標準，Q30>85%則表示可以進行后續的數據分析。3齡中華鱉精巢組織Q30≥92.7%，3齡中華鱉卵巢組織Q30≥90.3%(表1)。

表1 使用Illumina Hiseq X平臺測序過濾后的數據結果

2.2 測序數據與參考基因組比對結果

Mapping rate是指映射clean reads的比率，是轉錄組測序數據使用率最直接的表示(表2)。在此次的比對結果中，3齡中華鱉精巢和卵巢組織的Reads與參考基因組的比對效率在53.7%～67.5%之間。通過分析比對效率的結果，可以確定所選的中華鱉參考基因組的組裝滿足信息分析需求。

表2 Reads與參考基因組的比對情況

2.3 差異基因的表達水平分析

通過比較分析3齡中華鱉精巢和卵巢組織的差異表達基因，總的差異基因10 145個，其中在雌性卵巢中表達上調4 034個、表達下調6 111個。為直觀顯示兩組樣本間FDR和差異倍數FC的分布情況，繪制了MA plot和volcano plot(圖1)。

2.4 差異基因的GO分類及富集分析

以生物學過程、細胞組分和分子功能為主要內容的GO數據庫對轉錄組數據進行了系統分類與功能注釋，總共有10 145個差異基因注釋到10 291個GO分支中，每個類別列出10個以上注釋最多的子類別。

在GO-生物學過程中注釋到細胞過程(Cellular process)、單一生物過程(Single-organism process)、代謝過程(Metabolic process)、生物調節(Biological regulation)、生物過程調節(Regulation of biological process)的基因個數依次是1 999、1 701、1 486、1 424和1 358個；在GO-細胞組分數據庫中，注釋到細胞(Cell)和細胞組分(Cell part)的均為1 945個，注釋到細胞器(Organelle)為 1 593個，注釋到細胞膜(Membrane)的有1 206個；在GO-分子功能數據庫中，1 722個基因注釋到連接(Binding)，951個基因注釋到催化活性(Catalytic activity)，以上通路在各類別中注釋的基因較多(圖2)。

為進一步分析中華鱉雌雄性別差異方面的功能，對生殖過程，配子產生、性腺分化等通路進行分析，其中注釋到生殖過程(Reproductive process)的有105個；注釋到單一有機體生殖過程(Single organism reproductive process)的有97個；注釋到多有機體生殖過程(Multi-organism reproductive process)76個；注釋到有性生殖(Sexual reproduction)的有71個；注釋到胚胎發育結束于出生與卵孵化(Embryo development ending in birth or egg hatching)的有69個；注釋到多細胞生物的生殖過程(Multicellular organismal reproductive process)的有64個；注釋到配子產生(Gamete generation)的有56個；注釋到性別分化(Sex differentiation)的有21個；注釋到性腺發育(Gonad development)的有18個等。將上述通路篩選到差異基因進行匯總，得到122個精巢和卵巢差異表達的基因(表3)。

表3 中華鱉性別差異表達基因GO注釋結果

2.5 差異表達基因KEGG注釋及富集分析

為鑒別差異基因在中華鱉性腺中的具體功能，通過KEGG數據庫對富集的信號通路進行分析，在中華鱉雌雄差異基因表達研究時發現，富集通路最多的是代謝和有機體系的通路；而在KEGG富集基因通路中，富集基因數目上依次是信號轉導(Signal transduction)、免疫系統(Immune system)、內分泌系統(Endocrine system)、細胞過程中分類運輸和分解代謝(Transport and catabolism)、真核細胞社區(Cellular community-eukaryotes)，富集的基因分別為641、317、296、305、244(圖3)。這為之后中華鱉性別差異研究和性別決定基因篩選提供了參考。

此外，實驗通過富集最顯著的20條通路，繪制了差異表達基因的KEGG富集散點圖(圖4)，可以對差異表達基因的途徑進行整體分析。圖中Rich factor指通路中的富集程度，Q-value是校正之后的P-value，數值越小表示富集越顯著。在3齡中華鱉精巢與卵巢轉錄組差異的研究中，綜合分析通路富集程度大及富集顯著的主要有多能干細胞調節信號通路(Signaling pathways regulating pluripotency of stem cells)、催產素信號通路(Oxytocin signaling pathway)、精氨酸生物合成信號通路(Arginine biosynthesis)、磷酸鹽和磷酸鹽代謝(Phosphonate and phosphinate metabolism)、磷脂肌醇信號通路(Phosphatidy-linositol signaling system)和MAPK信號通路(MAPK signaling pathway)等。這些富集了不同差異表達基因的通路，影響中華鱉機體的個體機能，其中催產素信號通路(Oxytocin signaling pathway)和MAPK信號通路(MAPK signaling pathway)等通路含有部分性別相關基因，對中華鱉性別維持等方面具有重要的作用。

為進一步鑒定差異基因在精巢和卵巢中的功能。本研究將10 145個差異基因注釋到283個信號通路進行分析。差異基因注釋比例最多的信號通路是PI3K-Akt信號通路(上調85個、下調48個)，其次一些包括性別相關的信號通路如：Wnt信號通路(上調37個、下調14個)，催產素信號通路(上調31個、下調32個)，Foxo信號通路(上調34個、下調13個)，卵巢類固醇生成素信號通路(上調14個、下調15個)等差異基因也較多。

2.6 蛋白-蛋白互作網絡

將差異蛋白結果數據文件導入到Cytoscape軟件進行可視化編輯，構建了成熟期中華鱉精巢和卵巢組織差異蛋白互作網絡及其子網絡，發現Foxl2與中華鱉性別決定相關，該蛋白的子網絡蛋白互作網絡節點較為繁簡適中(圖5)；而在雄性蛋白互作子網絡中，SMC通路在雄性相關蛋白的子網絡中起著調控作用(圖6)。

注：紅色代表在卵巢中上調，藍色代表在卵巢中下調，顏色越深表示上調/下調越明顯；PO(positive)：代表兩個蛋白呈正相關；NE(negative)：兩個蛋白之間呈負相關。

2.7 熒光定量PCR(RT-qPCR)驗證

為驗證轉錄組數據的可靠性，隨機選取了12個基因(引物見表4)進行RT-qPCR驗證，結果證實RT-qPCR得到的數據和轉錄組數據基本一致(圖7)，轉錄組數據可靠。

表4 熒光定量PCR引物列表

3 討論

3.1 中華鱉雌性差異表達基因

本研究結果顯示，中華鱉雄性和雌性各3個轉錄組的clean data分別為7 546 271 100、6 996 473 400、6 951 651 900 bp和10 234 669 800、7 293 149 700、7 077 323 100 bp。以Q30作為堿基測量的標準，Q30均大于91.6%，在隨機打斷的cDNA片段測得的GC含量大于51.4%。將3齡中華鱉精巢和卵巢的轉錄組進行比較分析，共得到10 145個差異基因，其中在卵巢中4 034個基因上調、6 111個基因下調。基于對數據庫中的比對信息，對性別相關的GO富集前26個主要富集通路基因進行初步篩選，獲得117個性別表達差異的基因，結合KEGG富集及蛋白互作發現，對中華鱉性別決定有影響的性別相關基因，如在卵巢中高表達的Foxl2、Sox3、Sox7、Bmp4、Csf1、Wnk2、Mtmr4、Gsc、Nr2f2、Fgf9、Tgfβ2、Tgfβ3、Rara、Smd2、Wdr1、Smd3和Smad5等。

Foxl2是一個重要的雌性偏好基因，在模式生物斑馬魚中，Foxl2a和Foxl2b作用于雌性卵巢的分化和維持，Foxl2b還有抑制卵巢分化的作用[19]。Foxl2在虹鱒卵巢中也特異性表達[20]。Foxl2在泰國斗魚雌性性腺中的表達量要高于雄性，組織表達具有顯著的特異性和性差異，對泰國斗魚卵巢發育和性別調控有重要的作用[21]。高麗麗等對Foxl2基因的研究中，雌性中華鱉Foxl2基因對性激素的應答基本相同，但在雄性中Foxl2的應答則不同[22]。此外，Foxl2基因在香港牡蠣性腺中的表達量也要高于其他組織[23]。成纖維生長因子9(Fgf9)是性腺發育的重要信號因子，敲低Fgf9抑制山羊睪丸間質細胞的T分泌[24]；Fgf9敲除的小鼠還會發生性逆轉，嚴重時出現死亡[25]。Bmp4是生殖細胞發生分化和形成配子的重要信號因子[26]，Bmp4與靶細胞上的受體結合激活細胞內Smads信號轉導蛋白，將信號傳入細胞核調節原始生殖細胞的生成[27]。Smd2和Smad3是調節卵巢發育的重要因子[28]。Yao B等研究發現，在斜帶石斑魚中Sox3基因持續表達誘導原始生殖細胞向卵原與卵母細胞發育，若停止表達則反之向精原細胞發育[29]。敲除Sox3基因后的小鼠出生后的生長繁殖受到影響[30]。Sox3基因在金線魚性腺中的表達量最高，存在組織表達和性別表達差異[31]。Sox7是爪蟾生殖細胞發育所必須的轉錄因子，對爪蟾生殖質的定位、轉錄激活、生殖細胞維持有重要的作用[32]。肌維管束蛋白14(Mtmr14)在睪丸、大腦等組織廣泛表達[33]；過表達Mtmr4會減弱早期核體內Smad3的活性，從而影響其核定位[34]。在本實驗蛋白與蛋白相互作用中Bmp4、Mtmr4均與Smad2、Smad3存在正相關作用，而Smad2、Smad3與Foxl2存在正相關作用，由此表明，在中華鱉中Bmp4與Mtmr4可能通過作用于Smad2和Smad3來影響Foxl2的表達。

此外，Diakite Haidou研究發現，Nr2f2在尼羅羅非魚中有Nr2f2a和Nr2f2b兩個基因，Nr2f2a和Nr2f2b在卵原細胞和初級卵母細胞中均有表達，Nr2f2a在體細胞和雄性生殖細胞中均有表達，但是在尼羅羅非魚雄魚精巢中沒有檢測到Nr2f2b的表達[35]。此外，下調Sara的表達水平或者阻斷Tgfβ信號將改善Wdr81不足引起的海馬神經發生障礙[36]，即Wdr81可能對Tgfβ相關基因的表達會產生一定的影響。吳燕華研究發現，敲除小鼠卵巢中的Id4基因后，小鼠仍然具備繁殖能力，在哺乳動物卵巢發育中，敲除Id4促進增殖的作用可能是通過其他轉錄因子來補償的[37]。在中華鱉性別表達基因，Wdr1和Nr2f2在卵巢中的表達量顯著高于雄性，可能與卵巢的發育有關。

3.2 中華鱉雄性差異表達基因

本研究中華鱉雄性中表達量較高的基因有Dmrt1、Amh、Sox30、Cct4、Cct5、Cldn11、Crem、Smc4、Top2a、Hook1、Spag6、Spata16、Spata24、Hormad1、Dhcr24、Dnai2、Dynll1、Gabarap、Herc4、Herpud2、Hexb、Hpgds、Ica1l、Meioc、Mlh1、Nme5、Nphp1、Pacrg、Ran、Rnf17、Rnf212b、Ropn1l、Saxo1、Spa17、Spire1、Stra13、Syce3、Tekt2、Tekt3、Ttc5和Zmynd15等。

Dmrt1基因是最古老的性別相關基因，在精巢的分化中起著重要的作用[38]。黃洋等的研究發現Dmrt1基因在泰國斗魚雄性的性腺和脾臟中特異性表達，在雌性性腺和肝臟中不表達，可能只參與泰國斗魚雄性的生殖發育[39]。在黃喉擬水龜中Dmrt1在精巢中的表達高于卵巢，可能參與黃喉擬水龜的性別決定過程[40]。孫偉等的研究發現，Dmrt1基因過表達后，會引起ZW型胚胎向雄性分化，在中華鱉雄性性別分化的過程中發揮著重要的作用[41]。本研究對3齡中華鱉的研究結果顯示，Dmrt1、Dmrta2和Dmrtb1在精巢中表達顯著高于卵巢，表明Dmrt1、Dmrta2和Dmrtb1為中華鱉雄性偏好基因。谷偉等研究發現，激素誘導使偽雄虹鱒魚腦和性腺中的抗穆勒氏管激素Amh基因表達量顯著高于雌性虹鱒，在雄性性腺分化中發揮著重要的作用[42]。在斜帶石斑魚中，Amh-Smad-Cyp19ala有信號通路存在，Amh-Smad和Amh-Cyp19ala呈負反饋調節，Amh激活Smad信號通路，而Smad和Cyp19a1a啟動子Smad結合位點結合，降低Cyp19ala的表達量，使得雌激素表達量下降，雄激素表達量上升[43]。Amh基因在黃顙魚中表達水平最高的是性腺，其次是腦和肝臟，1齡黃顙魚雄性性腺中的表達量最高，Amh基因可能作用于黃顙魚雄性性腺分化[44]。Amh基因是中華鱉早期性腺分化的關鍵因子，Amh在雄性中華鱉成體性腺中高水平表達，在胚胎期，使用芳香化酶抑制劑后Amh表達量上升，使用雌二醇誘導性逆轉則Amh表達量下降[45]。本研究在3齡中華鱉中Amh基因在精巢中的表達量顯著高于卵巢，因此，Amh為中華鱉雄性性別決定的關鍵因子[11]。Cldn11缺失會影響雄性小鼠的精子發育，而對雌性的生育能力無影響[45]。鉤狀同系物蛋白基因hook1與精子的連接有關，hook1基因缺失將對精子的連接產生影響[11]。Herc4基因對精子的產生和功能發揮起作用[46]。三倍體泥鰍Hormad1基因在雄性性腺中的表達量最高，雌性則在鰭中的表達量最高，可能調控其早期性腺分化[47]。Ptgds在公雞弱精癥的睪丸組織中的表達量顯著高于對照組[47]。Spag6基因被敲除的小鼠精子數量減少，精子結構被破壞，Spag6在精子運動及發生中發揮著重要的作用[48]。Sox30為睪丸特異轉錄因子，在出生后7 d小鼠的睪丸組織中就可以檢測到Sox30 mRNA[49]；Sox30基因缺失后，小鼠一部分精母細胞發生缺陷，而精子變形過程停止，同時缺失Sox30基因的小鼠表型與敲除Crem轉錄因子的小鼠表型基本一致[50]；Sox30在尼羅羅非魚性腺中特異性表達，在精巢的生殖細胞和卵巢的體細胞中表達[51]。唐耀浩對羅非魚Sox30基因的研究發現，Sox30在尼羅羅非魚精巢中表達，在90～180 d尼羅羅非魚的卵巢中不表達，Dmrt1敲低后，尼羅羅非魚精巢中Sox30的表達量降低，兩者可能存在正相關作用，Dmrt1過表達后Sox30的轉錄活性升高，Dmrt1與Sox30啟動子上的順式作用元件結合正向調控Sox30基因的轉錄[52]。此外，環磷酸腺苷反應元件調節物基因Crem基因的缺失會引起小鼠不育[53]。

在中華鱉雄性蛋白-蛋白互作子網絡中，Smc4、Msh4、Msh5和Actl6b正向調控Top2a和Top2b的表達，Atf1與Crem兩者正相關。雌性相關蛋白Wdr1通過Wdr1-Acta1-Smarcd1-Baz1b負反饋作用于雄性相關蛋白Top2a和Top2b的表達；同時雌性相關蛋白Wdr1還通過Wdr1-Acta1-XLOC_043302-Egr1和Atf1負反饋作用于Crem。Smc4、Top2a、Top2b和Crem在中華鱉精巢的發育過程中發揮著重要的作用。

本研究通過GO、KEGG注釋及蛋白-蛋白互作，獲得了大量的差異基因的功能信息。在GO數據庫中對功能信息進行注釋并分類，包括生殖過程、配子產生、性別分化、性腺發育等。這些差異基因功能注釋信息表明其可能在中華鱉的生殖、配子產生、性別分化和性腺發育等方面發揮作用，在中華鱉性別形成的二態性的形成過程中產生作用。3齡中華鱉為性別成熟期，進入繁殖階段，該階段產生的性別差異較多，有利于從廣泛篩選性別相關基因方面，為后續中華鱉性別決定機制的研究提供重要的參考。

KEGG注釋結果同樣豐富了功能通路信息，如PI3K-Akt信號通路、Wnt信號通路、催產素信號通路、Foxo信號通路、卵巢類固醇生成素信號通路等。其中PI3K-Akt信號通路在雌性生殖細胞生長發育中產生作用[54]。此外，催產素信號通路和卵巢類固醇生成素信號通路等與性腺的分化和發育有重要的作用。同時，本實驗還對差異基因進行蛋白-蛋白互作網絡構建，結果表明Bmp4、Smad2和Smad3調控Foxl2基因的表達，對中華鱉卵巢的發育有著重要的作用。