罗氏沼虾(Macrobrachium rosenbergii)由于食性广、病害少、生长快、营养价值高, 因此具有重要的经济价值, 它是世界上许多地区内陆水产养殖中重要的经济甲壳动物之一(陈雪峰等, 2012)。我国于1976年从日本引进该虾, 经40多年的发展, 我国罗氏沼虾相关产业处于世界第一。据不完全统计, 至2016年, 我国苗种年生产量达到200亿尾, 养殖面积超过3万公顷, 总产量超过15万吨。近年来整个养殖业又呈现出向北方及内陆地区扩展的趋势, 养殖面积的迅速扩增对苗种的稳定供应提出了挑战。然而目前人工育苗过程中往往会出现雌虾产卵相对不集中的现象, 从而导致幼体捕获量参差不齐、育苗周期延长、生产成本加大。掌握与运用性腺发育规律是进行罗氏沼虾的人工繁殖、规模化育苗的基础。

卵巢是甲壳动物的生殖器官, 具有产生成熟卵泡、分泌雌性激素、促进受精等重要功能(陆洋等, 2018)。其发育是复杂的组织发生过程, 包括原始生殖细胞形成、性别决定分化、卵母细胞发育成熟等, 该过程涉及大量重要功能基因的转录调控(谢波波, 2017)。目前为止, 针对罗氏沼虾卵巢的研究多集中于其形状、解剖结构以及相关基因的功能分析(赵维信等, 1995;刘金磊等, 2018), 系统的卵巢发育分子机制和调控机理尚未被完全阐明。

转录组测序(RNA sequencing)是一种可以获得细胞在某一状态下所能转录出来的mRNA、small RNA和non-coding RNA (ncRNA)的高通量测序分析技术(王跃等, 2017)。与传统芯片技术相比, RNA-seq能在全基因组范围内以单碱基分辨率检测和量化转录本, 具有信噪比高、灵敏度高、分析范围广、成本低、重复性好等优势, 从而全面剖析基因表达谱, 发现未知基因, 并精确地识别可变剪切位点及单核苷酸多态性(SNP)位点和非翻译区(UTR)。该技术已经应用于多种经济虾蟹类的性腺发育相研究(Liao et al, 2013;刘金磊等, 2018)。本研究采用RNA-Seq测序对罗氏沼虾卵巢发育的四个代表性时期的卵巢组织进行转录组学研究, 期望通过对卵巢发育过程中细胞进程相关基因的分析, 从宏观上获取与卵巢发育和卵子发生相关的基因, 以全面了解在这四个关键时期卵巢组织中基因的表达变化情况, 从而为分子水平探索罗氏沼虾卵巢发育调控机制提供理论基础。

1 材料与方法 1.1 实验虾

在浙江国家级罗氏沼虾遗传育种中心育种基地的无特定病原(SPF)罗氏沼虾核心育种群体中采集大小规格基本一致的雌虾共30尾, 平均体重约20g。实验虾暂养于水温28℃、盐度0.5%的人工配置海水中。取下每尾虾的卵巢组织分成2份, 1份存放液氮中, 1份存放在组织固定液中备用。

1.2 方法 1.2.1 卵巢组织RNA提取

通过组织切片观察确定不同发育时期卵巢组织(颜素芬等, 2004)。不同发育期卵巢组织RNA的提取使用Invitrogen公司的TRIzol^® Reagent试剂盒进行。抽提RNA满足OD₂₆₀/OD₂₈₀为1.8—2.2, RNA总量≥5μg, 浓度≥200ng/μL, RIN≥8, 进行下一步的文库构建。

1.2.2 文库构建与测序

采用TruSeq RNA Sample Preparation Kit (Illumina)进行卵巢cDNA文库构建。取5μg总RNA, 以带有Oligo dT的磁珠富集mRNA, 加入fragmentation Buffer, 将得到的mRNA随机断裂成200bp左右的小片段。以mRNA小片段为模板, 利用随机引物, 进行反转录, 合成双链cDNA, 加入End Repair Mix补平末端, 在3’末端加A, 连接接头。连接产物纯化后, 进行PCR扩增(15个循环), 扩增产物用2%琼脂糖胶回收200bp左右的目的条带, 得到测序用双末端文库。通过TBS380 Picogreen对文库进行定量, 用Illumina公司的簇生成试剂盒(TruSeq PE Cluster Kit v3-cBot-HS)进行桥式PCR扩增, 簇生成反应后, 用高通量测序平台Illumina HiSeq^TM 4000进行2×151bp测序, 获得的测序数据进行质量控制。

