黄河鲤生长性状的转录组解析

李佩珍，许建，朱优秀，冯建新，江炎亮，张瀚元*

(1. 上海海洋大学， 水产科学国家级实验教学示范中心， 上海 201306；2. 中国水产科学研究院，农业农村部水生动物基因组学重点实验室， 北京 100141； 3. 河南省水产科学研究院， 河南 郑州 450044)

鲤(Cyprinuscarpio)是世界上淡水鱼类中养殖范围最广的物种之一，目前有100多个国家和地区养殖[1]。中国是世界最大的鲤养殖和消费国，2020年养殖产量289.7万t，在水产养殖中占有重要的地位[2]。黄河鲤具有金鳞赤尾、体型梭长、肉质细嫩、气味清香等特点，“豫选黄河鲤”新品种2005年通过全国水产原种和良种审定委员会审定后，迅速在沿黄各地推广开来，形成了年产数十万吨的养殖规模，在沿黄地区渔民增收、产业增效和优质水产品供给等方面贡献显著。然而，由于长期自繁自养造成种质退化，生长减慢，严重影响养殖效益。转录组学分析是筛选分子辅助育种标记的重要手段，对解析鱼类生长等经济性状的调控机制具有重要意义[3]，因此有必要借助鱼类转录组学研究，系统解密鱼类发育、生长等问题。目前已有大量研究通过分析多种不同鱼类的转录组数据，探究其生长性状相关基因，并取得了一定成果，Sun等[4]利用RNA-seq鉴定了珍珠龙胆石斑鱼(Epinephelusfuscoguttatus×E.lanceolatus)脑和肝脏组织的差异表达基因并发现GH / IGF轴和下游信号传导途径中的差异表达基因能够潜在提高杂交种的生长优势。Fu等[5]对高体重组和低体重组鳙(Hypophthalmichthysnobilis)的肝脏、下丘脑-垂体轴进行RNA-seq分析，分别鉴定出173和204个差异基因，其中9个基因位于生长相关的数量性状位点(QTL，quantitative trait loci)内。Guan等[6]通过对鳜(Sinipercachuatsi)生长差异的个体进行转录组测序，鉴定出与生长速率相关的32个差异表达基因，功能主要涵盖脂肪酸生物合成(fasn和acacb)、集合管酸分泌(atp6e和kcc4)、细胞周期(cdc20和ccnb)和胰岛素样生长因子(igfbp1)调控。然而关于鲤生长性状相关基因的挖掘，目前主要通过全基因组关联分析方法来筛选鲤生长性状相关QTL和SNP。Peng等[7]通过构建鲤高密度遗传连锁图谱，运用QTL定位和关联分析方法鉴定出22个生长相关QTL，并确定了候选基因，包括kiss2、igf1、smtlb、npffr1和cpe。Feng[8]等通过限制性位点相关的DNA标记和微卫星标记构建了长江鲤的高分辨率遗传连锁图谱，检测到21个生长相关性状的QTL。陈琳[9]基于全基因组关联分析和QTL方法，鉴定了parv、srpk2、fsrp5、igf1、igf3、grb10、igf1r、notch2、sfrp2等与鲤骨骼和个体生长相关的基因，鉴定了441个生长相关的SNP。因此有必要基于黄河鲤转录组数据，探究其生长性状相关基因，且研究发现黄河鲤生长性状相关基因主要在其脑、肝脏和肌肉组织中表达。研究表明在硬骨鱼类中，体重、体长等生长性状与神经激素调节密切相关，机体通过下丘脑-垂体-肝脏轴，调控生长激素和类胰岛素生长因子轴(GH/IGF轴)分泌重要的信号因子，如生长激素(GH)、类胰岛素生长因子1(IGF-1)和类胰岛素生长因子2(IGF-2)等[10]，并通过血液循环作用于肝脏细胞表面的受体，最终促进肌肉组织细胞的生长和分化[11]。本研究通过转录组学方法，挖掘体重差异显著的两组黄河鲤脑、肝脏和肌肉组织中的差异表达基因，通过功能富集探究影响其生长性状的功能基因和代谢通路，并进行实时荧光定量实验验证基因功能，为黄河鲤生长性状的分子标记辅助育种提供了基础数据和依据。

1 材料与方法

1.1 材料

