基于卵巢转录组中华绒螯蟹大规格亲本育苗优势的分子机理

宣富君 王坤 曹佳欣 付龙龙 张建光 潘建林 管卫兵 成永旭

摘要：种质退化是当下河蟹产业健康发展所面临的重要课题之一，目前科研单位和养殖户逐步认识到选用大规格亲本（雄性：≥ 200 g;雌性：≥ 150 g）所育苗种在后续成蟹养殖上的优势，纷纷推荐和选用大规格亲本苗种，但其潜在的分子机制尚不清楚。本研究首次基于卵巢比较转录组（大规格与一般规格）探讨中华绒螯蟹大规格亲本潜在的分子优势。结果表明：2组共检测到8 772个差异表达基因（5 307个上调，3 465个下调），其中43个常见差异基因出现在与生殖、免疫和生长相关的15条通路中。经qRT-PCR验证，对其中9个相关基因（包括TRINITY_DN13931_c0_g1、TRINITY_DN1908_c0_g2和TRINITY_DN6686_c0_g1）开展进一步研究。这些被证实的常见差异表达基因主要富集于细胞凋亡、胰岛素信号和mTOR信号等途径。为进一步探讨利用大规格亲本遗传育种奠定了重要的分子基础，并在生长、繁殖和免疫方面指出了利用大规格亲本育苗潜在的分子优势，不仅可以丰富甲壳动物繁育生物学的基本理论，而且可为今后河蟹大规格优质苗种的推广提供必要的理论支撑。

关键词：中华绒螯蟹;大规格亲本;育苗;转录组;分子机理

中图分类号： S968.25  文献标志码： A

文章编号：1002-1302（2021）23-0166-07

收稿日期：2021-08-01

基金项目：江苏省种业振兴“揭榜挂师”专项（中华绒螯蟹种质资源）;江苏现代农业产业技术体系建设专项（编号：JATS[2020]362）;宁夏自治区科技厅重点项目（编号：2020ZDYF0860）;江苏省科技厅省政策引导类计划-苏北科技专项（编号：SZ-YC202041）;江苏省产学研项目（编号：BY2020587）。

作者简介：宣富君（1981—），男，浙江绍兴人，博士，副研究员，主要从事甲壳动物生殖生物学研究。E-mail：swimming_crab@126.com。

中華绒螯蟹（Eriocheir sinensis）别称河蟹、大闸蟹，是我国重要的经济水产物种。得益于20世纪70年代河蟹人工育苗技术的突破[1]和土池育苗技术的推广应用[2]，目前河蟹养殖已经遍布全国31个省（市、自治区），包括台湾[3]和西藏自治区[4]。河蟹养殖产量从2000年的23万t增加到2018年的75万t[5]，产值近500亿元，被称为特种水产养殖业中的“奇葩”，超过世界其他经济蟹类养殖和捕捞总和[6-7]。

种质退化是当下河蟹产业健康发展所面临的重要课题之一。育苗单位曾盲目追求经济利益和降低生产成本，无序引进其他各大水系的苗种及采用小规格亲本用于苗种生产，造成商品蟹形状退化，规格小，养殖效益低下[8]。获取高质量的亲本是河蟹育苗成功的关键[6]。近几年在科研单位的积极引导下，越来越多的育苗单位和养殖户认识到选用大规格亲本（雄性：≥200 g;雌性：≥150 g）苗种在后续养成优质商品蟹上的优势。例如，从不同规格亲本的繁殖性能和育苗效果看，大规格亲本的抱卵量、生殖力和产卵量明显高于其他规格。虽然不同规格亲本的胚胎大小无明显差异，但研究表明大规格亲本所产子代可以遗传其生长速度快的优良性状，大规格亲本所产大眼幼体的售价也往往高于普通大眼幼体[8-9]。在大眼幼体养成扣蟹阶段，大规格亲本苗种明显生长快，营养状况好，性能上更适合成蟹养殖[10-11]。更为重要的是，实践证明大规格苗种在成蟹养殖阶段可以显著提高河蟹的规格，产出更高比例的优质商品蟹，效益增加明显[12]。然而，与其他水产经济物种相比，中华绒螯蟹在种质资源机理方面的研究相对滞后，大规格亲本育苗潜在优势分子机制方面的研究至今尚未见报道。

