化橘红为芸香科植物化州柚Citrus grandis ‘Tomentosa’或柚Citrus grandis (L.) Osbeck的未成熟或近成熟的干燥外层果皮，前者习称“毛橘红”，后者习称“光七爪”“光五爪”或“光橘红”，具有燥湿化痰、理气、消食的功效^[1].有文献报道将目标起始密码子多态性(SCoT)分子标记^[2]、简单重复序列(SSR)分子标记技术^[3]、DNA条形码^[4]等基因分析技术用于不同品种化橘红的亲缘关系分析.目前，已知柚皮苷、野漆树苷、柚皮素、新橙皮苷、芹菜素等黄酮类成分是化橘红的主要活性成分^[5-9]，但毛橘红与光橘红的黄酮类成分存在很大差异^[10-11]，由于其黄酮类成分的生物合成相关的基因序列还不清楚^[12]，所以不同品种化橘红的黄酮类成分差异形成的分子基础有待深入研究.

随着高通量测序技术的不断发展，测序成本不断降低，转录组研究已经成为挖掘药用植物功能基因的重要手段之一.近年来，国内外已经开展了三七^[13]、茴香^[14]、决明子^[15]、铁皮石斛^[16]、薯蓣^[17]、马蓝^[18]、当归^[19]、泽泻^[20]多种药用植物转录组的研究; 在黄酮类成分的生物合成方面，有文献报道采用二代测序技术研究了紫丁香^[21]、三叶青^[22]、甘草^[23]、杭菊^[24]、金银花^[25]、红花^[26]、雪莲^[27]等药用植物体内的黄酮类生物合成相关基因及其表达差异.目前，柚及其同属植物的转录组也有少量报道，例如：利用转录组研究低温对柑橘长时间贮藏后果实品质的影响^[28]、柚子分子标记的开发^[29]、柚子果皮颜色突变相关差异表达基因^[30]等.本文采用二代高通量测序技术，对化州柚和柚进行转录组测序、拼接组装，建立化州柚和柚叶的转录组，对黄酮类成分生物合成相关的差异表达基因进行了功能注释与比较分析，为阐明化州柚和柚的黄酮类有效成分含量差异的分子基础研究奠定基础.

1.   研究方法

1.1   材料与试剂

化州柚和柚的嫩叶样本采自广东食品药品职业学院药圃. Omega Plant RNA kit试剂盒(美国Omega)，NEB#7530试剂盒(新英格兰Biolabs)，DNA 1000 assay Kit (美国安捷伦科技有限公司).

1.2   仪器

PCR仪(东胜兴业科学仪器有限公司)，高速离心机(德国艾本德股份公司)，Nanodrop 2000 (美国赛默飞世尔科技公司)，Illumina HiSeq 2500 (美国Illumina公司).

1.3   样本采集、总RNA提取及转录组测定

采集化州柚和柚的嫩叶适量，在液氮中速冻15 min后，保存于-85 ℃超低温冰箱备用.按照Qiagen RNeasy Mini Kit提取试剂盒说明书的方法提取各样本总RNA，用Nanodrop 2000与Agilent 2100 Bioanalyzer检测所提取总RNA的浓度与纯度，A₂₆₀/A₂₈₀比值在1.8~2.0，RNA分子完整数(RIN，RNA Integrity Number)大于8，则符合转录组测序样本要求，转录组测定由广州基迪奥生物科技有限公司完成.

1.4   转录组数据组装

测序产生的原始数据(Raw Reads)，经过初步过滤后得到待分析数据(Clean Reads)，去除有带接头的Reads、N的比例大于10%的Reads以及质量很低的Reads后得到高质量的Clean Reads.使用短序列组装软件Trinity将高质量待分析数据从头组装得到单基因(Unigene)，用N₅₀数值来评估组装结果质量. N₅₀越长，数量越少，说明组装质量越好. Unigene表达量的计算使用RPKM法(Reads Per kb Per Million Reads)，其计算公式为：RPKM=(1 000 000×C)/(N×L/1 000)，假设RPKM为某Unigene的表达量，则C为比对到某Unigene的Reads数，N为比对到所有Unigene的总Reads数，L为某Unigene的碱基数.

1.5   Unigene基本功能注释

首先，通过比对软件Blastx将Unigene序列比对到美国国家生物信息中心(NCBI)非冗余蛋白序列数据库(Nr)、由欧洲生物信息学研究所(EBI)维护的经过注释的蛋白质序列数据库(Swiss-Prot)、京都基因与基因组百科全书数据库(KEGG)和蛋白相邻类的聚簇数据库(COG/KOG)，得到与给定Unigene具有最高序列相似性的蛋白，从而得到该Unigene的蛋白功能注释信息.

1.6   差异表达基因分析

