分享至
花生转录组denovo测序揭示果针发育相关的功能基因
近日，由广东省农业科学院作物研究所梁炫强研究员、陈小平博士，山东省农业科学院花生研究所所长禹山林研究员及国际半干旱热带作物研究所(ICRISAT)基因组中心负责人Rajeev Kumar Varshney等完成的花生转录组de novo测序的结果发表于Plant Biotechnology Journal。该杂志创刊于2003年，影响因子5.44，是专注于应用植物科学的期刊。本研究是中国花生基因组计划取得的又一阶段性成果，测序和部分信息分析工作均在Macrogen千年基因完成。
花生(Arachis hypogaea L.)是我国重要的油料和经济作物，因其“地上开花，地下结荚”的特性被称之为“落花生”。花生受精后子房柄(又称果针)开始伸长并启动向地性生长，从而将子房带入地下，到达一定深度后子房开始膨大并发育成荚果。研究表明，黑暗和机械刺激是果针发育的必要条件，而温度、水分、空气、营养等外界因素也影响果针发育的过程。因此，果针能否入土关系到花生能否正常结荚。但由于目前尚未完成花生全基因组图谱的绘制，关于果针发育的分子机制目前尚不明确，本研究将从转录组水平阐明影响果针发育的关键基因。
本研究应用Roche 454(GS FLX Titanium)测序平台，分别对花生的地上果针(AP)和2个不同发育时期的地下果针(SP1和SP2)进行1/2 run转录组de novo测序，得到274M(704,738条reads)、290M(711,496条reads)和238M(609,841条reads)原始数据(raw data)，测序reads的平均读长为396bp，去除接头等序列和低质量值reads后获得83%的有效数据(clean data)。相关测序结果已提交至NCBI的SRA数据库，登录号为SRA053198。
AP、SP1和SP2的测序数据合并后以 Newbler(v2.6)和TGICL2.0分别进行组装，并比较两个软件的组装性能。组装结果与现有的花生转录组数据合并去冗余后构建花生转录组参考序列，序列大小114.9M，约是花生基因组大小的4.1%。组装结果分别进行功能基因预测、GO(gene ontology)分类、KEGG(Kyoto Encyclopaedia of Genes and Genomes)代谢通路注释和基因差异表达分析，结果表明在2194个差异表达基因中，地上果针和地下果针分别有859和1068个转录本的表达显著性上调。在地上果针中，与光合作用和衰老相关的基因表达显著上调，这2个基因很可能阻止了果针的膨大从而导致败育。本研究的结果为探索果针发育的分子机制提供了充分的理论依据，并为花生转录组和功能基因组的研究提供了更多的数据资源。
文献解析：
1、 实验材料
本研究选取的花生栽培种为“粤油 7 号”，该品种是由广东省农业科学院作物研究所选育的高产花生品种，目前已在南方大面积推广种植。在花生受精完成后，通过人工干预的方法使部分果针未能伸入土壤中，并分别取开花后第8、10、12、16、20天的地上果针样品混合后成为一个AP样本。另外，对于正常发育的果针，分别取果针伸入土壤后第2、 4、6、8、10、12天的地下果针样品，并将第2、4、6天的样本混合成为SP1，第8、10、12天的样本混合成为SP2。所有样本均以改良的 CTAB法抽提总RNA，并以安捷伦2100检测RNA的质量。
2、 454测序
RNA样品的文库构建和测序均在Macrogen总部最高质量标准的基因组学实验室完成，实验试剂均采用进口原厂试剂，实验流程和质量控制均严格受Macrogen董事会监管。AP、SP1和SP2样本均进行1/2run的454测序，分别得到 274M(704,738条reads)、290M(711,496条reads)和238M(609,841条reads)原始数据，测序reads的平均读长为396bp，且大部分reads分布于500bp左右。将测序得到的原始数据首先通过SeqClean软件与UniVec数据库进行比对后，去除接头等序列，然后设置严格的过滤标准去除低质量reads。如果一条read中30%的碱基质量值小于Q20或者该read长度小于50bp，则会被过滤掉。经过以上数据处理，共得到1683,455条有效reads，有效数据(clean data)占原始数据(raw data)的比例为83%。相关测序结果已提交至NCBI的SRA数据库，登录号为SRA053198。
3、 转录组de novo组装
目前用于组装的软件种类较多，研究表明对于454测序数据的组装，Newbler和TGICL的组装性能高于其他软件。为了获得最佳的组装效果，本研究先应用RepeatMasker软件屏蔽简单重复等低复杂性序列，然后运用Newbler(v2.6)和TGICL2.0分别进行组装，并比较这两个软件对于花生转录组de novo组装的性能。其中，Newbler是Roche 454平台自带的组装软件，本研究开展时已更新至v2.6版本，现已更新至v2.8版本，TGICL是专门用于转录组组装的软件。Newbler组装得到近3万个contigs，其中63%的contigs大于1kb，contig N50为1730bp，平均contig长度为1447bp。TGICL组装得到近5万个contigs，其中28%的contigs大于 1kb，contig N50为974bp，平均contig长度为859bp。另外，为了进一步验证组装效果，本研究以大豆基因组序列预测得到的蛋白质序列作为参考序列进行同源性比对，结果表明TGICL的组装结果与大豆蛋白质序列的同源性更高。当单个contig与大豆蛋白序列比对的同源性设置为R50%时，Newbler 的组装结果可比对到13,000个大豆蛋白，而TGICL可比对到15,000个大豆蛋白。因此，Newbler和TGICL在花生转录组数据的组装过程中各有优势，TGICL能够得到更多的contigs数量，因而能够反应转录本的多态性，而Newbler能够得到更大的contigs长度。
4、 转录组参考序列的构建
将本研究的组装结果与现有的花生转录组数据合并去冗余后，得到151,533个转录本，共计114.9M，约是花生基因组大小的4.1%。其平均contig长度为758bp，contig N50为866bp，平均GC含量为40.5%。现有版本的花生转录组组装结果NCBI Peanut Unigene Build #3和PeanutDB分别是34.41和33.97M，将本研究的转录本分别与这两个版本比对得知，共有5,286个转录本无法比对至参考序列，这些很可能是本研究发现的新转录本。
| Transcriptome assembly | NCBI Peanut Unigene Build #3 | PeanutDB | |
| Contigs | 151,533 | 52,358 | 32,619 |
| Total Size (bp) | 114,867,135 | 34,412,752 | 33,971,914 |
| Contigs (≥100 bp) | 151,533 | 52,385 | 32,617 |
| Contigs (≥1 kb) | 30,044 | 9,960 | 12,225 |
| Maximum length (bp) | 8,119 | 7,622 | 29,281 |
| N50 (bp) | 866 | 922 | 1,167 |
| Average length (bp) | 758 | 656 | 1,041 |
| Contig with significant hit (%) | 98,731 | 36,627 | 26,416 |
| Contig with 80% or greater coverage (%) | 26,588 | 8,213 | 9,449 |
| Soybean protein hits (%) | 29,380 | 21,373 | 19,020 |
| Soybean proteins with 80% or greater coverage (%) | 15,330 | 8,713 | 9,990 |
表3. 现有版本的花生转录组组装结果比较
应用SSAHA2 v2.5.3将AP、SP1和SP2样本测序得到的reads分别与构建的转录组参考序列比对，约99%的reads能够匹配至参考序列，这表明该参考序列能够较全面的体现三个样本的转录情况。其中，三个样本中平均比对至每个转录本的reads数为10.5、10.9和10.1，分别有73.79%、 76.30%和75.17%的reads与参考序列有唯一匹配，约23%-26%的reads有多个匹配，这部分reads在后续的分析中均被过滤掉了。根据比对结果得知，共鉴定出74,974个转录本，其中AP、SP1和SP2样本中分别鉴定出49,632、49,952和50,494个转录本。虽然三个测序样本分别取自不同的发育时期，但27,732个转录本是三个样本共有的，且AP中约84%的转录本能够在SP1或者SP2中检测到。
通过RPKM值衡量三个样本的基因表达水平，结果表明AP、SP1和SP2的平均表达量分别为RPKM15.2、14.3和15.6，大部分基因是低表达和中等表达水平(RPKM值介于3和50之间)，只有1.2%-1.4%的基因是高表达水平的(RPKM值大于100)。
5、 基因功能的注释
为了鉴定相关的功能基因，将三个样本中共74,974个转录本利用BLASTX与UniProtKB进行序列同源性比较，E值设定为1E-5。其中，59,342个转录本表现出较高的同源性，共匹配了28,633个功能蛋白和大量未知功能的蛋白。另外，约20%的转录本未能与数据库中的任何蛋白匹配，很可能由于这部分转录本是花生特有的转录本。
对于已通过UniProtKB注释的转录本，利用GO对已注释基因的功能进行分类，共有42,995个转录本可比对至GO数据库，且38,280个转录本与分子功能相关，31,152个转录本与生物过程相关，17,881个转录本与细胞组成相关。另外，分子功能分类中与代谢活性和结合功能相关的转录本所占比例最大，约20%以上。在生物过程分类中，细胞过程和代谢过程分别占13.39%和6.33%。在细胞组成分类中，11.42%和4.53%的定位在细胞内和细胞膜，只有很少一部分位于细胞外基质、细胞连接区。
为了进一步了解转录本的生物学功能，进行KEGG代谢通路注释。利用KOBAS v2.0将蛋白数据库注释的转录本比对到KO数据库，结果表明共有2,968个转录本比对到174个KEGG代谢通路中，其中168个KEGG通路是地上果针(AP)和地下果针(SP)共有的，这表明两种样本表达的大部分基因是相同的。
6、 差异表达分析
利用DEGseq和 RPKM分别对差异表达基因进行分析，在地上果针和地下果针样本中分别鉴定出859和1068个表达上调的基因。其中，地上果针中与光合作用和衰老相关的基因表达显著上调，这2个基因很可能阻止了果针的膨大从而导致败育。另外，差异表达分析的结果也得到了RT-PCR的验证。