本研究拟通过对发育成熟的大规格雌性亲本（约 250 g）与一般规格（约100 g）的卵巢组织进行转录水平的高通量测序，并从生长发育、繁殖及免疫性能等方面进行相关基因的差异表达分析，进而探讨大规格亲本育苗潜在的分子优势。本工作不仅可以丰富甲壳动物繁育生物学的基本理论，而且可以为今后开展中华绒螯蟹的遗传育种及河蟹育苗产业推广普及大规格亲本提供必要的理论依据。

1 材料与方法

1.1 样本获取

2020年11月30日，从江苏大仁水产良种有限公司盐城射阳育苗基地（120.455°E、33.886°N）获取大规格雌性亲本3只[（248.60±3.12） g，F-O-D 组]，一般规格3只[（102.87±1.52） g，F-O-Z组]，分别暂养于盐城师范学院湿地学院配备的从大连汇新公司定制室内可控养殖系统的水桶内。3 d 后，待河蟹适应环境，个体经冰浴麻痹后迅速解剖;取出的卵巢组织分别放入2 mL冻存管，经液氮冷却，最终储存于-80 ℃超低温冰箱中。

1.2 RNA提取及测序

使用TRIzol 试剂盒（TaKaRa公司）提取不同卵巢样本中的总RNA，之后通过Agilent 2100 Bioanalyzer （美国）（D260 nm/D280 nm为1.8～2.2）和 15 g/L 琼脂糖凝胶电泳进行质检（28S rRNA/18S rRNA>1.0）。样品质检合格后，用干冰保护并送往上海元莘生物医药科技有限公司测序。过程简述如下：使用用于Illumina（美国马萨诸塞州新英格兰生物实验室）的TruSeqTM RNA样品制备试剂盒构建总RNA（5 μg/样品）的逆转录cDNA文库。然后，经修复末端修饰处理后，通过15个修复链式反应周期扩增cDNA。最后，采用配对末端法的Illumina HiSeq 4000平台用于cDNA簇测序[13]。

1.3 RNA-seq数据的预处理

Raw reads的质量控制是根据FastqStat.jar工具软件（1.0版本）实现的：（1）去除5′端含有非AGCT的条形码序列和碱基;（2）消除N>10%的成对读数;（3）排除Q<20的碱基的低质量读数;（4）排除长度小于25 bp的读数。Clean reads的质量控制是参照FastQC（0.11.4版本）软件执行的[14]。

1.4 转录组组装与基因表达分析

Trinity是一种基于Illumina短片段序列（Clean reads）的无参考基因组的从头组装软件[15]。在本研究中，利用Trinity Assembler（2.6.6版本）软件完成从Illumina平台获取序列的de novo组装，并在Trinity中使用ORF方法进行基因预测。然后，基于RNA-seq数据量化软件Kallisto（0.43.1版本），对不同样本中的基因表达水平进行定量分析[16]。

1.5 样品间的关联分析

为调查不同样本，尤其是重复样本之间的相关性，本研究进行了样本间的关联分析。此外，利用主成分分析（PCA）根据样本聚类的表达情况进一步识别大样本的影响。

1.6 差异表达基因（DEGs）分析

使用edgeR（3.24版本）对2组之间的基因进行差异表达分析[17]。选择错误发现率（FDR）<0.05和| log2差异倍数（FD）|>1作为DEGs调查的截止阈值。结果通过火山图和散点图可视化。然后，利用plot_cluster_exp（1.1.0版本）软件中的距离算法（样本间：Spearman相关系数;基因间：Pearson相关系数）进行DEGs的聚类分析，揭示具有相同或相似表达模式的基因。聚类结果通过热图可视化。

1.7 基于文献检索的常见差异表达基因（co-DEGs）探索

通过文献数据库检索筛选与生殖、生长和免疫相关的途径和富集基因。本研究检索的电子文献数据库，包括Embase、PubMed和Cochrane。此外，在本研究中，这些途径中富集的基因与DEGs交叉，可获得与生殖、生长和免疫反应相关的co-DEGs。

1.8 qRT-PCR验证