实验用黄河鲤采自河南省荥阳市王村镇河南省大宗淡水鱼产业技术体系黄河鲤育种基地。选取同一家系中体重最低的5尾和体重最高的5尾，记录体重(表1)。从高体重组和低体重组中，各选取3尾体重差异极显著的个体，作为研究生长性状遗传差异的样本。解剖采集其脑、肝脏及肌肉组织，装入无RNA酶管，立即在液氮中冷冻，并放置于-80 ℃超低温冰箱保存，用于转录组测序和基因表达分析。

1.2 方法

1.2.1 RNA提取、文库构建和转录组测序

采集所选样本的脑、肝脏及肌肉组织，分别在液氮中研磨，使用RNeasy Mini Kit试剂盒(Qiagen公司，中国)对上述样本进行RNA提取，琼脂糖凝胶电泳检测所提取RNA的质量。对质量检测合格的RNA样本进行转录组测序，每个样本测序数据不少于6Gb。使用FastQC、Trimmomatic软件去除转录组测序数据中低质量的序列，获得高质量的测序数据。

1.2.2 筛选差异表达基因

利用Bowtie2构建Bowtie索引，使用Tophat2将过滤后的序列比对到鲤参考基因组[12]。利用Samtools对bam文件建立索引，用Cufflinks和Cuffmerge计算每个基因单位长度内的表达片段数目(FPKM，Fragments per kilobase of exon model per million mapped fragments)并生成一个完整的转录本集合。将基因表达量标准化，设置筛选参数为错误发现率(FDR，false discovery rate)小于0.05和差异倍数(FC，fold change)大于2。用Cuffdiff鉴定差异表达的转录本，寻找各组织中的差异表达基因，用VennDiagram软件包绘制韦恩图。用ggplot2软件包绘制脑、肝脏及肌肉组织差异表达基因火山图，用pheatmap软件包绘制脑、肝脏及肌肉组织差异基因聚类热图。

1.2.3 功能注释和通路富集分析

得到候选基因列表之后，使用DAVID 6.8软件进行差异表达基因的功能注释(GO，Gene Ontoloty)和通路富集分析(KEGG，Kyoto Encyclopedia of Genes and Genomes)。使用clusterProfiler软件包构建GO和KEGG富集分析图；并使用Cytoscape软件绘制基因互作网络分析图，展示差异表达基因富集到的GO通路和KEGG代谢通路之间的关系。

1.2.4 定量PCR验证

利用定量PCR方法对黄河鲤脑、肝脏和肌肉组织中的差异候选基因进行验证。验证方法为：从黄河鲤群体中，选取3尾低体重和3尾高体重个体，提取其脑、肝脏和肌肉组织的RNA，对总RNA纯度和完整性进行检测，用反转录试剂盒制备合成cDNA(Qiagen，中国)。使用Primer Premier 5软件设计候选基因引物，内参基因选用经过验证的β-actin[13]，引物由Qiagen公司合成。定量PCR反应总体系为15 μL，其中cDNA为1 μL，上下引物各0.3 μL，SYBR Green Master Mix(TOYOBO，日本)为7.5 μL，ddH2O为5.9 μL。每个基因的程序设置为95 ℃持续2 min，接下来是95 ℃持续15 s，59 ℃持续30 s，72 ℃持续1 min进行40个循环，熔解曲线为默认条件，进而完成定量PCR检测。采用2-ΔΔCt法计算基因的相对表达量，以P值小于0.05作为显著差异水平，验证候选基因的表达趋势与转录组分析结果是否一致。

表2 脑组织候选基因的定量PCR引物Tab.2 Quantitative PCR primers of candidate genes for brain tissue

2 结果与分析

2.1 转录组数据比对

利用Illumina Hiseq 2500高通量测序系统(Illumina，美国)对所有样本进行转录组测序，质控和比对结果见表3。由表3可以看出，所有样本的测序数据有效率平均值为98.34%，比对率平均值为76.58%。由以上数据可以得出，该转录组测序数据质量和比对率都较高，可为后续数据分析和实验验证提供可靠的原始数据。

表3 转录组测序数据质控和比对结果Tab.3 Quality control and comparison results of transcriptome sequencing data

2.2 差异表达基因鉴定与分析结果

筛选高体重组和低体重组黄河鲤各组织的差异表达基因，结果显示，脑组织显著差异表达基因6 519个，肝脏组织显著差异表达基因499个，肌肉组织显著差异表达基因19个，其中三个组织共表达基因1个(图1)。分别对脑、肝脏、肌肉组织中所有差异表达基因和样品进行双向聚类分析，结果表明，相关性较高，均一性较好(图2)。对黄河鲤转录组数据进行生长性状候选基因挖掘，分析发现，脑组织中有6 519个差异表达基因，剔除错误发现率或q值缺失基因后，有效差异表达基因个数为6 385个，其中上调基因3 443个，下调基因2 942个(图3 A)。肝脏组织中有499个差异表达基因，有效差异表达基因个数为461个，其中上调基因275个，下调基因186个(图3B)。肌肉组织中有18个差异表达基因，有效差异表达基因11个，其中上调基因6个，下调基因5个(图3C)。