1.2.3 转录组de novo组装

Illumina测序得到的原始图像数据经过Base Calling转化为序列数据, 对原始数据(Raw Date)进行过滤, 去除含测序接头的reads、低质量reads(质量值Q≤10的碱基数占整个read的20%以上)、N率≥10%的reads及长度过短序列, 得到高质量的测序数据(Clean reads)。使用Trinity软件在默认设置下对测序数据(Clean reads)进行组装, 首先将reads拼接成长片段contigs, 然后去除冗余contigs, 将contigs进行分组, 连接成两端不能再延长的转录片段(Unigenes)。

1.2.4 转录组注释

注释前, 利用Trinity软件提供的ORF预测流程对组装得到的所有转录本序列进行基因预测。通过HMMER3 (Eddy, 1998)程序对组装的转录本在Pfam数据库中进行蛋白家族注释。设定期望值E Value < 1e-5, 使用BlastX (版本号: 2.2.25)软件对拼接所得的所有核苷酸序列分别与NR(邓泱泱等, 2006)、String(柳军涛, 2017)、Swiss-prot (Apweiler et al, 2004)、GO (Ashburner et al, 2000)、KEGG (Kanehisa et al, 2004)数据库进行比对获得相应的注释信息。

1.2.5 差异基因的表达分析

使用软件bowtie (Langmead et al, 2009)将clean reads比对到转录组序列上, 使用软件RSEM对bowtie比对结果进行表达量统计分析, 通过RPKM (FPKM)值(Reads/ (Fragments) Per Kilobase of exon model per Million mapped reads), 来衡量基因的表达量水平。使用edgeR软件进行基因表达差异分析。不同样本比较后FDR (False discovery rate)值≤0.05且基因表达量差异倍数以2为底数的对数值|log₂FC(Sample2/Sample1) |≥1, 则该基因为显著差异表达基因。

1.2.6 GO分类与Pathyway富集分析

利用Blast2go软件在GO数据库中, 进行Unigenes参与的生物学过程、细胞组分、分子功能分类注释。利用KEEG数据库获得Unigenes的代谢通路注释信息。

1.2.7 SSR与SNP分子标记鉴定

采用MISA软件对Unigenes进行SSR检测, 设置待检测SSR核心序列的碱基数分别为1、2、3、4、5、6, 对应的核心序列的最少重复次数分别为10、6、5、5、5、5。此外, 以组装好的转录本为模板序列, 利用Samtools和VarScan v.2.2.7软件寻找候选SNP。

2 结果 2.1 转录组测序结果 2.1.1 测序数据统计与组装

根据罗氏沼虾卵巢发育的4个不同时期的卵巢样品, 建立了4个单独的测序文库, 通过Illumina HiSeq^TM 4000高通量测序, 获得卵巢4个期Raw reads序列条数分别为47927784、70553888、57962050、54377814, 经原始数据过滤后, 得到Clean reads共221580604条, 4个时期分别为45860872、67775028、55624580、52320124条。原始数据过滤后, Q20% (Phred数值大于20的碱基占总体碱基的百分比)有明显提高, 达到98%左右。使用Trinity软件对转录组进行de novo组装, 结果表明, 组装得到的Unigene长度范围为201—14658bp, 4个时期样品总共组装得到平均长度为766.29bp的Unigene序列95379条, 平均GC含量为38.96%(表 1)。

所有Unigenes长度分布统计结果表明, Unigenes序列长度主要分布在1—4000bp, 其中长度范围在1—1000bp有77247个, 占80.98%, 长度范围在1001—2000bp的有9409个, 占9.87%, 长度范围在2001—4000bp的有6284个, 占6.59%。大于4000bp的Ungenes仅有2439个, 占2.56%(图 1)。

2.1.2 转录组注释

