翘嘴鳜(Siniperca chuatsi)是中国特有的重要经济鱼类, 分布于中国除青藏高原外的各大水系, 属于鲈形目(Perciformes), 鮨科(Serranidae), 鳜属。翘嘴鳜具有生长速度快、肉质好、个体大的特点, 是重要的养殖品种^[1]。与此同时, 与翘嘴鳜生长、应激相关基因的克隆和表达也受到关注^[2-3]。Wang等^[4]从GH (growth hormone)、IGF (Insulin-like growth factor)基因上筛选出6个SNP (single nucleotide polymorphisms)位点, 并验证其与翘嘴鳜生长具有显著的相关性。Zhu等^[5]通过饥饿再投喂实验对翘嘴鳜肌肉进行了microRNA分析, 并且筛选出4个在急性肌肉萎缩刺激条件下对肌肉中基因转录、翻译有重要作用的miRNAs位点。了解不同生长性状翘嘴鳜的分子表达模式差异, 可以帮助揭示翘嘴鳜生长调控网络, 筛选分子育种标记, 改良翘嘴鳜品质。

转录组测序可以帮助分析特定组织的基因表达情况, 寻找生物调控网络和机理。通过转录组测序可获得大量基因表达信息, 用于比较不同状态下基因表达量的差异, 获取与目的性状相关的基因表达情况, 从而了解调控目的性状的基因网络。随着测序技术的逐渐成熟, 转录组测序被应用于免疫、性别、生长、应激等多种表型相关的关键基因调控网络研究和重要性状分子标记鉴定, 水产动物中转录组测序研究也越来越多^[6]。Wang等^[7]用转录组测序技术筛选了抗肠道败血病候选基因和疾病相关表达差异显著的单核苷酸多态性位点; Sun等^[8]采用Illumina测序技术筛选斑点叉尾鮰(Ietalurus punetaus)性别相关基因, 并通过进一步的分析发现许多涉及性腺发生、精子发生、性别分化的基因; Jeukens等^[9]利用RNA-seq技术研究了短小型和正常体型白鲑(Salmonidae)的适应性转录组变化, 发现白鲑短小型高表达基因主要与免疫、DNA复制以及能量代谢有关, 而正常体型则与蛋白质合成有关。

鱼类的生长伴随鱼类整个生活周期, 受到基因和环境等多方面因素的共同影响, 由于其调控机制非常复杂, 受到了研究者广泛的关注。而转录组测序可以在没有研究物种基因信息的情况下进行分析, 对其功能基因组进行深入研究。本研究首次使用转录组测序对体重差异翘嘴鳜进行差异基因分析, 为翘嘴鳜生长相关基因的研究提供了完整的基因表达背景。

本研究根据全同胞家系翘嘴鳜出现体型显著差异现象, 通过对体型差异翘嘴鳜个体进行肌肉转录组测序, 对差异表达基因进行注释和通路富集分析, 研究同一家系在相同养殖条件下, 极大个体(Overweight)和极小个体(Underweight)翘嘴鳜肌肉基因表达模式差异, 了解大、小个体间差异基因表达特点, 为深入开发翘嘴鳜生长相关分子标记提供科学依据。

本研究中翘嘴鳜亲鱼来自于佛山市南海区百容水产良种股份有限公司, 经人工授精繁殖后, 获得翘嘴鳜全同胞家系。每天饱食投喂, 养殖水温在26~28℃, pH变化范围在6.5~7.0, 溶氧保持在6 mg/L以上。养殖3个月后, 对所有个体进行称重、测量, 从中选择体型极大个体和极小个体各5尾, 禁食24 h, 使用MS-222将实验鱼麻醉后快速对肌肉组织进行取样。

将样本在液氮中充分研磨, 按照E.Z.N.A.^® Total RNA Kit Ⅱ试剂盒操作规范进行RNA提取, 使用DNase I对提取到的总RNA消化1 h, 去除DNA残留。分别取1 μL RNA样品用核酸蛋白测定仪和Agilent2100对RNA样品的纯度、质量和浓度进行检测, 同时取800 ng RNA于1%琼脂糖凝胶进行电泳检测RNA样品完整性。