图1 两组黄河鲤差异表达基因的韦恩图Fig.1 Venn diagram of differentially expressed genes in two Cyprinus carpio groups

图2 差异表达基因聚类热图注：A-脑组织；B-肝脏组织；C-肌肉组织,下同。Fig.2 Clustering heatmap of differential expression genes in different tissueNote: A- Brain tissue; B-Liver tissue; C-Muscle tissue,the same below.

图3 黄河鲤生长性状差异表达基因的火山图注：红色圆点表示显著上调基因；蓝色圆点表示显著下调基因；灰色圆点表示差异不显著的基因。Fig.3 The volcano plots of differentially expressed genes related to growth traits in Cyprinus carpioNote: Red dot means significantly up-regulated gene; blue dot means significantly down-regulated gene; gray dot means no significant difference gene.

2.3 功能注释和富集分析结果

对脑、肝脏和肌肉组织中的差异表达基因进行GO和KEGG富集分析。图4 A展示了脑组织KEGG富集结果，筛选出与黄河鲤生长性状相关度较高的10个代谢通路，分别显著富集在丝裂原活化蛋白激酶(mitogen-activated protein kinase，MAPK)信号通路、钙(calcium)信号通路、Wnt信号通路、细胞外基质受体相互作用(extracellular matrix receptor interaction，ECM-receptor interaction)等信号通路。图4B展示了肝脏组织KEGG富集分析结果，筛选出与黄河鲤生长性状相关度较高的10个信号通路，发现胰岛素(insulin)信号转导通路、精氨酸和脯氨酸代谢(arginine and proline metabolism)通路、MAPK信号通路、ErbB信号通路、肌动蛋白细胞骨架调节(regulation of actin cytoskeleton)等通路与生长密切相关。由于肌肉组织显著差异基因较少，未富集到相应的KEGG通路。用Cytoscape绘制候选基因和通路之间的网络图(图5)，MAPK和Wnt信号通路均为富集到候选基因较多的代表性信号通路。根据文献中脊椎动物生长相关的经典信号通路调控模式和差异基因的功能，结合富集分析、网络图结果，在黄河鲤脑组织中初步筛选出13个参与调控生长性状通路的候选基因，分别为mapk1、mapk7、mapk10、mapk12、wnt16、wnt8b、wnt10b、fgf13、fgf19、tgfbr2、notch2、ednra和lpl(表4)。其中，参与MAPK信号通路的有4个基因，分别为mapk1、mapk7、mapk10和mapk12基因。在网络图结果中，mapk1和mapk12是MAPK通路的节点基因，且mapk1基因在4个代表性通路中都发挥作用，表明mapk1基因可能在黄河鲤生长性状中发挥重要的调控作用。参与Wnt信号通路的有3个基因，分别为wnt16、wnt8b和wnt10b基因。在网络图结果中，wnt8b和wnt10b基因是Wnt信号通路的节点基因，且参与mTOR通路的调控。此外，参与成纤维细胞生长因子/成纤维细胞生长因子受体(fibroblast growth factor/ fibroblast growth factor receptor，FGF/FGFR)信号通路的有3个基因，分别为fgf13、fgf19和tgfbr2，以及参与Notch信号通路的notch2基因，参与内皮素受体A(endothelin receptor type A，EDNRA)信号通路的ednra基因，和参与脂肪代谢的蛋白脂肪酶基因(lpl)。这些基因在高、低体重组黄河鲤中，表现出不同程度的上调或下调表达。在黄河鲤肝脏组织中筛选出9个与生长性状相关的基因，分别是fgf19、fasn、fgfr2、srebf1、acac、igf2、acod、fads6和acsl3；在黄河鲤肌肉组织中筛选出lpl基因与生长性状密切相关。

图4 差异表达基因的通路富集结果Fig.4 KEGG enrichment results of differentially expressed genes

2.4 qRT-PCR验证转录组数据结果

