基于转录组分析不同温度下池蝶蚌性别相关基因的表达

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

(南昌大学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 ℃以上的温度，适合雌性的生殖发育。本研究为探索池蝶蚌的性别分化和性腺发育做出一定的参考，关于温度对池蝶蚌性腺发育和性别分化的影响，还涉及诸多基因与其产物之间的相互作用，依然有待进一步的深入发掘和探索。