茎上承枝叶下接根部,在植物生长过程中主要起运输和支撑的作用。茎的生长包括初生生长和次生生长。初生生长是顶端分生组织,经细胞分裂分化形成表皮、皮层和维管柱的过程[1]。国外学者通过对毛白杨（Populus tomentosa）各节间进行转录组测序,发现其茎主枝伸长区含特殊的碳水化合物活性酶（carbohydrate active enzymes, CAZymes）及扩展蛋白（expansins）家族成员基因。其中CAZymes中的糖基水解酶（glycoside hydrolases, GHs）小家族基因在茎伸长区显著性表达,部分成员可能参与初级细胞壁代谢（primary cell wall metabolish）、细胞壁疏松（cell wall loosening）和细胞伸长（cell elongation）[2]。Expansins是植物细胞扩张和疏松的关键调控因子[3,4]。而在毛白杨茎发现数个Expansins基因在茎伸长区表达,可能参与初级细胞壁的合成和修饰[2]。三羧酸循环在毛白杨茎发育早期呈下降趋势,可能与丙酮酸脱氢酶（pyruvate dehydrogenase, PHD）、异柠檬酸脱氢酶（isocitrate dehydrogenase, ICDH）和琥珀酸脱氢酶（succinate dehydrogenase, SDH）蛋白丰度的降低有关[5]。
维管组织的出现,为植物适应陆地生活提供了结构基础[6]。光合作用是绿色植物维持正常生长所必需的能量和物质来源,而叶绿体的数目、形态与结构以及叶绿素含量等是影响此过程的重要因素[7]。叶绿素缺乏会导致植物发生白化,不能进行正常的光合作用[8]。白化现象在植物界普遍存在,目前在苦瓜（Momordica charantia）[9]、玉米（Zea mays）[10]、拟南芥（Arabidopsis thaliana）[11]、棉花（Gossypium spp）[12]及椰子（Cocos nucifera）[13]等植物上已有与白化现象相关的研究,并发现导致基因突变的过程主要涉及叶绿素合成与代谢、叶绿体发育以及其他生物色素合成[14]。然而这些研究主要在草本植物中开展,在木本植物中报道较少。
波罗蜜（Artocarpus heterophyllus）是桑科（Moraceae）木波罗属（Artocarpus）的常绿乔木,集医药、林木和果蔬鲜食于一体[15]。付影等[16]在野外发现了波罗蜜叶绿素缺失突变体（chlorophyll deficient mutant, CDM）并对其性状进行了相关研究,其中叶绿素含量明显降低、脯氨酸含量明显增加。谢柳青等[17]对波罗蜜CDM茎次生生长进行了转录组分析,并且近期波罗蜜基因组草图数据已公布,为其后续在分子生物学方面的研究奠定基础[18]。本研究通过转录组测序,分析叶绿素合成缺失对波罗蜜嫩茎发育的影响,为光合作用对木本植物嫩茎发育的研究提供数据参考。
1 材料与方法
1.1 材料
图1 波罗蜜正常苗与叶绿素缺失突变体茎节间表型性状CK：正常波罗蜜;CDM：波罗蜜叶绿素缺失突变体;虚线以上部分为嫩茎（CK为嫩绿色,CDM为白色）;标尺为10 cm。 Fig. 1 Stem internode phenotypic characters of normal seedlings and chlorophyll deficient mutants of A. heterophyllus CK: Normal A. heterophyllus; CDM: Chlorophyll deficient mutant of A. heterophyllus; Above the dotted line is primary stem (CK is light green, CDM is white); The scale is 10 cm. |
1.1.3 数据组装 参考波罗蜜茎叶三代全长转录组（PRJNA579273）将二代数据进行拼接。将原始测序数据进行过滤后得到有效数据。使用Trinity软件对有效数据进行组装,min_kmer_cov设置为3,长度短于150 bp的contigs将被舍弃。然后使用Tgicl对转录本进行聚类去冗余得到较高质量的Unigene序列[21]。以上流程均由华大基因公司完成。
本研究对5个样品进行分析,CK的2个生物学重复分别为：CK1和CK2,CDM的3个生物学重复分别为：CDM1、CDM2和CDM3。
1.2 差异基因筛选及分析
过滤平均表达量≤0.5的基因,通过TTEST函数计算P-value值（自由度为3）,以P-value< 0.001和log2|Fold change|≥1为标准筛选出全部的差异基因。
1.2.1 WGCNA分析 采用宋长新等[22]的方法在R软件中应用WGCNA软件包,对差异基因表达谱数据及表型数据,建立基因模块,并探索其与表型间的关联。本研究参数β范围设为1~30,最终取值为30,将距离较近的模块合并（mergeCutHeight=0.25）,据此构建加权基因共表达网络。
1.2.2 功能注释与分析 将差异基因进行GO、KEGG功能注释、分类及富集分析。GO分类：利用GO数据库对鉴定所得的差异基因按细胞组分、分子功能和生物学过程进行分类,生成相应的分组柱状图。同时利用KEGG数据库对鉴定所得的差异基因按代谢、遗传信息处理等过程进行分类,生成相应的分组柱状图。GO显著性富集：①将所有差异基因向GO数据库的每个Term映射,计算各Term的Unigenes数;②用超几何检验,与整个基因组相比,找出显著富集的GO条目;③P-value≤0.05的GO条目为显著性富集。同时将差异基因注释到KEGG数据库,其显著性富集分析和阈值与GO相同[23]。
2 结果与分析
2.1 测序质量分析
结果显示,CK1、CK2、CDM1、CDM2、CDM3测序原始数据分别为72.28 M、64.08 M、72.69 M、70.18 M和70.18 M,去接头和低质量序列后分别得到69.29 M、60.77 M、69.95 M、67.45 M和67.58 M高质量数据,且Q20均高于97%,Q30均高于88%,数据比对效率均大于94%。各样品的GC含量稳定在38.5%~42.5%范围内,且N50长度均高于1200 bp（表1）。同组样品基因表达模式趋于一致且Pearson相关系数均大于0.72,主成分分析（principal component analysis, PCA）结果表明两组处理之间是相互独立的（图2）。
表1 波罗蜜嫩茎测序结果的数据分析Tab. 1 RNA-Seq dates from primary stem tissues of A. heterophyllus |
| 样品 Sample | 原始数据 Total reads/M | 过滤后数据 Clean reads/M | 转录本 Transcripts | Q20质量值 Q20 quality value/% | Q30质量值 Q30 quality value/% | GC含量 GC content/% | N50值 N50 value/bp | 数据对比效率 Clean read ratio/% |
|---|---|---|---|---|---|---|---|---|
| CK1 | 72.28 | 69.29 | 212 772 | 97.48 | 89.30 | 38.89 | 1 353 | 95.86 |
| CK2 | 64.08 | 60.77 | 231 308 | 97.25 | 88.68 | 39.17 | 1 218 | 94.83 |
| CDM1 | 72.69 | 69.95 | 64 907 | 97.50 | 89.51 | 42.24 | 1 933 | 96.24 |
| CDM2 | 70.18 | 67.45 | 74 169 | 97.55 | 89.63 | 41.91 | 1 961 | 96.11 |
| CDM3 | 70.18 | 67.58 | 79 299 | 97.48 | 89.36 | 41.72 | 2 039 | 96.29 |
将全部295 869个Unigenes比对至各功能数据库进行注释,共有174 291个Unigenes获得注释,获得注释的数据库分别是NR（152 036,51.39%）、NT（156 922,53.04%）、Swissprot（111 271,37.61%）、KEGG（119 331,40.33%）、KOG（119 737,40.47%）、Pfam（110 505,37.35%）和GO（101 591,34.34%）（表2）。
表2 波罗蜜嫩茎转录组的7个数据库BLAST结果Tab. 2 Summary for the BLAST results of A. heterophyllus primary stem transcriptome against seven databases |
| 注释的数据库 Annotated database | 注释的基因数量 Annotated gene number | 基因比例 Gene percentage/% |
|---|---|---|
| NR | 152 036 | 51.39 |
| NT | 156 922 | 53.04 |
| Swissprot | 111 271 | 37.61 |
| KEGG | 119 331 | 40.33 |
| KOG | 119 737 | 40.47 |
| Pfam | 110 505 | 37.35 |
| GO | 101 591 | 34.34 |
| 总注释基因 | 174 291 | 58.91 |
2.2 差异基因的筛选与分析
图3 波罗蜜嫩茎差异表达基因统计与分析A：CK和CDM中不同表达量的基因数;B：CK和CDM差异表达基因的火山图;C：差异基因在各部分中表达韦恩图;D：差异基因长度的比较;E：基因聚类树以及其模块;F：各样品中差异基因表达热图及特征向量;Ⅰ：CK中特有基因;Ⅱ：CK与CDM共有基因在CDM中下调;Ⅲ：CK与CDM共有基因在CDM中上调;Ⅳ：CDM中特有基因。 Fig. 3 Statistics and analysis of DEGs in the primary stem of A. heterophyllus A: Number of genes with different expression levels in CK and CDM; B: Volcano map of DEGs between CK and CDM; C: Venn map of DEGs in each part; D: Comparison of DEGs length; E: Gene cluster tree and its module; F: Heat map and feature vector of DEGs in each sample;Ⅰ: Unique gene in CK; Ⅱ: Common genes of CK and CDM are down-regulated in CDM; Ⅲ: Common genes of CK and CDM are up-regulated in CDM; Ⅳ: Unique gene in CDM. |
2.3 差异基因GO、KEGG功能注释
2.3.1 GO功能注释 将22 988个差异基因进行GO功能注释和分析,可确定其重要的生物学功能。结果表明,共有3712个基因获得注释,并划分为48个功能组（图4）。因功能组的不同又可进一步分为三大主类：分子功能（molecular function）、细胞组分（cellular component）和生物学过程（biological process）。三类注释分别占79.85%、56.17%和46.82%,表示大部分差异基因与分子功能相关。分子功能注释的基因大部分集中在催化活性（catalytic activity）,且与CK相比,CDM中120个基因上调,1787个基因下调,其中CDM特有的有16个。在生物学过程中,以细胞过程（cellular process）的基因数最多,且与CK相比,CDM中大部分差异基因下调,少数上调,其中CDM特有的有8个;在细胞组分中,膜（membrane）功能组拥有的基因数最多,且CDM相对于CK发生大量下调,少数上调,其中CDM特有的有6个。
2.3.2 KEGG功能注释 对差异基因进行KEGG途径分析,根据对结果中的代谢（metabolism）、遗传信息处理（genetic information process）、环境信息处理（environmental information processing）、细胞过程（cellular processes）以及生物系统（organismal systems）5组类目进行归类,有助于深入了解CK与CDM间差异基因所在的通路情况（图5）。结果表明,2 080个Unigenes被注释到19个条目,其中10个在Ⅰ、Ⅱ、Ⅲ和Ⅳ 4部分中均有出现。代谢类目注释的Unigenes最多,有1152个（55.38%）,且以糖代谢（carbohydrate metabolism）注释的Unigenes最多,其中CK和CDM特有的分别有109和2个。而在遗传信息处理中,有698个（33.56%）Unigenes得到注释,翻译（translation）所含的Unigenes最多,且CDM相比于CK发生大量下调,少量上调,CDM特有的有2个。环境信息处理、细胞过程及组织系统所注释的Pathway较其他分类少,分别为2、1、1条。
2.4 GO富集分析
为进一步研究波罗蜜CDM嫩茎发育过程中基因的生物学功能,将差异基因进行GO富集分析（图6）。与CK相比,CDM有379个基因上调,22 609个基因下调。通过Q-value值分别筛选出20条最显著的条目。结果表明,在上调基因中,20个GO条目富集显著程度相同,而翻译（translation）和信号转导（signal transduction）是富集基因数较多的GO条目。在下调基因中,单加氧酶活性（monooxygenase activity）、多细胞生物发育（multicellular organism development）和作用于成对供体上的氧化还原酶活性,并结合或还原氧分子（oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen）3个GO条目最显著富集,而膜的组成成分（interal component of membrane）是富集基因最多的GO条目。
2.5 KEGG富集分析
通过KEGG Pathway富集分析,确定差异基因参与的最主要途径,在差异基因中找出显著富集的pathway（图7）。KEGG富集分析结果显示,脂肪酸降解（fatty acid degradation）、光合作用（photosynthesis）、苯丙氨酸代谢（phenylalanine metabolism）和色氨酸代谢（tryptophan metabolism）是上调基因和下调基因共有的富集代谢通路。在上调基因中,获得了5个氧化磷酸化（oxidative phosphorylation）途径的Unigenes,是富集基因最多的代谢通路。而在下调基因中,油菜素甾醇生物合成（brassinosteroid biosynthesis）和苯丙烷生物合成（phenylpropanoid biosynthesis）是最为显著富集的代谢通路。其次,二萜类生物合成（diterpenoid biosynthesis）、角质,琥珀和蜡的生物合成（cutin, suberine and wax biosynthesis）、黄酮类生物合成（flavonoid biosynthesis）和玉米素生物合成（zeatin biosynthesis）等是较为显著富集的通路,苯丙烷生物合成（phenylpropanoid biosynthesis）是富集基因最多的代谢通路。
3 讨论
近年来,随着转录组测序技术的发展,为波罗蜜分子生物学方面研究提供极大便利。学者通过对波罗蜜花被进行RNA-seq分析,探究其糖代谢相关基因差异表达的特点[24]。本研究通过转录组测序,分析了叶绿素合成缺失对波罗蜜嫩茎发育的影响。