为了筛选与测序结果趋势相同的DEG，对9个FC较高的co-DEG进行基于qRT-PCR的定量验证。过程简述如下：使用TRIzol试剂（Invitrogen 公司，美国）从每组样品中提取并定量总RNA（10 μL），并使用SYBR Green PCR Mastermix（Applied Biosystems公司，美国）进行反转录。然后，在ABi PrisMR 7500序列检测系统（Applied Biosystemg公司，美国）协助下进行qRT-PCR验证。以ACTB基因为参考。扩增-转录反应条件如下：95 ℃ 30 s;95 ℃ 5 s，40个循环，60 ℃下 34 s;在每个延长周期的终点，获得信号，然后研究放大曲线。表1列出了引物序列，而目标基因的相对表达计算过程则参考2-ΔΔCT方法[18]。

1.9 通路富集分析

本研究利用go_enrichment（2.1.0版本）软件，进一步对已得到验证的co-DEG开展京都基因和基因组百科全书（KEGG）通路富集分析[19-20]。其中，P<0.05被认为是KEGG通路显著富集的截止值。

2 结果与分析

2.1 样本间关联分析

在质量控制和数据预处理后，将测序数据纳入F-O-D组和 F-O-Z 组之间进行关联分析，以检验试验的可靠性和样本选择的合理性。PCA分析表明，2组间的样本以分散的方式分布，而组内的样本则聚集在一起（图1）。相关聚类热图显示，参与本研究的样本的相关性很高（图2）。

2.2 差异表达基因（DEGs）分析

根据DEGs测序数据，在F-O-D组和 F-O-Z 组之间共发现5 307个上调基因和3 465个下调基因。DEGs分析结果的火山图如图3-A所示，它直接指示了不同组中样本之间基因的相对表达值。所有这些DEGs的散点图如图3-B所示。此外，对样本中的基因表达模式进行聚类分析显示，样本中的基因表达模式存在明显差异（图4）。

2.3 常见差异表达基因（co-DEGs）调查及qRT-PCR验证

根据文献检索，共有7条生殖相关途径（视黄醇代谢、细胞周期、mTOR信号通路、Wnt信号通路、胰岛素信号通路、卵巢类固醇生成和孕酮介导的卵母细胞成熟），4条生长相关途径（mTOR信号通路、PI3K-Akt信号通路、Wnt信号通路、TGF-β信号通路）和7条免疫相关途径[细胞周期、溶酶体、内吞、吞噬体、凋亡、细胞粘附分子（CAMs）和自然杀伤细胞介导的细胞毒性]被筛选出来。这15条途径（去除重叠途径后）中共富集338个基因，结合本研究调查的DEGs，找到43个与这些生殖、免疫和生长途径相关的co-DEGs。最后， qRT-PCR验证了其中9个co-DEG（TRINITY_DN6686_c0_g1、TRINITY_DN6296_c2_g1、 TRINITY_DN2026_c1_g1、TRINITY_DN1908_c0_g2、TRINITY_DN605_c0_g2、TRINITY_DN2717_c0_g1、TRINITY_DN608_c5_g1、TRINITY_DN13931_c0_ g1和TRINITY_DN35647_c0_ g1），结果显示其均与测序结果趋势一致（图5）。

2.4 验证co-DEGs的富集分析

对上述9个经qRT-PCR验证的co-DEGs开展KEGG富集分析（表2），结果表明，这些DEGs主要富集于细胞凋亡（ko04210，基因：TRINITY_DN13931_c0_g1、TRINITY_DN1908_c0_g2和TRINITY_DN6686_c0_g1）、胰岛素信号通路（ko04910，基因：TRINITY_DN13931_c0_g1、TRINITY_DN2717_c0_g1和TRINITY_DN608_c5_g1）和mTOR信號通路（ko04150，基因：TRINITY_DN13931_c0_g1和TRINITY_DN2717_c0_g1）。

3 讨论

目前河蟹育苗使用大规格亲本已成为行业发展的新趋势，但其潜在的分子优势仍不清楚。本研究首次基于卵巢比较转录组分析，共检测到8 772个DEGs（其中5 307个上调，3 465个下调），而所有这些DEGs与文献报道的15条跟生殖、生长和免疫相关通路中的338个基因整合时，共检索出43个co-DEGs。经过qRT-PCR验证，最终对其中9个co-DEGs，包括TRINITY_DN13931_c0_g1、TRINITY_DN1908_c0_g2和TRINITY_DN6686_c0_g1进行了富集通路分析，结果表明这些co-DEGs主要富集在胰岛素信号通路（生殖相关）、mTOR信号通路（生长相关）和细胞凋亡（免疫相关）等途径中。