分别将翘嘴鳜极大个体和极小个体各5尾的肌肉组织RNA等量混合, 得到两个RNA pools, 用Oligo(dT)磁珠富集纯化mRNA, 以片段mRNA为模板, 逆转录酶和六碱基随机引物合成cDNA第一链, 再用RNaseH和DNA聚合酶合成第二链。使用T4 DNA聚合酶、Klenow DNA聚合酶和T4多聚核苷酸激酶修复DNA末端, T4 DNA连接酶连接测序接头, 通过纯化选择200 bp+25 bp的cDNA片段扩增得到DNA文库。采用双端测序(paired-end)技术, 将构建好的两个测序文库用Illumina HiSeqTM2000进行测序, 获得900 bp左右的双端测序的片段(paired-end reads)。

测序得到的原始图像数据经碱基识别(base calling)转化为序列数据称之为原始数据(raw data)或原始片段(raw reads), 经过对raw reads序列进行筛选, 去除含adaptor、N比例大于5%及低质量reads(质量值Q≤20的碱基数占整个read的50%以上), 得到可供后续分析的过滤后片段(clean reads)。使用短reads组装软件Trinity做转录组从头组装, 使用序列聚类软件TGICL做进一步序列拼接和去冗余处理, 得到不含N的非重复序列基因(Unigene)。用SOAPaligner/soap2将clean reads比对到参考基因组和参考基因序列中进行基因注释。

将Unigene序列与蛋白数据库Nr (nonredundant protein)、Swiss-Prot、KEGG (Kyoto Encyclopedia of Genes and Genomes)和COG (Clusters of Orthologous Groups)做Blastx比对(evalue < 0.00001), 取比对结果最好的蛋白确定Unigene的序列方向。根据Nr注释信息, 通过Blast2GO软件对得到的Unigene进行GO(Gene Ontology)富集。用WEGO软件对所有Unigene做GO分类统计, 将Unigene描述为分子功能、细胞组分、参与的生物过程3个功能分类, 认识该物种的基因功能分布特征。

使用RPKM法(Reads Per Kb per Million reads)计算基因表达量。RPKM法可以消除基因长度和测序量差异对计算基因表达的影响, 计算得到的基因表达量可直接用于比较不同样品间的基因表达差异。对差异检验的P值作多重假设检验校正, 通过控制FDR(False Discovery Rate)来决定P值的域值。FDR值越小, 差异倍数越大, 表明表达差异越显著。差异表达基因评判标准为FDR≤0.001且倍数差异不低于2倍的基因。分别用WEGO软件对差异表达基因进行GO功能分类统计和用KEGG数据库对差异表达基因进行Pathway注释。

提取肌肉中RNA进行逆转录, 以肌肉cDNA为模板, 对转录组结果中6个差异表达显著的部分基因进行荧光定量PCR验证, 包括生长轴的生长激素受体(GHR2)和胰岛素样生长因子受体(IGFR1)、肌肉纤维组成蛋白(Mhc)、涉及蛋白和折叠和运输的热休克蛋白Hsp90(heat shock 90)、涉及蛋白质生成的4ebp、肌肉纤维调节蛋白(MyoD)。以18S、β-actin为内参基因, 通过RT-PCR(realtime fluores-cence quantitative PCR)检测不同基因的表达水平。

本实验在相同养殖条件下, 从同一家系的60尾健康翘嘴鳜个体中挑选极大个体(Overweight)、极小个体(Underweight)各5尾, 极大个体组平均体重(60.2±5.7) g; 极小个体组平均体重(13.2±1.4) g, 极大个体组平均体重可达极小个体组4倍, 体长、体高也分别达到极小个体组两倍, 具有明显的体型差异(表 1)。

表 1(Tab.1)

表 1 翘嘴鳜极大个体和极小个体生长性状 Tab.1 Growth traits of overweight and underweight Siniperca chuatsi 个体individual 体重/g weight 体长/mm length 体高/mm height O-1 62.2 138 46 O-2 80.1 147 52 O-3 55.2 132 43 O-4 54.3 131 45 O-5 59.3 135 44 U-1 12.6 79 24 U-2 7.8 71 23 U-3 16 85 26 U-4 12.9 82 24 U-5 12.8 79 26 注: O表示极大个体, U表示极小个体. Note: O means overweight; U means underweight.

