基于RNA-seq 技术的罗氏沼虾(Macrobrachium rosenbergii)不同组织转录组比较分析*

李喜莲 顾志敏 慎佩晶 徐 洋 张宇飞 高 强 程海华 陈雪峰

(农业农村部淡水渔业健康养殖重点实验室 浙江省淡水水产遗传育种重点实验室 浙江省淡水水产研究所 湖州 313001)

罗氏沼虾(Macrobrachium rosenbergii), 又称淡水长臂大虾、马来西亚大虾, 动物学分类上属于节肢动物门(Arthropoda), 软甲纲(Malacostraca), 十足目(Decapoda), 长臂虾科(Palaemonidae), 沼虾属(Macrobrachium), 在各种类型的淡水和半咸水中都能生活, 易养殖且体型肥大, 肉质鲜美。因其适应能力强、生长速度快, 罗氏沼虾已成为我国重要的淡水养殖虾类之一。

近年来对罗氏沼虾的研究主要包括遗传多样性研究, 性腺发育、基因克隆及生长等方面的研究。其中在多样性研究中, 以微卫星为研究对象的如吕敏等(2019)研究异型雄性罗氏沼虾遗传多样性, 戴习林等(2017)对种群 SSR 分析中样本量及标记量对遗传多样性指标的影响进行分析; 或对各不同群体进行遗传多样性研究(董丁健等, 2020; 冯艺, 2018; 孙成飞等, 2015); 周晓敏等(2020)选取60 个SNP 位点对养殖群体和选育群体进行多样性研究。罗氏沼虾雌激素相关受体(Estrogen-related receptor, ERR)成为研究卵巢发育的一个重要基因(赵苗鑫等, 2017; 刘金磊,2018; 刘金磊等, 2018)。药物对性腺分化及发育的影响研究主要包括壬基酚(薛海波, 2010; 李郁娇, 2011;朱春华等, 2017; Guo et al, 2019)、十氯酮(Lafontaine et al, 2016, 2017)及三丁基锡(薛海波, 2010; 李郁娇,2011; 吴维福等, 2013)这三种药物, 以期研究对罗氏沼虾卵黄蛋白原基因表达以及性腺发育的影响。Stalin 等(2019)研究了钴-60 射线对淡水对虾生殖障碍的影响; Tan 等(2019)研究了罗氏沼虾性别逆转与雄激素腺(AG)的关系。罗氏沼虾基因研究集中在免疫相关(刘伟利等, 2017; 江婷佳, 2017; 邱庆庆等, 2019)、性别相关(俞炎琴, 2013; 姜建萍等, 2019; Abayed et al, 2019)、生长相关(叶成凯等, 2019; 邱庆庆, 2019;杨光等, 2020; Dong et al, 2020)和酶类(田荣等, 2016;张夏青等, 2016; 卢志杰等, 2019)。罗氏沼虾生长研究集中在药物对生长的影响(程安达等, 2019; 潘忠超等,2019; Tadese et al, 2020)、微生态制剂(朱光来等, 2019;赵臣泽, 2019)、生长环境(戴习林等, 2016; 陈建酬等,2017; 何竺柳, 2018; 朱其建等, 2019; 张俊功, 2019;Chen et al, 2019; Dong et al, 2020; Manickam et al,2020)、饲料中不同配比(杨树浩等, 2018; 杨景丰等,2019; 张剑伟等, 2019; 单凡等, 2019; 黄黎明等,2019; Feng et al, 2019)和不同品系生长对比试验(蒋飞等, 2013, 2014, 2016)。