胰島素信号与生殖之间的密切关系已在蠕虫和昆虫的身体发育中得到揭示[21-22]，而控制个体规格大小涉及包括胰岛素在内的多条生长调节通路[23]。先前的一项研究表明，类胰岛素性腺激素基因广泛参与雌性青蟹（Scylla paramamosain）的生殖过程[24]。通过转录组测序，目前已在甲壳动物物种中鉴定出多个参与胰岛素信号通路的类胰岛素肽[25]。重要的是，雌性蟹已被认为是研究蟹类胰岛素肽生物学功能的重要载体[26]。例如，胰岛素在雌性青蟹卵黄发生和卵母细胞成熟中的抑制作用是通过某些信号途径实现的，如自分泌和旁分泌[21]。有趣的是，刺激卵黄蛋白原的表达则是通过蟹体内的胰岛素实现的[25]。这些胰岛素相关蛋白，包括类胰岛素肽，不仅参与卵巢发育[26]，还逆向调节葡萄糖代谢并参与中华绒螯蟹的抗病原体感染过程[27]。

mTOR是磷脂酰肌醇3-激酶相关激酶家族的成员，调节包括生长控制在内的各种生物过程[28]。其中，mTOR通路在动物发育、寿命和幼虫发育的调节过程中，与胰岛素信号通路直接相关[29]。由于甲壳类动物必须定期脱落外骨骼才能发育和生长，因此调节生长和蜕皮的mTOR信号通路已被人们广泛研究[30]。已经证明，mTOR信号通路相关基因的表达有助于地蟹（Gecarcinus lateralis）的生长，参与幼体热适应过程中的蜕皮调节[31]。此外，mTOR通过整合内在信号（蜕皮激素）和外在信号（热应激）来调节十足甲壳动物的蜕皮和生长。蜕皮腺是十足目甲壳动物类固醇激素产生和蜕皮周期调节的来源，对蟹的外部环境和内部生理信号都有反应[30]。在蟹等甲壳类动物中，mTOR的活性直接或间接控制着蜕皮腺驱动相关基因的转录[32]。

凋亡通路已被证明参与免疫过程[33]。最近，在青蟹中发现了效应caspase，如Sp-caspase，它表现出免疫反应和细胞凋亡[34]。在蟹类中，细胞凋亡可以诱导血细胞的活性，而血细胞在受到病原菌刺激后，在抵御病原菌入侵方面起着至关重要的作用[35]。对于中华绒螯蟹之前的一项研究表明，通过某些蛋白（如Es IAP1）减少凋亡，可以调节细胞在生长过程中的活性[36]。最近的转录组分析发现，中华绒螯蟹在生长过程中，由于各种因素（如肝胰腺坏死疾病） 引起的生长抑制和存活率降低均与细胞凋亡显著相关[37]。重要的是，凋亡通路中的nm23和caspase等基因参与了中华绒螯蟹个体生长过程中对不利环境的抵抗[38-39]。因此，中华绒螯蟹大规格亲本可能通过细胞凋亡的基因表达表现出较强的环境抗逆性，进而有利于自身及繁育后代在养殖过程中表现出优势。在本研究中该通路共涉及3个已验证的co-DEGs，分别为TRINITY_DN13931_c0_g1、TRINITY_DN1908_c0_g2和TRINITY_DN6686_c0_g1。这些co-DEGs在细胞凋亡、胰岛素信号和mTOR信号通路中亦有富集。因此进一步推测，TRINITY_DN13931_c0_g1、TRINITY_DN1908_c0_g2和TRINITY_DN6686_c0_g1富集的细胞凋亡、胰岛素信号通路和mTOR信号通路可能与大规格亲本繁育优势有关。

4 结论