表 1 翘嘴鳜极大个体和极小个体生长性状 Tab.1 Growth traits of overweight and underweight Siniperca chuatsi

经Illumina测序后, 两个翘嘴鳜转录组样本共获得95581082条raw reads, 其中极大个体组48963598条, 极小个体组46617484条; 经过筛选去除低质量序列后, 共得到91264026条clean reads, 平均长度约为100 bp, 占总reads数的95.48%;测序结果中GC含量为51.31%, Q20百分比为98.57%, 具有良好的测序质量, 两个转录组本间无显著数据量差异, 可用于分析两个转录本表达模式差异。经Trinity拼接组装后, 两个转录本共得到73353条Unigene, 平均长度为703 bp(表 2)。在去除短序列和重复序列之后, 经过Nr、Swiss-Prot、KEGG和COG4个数据库共定义39005条Unigene。Nr、COG、Swissprot及KEGG分别定义了38833条、10926条、33749条和19791条Unigene, 确定Unigene序列方向后, 通过定义最多Unigene数量的Nr数据库进行注释(图 1)。

表 2(Tab.2)

表 2 翘嘴鳜极大个体和极小个体转录组测序数据产出及拼接结果汇总 Tab.2 The brief information of deep sequencing and assembly for overweight and underweight Siniperca chuatsi 类别category 数据data 总reads数number of reads 95581082 可用reads数clean reads 91264026 (95.48%) 总碱基数/nt total nucleotides 9126402600 Q20比例/% Q20 ratio 98.57 N比例/% N ratio 0.00 GC比例/% GC ratio 51.31 Contigs总数number of contigs 80314 Contigs平均长度average length of contigs 621 Unigene总数number of unigenes 73353 Unigene平均长度/bp average length of unigenes 703

表 2 翘嘴鳜极大个体和极小个体转录组测序数据产出及拼接结果汇总 Tab.2 The brief information of deep sequencing and assembly for overweight and underweight Siniperca chuatsi

图 1(Fig. 1)

图 1 Nr、COG、KEGG、Swissprot数据库比对结果 Fig.1 Detection of homologous genes in public databases (Nr, COG, KEGG, Swissprot)

计算RPKM的大小衡量基因表达水平的高低, 经过RPKM标准化计算处理后, 共获得8942个差异表达基因, 其中7858个Unigene在极小个体中表达量较高。GO富集将5150个差异表达基因富集到3个分集, 分别是生物过程1530条、细胞组分1737条和分子功能1883条。

在生物过程分集中, 细胞(cell)和细胞成分(cell part)是主要的功能分组; 分子功能分集中, 连接(binding)是主要功能分组; 细胞代谢过程(cellular process)是生物过程分集中的主要功能分组(图 2)。

图 1(Fig. 1)

图 2 翘嘴鳜极大个体和极小个体转录组测序差异表达基因GO数据库比对结果 横轴表示GO注释的基因数目以10为底数的对数. Fig.2 The GO analysis of different expression genes between overweight and underweight Siniperca chuatsi The lateral axis indicates the number of genes annotated by GO after log transformation.

通过翘嘴鳜转录组数据分析得到了大量的差异表达基因, 其中大部分基因在极小个体组表达量更高, 如GHR2 (growth hormone receptor 2)、IGFR1 (type 1 insulin-like growth factor receptor)等(表 3)。通过KEGG分析, 以P < 0.05判定为有显著意义, 共获得25个通路。富集差异表达基因最多的通路是内质网蛋白加工(protein processing in endoplasmic reticulum)、泛素调节的蛋白质降解(ubiquitin mediated proteolysis)以及RNA转运(RNA transport), 分别是119条(3.77%)、101条(3.2%)和91条(2.88%)差异表达基因(表 4)。

表 3(Tab.3)

