1. 引言
猕猴桃溃疡病是由丁香假单胞菌猕猴桃致病变种(Pseudomonas syringae pv. actinidiae, Psa)引起的病害，严重影响猕猴桃产业的发展 [1]。该病害对各国猕猴桃生产造成了巨大的经济损失 [2] [3] [4] [5]。近年来针对猕猴桃溃疡病的防治开展了一系列研究，在生物防治上，已经筛选出多种细菌、真菌、噬菌体等对Psa有显著抑制效果的拮抗菌株 [6] [7] [8] [9] [10]，在化学防治上针对新型药剂进行研发，研究发现多羟基双萘醛提取物对Psa具有抑制效果 [11]。为了开发出更高效防治药剂，针对猕猴桃溃疡病菌的致病机理进行了一系列研究，其中勘探寄主致病相关基因显得尤其重要。鞭毛作为细菌的运动器官，与病菌的定殖﹑侵染能力、粘附、诱导防卫反应和生物膜形成能力相关 [12]，是重要的致病因子 [13]。目前，利用RNA-Seq技术在转录组水平上已开展了大量植物相关基因及抗病基因挖掘研究 [14] [15]，这些研究为发掘参与抗病性的转录因子及基因提供了理论依据和研究思路。本研究以贵长猕猴桃为材料，以接种溃疡病菌后不同时间点的贵长猕猴桃叶片为供试材料，利用RNA-Seq技术进行测序，在获得大量差异表达基因的基础上，对其中的鞭毛相关差异表达基因进行GO功能注释和KEGG分析，以期为猕猴桃溃疡病菌鞭毛致病基因挖掘提供基础。
2. 材料与方法
2.1. 试验材料及处理
供试材料为一年生健康贵长猕猴桃，取自贵州省植物保护所植病课题组基地。猕猴桃溃疡病菌菌株丁香假单胞菌猕猴桃致病变种(Pseudomonas syringae pv. actinidiae, Psa)，强致病力菌株(菌株保存号：GZCC7520201)，分离自贵州省修文县洒坪镇红星村，保存于本实验室(数据未发表)。分别于接菌0 h、1 h、2 d、20 d后对猕猴桃进行取样，样本重量50 mg~100 mg之间。重复5次，采用液氮速冻，用于转录组测序。
2.2. RNA-Seq测序
待全部样本采集完成后，通过干冰将样品寄到上海美吉生物医药科技有限公司，委托其对样品进行RNA提取以及后续RNA-Seq测序。公司利用TruSeqTMStranded Total RNA Library Prep Kit试剂构建文库，文库利用Illumina HiSeq平台测序构建。
2.3. 差异表达基因鉴定
根据样品的测序数据，采用FPKM算法(fragments per kilobase of transcript per million mapped reads)和控制阈值来筛选差异表达基因。Fold Change表示不同处理之间各个基因表达量的比值，对差异表达基因的筛选阈值为|log2FoldChange| ≥ 1，P ≤ 0.05且FDR < 0.01。
2.4. 鞭毛相关基因的GO和KEGG通路分析
将筛选出的鞭毛相关差异表达基因注释到GO数据库和KEGG数据库，通过GO注释或KEGG注释对差异基因功能进行预测。
2.5. 实时荧光定量PCR (qRT-PCR)
从转录组测序分析结果中随机挑取5个差异基因，以验证转录组测序的可靠性。分别以康为世纪生物科技有限公司的RNA提取试剂盒(DNase I)提取预留样品RNA，采用TOYOBO公司的Rever TraAceq PCR-RTKit试剂盒反转录合成cDNA第一链，作为实时荧光PCR反应的模板。利用primerpremier5.0软件设计qRT-PCR引物，以gyrb为内参基因，根据Go Taqq PCR Master Mix，SYBR Green I步骤分析表达水平，在BioRad实时荧光定量PCR仪上运行，设置反应参数95℃~5 min，1循环；95℃~10 s、58℃~30 s，40个循环；熔解曲线(55℃，0.5℃/s，至95℃)。各样本各基因均作3次重复，参照Schmittgen的方法，使用2−∆∆Ct计算相对表达量 [16]。引物详情如表1。
Table 1. Real-time quantitative PCR detection primers
表1. 实时荧光定量检测部分基因引物
3. 结果及分析
3.1. 测序数据质量评估
将Psa分别接种贵长0 h、1 h、2 d和20 d后，共8个样品在Illumina Hiseq测序平台采用双端测序，共获得101.01 Gb Clean Data，各样品Clean Data均达到4.85Gb以上，Q30碱基百分比在96.77%以上。0h样品的原始测序序列raw reads都在3900万~1.37亿之间，随后将raw reads中含有带接头的、低质量的reads去掉，最终得到的clean reads数在3800万~1.26亿之间。8个样品经过测序得到的总reads数再过滤后得到的clean reads数比例较高，均在97.66%以上。Q20在99.05~99.28之间，表明此次测序质量较好，能满足后续分析得基本要求。其余各项指标结果详见表2。
Table 2. Summary of data quality of RNA-Seq of each sample
表2. 各样本RNA-Seq的数据质量情况汇总
注：Psa代表猕猴桃溃疡病菌，0 h、1 h、2 d、20 d代表接种Psa的不同时间点，1和2代表两次重复。
3.2. 鞭毛相关差异基因鉴定及分析
根据FPKM值筛选标准，在8个样本中筛选出鞭毛相关差异表达基因(differentially expressed genes, DEGs) 31个(表3)。31个鞭毛相关差异表达基因中，溃疡病菌1 hpi相对0 h有26个基因上调，4个基因下调；2 dpi对于0 h有23个基因上调，7个基因下调；20 dpi对于0 h有12个基因上调、19个基因下调。并且上调基因的个数多于下调基因个数。
Table 3. FPKM expression analysis
表3. FPKM表达量分析
3.3. 鞭毛相关基因的GO注释
对31个差异基因进行GO注释。GO注释结果表明：差异基因主要分布于生物过程(Biological process)、细胞成分(Cellular component)及分子功能(Molecular function) 3大类，其中参与生物过程方面的基因数目最多，其次是细胞组分和分子功能。在三大类一级分类下有很多的二级分类，其中在生物过程中有4个二级分类，细胞组分中有2个二级分类，分子功能中有3个二级分类。在生物过程中的功能参与细胞过程(cellular process)、移动(locomotion)、代谢过程(metabolic process)等，其中参与细胞过程的基因数目最多；在细胞组分方面，主要作为细胞解剖实体(cellular anatomical entity)发挥作用。在分子功能方面，具有催化活性(catalytic activity)、结合性(binding)等，其中注释催化活性(catalytic activity)的数目最多(图1)。
Figure 1. GO function classification of differentially expressed genes. Note: The abscissa in the figure indicates the secondary classification terminology of GO, the left ordinate indicates the percentage of genes contained in the secondary classification in the total, and the right ordinate indicates the number of genes compared with the secondary classification
图1. 差异表达基因的GO功能分类。注：图中横坐标表示GO的二级分类术语，左边纵坐标表示包含在该二级分类中的基因占总数的百分比，右边纵坐标表示比对上该二级分类的基因数量
3.4. 鞭毛相关基因的KEGG注释
对31个差异基因进行KEGG注释。KEGG注释结果表明：鞭毛相关基因KEGG代谢通路可分为4个主要通路：环境信息处理(Environmental Information Processing)，细胞过程(Cellular Processes)，生物体系统(Organismal Systems)等。4个主要通路下面被分为6个次级通路，包括细菌性传染病(Infectious disease: bacterial)、免疫系统(Immune system)、环境适应(Environmental adaptation)、细胞运动(Cell motility)、信号传导(Signal transduction)、膜运输(Membrane transport)，其中注释在细胞运动(Cell motility)中的基因数目最多，为20个。其次是信号传导(Signal transduction)，有3个。其余的功能中注释的数目均为1个(图2)。
Figure 2. KEGG pathway of differentially expressed genes. Note: The ordinate is the name of KEGG metabolic pathway, and the abscissa is the number of genes annotated to this pathway
图2. 差异表达基因的KEGG代谢通路。 注：纵坐标为KEGG代谢通路的名称，横坐标为注释到该通路的基因数目
3.5. 差异表达基因的实时荧光定量PCR (qRT-PCR)验证
为了验证转录组测序的可靠性，本研究选取了5个差异基因，使用qRT-PCR来验证转录组分析的可靠性。这些基因包括III型分泌系统蛋白(NZ708_RS06865)、III型分泌蛋白HrpF (NZ708_RS06815)、鞭毛钩长控制蛋白FliK (NZ708_RS06870)、III型效应蛋白HrpK (NZ708_RS06745)、丁香内酯生物合成蛋白AvrD1 (NZ708_RS03570)，以gyrb基因为内参基因。结果表明，上述5个基因在接种不同时间点的表达量不同，5个基因在接种2 d后的表达量均达到上调倍数最高，总体来说qRT-PCR与RNA-Seq的表达趋势基本一致。基于RNA-Seq分析的数据可靠性，差异表达基因可以用于下一步的分析(图3)。
4. 讨论
建立在新一代高通量测序平 Illumina基础上的RNA-Seq技术成为研究、挖掘新基因的重要手段。越来越多的人针对不同生物体转录组测序展开了一系列研究，其中包括分析苹果树腐烂病真菌与植物互作机制 [17]，检测石榴中的miRNA [18]，井赵斌针对溃疡病菌与毛花猕猴桃互作转录组分析的研究 [15]。目前，利用RNA-Seq技术在转录组水平上已开展了大量相关植物的抗病基因挖掘研究，但关于猕猴桃与溃疡病互作在转录组水平上的研究鲜有报道。
本研究利用贵长猕猴桃为材料，利用RNA-Seq技术，以接种溃疡病菌0 h、1 h、2 d和20 d的猕猴桃枝条为供试材料，从RNA转录组水平上分析溃疡病菌侵染后不同时间点鞭毛相关基因表达情况，31个鞭毛相关差异表达基因中，溃疡病菌1 hpi相对0 h有26个基因上调，4个基因下调；2 dpi对于0 h有23个基因上调，7个基因下调；20 dpi对于0 h有12个基因上调、19个基因下调。并且上调基因的个数多于下调基因个数。并随机选取了5个差异基因，使用qRT-PCR来验证转录组分析的可靠性。总体来说 qRT-PCR与RNA-Seq的表达趋势基本一致。
鞭毛相关基因GO生物学分析过程表明，参与生物过程的差异基因表达的数目最多，细胞组分和分子功能次之。参与细胞组分的差异表达基因主要作为细胞解剖实体(cellular anatomical entity)发挥作用。在分子功能方面，参与催化活性(catalytic activity)的数目最多。鞭毛相关基因KEGG注释结果表明：鞭毛相关基因KEGG代谢通路可分为4个主要通路，其次4个主要通路下面被分为6个次级通路，其中注释在细胞运动(Cell motility)中的基因数目最多，为20个。其次是信号传导(Signal transduction)，有3个。这与井赵斌针对溃疡病菌与毛花猕猴桃互作转录组分析的研究结果类似，表明在受到溃疡病菌侵染后，贵长猕猴桃体内多个生物学过程进行了响应 [15]。
基于上述研究，根据31鞭毛相关差异表达基因的GO注释以及KEGG注释情况，结合鞭毛在细菌中的作用，可用于后续克隆功能分析，从而对贵长猕猴桃与溃疡病菌之间的互作机制进行进一步了解，更好地理解Psa的致病机制，为猕猴桃溃疡病的科学防控提供理论依据。
参考文献
NOTES
*通讯作者。