孙永峰，武惠岩，王子遒，路洪涛，李书哲，杨 威，隋玉健

（吉林农业大学动物科学技术学院，吉林长春 130118）

鹅等级前卵泡组织转录组测序分析

孙永峰*，武惠岩，王子遒，路洪涛，李书哲，杨 威，隋玉健

（吉林农业大学动物科学技术学院，吉林长春 130118）

为研究鹅等级前卵泡组织的转录组差异，丰富鹅卵泡转录组数据信息，实验选用3只经产处于高峰期的籽鹅等级前卵泡，通过Illumina HiSeq 2500 测序平台进行测序，并进行序列分析和注释。经测序质检后共得到224 929 214条序列，31 722 729 162 bp碱基，且质量大于20的碱基比例高于95 %，从头拼装后得到114 486 条Unigenes序列和145 020条Transcripts序列，通过比对分析共有42 453条序列注释到NR数据库，并进行物种同源性比对，发现绿头鸭与其同源性最高，通过GO数据库比对，发现共有1 045条Unigenes注释到GO数据上，主要涉及生物过程、分子功能和细胞组分。结果表明，共有 11 949 条Unigenes 比对到KEGG数据库，共涉及到341条通路，其中癌症通路、PI3K-AKT信号通路和MAPK信号通路显著富集。为验证基因表达差异，以FPKM值作为基因的表达量，FDR＜0.05，︱logFC︱≥1为标准筛选差异基因，其中以SY 等级卵泡为参考，在 PE等级卵泡组织间中得到15 516条差异基因，其中上调基因8 293条，下调基因7 223条。并对筛选得到的所有差异基因进行GO和KEGG比对。本实验条件下，测得各等级卵泡之间基因表达差异较大，且癌症通路、PI3K-AKT信号通路和MAPK信号通路较为活跃，可间接证明这3条通路对鹅等级卵泡发育影响较大，这对进一步分析鹅等级前卵泡发育机制提供基础。

鹅；等级前卵泡；转录组测序；差异表达基因；功能分类

转录组测序技术（Transcriptome Sequencing）是指利用第2代高能高通量测序技术记性cDNA测序技术，可全面快速地捕获某一物种特定器官或组织在某一状态下所有转录本。可从整体水平研究基因功能以及基因结构，发现不同生理或病理状态下细胞、组织或个体间差异表达基因。目前随着高通量测序技术的发展，基因研究速度加快，对畜禽卵泡、卵巢组织研究逐渐多了起来。吴阳升等[1]对绵羊小卵泡与中卵泡转录组差异特征分析，筛选出1 073条差异基因，其中上调基因744个，下调基因959个。朱志明等[2]对山麻鸭开产期和产蛋高峰期卵巢组织转录分析，得到1 929条差基因，通过GO和KEGG富集分析GnHR信号通路和TGF-β信号通路在卵巢的生理活动中发挥重要作用。韩坤鹏等[3]对京海黄鸡卵巢转录组研究后与所选参考基因组序列对比发现4 431个新基因，其中1 809个基因通过与各个数据库比对得到功能注释。Tariq等[4]对鹅全血进行Denovo 转录组分析共得到211 198条功能基因（Unigenes）。

鹅（Goose）属于鸟纲雁形目鸭科动物，杂食性水禽，是一种重要的经济动物，品种约150多种，广泛分布于世界各地。中国作为鹅品种资源最为丰富的国家，饲养着世界约90%的家鹅[5]。但鹅业发展一直落后于其他家禽业，这主要是由于鹅的卵泡发育与其他家禽不同，具有严格的等级性和选择性[6]，产蛋性能低下。为研究鹅卵泡的发育机制，本实验利用高通量RNA-seq，以鹅等级前卵泡组织为研究对象，对其转录组进行测序分析。通过分析鹅卵泡发育过程中的转录组基本信息，来探讨研究鹅卵泡的发育机制，并进一步研究和发掘卵泡发育过程中的差异基因及通路，为鹅品种优良性状的改良提供理论依据。

1 材料与方法

1.1 实验材料 实验选择经产处于产蛋高峰期籽鹅3只（由吉林农业大学鹅研中心种鹅繁殖基地提供）。实验药品及仪器有TRIzol®Reagent，RNA Purification Reagent，磁力架，TBS380 Picogreen，由Invitrogen公司生产；TruseqTMRNA sample prep Kit，cBot Truseq PE Cluster Kit v3-cBot-HS，Hisseq2000 Truseq SBS Kit v3-HS，Miseq Reagent Kit V2，由Illumina生产；Certifed Low Range UItra Agarose 由Bio-Rad生产。