通过与Pfam、NR、String、Swissprot、KEGG数据库比对, 有52029个Unigenes与已知基因同源, 占整个组装转录本的54.55%, Pfam、NR、String、Swissprot、KEGG数据库中同源序列条数分别为11308(11.86%)、16417(17.21%)、4281(4.49%)、11600(12.16%)、8423(8.83%)(表 2)。物种同源性分析表明, 本研究罗氏沼虾unigenes与NR库中无脊椎动物的基因具有同源性, 与NR数据库比对的16417个unigenes中, 与内华达古白蚁(Zootermopsis nevadensis)同源序列比例最高, 达到10.99%(1804个), 其次是与淡水枝角水蚤(Daphnia pulex)的同源序列比例, 达到6.62%(1087个), 而与Danaus plexippus的同源序列比例最低, 为0.98%(161个)(图 2)。

2.2 差异表达基因

对罗氏沼虾卵巢发育4个时期转录组测序文库进行两两比较分析, 根据不同样本比较后FDR(False discovery rate)值≤0.05且基因表达量差异倍数以2为底数的对数值|log₂FC (Sample2/Sample1) |≥1获得各个时期的基因差异表达情况, 共检测到6605个差异表达基因, 其中显著上调基因2410个, 显著下调基因4195个。Ⅱ期与Ⅲ、Ⅳ、Ⅴ期比较, 表达量差异显著的基因数目分别为1217(上调329, 下调888)、1435(上调1080, 下调355)、207(上调87, 下调120)个; Ⅲ期与Ⅳ、Ⅴ期比较, 表达量差异显著的基因数目分别为350(上调179, 下调171)、1508(上调328, 下调1180)个; Ⅳ、Ⅴ期比较, 表达量差异显著的基因数目为1888(上调407, 下调1481)个。

2.3 GO功能分类

GO是一个标准化的基因功能分类体系, 对基因和蛋白功能进行统一的限定和描述。利用GO数据库, 可以对1个或1组基因按照其参与的生物学过程(Biological process, BP), 基因的分子功能(Molecular function, MF)和细胞组分(Cellular component, CC)三方面进行分类注释。对全部序列进行GO功能分类, 共有6422条unigenes被归类到59个功能分类中(图 4), 这些功能主要涉及代谢相关、生长发育调控、凋亡、生物合成、免疫防御、分子加工、信号转导及转录调控等。

按生物学过程(Bilogical processes, BP)进行差异基因分类注释, 结果6个比较组分别有82、117、51、15、94、113个GO类别出现富集, 按细胞组分(Cellular components, CC)进行分类注释, 结果6个比较组分别有13、14、1、1、19、25个GO类别出现富集, 按分子功能(Molecular functions, MF)进行差异基因分类注释, 结果6个比较组, 分别有55、79、15、5、73、99个GO类别出现富集。差异基因在糖类转运、代谢、氨基酸合成等生物学过程有显著变化, 影响分子功能主要体现在酶活性, 影响细胞组分主要涉及细胞膜的组成(表 3)。

2.4 KEGG通路富集分析

KEGG是基因组破译方面的公共数据库, 该数据库从分子水平信息, 尤其是大型分子数据集生成的基因组测序和其他高通量实验技术的实用程序数据库资源中系统分析基因功能、联系基因组信息和功能信息的大型知识库。通过KEGG通路注释, 罗氏沼虾卵巢发育4个时期共有8423(8.83%)条unigenes被注释, 共涉及184种代谢途径, 基因数量排名前20位的KEGG途径如图 5所显示, 其中unigenes数量最多的为核糖体(ribosome)(308条, 3.66%), 其余依次为RNA转运(RNA transport)(171条, 2.03%)、帕金森氏病(Parkinson’s disease)(165条, 1.96%)、溶酶体(Lysosome) (156条, 1.85%)、吞噬体(Phagosome)(151条, 1.79%)等, 参与基因数最少的为钙离子信号通路(calcium signaling pathway)(94条, 1.12%)。在184个代谢通路中, 有7个代谢通路可能与卵巢发育相关(表 4), 分别为视黄醇代谢(Retinol metabolism)、孕酮介导卵母细胞成熟(Progesterone-mediated oocyte maturation)、胰岛素信号转导(Insulin signaling pathway)、细胞周期(Cell cycle)、Wnt信号通路(Wnt signaling pathway)、卵巢类固醇激素分泌(Ovarian steroidogenesis)、mTOR信号通路(mTOR signaling pathway)。

