中华鳖性别差异表达基因比较研究

周先文，王晓清，王 佩，江 辉，曾 丹，彭英海

(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基因的表达，对中华鳖卵巢的发育有着重要的作用。