本研究为进一步探索大规格亲本在中华绒螯蟹育种中的优势奠定了重要的分子基础，相关数据已上传NCBI（登入号：PRJNA669606）。通过基于卵巢比较转录组差异表达基因分析，首次探讨了利用大规格亲本在生长、繁殖、免疫等方面潜在的分子优势。其中，凋亡通路、胰岛素信号通路和mTOR信号通路可能在中华绒螯蟹大规格亲本的优势育种中发挥重要作用。

参考文献：

[1]赵乃刚. 河蟹育苗工厂化的进展[J]. 渔业现代化，1979，6（1）：34-35.

[2]韩炳炎. 河蟹养殖高产技术问答[M]. 北京：中国农业出版社，1996.

[3]虞丽娟，杨劲松，凌培亮，等. 基于物联网智慧服务的中华绒螯蟹蟹种质量动态追溯系统研究[J]. 水产学报，2013，37（8）：1262-1269.

[4]陈 伟，王 春，杨印蹼，等. 中华绒螯蟹在西藏高原条件下越冬期生化组分的变化[J]. 上海海洋大学学报，2014，23（5）：733-740.

[5]张显良. 中国渔业统计年鉴[M]. 北京：中国农业出版社，2019.

[6]Cheng Y X，Wu X G，Yang X Z，et al. Current trends in hatchery techniques and stock enhancement for Chinese mitten crab，Eriocheir japonica sinensis[J]. Reviews in Fisheries Science，2008，16（1/2/3）：377-384.

[7]Wang Q D，Liu J S，Zhang S Y，et al. Sustainable farming practices of the Chinese mitten crab （Eriocheir sinensis） around Hongze Lake，Lower Yangtze River Basin，China[J]. Ambio，2016，45（3）：361-373.

[8]王少兵，劉乃更，姜晓东，等. 三种亲本规格河蟹土池育苗的试验总结[J]. 科学养鱼，2019（6）：12-13.

[9]茅海成，王高龙，杨永超，等. 中华绒螯蟹不同规格亲蟹池塘生态育苗效果的生产性评估[J]. 水产科技情报，2014，41（5）：233-236.

[10]陈军伟，马旭洲，王 武，等. 不同规格中华绒螯蟹母本子代的生长特性比较[J]. 动物学杂志，2016，51（5）：895-906.

[11]何先林，姜晓东，王海宁，等. 池塘养殖不同规格中华绒螯蟹扣蟹生化组成的比较研究[J]. 水产科技情报，2019，46（6）：316-319，323.

[12]张红水. 大亲本蟹苗养殖试验研究[J]. 乡村科技，2018（4）：110-111.

[13]Ran M L，Chen B，Li Z，et al. Systematic identification of long noncoding RNAs in immature and mature porcine testes[J]. Biology of Reproduction，2016，94（4）：77，1-9.

[14]Andrews S. Babraham bioinformatics-FastQC a quality control tool for high-throughput sequence data[EB/OL]. [2021-06-20]. https：//www. bioinformatics. babraham.ac.uk/projects/fastqc/.

[15]Bankar K G，Todur V N，Shukla R N，et al. Ameliorated de novo transcriptome assembly using Illumina paired end sequence data with Trinity Assembler[J]. Genomics Data，2015，5：352-359.

[16]Du Y H，Huang Q H，Arisdakessian C，et al. Evaluation of STAR and kallisto on single cell RNA-seq data alignment[J]. Genes，Genomes，Genetics，2020，10（5）：1775-1783.

[17]Robinson M D，McCarthy D J，Smyth G K.edgeR：a bioconductor package for differential expression analysis of digital gene expression data[J]. Bioinformatics，2009，26（1）：139-140.

[18]Livak K J，Schmittgen T D.Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCT method[J]. Methods，2001，25（4）：402-408.

[19]Consortium G O.The Gene Ontology （GO） database and informatics resource[J]. Nucleic Acids Research，2004，32（S1）：D258-D261.

[20]Kanehisa M，Goto S. KEGG：kyoto encyclopedia of genes and genomes[J]. Nucleic Acids Research，2000，28（1）：27-30.

[21]Bao C C，Yang Y N，Huang H Y，et al. Inhibitory role of the mud crab short neuropeptide F in vitellogenesis and oocyte maturation via autocrine/paracrine signaling[J]. Frontiers in Endocrinology，2018，9：390.