表 3 翘嘴鳜极大个体和极小个体差异表达基因 Tab.3 The different expression genes between underweight and overweight Siniperca chuatsi 基因编号gene ID 定义annotation 小个体中RPKM值 underweight RPKM 大个体中RPKM值overweight RPKM log2比值log2 ratio P FDR Unigene0032360 生长激素受体-1前体growth hormone receptor isoform 1 precursor 5.014 0.001 -12.292 5.24E-51 3.03E-49 Unigene0032361 生长激素受体-1 growth hormone receptor 1 189.052 37.474 -2.335 3.22E-05 3.10E-05 Unigene0026866 胰岛素样生长因子受体-1 type 1 insulin-like growth factor receptor 1.517 0.001 -10.567 6.09E-07 5.42E-06 Unigene0033081 心肌样肌钙蛋白troponin T, cardiac muscle-like 4.009 0.001 -11.969 1.15E-19 2.87E-18 Unigene0068992 心肌样肌球蛋白轻链-1 myosin light chain 1, cardiac muscle-like 3.163 0.253 -3.646 3.44E-06 2.74E-05 Unigene0032402 原肌球蛋白α链tropomyosin alpha-1 chain 1.783 0.237 -2.909 2.22E-07 2.10E-06 Unigene0021259 快肌样肌球蛋白重链-1 myosin heavy chain, fast skeletal muscle-like 4.05 0.578 -2.808 5.39E-09 6.18E-08 Unigene0002916 肌动蛋白actin 5.893 0.741 -2.992 5.70E-12 8.74E-11 Unigene0068993 肌球蛋白轻链-4亚型X1 myosin light chain 4-like isoform X1 6.423 0.058 -6.789 6.95E-32 2.66E-30 Unigene0072056 热休克蛋白90β heat shock protein hsp90 beta 8.095 1.88 -2.106 6.28E-95 1.35E-05 Unigene0022967 真核翻译起始因子4E结合蛋白1 eukaryotic translation initiation factor 4E-binding protein 1 45.133 5.536 -3.027 3.95E-216 7.09E-07 Unigene0056264 丝/苏氨酸激酶D3激酶X2 serine/threonine-protein kinase D3 isoform X2 0.001 0.736 9.523 2.18E-05 5.91E-214 Unigene0026704 肠道样脂肪酸结合蛋白fatty acid-binding protein, intestinal-like protein 0.001 0.985 9.944 2.90E-06 0.0005445 Unigene0027806 免疫球蛋白超家族成员22 immunoglobulin superfamily member 22 0.176 1.609 3.193 1.97E-07 3.76E-21 Unigene0075065 肌源性调节因子myogenic regulatory factor MyoD1 3.248 6.625 1.029 8.51E-08 0.000379

表 3 翘嘴鳜极大个体和极小个体差异表达基因 Tab.3 The different expression genes between underweight and overweight Siniperca chuatsi

表 4(Tab.4)

表 4 翘嘴鳜极大个体和极小个体差异表达基因通路分析 Tab.4 The pathway analysis of different expression genes between underweight and overweight Siniperca chuatsi 通路编号 pathway ID 通路pathway 差异表达定义基因数量 DEGs with pathway annotation 所有定义基因数量 all genes with pathway annotation P ko03050 蛋白酶体proteasome 42 (1.33%) 73 (0.37%) 5.69E-16 ko03008 真核生物核糖体生物合成ribosome biogenesis in eukaryotes 62 (1.96%) 168 (0.85%) 3.04E-11 ko04141 内质网蛋白质加工protein processing in endoplasmic reticulum 119 (3.77%) 435 (2.2%) 7.21E-10 ko04120 泛素介导蛋白质水解ubiquitin mediated proteolysis 101 (3.2%) 357 (1.8%) 2.01E-09 ko03013 RNA转运RNA transport 91 (2.88%) 383 (1.94%) 4.08E-05 ko04974 蛋白质消化和吸收protein digestion and absorption 44 (1.39%) 177 (0.89%) 0.001422099 ko00230 嘌呤代谢purine metabolism 71 (2.25%) 341 (1.72%) 0.009785118 ko05414 扩张型心肌病dilated cardiomyopathy 80 (2.53%) 398 (2.01%) 0.01530128 ko03040 剪接体spliceosome 71 (2.25%) 353 (1.78%) 0.02097551 ko05410 肥厚型心肌病hypertrophic cardiomyopathy 77 (2.44%) 393 (1.99%) 0.02970835