通过qRT-PCR实验检测两组生长差异显著的黄河鲤脑、肝脏和肌肉组织的候选基因相对表达水平，并将实验结果与转录组分析结果相比较，结果发现脑组织部位的13个候选基因的表达水平与转录组预测趋势相吻合，然而肝脏组织候选基因的表达水平与转录组预测趋势差异较大，且肌肉组织候选基因数量太少，因此重点分析脑组织部位候选基因的相对表达量(图6)。结果表明，脑组织中mapk1、mapk7、mapk10、mapk12、wnt16、wnt8b、wnt10b基因的表达量均显著上调，fgf13、ednra、notch2、tgfbr2、fgf19、lpl基因的表达量均显著下调，这些基因qRT-PCR实验结果与转录组分析结果表达趋势相同，然而上调和下调程度有所差异。

图5 生长性状关键通路的基因互作网络分析图Fig.5 Gene interaction network analysis diagram among key pathways related to growth traits

表4 黄河鲤生长性状相关的13个候选基因信息Tab.4 Information of 13 candidate genes related to growth traits of Cyprinus carpio

图6 黄河鲤脑组织候选基因qRT-PCR结果Fig.6 qRT-PCR results of candidate genes in brain tissue of Cyprinus carpio

3 讨论

本研究基于两组高、低体重的黄河鲤群体脑、肝脏、肌肉组织的转录组数据，对生长性状相关基因进行转录组分析，研究发现双向聚类分析可将表达谱相似的基因归纳成簇，不仅有助于推断基因的功能，并且可以有效的发现基因之间存在的调控关系，如果高、低体重组中某样本归类错误，说明该样本数据出现问题，可信度不高，需要剔除。本研究双向聚类热图分析结果显示，高、低两组黄河鲤群体的差异表达基因相关性较高，均一性较好。火山图分析结果可直观、合理的筛选出两样本间的差异表达基因及其上调、下调水平。通过剔除log2FC或q值缺失基因后，火山图分析结果发现了脑、肝脏、肌肉组织中的有效差异表达基因。

在对脑、肝脏组织中差异表达基因的GO和KEGG富集分析，基因主要富集在MAPK、Wnt、胰岛素信号转导、ErbB、肌动蛋白细胞骨架调节等与机体生长、发育、代谢、蛋白质合成相关的信号通路。脊椎动物中生长发育过程主要涉及MAPK、FGF/FGFR、Wnt、Notch、TGF、EDNRA等经典信号通路。Dai等[14]通过研究青鱼(Mylopharyngodonpiceus)在饥饿条件下的基因表达，发现差异表达基因主要富集在胰岛素信号通路、cAMP信号通路、FoxO信号通路、AMPK信号通路、内吞作用和凋亡等信号通路。由于肌肉组织显著差异基因较少，未富集到相应的KEGG通路。候选基因和通路的网络图结果显示，MAPK和Wnt信号通路为富集到候选基因较多的代表性信号通路。

根据文献中脊椎动物生长相关的经典信号通路调控模式和差异基因的功能，结合富集分析、网络图结果，初步筛选出13个参与黄河鲤生长调控通路的候选基因。荧光定量实验结果发现mapk1、mapk7、mapk10、mapk12表达水平上调，这与肿瘤细胞增殖[15-16]和斑马鱼(Daniorerio)[17]发育过程中MAPK家族基因表达趋势一致，证明了本实验结果的可靠性，说明MAPK家族基因在黄河鲤生长过程中发挥了重要调节作用。MAPK家族基因在不同物种中相对保守，是调节生长和新陈代谢的关键基因，参与水生动物细胞的生长、增殖、分化、迁移和死亡等多种生理过程。其通过保守的三级激酶级联反应将细胞外信号传递到细胞内，磷酸化底物蛋白或转录因子发挥作用[18]，且MAPK信号通路在胰岛素样生长因子1(Insulin-like growth factor 1，Igf1)受体的下游。通过Tian等[19]对鲈鱼(Lateolabraxjaponicus)MAPK信号通路的研究可知，根据基因的磷酸化作用位点，可将MAPK信号通路分为三个亚家族，其中mapk1、mapk7基因属于细胞外信号调节激酶(extracellular-signalregulated protein kinase，ERK)亚家族，mapk10基因属于c-Jun氨基末端激酶(JNK)亚家族，mapk12基因属于p38丝裂原活化蛋白激酶(p38 MAPK)亚家族。每个MAPK亚家族磷酸化靶蛋白底物的特定丝氨酸和苏氨酸，并介导生物化学上不同的信号级联，控制细胞对环境的反应，并调节基因表达、细胞生长和凋亡。