将6个比较组的差异基因进行KEGG Pathway通路富集分析, 结果显示, 6个比较组分别有78、210、67、34、220、223个通路出现富集, 6个比较组中富集水平最高的通路分别为药物代谢细胞色素P450通路(Drug metabolism-cytochrome P450)、药物代谢细胞色素P450通路(Drug metabolism-cytochrome P450)、氨基丁酸神经元突触(GABAergic synapse)、鞘糖脂生物合成(Glycosphingolipid biosynthesis-globo series)、氮素代谢(Nitrogen metabolism)、药物代谢细胞色素P450通路(Drug metabolism-cytochrome P450), 此外, 还有细胞色素P450异物代谢(Metabolism of xenobiotics by cytochrome P450)、视黄醇代谢(Retinol metabolism)、谷氨酸能突触(Glutamatergic synapse)、缝管隙连接(Gap junction)等通路出现在富集前5的通路中(表 5)。

2.5 SSR与SNP多态性分析

见图 6和图 7。共查找到18592个由1—6个核苷酸重复序列组成的SSR位点, 占Unigenes总数的比例为19.49%。单核苷酸重复类型比例最高, 达52.71%, 随后依次为二、三、四、五、六核苷酸重复类型, 占比分别为32.61%、15.90%、0.80%、0.05%、0.04%。在单、二、三核苷酸重复类型中, 数量最多的重复基元分别为A/T (9407个)、AG/GA (2244个)、AAT/ATA/TAA (344个), 数量最少的重复基元分别为C/G (392个)、CG/GC (12个)、CGG/GGC/GCG (15个)。

在整个转录组序列中, 共寻找到111163个候选SNP位点, 碱基转换位点80379个, 占72.31%, 其中A/G转换40126个, A/G转换40253个; 碱基颠换位点30784个, 占27.69%, A/T、A/C、T/G、C/G颠换位点分别为8149、6661、6397、9577个。

3 讨论 3.1 罗氏沼虾转录组测序

转录组测序技术作为一种快捷、高效的分子生物学研究手段, 已经被广泛应用于重要经济水产甲壳动物的研究中(Liao et al, 2013;刘金磊等, 2018)。有关罗氏沼虾的转录组测序研究工作也已陆续开展, Jung等(2016)对罗氏沼虾肌肉、卵巢、精巢等组织制成混合样品进行了转录组测序, 获得了787731个reads, 经过拼接组装得到123534个Unigenes, 并对参与肌肉生长的相关功能基因进行了挖掘; Rao等(2015)对弧菌感染与未感染的罗氏沼虾肝胰腺组织进行了转录组测序分析, 组装后, 对照组得到59050个Unigenes, 感染组得到73946个Unigenes, 显著表达差异基因14569个, 其中显著上调基因11446个, 显著下调基因3103个; Thanh等(2009)对养成规格大小具有明显差异的2个养殖场的罗氏沼虾的肝胰腺进行了转录组测序, 组装获得66953个isogene, 2组样本鉴定得到了147个差异表达基因。本研究对4个时期卵巢样品进行转录组测序, 总共组装得到平均长度为766.29bp的Unigenes序列95379条, 长度范围在1—1000bp占80.98%, 长度范围在1001—2000bp的占9.87%, 长度范围在2001—4000bp的占6.59%, 长度大于4000bp的占2.56%, 随着序列长度的增加, Unigenes分布数量呈逐级递减, 没有明显中断, 表明本次转录组测序质量较高。本研究首次对罗氏沼虾卵巢发育4个时期的卵巢样品进行了转录组测序, 组装后, 进行了GO、KEEG、差异基因筛选等分析, 进一步丰富了罗氏沼虾转录组序列数据库, 为研究罗氏沼虾卵巢发育的分子机理提供了非常重要的信息。

3.2 转录组文库注释分析