1.2 实验方法

1..2.1 实验样品的采集及等级划分 将籽鹅用乙醚麻醉后放血处死，迅速将卵泡取出并用预冷生理盐水清洗干净，同时根据卵泡直径大小对其进行等级划分。同时去除卵泡外层血管膜，并释放卵泡液。清洗干净后迅速放入液氮中保存。

1.2.2 卵泡等级划分 鹅卵泡可根据卵泡的发育将其分为等级卵泡和等级前卵泡。本实验所需卵泡为等级卵泡，等级卵泡一般由大到小分为F1、F2、F3、F4、F5、F6共6 个等级。等级前卵泡根据其直径大小可分为初级卵泡（Primary Follicles，PE）直径小于2 mm、小白卵泡（Small White Follicles，S）直径在2～4 mm、中白卵泡（Middle White Follicles，M）直径在4～6 mm、大白卵泡（Large White Follicles，L）直径在6～8 mm、小黄卵泡（Small Yellow Follicles，SY）直径在8～10 mm[7]。

1.2.3 总RNA提取及cDNA文库的建立 利用Trizol法将实验所采集的卵泡组织进行总RNA提取，通过Nanodrop 2000 对所提RNA 的浓度和纯度进行检测，琼脂糖凝胶电泳检测RNA 完整性，Agilent2100 测定RIN 值。将所提取的RNA通过Oligo（dT）磁珠富集mRNA用于分析转录组信息，利用金属离子将富集mRNA随机断裂成200 bp左右的小片段。在逆转录酶的作用下，以mRNA 为模板，利用随机六碱基引物，反转录合成cDNA。加入End Repair Mix对双链cDNA粘性末端进行修复，PCR扩增15个循环后，进行目的条带回收，TBS380（Picogreen）定量，在cBot上进行桥式PCR扩增，生成Clusters后，利用Hiseq2500

1.3 数据分析

1.3.1 原始数据质控及转录组拼接组装 通过Illumina 测序平台对卵泡组织进行测序后，因得到的数据量十分庞大，为宏观上直接反映样本的测序质量和文库的构建质量，对所有的测序reads 进行碱基分布和质量波动统计，主要包括碱基含量分布统计、碱基质量分布统计和碱基错误率分布统计。同时由于鹅属于无参考基因组动物，针对无参考基因组的转录组研究，将测序数据进行质量剪切，将测序数据中测序接头序列末端质量低于20 的碱基、含N比率超过10 %和质量修剪后长度小于20 bp的序列去除，获得RNA-seq高质量测序数据。得到高质量的测序数据后将所有测序读段通过Trinity软件进行从头组装，并对组装后得到的数据进行统计分析[8]。

1.3.2 序列分析 序列分析可对组装后得到的转录本进行功能注释，来探究这些转录本的生物功能，在注释前需利用Trinity软件提供的ORF预测流程对组装得到所有转录本序列进行核苷酸序列和蛋白质序列预测，根据真核生物默认序列含有非编码区原理，利用Markov模型预测该序列中的最佳ORF区域，并使用Pfam数据库对预测结果进行校正，保留比对到数据可预测出ORF的核苷酸序列和蛋白序列，之后将拼接所得到的所有核苷酸序列别于NR、GO、COG、KEGG数据库进行比对，获得相应的注释信息对其进行功能注释[9-10]。

1.3.3 差异因的筛选 为实现所得到的基因表达量可视化，通过RSEM对不同样品基因的表达量进行计算，并以FPKM值作为衡量基因表达水平的标准，使用EdgeR软件进行差异基因计算[11]，本次实验中，差异表达基因的筛选标准为FDR＜0.05，︱logFC︱≥1，同时将所得到的差异基因GO数据库和KEGG数据库富集进行注释与统计分析。

2 结果与分析

2.1 原始数据分析 实验采用Illumina Hiseq 测序平台对鹅等级前卵泡组织完成转录组测序，并构建Illumina PE文库进行2×101 bp测序，为了保证后续试验质量，对获得的测序数据进行RNA质量验证并对原始测序数据质控。

2.1.1 RNA质量验证 将从鹅等级前卵泡组织中提取的总RNA进行质量验证，其中由SY到PE各等级RNA 量为17.7、23.3、19.6、17.6、8.2 μg，且RNA浓度大于200 ng/μL，符合单次建库RNA质量要求。经琼脂糖电泳检验，5S、18S、28S条带清晰（图 1）。验证结果表明RNA质量良好，未发生降解，符合下一步建库要求。

图1 鹅5个等级前卵泡电泳结果图