表 4 翘嘴鳜极大个体和极小个体差异表达基因通路分析 Tab.4 The pathway analysis of different expression genes between underweight and overweight Siniperca chuatsi

RT-PCR扩增产物结果显示, 熔解曲线单一; 实验有限循环值在15~30之间, 且内参基因18S的C_t比其他待测基因C_t值小, β-actin在不同个体间表达量稳定。通过RT-PCR, 基因表达差异趋势与转录组结果相同(图 3)。

图 1(Fig. 1)

图 3 翘嘴鳜极大个体和极小个体差异表达基因荧光定量PCR验证 y轴表示相对表达量, 误差线表示标准差, *表示显著差异(P < 0.05). Fig.3 The RT-PCR analysis of gene expression in overweight and underweight Siniperca chuatsi The y-axis indicates relative expression level. Error bars mean standard deviation. * indicate the differences between the two groups are significant (P < 0.05).

鱼类的生长伴随于鱼类整个生活周期, 与基因、环境等多方面因素有关, 在体内受到生长与生殖调控轴相互作用的影响^[10]。过去的研究中发现, 许多鱼类中都存在不同家系个体具有显著生长差异的特性, 甚至来自同一家系的个体也会出现明显生长差异的现象, 具有显著生长差异的个体成为研究生长差异调控机制的重要对象^[11-13]。Kocmarek等^[13]对同一品系虹鳟在不同产卵期的大个体和小个体肌肉和肝脏分别进行mRNA微阵列分析, 评估不同季节(晚秋vs春初)和体型(大vs小)间转录组水平差异。王博等^[14]通过对生长存在显著差异的翘嘴鳜个体进行肌肉microRNA转录组测序, 发现miR-122、miR-192、miR-451、let-7j-5p共4个表达显著差异的miRNA, 差异表达miRNAs参与生长、发育、代谢等重要途径。与其他鱼类类似, 本研究发现, 来自同一家系的翘嘴鳜在相同养殖条件下生长至3月龄, 出现明显的生长差异现象。通过挑选出极大个体和极小个体各5尾, 极大个体组体重可达极小个体组4倍以上, 体长、提高均可达其2倍, 具有明显的体形差异。

在活鱼体内, 鱼体体重的30%~80%为骨骼肌重量, 是鱼类经济价值的主要体现。鱼类出生后肌肉的生长主要通过纤维细胞的增殖(hyperplasia)和肥大(hypertrophy)达到肌纤维数量、直径和长度的增加, 最终达到生长的目的, 而这一过程受到一系列基因在转录水平上的综合调控^[15]。生长因子、调控蛋白、转录因子都被证实参与了鱼类肌肉形成的调控过程^[16]。本研究通过对极大个体和极小个体肌肉组织进行转录组测序, 共获得8942个差异表达基因, 其中大部分差异表达基因在极小个体中表达量增加。极小个体组包含更多的表达量上调的基因说明极小个体组具有更高的基因表达活性。相似的结果在其他研究中也观察到。Kocmarek等^[13]对不同大小的虹鳟微阵列分析结果显示, 无论在春季或秋季, 小鱼组肝脏和肌肉中包含更多的表达上调的基因, 表达量上调基因涉及蛋白转录、翻译以及细胞凋亡, 造成了小鱼肌肉中较高的细胞周转率, 从而导致体型差异。在具有不同生长率的珍珠贝中, 通过分别对大小个体进行转录组分析发现, 小个体组具有更多表达量上调的基因, 包括与贝壳形成、能量代谢、肌肉生成等有关的基因, 具有小个体追赶生长的趋势^[17]。