本研究对罗氏沼虾卵巢组织Unigenes进行生物学功能注释, 有52029个Unigenes与已知基因同源, 注释率超过一半, 占整个组装转录本的54.55%, 与NR数据库比对结果显示, 有接近20%(17.21%)的序列同源性, 相对而言, 本文测定的转录组序列注释率较高, 原因可能是由于相近物种的基因注释信息在现有的国际公共数据库中收录较多, 尤其是罗氏沼虾自身的转录组测序研究已经得到开展, 从而导致了较高的注释率。同样, 尚有45.45%的序列没有注释成功, 这些基因序列与已有数据库中基因序列相似性较低, 可能为非编码RNA, 或是无功能结构域序列, 这些基因可能参与重要的生物学过程, 有待进一步研究。GO、KEEG等数据库是为系统、深入了解基因的功能、基因序列同源性及功能基因的代谢通路信息而开发的。本研究测定的转录组序列经COG功能分类, 共有2719个unigenes被归类到预测功能、蛋白质翻译、核糖体结构和生物合成相关基因等25个功能类别中; 通过GO功能分类, 共有6422条unigenes被归类到代谢相关、生长发育调控、凋亡、生物合成等59个功能分类中; KEGG通路注释, 共有8423(8.83%)条unigenes被注释, 共涉及核糖体、RNA转运、帕金森氏病等184种代谢途径。这些功能基因及重要代谢通路途径的发现为后期在分子水平上对罗氏沼虾重要生物学功能的机制机理研究提供了较为全面的基础数据。

3.3 差异表达基因分析

经济甲壳动物的卵巢发育调控机制一直是研究热点, 已有研究证实, 罗氏沼虾的卵巢发育受神经多肽、类固醇激素、神经递质以及外源环境因子的综合调控(姚泊, 1991;张艳, 1999;黎东, 2013;吴维福等, 2014), 本研究共检测到6605个差异表达基因, 其中显著上调基因2410个, 显著下调基因4195个, Ⅱ期与Ⅲ期比较有1217个差异表达基因, Ⅲ期与Ⅳ期比较有350个差异基因, Ⅳ期与Ⅴ期比较有1888个差异表达基因, 这表明参与卵巢发育过程的调控基因较多, 多种因素的共同作用影响控制卵巢发育的相关基因的表达与抑制, 从而推动整个卵巢发育的生物学过程(陈娈娈, 2014)。差异基因GO功能分析表明, 糖类转运、代谢、氨基酸合成等生物学过程有显著变化, 影响分子功能主要体现在酶活性, 影响细胞组分主要涉及细胞膜的组成。差异基因KEGG富集结果显示药物代谢细胞色素P450通路、氨基丁酸神经元突触、鞘糖脂生物合成、氮素代谢等通路富集水平较高, 这些代谢通路与蛋白质合成、脂类合成及神经递质转导有关, 在卵巢发育过程中的重要作用有待进一步研究。

3.4 卵巢发育相关通路

本次转录组测序中, 视黄醇代谢、孕酮介导卵母细胞成熟、胰岛素信号转导、细胞周期、Wnt信号通路、卵巢类固醇激素分泌、mTOR信号通路等7个与卵巢发育相关代谢通路被发现, 其中细胞周期通路中的一些关键基因家族如Cdk基因家族、Cyclin基因家族调控卵巢发育在其他虾蟹类的研究中已有报道, Cyclin基因家族通过磷酸化Cdk基因家族的基因, 从而调节细胞周期, 郭松等(2014)研究表明, Cdk8基因对斑节对虾的卵巢发育启动具有重要作用, 贾复龙等(2016)研究发现在三疣梭子蟹卵巢发育过程中, Cyclin H在卵巢发育Ⅰ期与Ⅱ期的表达量显著高于其他各期, 房君江(2008)分析了罗氏沼虾cyclin B基因在卵子发生过程中的表达情况, 表明了cyclin B基因可能与卵原细胞增殖和卵母细胞减数分裂成熟存在密切关系。此外, 孕酮介导卵母细胞成熟通路中的cdc2、cdc25等基因在虾蟹类卵巢发育中的作用研究也有报道(周俊等, 2012; Chen et al, 2013)。视黄醇代谢、胰岛素信号转导、Wnt信号通路、卵巢类固醇激素分泌、mTOR信号通路等的关键基因对卵巢发育的调控在其他生物中已有较多研究报道, 这些通路中的基因在罗氏沼虾卵巢发育过程中如何调控关键生物学过程有待进一步的研究。

4 结论

采用Illumina HiSeqTM 4000高通量测序平台对罗氏沼虾卵巢发育四个时期卵巢组织进行转录组测序, 获得95379个Unigenes, 分别对四个样品测序文库进行两两比较, 共检测到6605个差异表达基因, 对差异表达基因进行了GO功能分类及KEGG通路注释, 同时对测序序列中的SSR与SNP分子标记进行了鉴定。研究结果为进一步探索罗氏沼虾卵巢发育机制提供了基础信息。