2.1.2 原始测序数据质控 经测序共得到2亿多条序列，为宏观反映出每条序列的质量，通过统计学的方法对A/T/G/C碱基含量、碱基质量、碱基错误率分布统计。经碱基含量分布统计结果表明，在整个测序过程中，AT、GC含量相同，且除序列起始和测序接头位置由于6 bp随机引物引起A、T、C、G 出现正常波动外，其余呈现稳定状态，且模糊碱基N比率在合理范围内。碱基质量检测是对所有测序循环下的测序读段的质量波动进行统计分析，对于RNA-seq技术，根据下述公式进行计算质量值Q。通过统计分析发现随测序进行，Q值随着测序试剂不断消耗而下降，这属正常现象。且整体测序错误率较低，读段碱基质量较高，所获得测序数据达到后续分析要求。

通过碱基错误率分布分析，发现碱基错误率会随着测序序列长度增加和测序化学试剂消耗而增加，且建库过程中反转录所使用随机6碱基引物与RNA模板结合不完全导致测序前6个碱基位置错误率较高。但所有测序单个碱基位置平均错误率低于0.1%，符合实验要求，可用于后续实验。为保证后续的信息分析的准确性，对进行测序的数据进行质量剪切，原始测序数据过滤，经质量剪修后统计共得到224 929 214条序列。具体结果见表 1。

表1 不同胎次繁殖母猪预定的淘汰率和产活仔数

2.2 测序数据分析 鹅属于无参考基因组动物，通过Trinity软件将所有测序读段从头组装成重叠峰（Conting）和单一序列（Singleton），得到 277 077 条Unigenes序列和314 320条Transcripts序列，并对其长度分布进行统计，其中长度在1～600的Unigenes序列和Transcripts序列所占的总体比例较高，具体见表2。

表2 拼装结果统计

2.3 序列分析

2.3.1 NR库注释结果 为对所组装的转录本进行功能注释，注释前需先利用Trinity软件对组装得到的所有转录本进行基因预测，经预测共得到74 405条核苷酸序列。将预测所得的核苷酸序列与NR数据库进行比对，可得到本物种转录组序列与相近物种的相似情况及同源序列的功能信息。NR数据库为NCBI非冗余蛋白库，是一包含SwissProt、PIR、PRF、PDB蛋白质数据库中非冗余的数据以及从GenBank和RefSeq的CDS数据库中翻译蛋白质数据的综合数据库。通过比对后统计共有42 453条序列注释到NR数据库，且进行物种同源性比对，具体结果见图 2。图2反映了测序样本比对近源物种的情况，图中每一扇形表示一个物种，扇形面积越大表示比对上序列越多，同源性越高。其中绿头鸭比对序列的数目最多，因此可认定绿头鸭与鹅的同源性最高。

图2 对比物种分布

图3 GO二级分类统计图

2.3.2 GO 注释结果 GO数据库是基因本体论联合会对不同数据库中关于基因和基因产物的生物学术语进行标准化及基因和蛋白功能进行统一限定和描述的综合数据库。将所得Unigenes通过GO数据库按照其参与的生物过程（Biological Process，BP），分子功能（Molecular Function，MF）和细胞组分（Cellular Component，CC）进行分类注释，其三大分支下又有很多小的分级，级别数越多，功能注释越详细。经比对共有1 045条Unigenes对比到GO数据上，其中在BP分类中细胞过程、代谢过程和单有机体过程类别所占比例最多；在CC分类中细胞、细胞部分分类和细胞器所占比例最多；在MF分类中，分子绑定功能和催化活性功能类别所占比例最多，具体结果见图3。

2.3.3 GOG注释结果 GOG数据库是蛋白直系同源簇数据库，通过比对可以进行功能注释，归类以及蛋白进化分析，从而可对转录本进行功能分类，其中功能预测所占的转录本数目最多，其次为信号转导机制和翻译后修饰。具体见表3。

2.3.4 KEGG通路注释结果 与基因组破译公共数据库（KEGG）进行比对分析，可通过KEGG数据库中丰富的通路信息系统的了解与分析基因的生物学功能。通过分析结果发现，共有11 949 条Unigenes 比对到KEGG数据库，共涉及到341条通路，其中富集最多的通路为癌症通路，其次为PI3K-AKT信号通路和MAPK信号通路，说明这3条通路在鹅等级卵泡组织中最为活跃，具体见图4。与KEGG数据对比可获得转录本对应的KO编号，根据KO编号可知其可能参与的具体生物学通路，可根据其KO注释对其所参与的KEGG通路分类为细胞过程、环境信息处理、遗传信息处理、代谢和有机系统5个分支，具体见图 5。