[22]Dissous C.Venus kinase receptors at the crossroads of insulin signaling：their role in reproduction for helminths and insects[J]. Frontiers in Endocrinology，2015，6：118.

[23]Nijhout H F，Riddiford L M，Mirth C，et al. The developmental control of size in insects[J]. Wiley Interdisciplinary Reviews：Developmental Biology，2014，3（1）：113-134.

[24]Huang X S，Ye H H，Huang H Y，et al. An insulin-like androgenic gland hormone gene in the mud crab，Scylla paramamosain，extensively expressed and involved in the processes of growth and female reproduction[J]. General and Comparative Endocrinology，2014，204：229-238.

[25]Huang X S，Feng B Y，Huang H Y，et al. In vitro stimulation of vitellogenin expression by insulin in the mud crab，Scylla paramamosain，mediated through PI3K/Akt/TOR pathway[J]. General and Comparative Endocrinology，2017，250：175-180.

[26]Huang X S，Ye H H，Feng B Y，et al. Insights into insulin-like peptide system in invertebrates from studies on IGF binding domain-containing proteins in the female mud crab，Scylla paramamosain[J]. Molecular and Cellular Endocrinology，2015，416：36-45.

[27]Wang L，Chen H，Wang L L，et al. An insulin-like peptide serves as a regulator of glucose metabolism in the immune response of Chinese mitten crab Eriocheir sinensis[J]. Developmental & Comparative Immunology，2020，108：103686.

[28]Sarbassov D D，Ali S M，Sabatini D M. Growing roles for the mTOR pathway[J]. Current Opinion in Cell Biology，2005，17（6）：596-603.

[29]Jia K L，Chen D，Riddle D L.The TOR pathway interacts with the insulin signaling pathway to regulate C. elegans larval development，metabolism and life span[J]. Development，2004，131（16）：3897-3906.

[30]Das S，Pitts N L，Mudron M R，et al. Transcriptome analysis of the molting gland （Y-organ） from the blackback land crab，Gecarcinus lateralis[J]. Comparative Biochemistry and Physiology Part D：Genomics & Proteomics，2016，17：26-40.

[31]Cosenza K S. Role of ecdysteroids on Myostatin and mTOR signaling gene expression in molt-dependent growth and atrophy of skeletal muscle in Gecarcinus lateralis and Carcinus maenas[D]. The Colorado State University Libraries，2016.

[32]Shyamal S，Das S，Guruacharya A，et al. Transcriptomic analysis of crustacean molting gland （Y-organ） regulation via the mTOR signaling pathway[J]. Scientific Reports，2018，8（1）：7307.

[33]Lum J J，Bauer D E，Kong M，et al. Growth factor regulation of autophagy and cell survival in the absence of apoptosis[J]. Cell，2005，120（2）：237-248.

[34]Li J S，Dong L X，Zhu D P，et al. An effector caspase Sp-caspase first identified in mud crab Scylla paramamosain exhibiting immune response and cell apoptosis[J]. Fish & Shellfish Immunology，2020，103：442-453.

[35]Johansson M W，Keyser P，Sritunyalucksana K，et al. Crustacean haemocytes and haematopoiesis[J]. Aquaculture，2000，191（1/2/3）：45-52.

[36]Qu C，Sun J J，Xu Q S，et al. An inhibitor of apoptosis protein （EsIAP1） from Chinese mitten crab Eriocheir sinensis regulates apoptosis through inhibiting the activity of EsCaspase-3/7-1[J]. Scientific Reports，2019，9（1）：20421.

[37]Yan B Y，Liu X M，Zhou Y Q，et al. Transcriptomic analysis reveals that hepatopancreatic necrosis disease in Eriocheir sinensis （Chinese mitten crabs） may be the result of autophagy and apoptosis[J]. Aquaculture，2020，515：734579.

[38]Jin X K，Li W W，He L，et al. Molecular cloning，characterization and expression analysis of two apoptosis genes，caspase and nm23，involved in the antibacterial response in Chinese mitten crab，Eriocheir sinensis[J]. Fish & Shellfish Immunology，2011，30（1）：263-272.

[39]Jin X K，Li W W，Wang Q.Caspase and nm23：apoptosis genes linked to the antibacterial response of the Chinese mitten crab[J]. Bioengineered Bugs，2011，2（3）：174-177.