随着高通量测序技术的发展, 转录组技术也应用于罗氏沼虾的研究中。郭梁等(2018)利用高通量测序技术和数字基因表谱等技术对感染螺原体的罗氏沼虾进行免疫通路相关基因及其差异表达分析, 获得转录本43405 个; 严赛峰等(2018)、李俊杰等(2018)和邓泽森等(2018)对感染螺原体的罗氏沼虾高通量测序结果开展了 SSR 位点、SNP 位点和重要免疫通路相关 microRNA 研究和分析; 李威霖(2018)对罗氏沼虾肢体再生组织进行了转录组及其差异表达基因研究, 得到 Unigenes 总数目为 87783; 王传聪等(2018)对罗氏沼虾肝胰腺组织转录组测序结果进行 SSR 检测与分析, 获得15356 个SSR 位点; 陈雪峰等(2019)采用Illumina HiSeqTM4000 高通量测序研究罗氏沼虾卵巢发育四个时期卵巢组织的差异, 卵巢发育四个时期共获得 95379 个 Unigenes。Jiang 等(2019a)通过高通量测序技术研究了罗氏沼虾雌性(ZW)、雄性(ZZ)和超雌个体(WW)的性腺差异。Pasookhush 等(2019)采用转录组技术分析了罗氏沼虾幼虾感染新型冠状病毒的反应。Jiang 等(2019b)对未成年雌虾和成年雌虾的眼柄组织进行转录组分析, 共获得 53878 个Unigenes; Cao 等(2017)研究罗氏沼虾正常样本和WSSV 感染样本的淋巴组织转录组, 分别获得73658和72374 个Unigenes; Rao 等(2016)分别研究了罗氏沼虾正常样本和 WSSV 感染样本的肝胰腺转录组, 共获得 63584 个 Unigenes。

本研究以罗氏沼虾7 个组织(眼柄、肝脏、卵巢、鳃、心脏、肌肉、精巢)为研究对象, 采用新一代高通量测序技术进行转录组测序分析, 通过 Trinity 软件组装、数据库功能注释、基因表达差异分析和SSR、SNP 位点筛选, 获得罗氏沼虾分子遗传信息, 以期为进一步研究罗氏沼虾遗传多样性、功能基因及基因表达差异提供理论数据, 同时为深入研究罗氏沼虾生长阻滞现象提供一定的数据基础。

1 材料与方法

1.1 实验材料

实验用罗氏沼虾(Macrobrachium rosenbergii)样品取自浙江省淡水水产研究所八里店综合试验基地,体重为(5.34±1.22) g, 体长为(6.13±0.42) cm。样品经解剖后取各组织(眼柄、肝脏、卵巢、肝脏、心脏、肌肉、精巢), 共7 个组织。每个组织取自3 个个体(见表 1)。将分装好的样品迅速投入液氮中速冻, 置于-80°C 冰箱中保存, 干冰运输。

表1 组织及样品名称列表Tab.1 List of tissue and sample names

1.2 方法

1.2.1 RNA 提取 各组织样品均选用3 个个体的混合样组织, 使用 TRIzol®试剂按照制造商的说明从各组织中提取总 RNA (Invitrogen), 使用 DNase I(TaKaRa)去除基因组DNA。使用Nanodrop ND-2000分光光度计(美国赛默飞)、Aglient 2100 分析仪器对总RNA 的纯度、浓度和完整性进行检测。RIN 值＞7的RNA 用于下游实验。

1.2.2 mRNA-seq 文库构建和Illumina 测序 使用mRNA-seq 样品制备试剂盒(Illumina, San Diego, CA)按照试剂说明书步骤构建了 mRNA-seq 文库和Illumina 测序文库。