2.4 基因的表达 通过计算可得到鹅等级前卵泡组织的FPKM值，其中各等级前卵泡中FPKM值大于1的基因占20 % 以上，大于10的基因占3 % 以上（表4）。同时通过FPKM值对所有基因的表达量概率密度做图，其中横坐标为log10FPKM，此数值越高说明基因的表达量越高，纵坐标对应基因的密度，密度曲线的峰值为整个样本中基因表达量最集中的位置，具体结果见图 6。

2.5 差异表达基因的分析

2.5.1 差异基因的筛选与注释分析 通过EdgeR软件对鹅等级前卵泡中的差异基因计算，筛选标准为FDR＜0.05，︱logFC︱≥1，即基因在2个样本间的差异显著性结果校验结果小与0.05，且样本间差异倍数以2为底的对数绝对数值大于1可认为此基因在这两样本间差异表达。SY与PE等级卵泡均为鹅等级卵泡，SY等级卵泡为等级前卵泡发育的等级卵泡中最后1个等级，在众多的SY等级卵泡中只有1个可以发育为等级前卵泡，而PE等级卵泡是等级卵泡中的初始等级卵泡，通过研究分析发现这两等级卵泡中所存在的差异基因数量最多，其中以SY等级卵泡为参照，在PE 等级卵泡组织中共获得15 516条差异基因，其中上调基因8 293条，下调基因7 223条，具体见图 7 。

表3 GOG数据库注释结果

图4 KEGG pathway中前20个信号通路

图5 KEGG注释统计图

表4 鹅等级前卵泡组织FPKM值

图6 表达量分布

图7 PE-SY 卵泡组织存在的差异表达基因可视化图

图8 PE-SY 卵泡组织、上下调基因注释柱形图

通过与GO数据库比对筛选得到的15 516条差异基因进行分子功能分类，共比对注释到3 626条基因，其中上调基因1 450条，下调节基因2 176条，共注释到参与的生物过程，细胞组分和实现分子功能3大类的22个小类上，其中生物过程24类，细胞组分17类，分子功能14类，具体结果见图 8。

通过KEGG数据库比对分析，共筛选出4 950条差异基因注释到299条KEGG通路上。为了能够识别与生物现象最相关的生物过程，提高研究的可靠性，使用KOBAS（http://kobas.cbi.puk.edu.cn/ home.do）对所得到的KEGG PATHWAY富集分析，经计算校正后P≥0.05的KEGG通路定义为在差异表达基因中显著富集的KEGG通路。其中显著富集通路有174条。其中显著富集的通路有PI3K-AKT通路，癌症通路和黏着斑通路，具体结果见图 9。

图9 PE-SY 卵泡组织KEGG通路富集图

3 讨 论

鹅卵泡发育过程是较为复杂的生物过程[12]。为探究这一过程，通过RNA-Seq技术对产蛋高峰期鹅各等级前卵泡组织进行转录组测序分析，RNASeq技术作为第2代高通量测序技术，对揭示研究动植物差异基因表达及调控发挥了重要作用[13]。本实验通过Illumina HiSeq 2500测序平台对鹅等级前卵泡组织进行转录组文库的构建，并获得大量转录组信息，同时通过NR、GO和KEGG数据库进行比对，进行功能注释及分类。共得到224 929 214条序列，Q20比例为95.5%以上，GC含量为47%以上，说明文库构建成功且质量较好。经dnovo组装拼接后共得到114 486 条Unigenes序列，并预测出74 405条核苷酸序列，说明测得到基因74 405个，经NR数据库对比发现，绿头鸭、原鸡和原鸽与其同源性较高。通过KEGG通路比对分析，其中最为活跃的通路为癌症通路，PI3K-AKT通路和MAPK通路，而PI3K-AKT通路和MAPK通路在卵泡发育过程中的作用及功能以广泛的被验证及证实。PI3K-AKT通路为动物体内一基本信号通路，可灵活的参与体内卵母细胞增长和原始卵泡发育过程，同时还涉及颗粒细胞增殖、分化和凋亡，继而调节卵泡发育。MAPK通路在动物的胚胎发育和细胞增殖、分化、凋亡中同样也具有重要的调节作用[14-15]。由此可以说明，卵泡的发育是一个动态的空间调节过程，并不受单一通路影响。同时对SY与PE等级卵泡组织间对比共获得15 516条差异基因，其中上调基因8 293条，下调基因7 223条。通过GO数据库比对发现差异基因主要集中在生物过程类中。通过KEGG富集分析，富集最多的通路为PI3K-AKT通路、癌症通路和黏着斑通路，这与兰道亮等[16]在对牦牛卵巢组织做转录组测序数据分析时所得结果相同，可表明这3条通路在卵泡发育过程中的活跃状态和对卵泡发育中所起到的作用。目前研究还尚未发现有对鹅等级前卵泡进行等级分类测序研究实例，本实验极大丰富了鹅等级前卵泡发育时期组织的基因数据，同时通过各等级间卵泡组织所得差异表达基因数据分析，可进一步对鹅等级前卵泡发育关键基因进行发掘。