在本研究中，fgf13基因表达上调，这与Li等[20]在轴突再生中的研究结果一致，说明fgf13基因上调有助于发育。成纤维细胞生长因子13(fibroblast growth factor 13，fgf13)、成纤维细胞生长因子19(fibroblast growth factor 19，fgf19)基因属于FGF家族，FGF家族基因在发育个体和成体组织中广泛表达，在体内和体外都具有各种生物学活性，包括血管生成、促有丝分裂、细胞凋亡、细胞生存、细胞趋性、细胞黏附、细胞迁移、细胞分化、细胞增殖和组织损伤修复等[21]。FGF13是FGF家族的一种非分泌性蛋白，通过稳定微管在皮质神经元的发育中起着至关重要的作用。然而，却发现同属于FGF家族的fgf19基因在高体重组黄河鲤中表达水平下调，这与Tomlinson等[22]研究结果相似，fgf19转基因小鼠(Musmusculus)通过增加能量消耗降低体质量。fgf19是一种由远端小肠分泌的激素，作用类同胰岛素，可刺激肝脏蛋白质和糖原合成，但不能促进脂肪合成，可调节人体内葡萄糖、脂质和能量稳态，具有减轻体重、增强胰岛素敏感性等作用[23]。因此推测该基因表达下调，将有助于黄河鲤体重增长。

Wnt通路中的wnt8b、wnt10b、wnt16基因表达均上调，这与之前相关研究的结果吻合，Wnt信号通路在所有物种中都很保守[24]。研究发现过表达这些基因可促进癌细胞增殖[25]、成骨细胞分化[26]，说明Wnt家族基因表达的上调，有助于体重增长。此外，发现生长快的黄河鲤群体，tgfbr2基因表达下调。tgfbr2基因即TGF-β II 型受体(transforming growth factor beta receptor II)，属于转化生长因子(transforming growth factor beta，TGF-β)家族，是TGF-β信号传导的关键调节因子。已知tgfbr2基因在各种类型的癌症中起到肿瘤抑制作用，如胃癌[27]和食管癌[28]，说明该基因过表达，将阻碍细胞增殖和分化。本研究中发现，ednra基因在高体重组黄河鲤中表达下调，这与Jacobs等[29]研究结果一致，该研究发现miR-200c可通过下调ednra的表达调控胃癌细胞的增殖、凋亡和侵袭，说明ednra基因的过表达可阻碍细胞增殖。ednra基因广泛地表达于动物的心血管系统、神经系统和胃肠道之中，与细胞增殖、血管收缩舒张、胃肠道动力和腺体分泌等有密切的关系，在小鼠和斑马鱼之间高度保守[30]。Notch信号通路作为进化上保守的信号通路，参与多种细胞发育过程[31]。根据肿瘤类型和细胞生长环境不同，Notch受体的激活可以致癌或抑癌[32]，本研究发现notch2基因的相对表达量下调与黄河鲤体重增长有一定联系。脂蛋白脂酶(lipoprotein lipase，lpl)基因是脂肪酶家族的一员，主要存在于脂肪和肌肉组织中，可在组织损伤、炎症和肿瘤等细胞中可大量表达。虽然LPL蛋白和多肽生长因子结构不同，但均可诱导转录活性蛋白的合成，从而促进转录活性蛋白与多个生长相关基因启动子的关键上调元件结合，介导细胞生长、分化和活动。Goetzl等[33]研究发现lpl可介导溶血磷脂酸(lysophosphatidic acid，LPA)和1-磷酸鞘氨醇(sphingosine 1-phosphate，S1P)表达，而LPA和S1P的功能主要与生长有关，例如诱导细胞增殖、改变细胞分化和存活状态以及抑制细胞凋亡。本研究发现lpl基因在生长较快的黄河鲤脑组织中下调表达，这与以下研究结果相吻合，Richelsen等[34]研究发现生长激素可显著减少肥胖女性腹腔内脂肪组织的含量，同时使lpl活性降低约50%(P<0.01)。

综上所述，本研究对体重差异显著的两组黄河鲤的脑、肝脏、肌肉组织进行转录组测序，挖掘差异表达基因，鉴定与生长相关的候选基因和通路。最终在黄河鲤脑组织中筛选出13个与生长相关的候选基因和MAPK、Wnt、ErbB等信号通路。研究结果为深入挖掘黄河鲤生长性状的分子调控机制提供了研究基础，为加速黄河鲤生长相关的分子辅助育种提供了研究依据。