甘蔗主要分布于热带和亚热带地区[1],在我国主要分布于广西、云南、广东、海南等省（区）,是我国主要的糖料作物。脱叶是甘蔗本身所具有的不同于一般落叶现象的性状,难脱叶的甘蔗在制糖过程中会含有大量杂质,对蔗糖的加工造成极大的干扰,影响制糖[2,3],我国糖厂规定甘蔗人工收获的含杂率不能超过1%,而对于机械收获则通常扣杂较高,可见,难脱叶的甘蔗极大影响甘蔗机械收获的质量和人工收获的效率。国外有在成熟期焚烧甘蔗老叶以助脱叶的处理方法,但会产生有毒气体,污染环境[4],同时会影响土壤微生物的活动,对土壤肥力产生不良影响[5,6];另一方面,甘蔗在成熟期叶片的脱落不但不影响产量,适当剥叶还可以提高产量和品质[7]。因此,在甘蔗新品种选育过程中把成熟叶片易脱落的性状作为选育的重要指标具有重要的生产意义。
甘蔗脱叶性状是数量性状,受多基因控制,植物脱落发生于离区[8],对于甘蔗而言,就是叶鞘和茎秆连接的部位。一般认为脱落的过程分成4个阶段：（1）离区的形成;（2）感受到脱落信号;（3）脱落信号的传导;（4）脱落后组织愈合。这4个过程都需要特定基因的表达,其中第2、3阶段在其他植物的叶、花、果等部位的脱落上有不少的研究[9,10]。集群分离分析（bulk segregant analysis, BSA）是QTL技术的一个分支,可以快速地对数量性状的主效基因进行定位,转录组学是研究基因的转录水平以及调控机制规律的技术,Bulked segrant RNA-Seq（BSR-Seq）将BSA技术和转录组学结合在一起,已经应用于小麦抗条锈病基因[11],玉米黄花突变基因的定位[12],大白菜黄化突变基因Brpem2的定位[13],猕猴桃抗寒机制的研究[14],小麦抗白粉病基因PmSESY的定位[15]等。
Li等[16]利用RNA-Seq技术发现了甘蔗脱叶相关基因在植物-病原菌相互作用、应激反应及与脱落酸相关通路中富集,关于甘蔗脱叶方面的研究报道还很少。有报道将甘蔗的脱叶性状分为三级,研究了不同组合选配方式的脱叶遗传情况,发现三级脱叶更易遗传给后代,其次是一级脱叶,最后是二级脱叶,也就是说中间类型的脱叶性状不容易稳定遗传[17]。本研究基于BSR-Seq技术对甘蔗脱叶相关的主效、稳定遗传的候选基因进行了筛选,并利用荧光定量分析对极端材料的+ 1~+7叶叶鞘离区部位进行了不同群体的分析,验证基因是否在不同群体中都具有稳定遗传的效应。通过对稳定遗传基因的生物信息学分析了解该基因的特点,为甘蔗脱叶的分子研究提供理论基础。
1 材料与方法
1.1 材料
甘蔗亲本（‘ROC25’ב云蔗89-7’）及其杂交后代F1（17 000多个品系）种植于广西大学扶绥基地,‘ROC25’表现出相对难脱叶的特点,‘云蔗89-7’则表现出相对易脱叶的特点,从F1中挑选极端难脱叶个体（叶鞘与茎秆包裹紧密,只有基部3~4片叶脱落,其他成熟叶鞘难以从茎秆上剥离）113个,易脱叶个体（+7~+10叶位以下全部脱落）116个,取其+6叶叶鞘离区暂时保存于液氮,带回实验室-80 ℃保存,同时将取样品系扩繁,观察1 a,最后筛选出难脱叶品系43个,易脱叶品系48个,选取其中的极端难脱叶品系5-94和极端易脱叶品系40-159的+1~+7叶位的叶鞘保存。取另一杂交群体（‘35-9’בCP08-1553’）的杂交后代个体中的极端难脱叶品系16-224和极端易脱叶品系16-226的+1~+7叶位的叶鞘用于验证。
主要试剂：RNA-easyTM Isolation Reagent Vazyme Cat（南京诺维赞生物科技股份有限公司）,反转录试剂盒TransScript® One-Step gDNA Removal and cDNA Synthesis SuperMix（全式金生物技术有限公司）,荧光定量试剂盒ChamQTM Universal SYBR qPCR Master Mix（南京诺维赞生物科技股份有限公司）。
1.2 方法
1.2.1 BSR-Seq建池 RNA提取采用TriZol法,提取后的RNA用Thermo的微量核酸纯度浓度测定仪检测浓度和纯度,利用琼脂糖凝胶电泳检测RNA完整性,各项指标合格后,将易脱叶的43个样品等量混合构建混池T03,难脱叶的48个样品等量混合构建混池T04,母本‘ROC25’单个样品存在记为T01,父本‘云蔗89-7’单独存在记做T02,最后送北京百迈客生物科技有限公司进行上机分析。
1.2.3 与参考基因组比对和候选区域定位 以甘蔗栽培品种‘R570’基因组为参考,利用STAR将过滤后的片段比对到‘R570’基因组上,将样品基因组组装到染色体水平。使用GATK[20]软件工具包对样品与参考基因组间及样品之间的单核苷酸多态性（SNP）进行检测,只有比对到参考基因组的唯一unigene才用于SNP挖掘。依靠欧式距离（Euclidean distance, ED）和SNP-index关联分析的方法分别筛选脱叶性状关联的区域,取2种关联方法的交集定位脱叶候选区域。
1.2.4 基因表达量和差异表达 基于负二项分布的基因模型,使用Cufflinks软件的Cuffquant和Cuffnorm组件,对基因和转录本的表达水平进行定量,由于抽取的转录本片段数目与测序数据量,转录本长度,转录表达水平有关,因此需要采用FPKM进行归一化[21],FPKM的计算公式为：
$FPKM =\frac{cDNA\ \ Fragments}{Mapped\ \ Fragments(Mllions) \times Transcript\ \ Length(kb)}$
公式中,cDNA Framents表示双端reads数目;Mapped Fragments表示mapped reads总数;Transcript Length表示转录本长度。
借助软件DEGseq检测差异表达基因,以q-value < 0.005且log2 Fold Change > 1作为筛选标准,筛选各个样品中的差异表达基因[22],通过聚类热图进行差异表达基因的可视化分析。利用COG、GO、KEGG、KOG、NR、Pfam等注释数据库,对样本得到的各类基因进行注释。
1.2.5 引物设计和qPCR 荧光定量引物的设计利用Primer Premier 5完成,使用TriZol法提取甘蔗叶鞘RNA。通过qRT-PCR技术检测不同甘蔗叶位的叶鞘中基因的相对表达水平,qRT-PCR反应溶液的配置参照定量试剂盒的说明书操作[23]。以25s基因作为内参,采用两步法在LightCycler®96仪器上运行,为了提高扩增效率,程序设为：95 ℃,30 s;95 ℃,10 s;60 ℃,60 s。记录CT值,以2‒ΔΔCT表示样品中基因的表达量。
1.2.6 候选基因的生物信息学分析 使用ORF finder（https://www.ncbi.nlm.nih.gov/orffinder/）对甘蔗脱叶相关的候选序列进行编码区查找,通过NCBI找到与该基因亲缘关系相近的其他物种的同类基因,利用MEGA-X构建进化树,并利用在线工具ITOL（https://itol.embl.de/）对构建的进化树进行处理;运用在线工具ExPASy（https://web.expasy.org/protparam/）分析其理化性质比如相对分子质量、等电点、脂肪族化、不稳定系数等指标;使用在线软件SOPMA（https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?page=npsa_sopma.html）分析其二级结构、SWISS-MODELR（https://swissmodel.expasy.org/interactive）分析三级结构、SignalP（http://www.cbs.dtu.dk/services/SignalP/）分析蛋白的信号肽、Server（http://www.cbs.dtu.dk/services/TMHMM/）预测跨膜结构、Plant-mPloc（http://www.csbio.sjtu.edu.cn/bioinf/plant-multi/）预测编码蛋白的亚细胞定位情况。
2 结果与分析
2.1 组装和比对
通过Illumina Hiseq 2000测序平台对4个样本进行测序,共产生60.46 Gb clean data,T01、T02、T03、T04的clean data分别为8.64、8.74、21.66、21.40 Gb,Q30值分别为93.67%、93.80%、93.97%、94.10%。GC含量分别为52.19%、53.15%、52.24%、51.71%,共产生的片段数为403 279 392。与参考基因组R570比对,比对上的片段比率分别为41.52%、40.28%、43.64%和43.34%（表1）。
表1 测序组装和比对情况统计表Tab. 1 Statistical table of sequencing assembly and comparison |
| 指标Index | T01 | T02 | T03 | T04 |
|---|---|---|---|---|
| 有效片段 | 28 845 092 | 29 167 652 | 72 248 183 | 71 378 769 |
| 高质量碱基 | 8 643 893 628 | 8 745 375 900 | 21 666 853 410 | 21 402 859 902 |
| Q30/% | 93.67 | 93.80 | 93.97 | 94.10 |
| GC/% | 52.19 | 53.15 | 52.24 | 51.71 |
| 总片段 | 57 690 184 | 58 335 304 | 144 496 366 | 142 757 538 |
| 比对片段 | 23 954 010 | 23 495 340 | 63 055 434 | 61 878 087 |
| 比对率/% | 41.52 | 40.28 | 43.64 | 43.34 |
2.2 差异基因表达和候选区域关联分析
对于所检测到的差异基因,样品T01 vs T02中,共有3223个差异基因,其中1621个基因上调,1602个基因下调。样品T03 vs T04中,共有842个差异基因,其中574个基因上调,268个基因下调。欧式距离关联分析得到的候选区域为14个,总长度为4.94 Mb,共包含328个基因,非同义突变SNP位点的基因42个。SNP-index关联分析共得到14个候选区域,总长度2.71 Mb,包含183个基因,非同义突变的基因20个。2种方法取交集得到4个与性状相关的侯选区域,总长度为1.40 Mb,候选区域内亲本之间非同义突变SNP共22个,混池间存在的非同义SNP共57个。候选区域共注释到基因86个,亲本间注释到的非同义突变15个。图1显示,在T01 vs T02差异基因、T03 vs T04差异基因,候选区间的交集中,有2个候选基因在亲本之间和子代混池中都表现出了差异,4个候选基因只在亲本间显示出了差异。
2.3 候选基因定量分析
在亲本之间和子代之间都存在差异的候选基因说明存在遗传效应,很可能与甘蔗脱叶有极其重要的关系。维恩图显示出了Sh09_g020080和Sh09_g020620这2个候选基因满足了各种筛选条件（图2）,特别是Sh09_g020080,无论是亲本之间还是子代混池之间都表现出易脱叶的表达量高于难脱叶,很可能存在稳定的遗传差异。Sh09_g020080的注释结果显示,具有σ因子活性,参与红光、蓝光、远红光反应,Sh09_g020620具有半胱氨酸内肽酶的活性,参与蛋白质水解的过程。Sh09_g020080和Sh09_g020620荧光定量引物设计如表2。BSR亚群体中选取难易脱叶品系的+1~+7叶进行荧光定量分析,结果如图3A显示,在+3、+4叶位,易脱叶品系中Sh09_g020620的表达量高于难脱叶品系,而在+1、+2、+5、+6、+7叶位,难脱叶品系的Sh09_g020620表达量高于易脱叶品系,在+7叶位中达到最大值。而对于Sh09_g020080来说,如图3B显示,在+2、+7叶位中,易脱叶品系的表达量相对较高,而在+1、+3、+4、+5、+6叶位中,则易脱叶品系的表达量较高,考虑到Sh09_g020080可能存在稳定的遗传效应,因此选择遗传基础一致的其他亚群体进行定量验证,其结果如图3C,在+1、+2、+3叶位中,难脱叶品系的Sh09_g020080表达量高于易脱叶品系。在+4、+5、+6、+7叶位中,易脱叶品系的Sh09_g020080表达量高于难脱叶品系。综合考虑,Sh09_g020080在难、易脱叶品系中存在稳定的遗传差异,推测很可能在甘蔗脱叶过程中是起着重要的作用。
图2 差异表达基因和候选区域基因的韦恩图T01 vs T02和T03 vs T04都表示差异基因集,CG（candidate gene）表示候选区域内的基因。 Fig. 2 Wayne map of differential expression genes and genes in candidate regions Both T01 vs T02 and T03 vs T04 represent differential gene sets, and CG (candidate gene) represents genes within candidate region. |
表2 候选基因和内参基因的引物序列Tab. 2 Primer sequences of candidate genes and internal reference genes |
| 基因ID Unigene ID | 正向引物（5'-3'）Former Primer (5'‒3') | 反向引物（5'-3'）Reverse Primer (5'‒3') |
|---|---|---|
| Sh09_g020080 | CGTTGGTCTGGTCTTTGGC | AATGGCGTGGACGGAATC |
| Sh09_g020620 | CGTCGTCATCAACCACCCA | TGTTCACGTTCACGTTCTTGC |
| 25s | ATAACCGCATCAGGTCTCCAAG3' | CCTCAGAGCCAATCCTTTTCC |
图3 Sh09_g020620和Sh09_g020080基因在难、易脱叶品系不同叶位中的相对表达量A：易脱叶品系为对照组的条件下,Sh09_g020620在+1~+7叶位的叶鞘离区的相对表达量;B：难脱叶品系为对照组的条件下,Sh09_g020080在+1~+7叶位的叶鞘离区的相对表达量（BSR亚群体）;C：难脱叶品系为对照组的条件下,Sh09_g020080在+1~+7叶位的叶鞘离区的相对表达量（验证群体）。 Fig. 3 Relative expression of Sh09_g020620 and Sh09_g020080 in different leaf positions of difficult and easy defoliation lines A: Relative expression of Sh09_g020620 in delamination region of +1-+7 leaf position under condition of the control group; B: Relative expression of Sh09_g020080 in extracellular zone of leaf sheath of +1-+7 leaf position (BSR subpopulation) under condition of control group; C: Relative expression of Sh09_g020080 in isolated area of leaf sheath of +1-+7 leaf position under the condition of the control group (verification population). |
2.4 甘蔗脱叶相关候选基因的理化性质分析
进化树分析发现Sh09_g020080与编码高粱的σ因子的亲缘关系最近（图4）,NCBI中显示,与XP_002440351.1、XP_021302430.1的一致性分别为85.53%和84.43%,可以通过与高粱的σ因子的特性相比较,了解甘蔗中该基因的特性和作用。利用软件分析发现Sh09_g020080 Unigene全长1866 bp,编码621个氨基酸,编码蛋白的相对分子质量为69 533.31 Da。该基因编码蛋白的理论等电点为10.01,据此可以判断为碱性蛋白质,不稳定指数显示为60.05,归为不稳定蛋白[24],平均疏水性为-0.665,为亲水蛋白。信号肽的分析发现有信号肽的可能性仅为0.0009,可以推断出蛋白不具备信号肽的功能,跨膜结构分析显示无跨膜结构,σ因子的作用是与RNA的转录有关,起到和核心酶启动子结合的作用[25],在甘蔗中可能会借助其他物质辅助。亚细胞定位显示该蛋白存在于叶绿体,这和注释结果的情况是一致的,该蛋白的作用是参与红光、远红光、蓝光反应,所在的位置决定了它的功能。二级结构分析发现,该蛋白具有α-螺旋、延伸链、β转角和无规则卷曲,比例分别为47.18%、12.24%、4.83%和35.75%。图5显示,三级结构显示有大量的α-螺旋结构组成一些单独的区域,中心部分也是以α-螺旋为主,这些部位很可能就是Sh09_g020080的保守区域。
3 讨论
无论是机械收获还是人工收获,成熟期甘蔗叶片的自然脱落都可以很大程度地减少成本,降低含杂率。BSA技术广泛使用于各种作物研究,如白菜头型[26]、向日葵耐盐性[27]等。BSR技术已经在小麦[11]、玉米[12]、白菜薹[28]、大白菜[29]等作物上得以应用。本研究利用BSR技术筛选出2个主效候选基因（Sh09_g020620,Sh09_g020080）,其中Sh09_g020620在亲本之间的差异状况和子代混池中的差异是异向的,具体表现为：亲本中为易脱叶品种的表达量高于难脱叶,子代混池中则相反。选取BSR亚群体进行单一品系的荧光定量分析,结果显示+6叶位叶鞘的表达和子代混池一致,推测该基因的表达可能受环境或基因互作影响,在遗传过程中其表达量发生了变化。Sh09_g020080在子代混池之间和亲本内却表现出了同向的差异（易脱叶品系的表达量比难脱叶的高）,选取BSR亚群体单一品系进行分析,与子代混池的差异状况是一致的。为了进一步验证Sh09_g020080是否具有不同群体的稳定差异遗传效应,选取‘B35-9’בCP08-1553’的杂交后代个体中的极端易脱叶品系16-226和极端难脱叶品系16-224进行荧光定量分析,整体的表达趋势为易脱叶品系的表达量高于难脱叶,和BSR群体中单一品系的不同叶位整体表达趋势一致。据此,可以初步判断出Sh09_g020080在难易脱叶品系中的差异是稳定存在的。
Sh09_g020080参与红光、蓝光、远红光反应,σ因子是和光合作用相关的,σ因子具有介导叶绿素积累的作用[31],推测Sh09_g020080可能也具有积累叶绿素的作用。根据Sh09_g020080的定量结果,该基因在难易脱叶品系中的表达是有差异的,在易脱叶品系中的整体表达趋势要高于难脱叶,可以推测易脱叶品系有更强的光合作用,相对少量的叶片足够满足光合所需,因此,在成熟期易脱叶品系下部有较多的叶片自动脱落。由于关于甘蔗脱叶的研究较少,而且σ因子的研究大多针对菌类,所以只能推测Sh09_g020080和甘蔗脱叶相关。甘蔗脱叶性状复杂,受较多基因控制,通过本研究可以为甘蔗脱叶研究提供新的见解。