4 结 论

本实验通过RNA-Seq技术首次对鹅等级前卵泡进行高通量转录组测序和分析，研究了鹅不同等级下的卵泡组织的差异表达基因，同时对所获得差异表达基因进行功能分类和代谢通路划分，不仅丰富了鹅卵泡组织的转录信息，还为以后开展鹅遗传育种研究及分子调控研究提供理论基础。

[1] 吴阳升, 林嘉鹏, 汪立芹, 等. 绵羊小卵泡与中卵泡转录组差异特征分析[J]. 江苏农业学报, 2016, (4):832-842.

[2] 朱志明, 陈红萍, 林如龙, 等. 山麻鸭开产期和产蛋高峰期卵巢组织转录组分析[J]. 中国农业科学, 2016, (5):998-1007.

[3] 韩昆鹏, 段炼, 李婷婷, 等. 京海黄鸡卵巢转录组研究:基因结构分析与新基因发掘注释[J]. 中国畜牧兽医, 2016, (4):854-861.

[4] Tariq M, Chen R, Yuan H Y, et al. De novo transcriptomic analysis of peripheral blood lymphocytes from the Chinese goose:gene discovery and immune system pathway description[J]. PLoS One, 2015, 10(3): e0121015.

[5] 王丹. 维生素A对朗德鹅肥肝抗氧化性能和血液生化指标影响的研究[D]. 石河子:石河子大学, 2011.

[6] 吴慧英, 贾汝敏, 刘铀, 等. 鹅繁殖就巢机理研究进展[J].中国草食动物, 2009, (5):56-58.

[7] Johnson A. Ovarian follicle selection and granulosa cell differentiation [J]. Poult Sci, 2015, 94(4): 781-785.

[8] Grabherr M G, Haas B J, Yassour M, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome[J]. Nat Biotechnol, 2011, 29(7):644-652.

[9] Camacho C, Coulouris G, Avagyan V, et al. Blast+: architecture and applications[J]. Bmc Bioinformatics, 2009, 10(1):421.

[10] Conesa A, Cotz S, Terol J, et al. A universal tool for annotation, visualization and analysis in functional genomicsresearch.Bioinformatics[J]. 2005, 21:3674-3676.

[11] Robinson M D, Mccarthy D J, Smyth G K. A Bioconductor package for differential expression analysis of digital gene expression data[J]. Bioinformatics, 2010, 26(1):139-140.

[12] 张林媛. 鸿雁卵泡的组织学研究[D].哈尔滨:东北林业学, 2007.

[13] Zeng T, Zhang L P, Li J J. De novo assembly and characterization of muscovy duck liver transcriptome and analysis of differentilly regulated genes in response to beat stress[J].Cell Stress Chaperon, 2015,( 2):483-493.

[14] 甄艳红. MAPK-CEBPβ对猪颗粒细胞凋亡、周期及类固醇激素分泌的调控作用和机制的研究[D]. 武汉:华中农业大学, 2014.

[15] 吴艳青, 陈丽云, 张正红, 等. 磷脂酰肌醇-3激酶/蛋白激酶B/骨形态发生蛋白-15通路与哺乳动物卵巢卵泡发育[J]. 中国医学科学院学报, 2013, (2):224-228.

[16] 兰道亮, 熊显荣, 位艳丽, 等. 基于RNA-Seq高通量测序技术的牦牛卵巢转录组研究:进一步完善牦牛基因结构及挖掘与繁殖相关新基因[J]. 中国科学:生命科学, 2014, (3):307-317.

S835.2

：A

：10.19556/j.0258-7033.2017-06-063

2016-11-29；

2016-12-18

吉林农业大学青年拔尖人才支持计划项目（2016）；吉林省留学回国人员择优资助项目（2015）；吉林省科技厅重点科技攻关项目（20150204025NY）；国家科技部星火计划项目（2015GA660001）；吉林省现代家禽产业技术体系项目（2016）作者简介：孙永峰（1977-），男，吉林长春市人，博士，副教授，主要从事家禽遗传育种工作研究，E-mail：sunyongfeng1977@

126.com

* 通讯作者