1.2.3 质量控制和从头转录组组装 使用Fast QC程 序 (http://www.bioinformatics. babraham.ac.uk/projects/fastqc/)检查Illumina Hiseq 2500 平台在校准前产生的原始序列读取的质量。低质量读数低于阈值质量20; 将长度小于50 bp 的reads 以及包含适配序列、ploy-N 和来自原始数据的测序引物的reads 去除,得到干净的reads同时, 对质控数据的误差率%、Q30、GC-含量%和sequence 重复水平进行了评价。所有的后续分析和注释都依赖于高质量的clean reads。

使用Trinity (v 2.8.5)软件(https://sourceforge.net/projects/trinityrnaseq/)对质控后的数据进行从头组装,获得转录本序列和Unigenes 序列。最后, 只有长度≥300 bp 的unique contigs 才能用于组装后的下游研究。

1.2.4 生物信息学分析 功能注释。拼接得到的Unigenes 序列, 使用 BLASTX 比对(BLAST+2.7.1, 比对标准: E 值不大于 1e-5)与非冗余蛋白数据库(Non-Redundant Protein Sequence Database, NR)、SWISS-PROT 蛋白质序列数据库(Manually Annotated and Reviewed Protein Sequence Database, SWISS-PROT)、基因本体数据库(Gene Ontology, GO)、直系同源蛋白数据库(Cluster of Orthologous Groups of Proteins,COG), 京都基因和基因组百科全书 KEGG 数据库比对。使用Trinity 软件自带的ORF 预测模块进行开发阅读框(ORF)预测。

1.2.5 差异表达分析 通过 Bowtie2 软件将 reads映射到组装的 Unigenes 序列, 使用 RSEM 软件根据比对结果来计算特定转录本的表达量水平。衡量基因表达水平的标准为 RPKM 值(Reads per kilobase of exon model per million mapped reads), 即每一百万条序列中, 每个基因以一千个碱基为单位, 比对上的reads 个数。

使用DE-Seq 软件分析各个组织之间的差异表达,从而找到差异基因组。以差异倍数(Foldchange)＞ 2 和假发现率(FDR)调整显著性值≤0.05 为判断Unigenes表达显著性的阈值。利用GO、egg NOG、KEGG 或thology (KO)和KEGG 通路富集分析对DEGs 进行分类, 并对潜在的生物学途径中的DEGs 进行评价。GO和KEGG 途径富集分析中P 值小于0.05 的过程、功能或成分在DEGs 中被认为存在显著差异。

2 结果

2.1 原始测序数据质控及统计

罗氏沼虾7 个不同组织测序共得到344151324 条原始序列, 质控后得到有效RNA-seq 311475706 条。每个个体的测序量为6.65—9.04 Gb, 平均7.38 Gb。去除接头序列, 截去连续4 个碱基平均质量值低于20的部分, 舍去长度少于 50 bp 的 reads。共得到311475706 条reads。过滤后每个个体的测序量为5.3—8.45 Gb, 平均6.60 Gb(见表2)。每个个体的测序量为每4 个碱基的平均质量值均大于20。

表2 质控后数据统计表Tab.2 Statistics of data after quality control

2.2 转录本拼接

转录本拼接后共获得 95220 个 Unigenes, 总Unigenes 长度为 101401098 bp。总的转录本数量为145717 个, 总的转录本长度为 207379988 bp。每个Unigenes 的平均长度为 1064.9 bp, 最长的 Unigenes长为36137 bp, N50 值为1553。大部分功能基因长度分布在401—600 bp, 占比达到34.27%(见图1)。

2.3 功能注释

在NR、GO、COG、KEGG、SWSS-PROT 五个公共数据库中对获得的 95220 个 Unigenes 进行功能注释, 其余未在NR 数据库中找到的Unigenes 可能为新的蛋白。

根据GO 数据库, 总共有18485 个基因被归类到三个主要的GO 类别中: 生物过程、分子功能和细胞成分。其中, 以“细胞”(14938)、“细胞部分”(14925)、“细胞器”(10827)、“生物调节”(9861)为主(见图 2)。

根据COG 功能分类分成26 类, 其中数量前三位的为功能未知 2172 个(13.56%), 信号转导机制 2070个(12.92%), 翻译后修饰, 蛋白质转换, 伴侣1660 个(10.36%)(见图 3)。

图1 组装序列长度分布图Fig.1 Assembly sequence length distribution diagram

图2 GO 统计二级统计图Fig.2 Secondary chart of GO level 2

五个数据库 NR、GO、COG、KEGG、SWSS-PROT分别注释到 19881、18485、15798、9147、13684 个Unigenes, 在 NR 数据库中比对注释的信息最多, 注释19881 个Unigenes, 占比达20.88%。在NR、GO、COG、KEGG、SWSS-PROT 五个数据库中都注释到的 Unigenes 有 7848 个(见图 4)。

将基因根据参与的 KEGG 代谢通路分为 5 个分支(见图5): 细胞过程(A, Cellular Processes), 环境信息处理(B, Environmental Information Processing), 遗传信息处理(C, Genetic Information Processing), 代谢(D, Metabolism), 有机系统(E, Organismal Systems)。以上5 个分支中含量的最多的类型分别为: 全局和概率地图(Global and overview maps), 转化(Translation),单组织过程(Signal transduction), 运输和分解代谢(Transport and catabolism)、内分泌系统(Endocrine system)。

2.4 差异基因表达分析

图3 COG 分类统计图Fig.3 COG classification statistical chart

图4 注释信息统计韦恩图Fig.4 Annotate the information with statistical Venn diagrams

7 个不同组织鉴定到的mRNA 数量由多到少排序为: G(鳃) ＞ T(精巢) ＞ H(心脏) ＞ L(肝脏) ＞ E(眼柄)＞ O(卵巢) ＞ M(肌肉)(见图 6)。在 7 个组织中共同表达的基因数为15260 个。

2.5 KEGG 通路分析

对罗氏沼虾转录组的 KEGG 分析显示, 9148 个Unigenes 被注释到330 条KEGG 通路中, 其中代谢途径(1578 个)、核糖体(540 个)、次生代谢物的生物合成(528 个)通路数量居前三(见图7)。

330 条通路中, 其中信号通路包括: PI3K-Akt 信号通路、Rap1 信号通路、催产素信号通路、Hippo信号通路、MAPK 信号通路、胰岛素信号通路、雌激素信号通路、Ras 信号通路、甲状腺激素信号通路、钙信号通路、HIF-1 信号通路、FoxO 信号通路、生成信号通路、Wnt 信号通路、PPAR 信号通路、促性腺激素信号通路、趋化因子信号通路、Adipocytokine信号通路、点样受体信号通路、ErbB 信号通路、mTOR信号通路、Jak-STAT 信号通路、T 细胞受体信号通路、p53 信号通路、肿瘤坏死因子信号通路、VEGF 信号通路、B 细胞受体信号通路、Notch 信号通路、催乳激素信号通路、NF-kappa B 信号通路、Fc epsilon RI信号通路等。

2.6 SSR 分析

在得到的95220 条序列中筛选SSR 位点共找到SSR位点37751 个, 这些位点存在于25520 条序列中, SSR 发生频率为 26.80%。其中单条序列中含多个SSR 位点的序列有7962 条, 复合型的SSR 位点3384 个。SSR 位点中单碱基重复、二碱基重复、三碱基重复、四碱基重复、五碱基重复、六碱基重复个数分别为 14919(39.52%)、14715(38.98%)、7577、488、32 和 20, 其中单碱基和二碱基重复含量居第一和第二位(见图8)。

2.7 SNP 分析

对获得的序列进行 SNP 位点分析, 共发现3228575 个 SNP 标记(见表 3), 其中包括 C:G-＞A:T、C:G-＞G:C、C:G-＞T:A、T:A-＞A:T、T:A-＞C:G、T:A-＞G:C6中类型, 其中T:A-＞C:G 和C:G-＞T:A 这两种碱基替换占比较高, 分别占总数的 33.45%和 33.36%。7 个不同组织中发现的SNP位点数量在鳃(23.12%)上最高, 其次是心脏(17.89%)、精巢(17.51%)和肝脏(16.28%), 在肌肉组织中的SNP 位点数量最少(5.43%)。

图 5 KEGG 注释Fig.5 The KEGG annotation

图6 每个样本鉴定得到的mRNA Upset 图Fig.6 The upstate map of mRNA from each sample

图 7 KEGG 通路列表(前20)Fig.7 The KEGG pathway list (Top 20)

图8 SSR 位点的分布情况Fig.8 The distribution of SSR sites

3 讨论

3.1 转录组技术在其他虾上的应用

随着转录组技术的迅速发展, 转录组技术被应用在水产动物的研究中。其中虾类转录组研究主要集中在差异基因的筛选和候选基因的发掘上。

红、黄和透明3 种纯色米虾皮肤组织的转录组获得 45434 条 Unigenes(林师, 2017)。波纹龙虾肝胰腺和卵巢组织进行了转录组测序获得Unigenes 共74124个(李斌等, 2016)。脊尾白虾成虾样品及胚胎样品共获得 47574 条 Unigenes(孙政, 2012)。

表3 不同SNP 类型在各组织之间的分布情况Tab.3 Distribution of different SNP types among different tissues

红螯螯虾肝脏、精巢和卵巢组织共获得了67369个 Unigenes(李喜莲等, 2019); 次级卵黄发生时期卵巢组织共得到69261 条Unigenes(康鹏飞, 2017)。

日本沼虾正常性成熟的和性早熟的卵巢组织中共获得63336 个Unigenes(江红霞, 2017); 感染白斑综合征病毒(WSSV)个体的肝胰腺转录组共获得 64049个 Unigenes(赵才源等, 2018); 亚硝酸盐胁迫下肝脏共获得19022 个Unigenes, 氨氮胁迫下肝脏组织共获得 63453 个 Unigenes (于杰伦, 2019)。

凡纳滨对虾(Litopenaeus vannamei)五个早期发育时期(卵裂期、囊胚期、原肠期、肢芽幼体期发育至膜内无节幼体期)转录组测序, 共得到 66815 条Unigenes(隗健凯, 2015); 性腺组织获得 Unigenes 65218 个(韦嫔媛, 2017); 低温胁迫下肝胰腺组织获得50921 条 Unigenes (董丽君等, 2019)。

本研究从正常个体 7 个组织转录组文库中共获得99520 个Unigenes, 这个数量较前人对于虾类转录组研究得到的 Unigenes 数量都高, 这可能与分析的组织数量大、覆盖面广有一定的关系。

3.2 功能注释

本研究从七个组织共获得 Unigenes 99520 个,N501553, 平均长度为1064.9 bp。在五个数据库中注释到的Unigenes 共有20368 个, 占到总数的21.39%,这比前人研究的转录组结果Unigenes 的注释率都低。肝胰腺注释率31%和35.31%, 卵巢31.46%和54.44%,眼柄29.3%, 淋巴器官29.46%, 再生肢体基部注释率37.23%。估计与本研究中首次对罗氏沼虾的心脏、肌肉、精巢和鳃等组织进行转录组高通量测序, 引入了较多在这几个组织中特异表达的基因, 这些新获得的序列在以上五个数据库中得不到注释, 从而降低了注释率。还有部分序列未能被注释, 可能与序列的长度有关, 过短的序列也会造成无法注释和分类; 无法注释的另一种可能是近缘物种序列信息的缺乏,导致无法通过同源序列比对得到注释。

3.3 KEGG 通路

对罗氏沼虾转录组的 KEGG 分析显示, 9148 个Unigeness 被注释到 KEGG 数据库中并分布在330 个已知途径中, 与免疫相关的通路如FoxO 信号通路和Rap1 信号通路。FoxO 转录因子家族在细胞生理事件中调控基因的表达, 包括凋亡、细胞周期控制、葡萄糖代谢、氧化应激抵抗和寿命。FoxO 蛋白的一个主要调控机制是对磷脂酰肌醇3 激酶(PI3K)下游的丝氨酸苏氨酸激酶Akt/蛋白激酶B(Akt/PKB)进行磷酸化,这是对胰岛素或几种生长因子的反应。FoxO 信号通路在本研究结果中涉及 93 个 Unigeness。Rap1 是一种小型GTPase, 它控制多种过程, 如细胞黏附, 细胞-细胞连接的形成和细胞极性。Rap1 通过调节多种细胞类型中整合素等黏附分子的功能, 在细胞-细胞和细胞-基质相互作用的控制中发挥主导作用。Rap1 还以高度依赖于细胞类型的方式调控 MAPK 活性。Rap1 信号通路在本研究中共涉及201 个Unigeness。这些结果的获得都将为进一步研究罗氏沼虾抗性相关基因提供科学依据。

4 结论

本研究对罗氏沼虾 7 个不同组织的转录组数据进行分析, 获得99520 个Unigenes, 这比单一组织或较少组织受到环境胁迫或细菌、病毒感染获得数据具有较高的可信度, 这一结果将大大丰富罗氏沼虾的基因数据库资源。与此同时, 20368 个Unigenes 在五大数据库中得到注释。各个组织间差异基因也得到进一步的分析, 本研究还筛选得到大量的 SSR 位点和SNP 位点, 这些分子标记也将在后续罗氏沼虾分析标记的开发中起到重要的作用, 为进一步深入挖掘和开发利用罗氏沼虾功能基因提供参考。