先用edgeR软件包进行化州柚与柚的转录组差异统计分析，利用错误发现率(FDR)与log 2FC(FC为两样品间表达量的比值)来筛选差异基因，筛选条件为FDR < 0.05且|log 2FC|>1，得到差异表达基因之后，对差异表达基因再进一步做基因本体(GO)功能分析和KEGG Pathway分析. GO功能分析给出差异表达基因的GO功能分类注释和GO功能显著性富集分析.通过Pathway显著性富集可以确定化州柚和柚转录组的差异表达基因参与的最主要生化代谢途径和信号转导途径.

2.   结果

2.1   测序与组装结果

本次测定共获得86 298 757 bp碱基，组装得到116 202个单基因(Unigene)，N₅₀为1 543，GC%为42.077 5，最长基因长度20 219 nt，最短基因长度201 nt，基因平均长度742 nt. Unigene长度分布见图 1.

图  1  化州柚和柚转录组Unigene基因的长度分布图

Figure  1.  The Unigene length distribution in the transcriptome of CGT and CGO

下载: 全尺寸图片 幻灯片

2.2   Unigene基本注释

运用BLAST比对软件将组装得到的116 202个Unigene与已知的4个公共数据库Nr、Swiss-Prot、KOG、KEGG进行比对，其中Nr注释了53 738个基因，Swiss-Prot注释了52 965个基因，KOG注释了40 293个基因，KEGG注释了24 891个基因，总共获得同源比对信息的Unigene有68 923条，占59.31%，没有获得同源比对信息的基因有47 279个，占40.69%(图 2).

图  2  化州柚与柚转录组在4个数据库基因注释的维恩图

Figure  2.  The Venn diagram of the transcriptome annotated in the four databases of CGT and CGO

下载: 全尺寸图片 幻灯片

KOG注释结果显示：化州柚和柚的40 293个Unigene分布于25个KOG分类中，其中聚类到一般性功能预测、信号转导机制、翻译后的修饰、蛋白转换、分子伴侣等功能基因分布较多(图 3).

图  3  化州柚与柚转录组Unigene基因的KOG功能分类图

Figure  3.  The functional classification of KOG in the transcriptome of CGT and CGO

下载: 全尺寸图片 幻灯片

化州柚与柚转录组GO功能注释的具体亚类与基因数量见图 4. 图 4显示：分子功能(Molecular Function)、细胞组分(Cellular Component)、生物过程(Biological Process)又可细分为50个功能亚类，其中生物学过程分为23个亚类，细胞组分分为16个亚类，分子功能分为11个亚类.在生物学过程中，代谢过程(Metabolic Process)、细胞过程(Cellular Process)和单有机体过程(Single-Organism Process)所占基因数较多; 在细胞组分中，涉及细胞(Cell)、细胞组分(Cell Part)和细胞器(Organelle)的基因较多; 在分子功能中，涉及催化活性(Catalytic Activity)和结合(Binding)基因较多.

图  4  化州柚与柚转录组Unigene基因的GO功能分类图

Figure  4.  The functional classification of GO in the transcriptome of CGT and CGO

下载: 全尺寸图片 幻灯片

KEGG数据库注释结果(图 5)显示：具有通路注释的基因共14 833个，涉及136个代谢通路.其中代谢途径相关基因6 430个(43.35%); 次级代谢产物的生物合成相关基因3 529个(23.79%); 核糖体相关基因1 665个(11.22%)，碳代谢相关基因1 423个(9.59%)，氨基酸的生物合成相关基因1 032个(6.96%)，植物-病原体相互作用相关基因743个(5.01%)，内质网中的蛋白质加工相关基因684个(4.61%)，氧化磷酸化相关基因671个(4.52%)，糖酵解/糖异生相关基因661个(4.46%)，嘌呤代谢相关基因632个(4.26%).

图  5  化州柚与柚转录组Unigene基因的KEGG通路注释

Figure  5.  The KEGG pathway annotation in the transcriptome of CGT and CGO

下载: 全尺寸图片 幻灯片

2.3   基因差异性表达分析

本研究共筛选得到化州柚与柚叶转录组的差异表达基因共6 419个，其中化州柚相对于柚上调基因有3 799个，占差异表达基因的59.18%，下调基因2 620个，占40.82%.获得GO注释的差异表达基因有4 826个，1 346条GO注释中，937个注释为生物学过程相关，111个注释为细胞组分相关，298个注释为分子功能相关.

与KEGG数据库进行比对，共有771个基因获得注释，涉及125条代谢通路，其中显著差异表达富集通路的基因有：植物-病原体相互作用相关基因134个(17.38%); MAPK信号传导途径-植物相关基因59个(7.65%); 植物激素信号转导40个(5.19%); 苯丙烷生物合成相关基因37个(4.80%)、倍半萜类化合物和三萜类化合物生物合成相关基因10个(1.30%); 胡萝卜素生物合成相关基因17个(2.20%); 油菜素类固醇生物合成相关基因7个(0.91%); 泛醌和其他萜类化合物-醌生物合成相关基因相关基因14个(1.82%); 类黄酮生物合成相关基因11个(1.43%); 氨基糖和核苷酸糖代谢相关基因36(4.67%); 次级代谢产物的生物合成相关基因223个(28.92%); 芪类化合物、二苯基庚烷类和姜辣素生物合成相关基因8个(1.04%); 黄酮和黄酮醇生物合成相关基因3个(0.39%); 谷胱甘肽代谢相关基因29个(3.76%)等.