本实验中发现大量在极大个体和极小个体间表达量差异显著的基因, GH-IGF通路受体基因GHR2、IGFR1, GH-IGF通路下游基因4ebp (initiation factor 4E-binding protein 1)、转录延长因子, 肌肉纤维组成基因Mhc (myosin heavy chain)、Mlc (myosin light chain 1)、Troponin等在极小个体的翘嘴鳜中表达量升高, 但MyoD1在极大个体和极小个体间表达水平并无显著差异。通过对差异表达基因进行KEGG富集通路分析结果显示, 大量的基因富集在与蛋白生成、消化有关的通路, 与差异表达基因分析趋势一致。GH-IGF是硬骨鱼中重要的生长调控通路, GH促进肝脏分泌IGFs, 与IGFR结合, 通过PI3K/AKT和ERK/MAPK两条通路调控生长, 前者促进蛋白生成, 后者促进细胞增生及增殖^[18]。PI3K通路磷酸化, 调节下游分子4ebp, 进而参与调节转录延长因子、转录起始因子以调控生长。通过在牙鲆饥饿再投喂实验中发现, 硬骨鱼PI3K通路和ERK/MAPK通路在肌肉生长中具有重要作用^[19]。然而GHR、IGFR表达量升高与生长率并非正相关, 4月龄的虹鳟用GH处理后出现生长抑制, 而其GHR1和GHR2表达量则升高; 虹鳟、翘嘴鳜等鱼类在经过饥饿后, GHR2、IGFR1表达量升高。与本研究结果一致, 在受到饥饿或生长抑制时, 虽然生长率与GHR2、IGFR1基因的表达成反比趋势, 可能参与了不利条件下肌肉能量储存代谢通路的激活^[20], 同时有观点认为GHR2可能促进了tGHR (truncted GH receptor)的生成, 而tGHR诱导了某些分解反应, 从而限制了GH信号通路^[21], 小个体翘嘴鳜中, GHR2表达量升高符合GH信号通路被抑制的趋势。Mhc和Mlc是两种重要的组成肌肉纤维的蛋白链, 肌钙蛋白(Troponin)是间质细胞决定分化为肌细胞后, 所表达的肌肉特异性蛋白, 参与调节肌肉收缩和舒张, Troponin表达量增加利于增强肌肉收缩, 进而通过促进肌纤维的合成^[22]。在本研究中, Troponin在小个体中具有更高的表达量。与凡纳滨对虾研究结果一致, 李喜莲通过对不同个体凡纳滨对虾肌肉中表达基因进行抑制性消减杂交发现, Troponin在小个体中的表达量明显高于大个体^[23]。

肌肉中大量组成成分是蛋白质, 鱼类的生长与蛋白合成及消化紧密相关。KEGG聚类结果显示, 富集差异表达基因最多的通路是内质网蛋白加工(protein processing in endoplasmic reticulum)、泛素调节的蛋白质降解(ubiquitin mediated proteolysis)以及RNA转运(RNA transport), 分别是119条(3.77%)、101条(3.2%)和91条(2.88%)差异表达基因(表 4)。RNA转运是细胞基因表达的基础, 而内质网蛋白质加工则帮助蛋白质的正确折叠, 并对错误折叠的蛋白质进行降解, 是蛋白生成的基础。与差异表达基因趋势一致, 在小个体翘嘴鳜中大部分基因表达量更高, 具有更活跃的表达模式。泛素调节的蛋白质降解通路是蛋白质降解的重要途径, 通过与酶结合后产生级联反应达到蛋白质降解的目的^[24]。本研究中小个体翘嘴鳜蛋白质降解通路表达量增加, 蛋白质积累减少, 具有更低的生长率。其他鱼类研究显示, 智利羽鼬鳚(Cenypterus chilensrs)在应激条件下, 泛素调节的蛋白质降解具有更高的基因表达, 也证实了在应激条件下, 鱼类的低生长率与蛋白质合成代谢有关^[25]。

综上所述, 来自同一家系翘嘴鳜个体在相同养殖条件下仍然产生明显的体型差异, 通过对极大个体、极小个体肌肉组织进行转录组测序, 比较基因表达差异发现, GHR2、IGFR1、4ebp、Mhc、Mlc、Troponin等多个基因在极小个体中表达量增加, 而大量差异表达基因富集在蛋白生成和消化相关通路, 可能参与了翘嘴鳜生长调控。本实验通过转录组测序获得了大量翘嘴鳜转录本信息, 研究结果为翘嘴鳜更深入的基因调控网络研究提供丰富的基因资源。