3个与黄酮和黄酮醇生物合成相关基因注释为黄烷酮7-O-葡萄糖苷2″-O-β-L-鼠李糖基转移酶(C12RT1)、黄酮醇-3-O-葡萄糖苷L-鼠李糖基转移酶(RT)、黄酮醇-3-O-葡糖苷/半乳糖苷葡糖基转移酶(3GGT); 11个与类黄酮生物合成相关的差异表达基因注释为反式肉桂酸4-单加氧酶(CYP73A12)、咖啡酰-CoAO-甲基转移酶(TSM1)、查尔酮合成酶(CHS1)、黄酮醇合成酶(FLS)、莽草酸O-羟基肉桂酰转移酶(HCT)(5个)、黄烷酮7-O-葡萄糖苷2″-O-β-L-鼠李糖基转移酶(C12RT1)、白细胞花青素还原酶(LAR)，这14个黄酮类生物合成相关差异表达基因表达量rpkm值见表 1.

表  1  化州柚与柚转录组中与黄酮类成分生物合成相关的差异表达基因

Table  1.  The differentially expressed genes related to flavonoids biosynthesis in the transcriptome of CGT and CGO

基因编号	基因	CGO_rpkm	CG_rpkm	P	错误发现率(FDR)
Unigene0018220	C12RT1	27.22±3.76	9.51±2.18	1.81E-08	9.79E-07
Unigene0001546	RT	9.36±2.44	23.75±3.62	3.20E-06	1.21E-04
Unigene0062842	3GGT	18.53±5.34	0.05±0.072	7.37E-51	5.51E-48
Unigene0102621	CYP73A12	0.054±0.004 0	0.40±0.34	1.72E-03	3.11E-02
Unigene0032665	TSM1	0.078 5±0.07	2.47±0.47	2.16E-22	3.98E-20
Unigene0000397	CHS1	0.46±0.28	1.29±0.07	8.21E-04	1.68E-02
Unigene0024628	FLS	15.06±4.48	4.97±0.48	1.50E-07	7.08E-06
Unigene0011342	HCT	8.50±4.28	0.06±0.04	2.30E-34	8.25E-32
Unigene0029908	HCT	0.32±0.13	2.65±0.53	1.32E-14	1.39E-12
Unigene0037769	HCT	53.54±1.68	25.88±1.19	2.66E-07	1.21E-05
Unigene0093083	HCT	0.77±0.24	0.09±0.06	8.08E-08	3.98E-06
Unigene0115527	HCT	0.005 6±0.009 7	2.01±0.74	1.35E-27	3.41E-25
Unigene0018220	C12RT1	27.22±3.76	9.51±2.18	1.81E-08	9.79E-07
Unigene0015054	LAR	0.90±0.38	0.08±0.05	3.19E-08	1.67E-06

下载: 导出CSV 

| 显示表格

3.   讨论

本次研究发现化州柚和柚的转录组有6 419个差异性表达基因，Pathway富集分析结果显示：与次级代谢产物的生物合成相关基因最多(占28.92%)，这为化州柚和柚的化学成分及含量差异分子基础研究奠定了基础; 与植物-病原体相互作用相关基因也较多(占17.38%)，为阐明二者植物生态适应性差异奠定了基础.另外，本研究结果显示：与黄酮和黄酮醇生物合成相关差异表达基因仅3个，与类黄酮生物合成相关的差异表达基因有11个，说明化州柚和柚在黄酮和黄酮醇生物合成方面很类似，但在类黄酮生物合成方面差异较大.由于本次研究采集的是化州柚和柚的嫩叶样本，而化橘红药用部位是未成熟或近成熟的干燥外层果皮，所以本次研究为不同品种来源化橘红黄酮类成分差异奠定了基因数据基础，其生物合成的相关合成酶差异表达还有待于下一步深入研究.尚有47 279个(约40%)Unigene未获得功能注释，与二代测序技术所获得的基因序列较短，不易获得同源基因信息有关; 本次研究是无参考基因组的转录组测序，由于目前还未见化州柚和柚的全基因组信息，只能参考相关物种的基因组数据进行分析，故所得注释信息有所欠缺，这也是目前药用植物转录组数据分析存在的共性问题^[31-35].

4.   结论

本研究首次获得了化州柚和柚的转录组数据，共得到116 202个Unigene，平均长度为742 nt，组装效果较好.通过与Nr、Swiss-Port、KOG、KEGG等数据库比对获得的生物信息学注释结果，揭示了化州柚和柚转录组主要的整体表达特征(约60%).