黄萎病胁迫下抗病陆地棉茎组织lncRNA的鉴定与分析

王国宁, 张艳, 宋俊丽, 杨君, 王省芬, 吴立强, 张桂寅*

(1.河北农业大学，华北作物改良与调控国家重点实验室，河北 保定 071001; 2.涿州市农业农村局，河北 涿州 072750)

棉花是重要的多用途经济作物，在世界范围被广泛种植。棉花生长发育过程中易遭受多种病菌危害而降低产量和纤维品质，尤其黄萎病菌(Verticilliumdahliae)危害最为严重。黄萎病菌是一种土传半营养型真菌，导致严重的黄萎病(Verticillium wilt)[1-2]。

长链非编码RNA(long non-coding RNA, lncRNA)是一类长度大于200 bp的非编码RNA[3]，能够在表观遗传水平、转录水平和转录后水平调控基因表达，在生长发育、病原菌防御等过程中起重要作用。研究发现，TalncRNA18、TalncRNA73、TalncRNA106和TalncRNA1084对小麦条锈病胁迫应答起重要调控作用[4]。油菜中发现一些lncRNA作为miRNA前体参与核盘菌胁迫应答[5]。拟南芥lncRNA ELENA1通过调控PR1提高其对假单胞菌的抗性[6]。Jiang等[7]揭示番茄lncRNA23468通过抑制miR482b而影响NBS-LRRs的积累，最终提高番茄对晚疫病菌的抗性。

近些年来，棉花根组织中发现大量lncRNA响应黄萎病菌侵染。Zhang等[2]通过转录组分析发现，抗病海岛棉和感病陆地棉中分别有1 236和1 907条lncRNA参与根组织对黄萎病菌胁迫应答。本课题组前期对抗病陆地棉转录组分析发现，4 277条lncRNA参与根组织拮抗黄萎病菌侵染[8]。黄萎病菌侵染根组织后，通过茎组织扩展至叶，并伴随病症的发生。Bae等[9]研究发现，黄萎病菌在不同抗性马铃薯品种茎组织中的扩展速度不同，表明茎组织在寄主拮抗黄萎病菌扩展中起重要作用。Zhang等[10]发现，与感病陆地棉相比，抗病海岛棉的茎部组织对黄萎病菌具有显著的抑制作用。然而，棉花茎组织抗黄萎病相关lncRNA还未见报道。本研究通过转录组分析发现971条lncRNA参与茎组织拮抗黄萎病菌侵染，并对其功能机制进行了初步解析，为抗黄萎病相关lncRNA功能研究提供了理论依据，为棉花抗病性遗传改良提供了新的见解。

1 材料与方法

1.1 植物材料与病原菌培养

供试棉花品种为抗病陆地棉农大601，由河北农业大学棉花品种创新与产业化团队培育[11]。通过浸种、催芽处理的种子于10 h光照/14 h黑暗和25 ℃条件下的温室进行培养。生长至4片真叶时，选取健康一致的幼苗进行黄萎病菌接种处理，同时以接水处理为对照组。本课题组前期鉴定的强致病力菌系临西2-1[12]于PDA平板培养基上25 ℃活化培养15 d，转接到Czapek培养液中，25 ℃震荡(130 r·min-1)培养10 d左右，将孢子悬浮液浓度稀释到1.0×107个孢子·mL-1，现配现用。

Czapek培养液配方：NaNO33 mg·mL-1，KH2PO41 mg·mL-1，MgSO41 mg·mL-1，FeSO4·7H2O 0.02 mg·mL-1，KCl 1 mg·mL-1，蔗糖30 mg·mL-1。

1.2 RNA提取与测序

根据前期对黄萎病菌侵染过程的研究结果[10]，分别于接种处理后2、6、12、24和48 h(hpi, hours post infection)取棉花茎组织，分别命名为VD2hpi、VD6hpi、VD12hpi、VD24hpi和VD48hpi；同时接水处理组相应时间点棉花茎组织也被收取(MT2hpi、MT6hpi、MT12hpi、MT24hpi和 MT48hpi)。利用RNA prep Pure Plant Kit试剂盒(北京天根生化科技有限公司)提取棉花茎组织总RNA，使用Ribo-zero rRNA Removal Kit试剂盒(美国Epicentre公司)去除rRNA，应用NEBNext Ultra Directional RNA Library Prep Kit试剂盒(美国NEB公司)构建cDNA文库，采用Illumina Hiseq 4000平台进行转录组测序。

1.3 lncRNA鉴定

为保证后续组装质量，提高生物信息分析结果的科学性，将原始数据去除含接头、含ploy-N reads和低质量的reads，获得高质量的clean reads。采用TopHat2[13]将clean reads对齐到棉花参考基因组[14]，已比对上的clean reads用Cufflinks[15]进行组装和FPKM(Fragments Per Kilobase of exon per Million fragments mapped)分析。依据lncRNA特征，对完成组装的转录本进行筛选和质量控制：①去除长度<200 bp转录本；②通过Cuffcompare软件筛除与数据库注释外显子区域有重叠的转录本；③去除FPKM<0.5的转录本；④去除具有编码潜能的转录本。

1.4 lncRNA差异表达及功能分析

使用Cuffdiff软件对lncRNA进行差异表达分析，筛选阈值设定为q值<0.05，且差异倍数≥2。q值为校正后的p值。为了阐明lncRNA的生物学功能，对差异表达lncRNA进行顺式和反式靶基因预测，采用KOBAS 3.00(http://kobas.cbi.pku.edu.cn/)对靶基因进行GO富集分析。lncRNA上、下游100 kb内的mRNA通过STEM软件进行筛选，得到顺式靶基因。采用Pearson相关系数法分析lncRNA与mRNA表达水平的相关性，以相关系数绝对值>0.95为阈值，筛选反式靶基因。利用UniProt数据库(https://www.uniprot.org/)搜索候选lncRNA靶基因在其他植物中的同源基因的功能，以预测靶基因在棉花植株中的功能。

1.5 qPCR

使用RNA prep Pure Plant Kit试剂盒提取高质量总RNA。各样本分别取0.5 μg总RNA，利用PrimeScript RT-PCR kit反转录试剂盒获得cDNA。利用Oligo 7设计特异性引物(表1)，采用SYBR©PremixExTaqTMⅡ实时定量试剂盒，在ABI 7500荧光定量PCR仪上进行扩增与基因表达定量分析。设置反应程序：95 ℃预变性2 min；95 ℃变性10 s，55～60 ℃退火20 s，72 ℃延伸20 s，循环40次。每个样品检测3次生物学重复，基因表达水平以GhUBQ14为内参进行归一化处理[16]。

表1 qPCR引物

2 结果与分析

2.1 获得971条棉花茎lncRNA响应黄萎病菌胁迫

不同处理不同时间点的10个样本转录组测序共获得1 164 310 812条reads，从中筛选得到1 126 308 443条clean reads。84.0%～90.5% clean reads比对到参考基因组(表2)。通过组装、筛选共得到12 772条lncRNA。使用Cuffdiff软件对lncRNA进行差异表达分析，共筛选出971条差异表达lncRNA(图1)。

图1 差异表达lncRNA诱导表达谱

表2 棉花茎组织测序数据

2.2 差异表达lncRNA响应黄萎病菌的时序性

从图2可以看出，差异表达lncRNA随黄萎病菌侵染进程而变化，下调表达lncRNA较上调表达lncRNA变化更为显著。整体趋势为响应黄萎病菌侵染的lncRNA数量随时间逐渐增多，至24 h达到高峰，48 h开始减少。进一步分析发现，接菌后不同时间点间存在高比例特异性差异表达lncRNA，尤其24 h达到67.18%。表明棉花茎组织lncRNA响应黄萎病菌侵染具有明显的时序性。

图2 差异表达lncRNA特征

2.3 lncRNA对靶基因的调控方式

lncRNA能够通过顺式作用和反式作用调控靶基因表达[17]，因此预测了差异表达lncRNA的顺式靶基因和反式靶基因。共鉴定108条lncRNA靶调控314条mRNA，构成580对lncRNA-mRNA调控关系。其中23对为顺式作用，557对为反式作用，表明抗病相关lncRNA对靶mRNA的调控方式以反式作用为主。通过对lncRNA和mRNA之间调控关系的研究，检测到40条lncRNA调控作用较为简单，每条lncRNA只调控1条mRNA；其余68条lncRNA调控作用较复杂，每条lncRNA能够靶向调控多条mRNA(图3)。相关性分析显示，17个顺式lncRNA正调控靶基因，5个顺式lncRNA负调控靶基因，而反式lncRNA均正调控靶mRNA。这些结果表明，lncRNA主要通过反式正调控作用影响靶mRNA表达，参与棉花茎组织抵御黄萎病菌侵染。

图3 差异表达lncRNA对mRNA的调控网络

2.4 GO富集分析

为解析差异表达lncRNA的潜在功能，对其靶mRNA进行了GO富集分析。GO注释可分为生物过程、细胞组分和分子功能三类。生物过程类别中，主要富集到对几丁质、缺水、缺氧等刺激应答反应，JA、ABA等植物激素应答反应，对真菌的防御反应和植物型超敏反应；细胞组分中，细胞膜、细胞壁、胞间连丝、叶绿体、过氧化物酶体、细胞胞液、内质网、高尔基体等被显著富集；分子功能主要包括脱氢酶、脱水酶、激酶、转移酶等酶活性，DNA结合、血红素结合、钙调素结合、蛋白结合等结合功能(图4)。以上结果表明，棉花茎组织差异表达lncRNA通过多种分子功能响应由黄萎病菌诱导的各种胁迫和激素信号。

2.5 候选lncRNA功能分析

大量研究表明，转录因子、激酶和细胞色素P450家族蛋白在植物生长发育和抗逆反应中起重要作用，本研究结果显示，激酶、转录因子和血红素结合GO条目被显著富集(图4)。为此进一步对靶调控这三类基因的18条lncRNA进行了深入分析。转录组数据显示，除了lnc_004530、lnc_004943和lnc_013707外，这些候选lncRNA受黄萎病菌诱导表达模式为：2～24 hpi上调表达或变化不显著，48 hpi下调表达或变化不显著(图5)。

图4 差异表达lncRNA靶基因GO富集分析

注：“In”表示诱导表达；“Lo”表示丧失表达。

2.5.1对激酶基因的调控作用 激酶在植物抗逆信号传递中起重要作用。通过对差异表达lncRNA靶基因注释分析，鉴定出15条lncRNA正调控CDPK、CRK、LECRK和WAK等32个激酶基因(图3，表3)。利用UniProt数据库搜索这些激酶基因的功能，结果发现，LecRK-IX.1、CDPK4、CRK5和CRK6在拟南芥等植物中的同源基因，通过调控H2O2积累、细胞程序性死亡、病程相关基因表达、气孔运动等过程参与植物抗病。而CRK10、CRK15与CRK4、CRK5、CRK6属于同一CRK亚家族[18]。基于以上结果，推测lnc_010753、lnc_004943、lnc_001364、lnc_003725、lnc_001285通过调节LecRK-IX.1、CRK10、CRK15和CDPK4等激酶基因的表达，激活下游一系列抗逆反应，最终增强棉花茎组织对黄萎病菌的抗性(图6)。

图6 候选lncRNA通过调节激酶基因表达介导棉花抗黄萎病的潜在机制

表3 部分候选lncRNA对激酶基因表达的调控作用

2.5.2对转录因子基因的调控作用 本研究共鉴定9条lncRNA对ERF、NAC和WRKY等25个转录因子基因具有潜在的反式正调控作用(图3，表4)。利用UniProt数据库搜索这些转录因子的功能发现，WRKY6、WRKY33、ERF1A和NAC072转录因子，在拟南芥等植物中通过激活抗病相关基因的转录提高植物抗病性。依据这些结果，推测lnc_004653、lnc_004530、lnc_007238、lnc_001285和lnc_013701通过调节WRKY6、WRKY33、ERF1A和NAC072转录因子基因的表达，激活下游抗逆基因的转录，最终增强棉花茎组织对黄萎病菌的抗性(图7)。

图7 候选lncRNA通过调节转录因子基因表达介导棉花抗黄萎病的潜在机制

表4 部分候选lncRNA对转录因子基因表达的调控作用

2.5.3对细胞色素P450基因的调控作用 本研究发现，lnc_003725、lnc_004943、lnc_001285、lnc_007238、lnc_011764、lnc_004530通过反式作用正调控CYP710、CYP71、CYP82和CYP75等9个细胞色素P450家族基因(图3，表5)。利用UniProt数据库查询这些细胞色素P450的功能，发现CYP710A1和CYP75B1分别催化合成抗病相关物质豆甾醇和槲皮素。大豆CYP82A3正调控病程相关基因PR3和PR4的表达[19]。根据以上结果推测，lnc_003725、lnc_004943、lnc_001285、lnc_007238、lnc_011764和lnc_004530正调控CYP82A3、CYP710A1和CYP75B1基因表达，进而上调抗病相关基因PR3和PR4表达或催化抗病物质的合成，提高棉花对黄萎病菌的抗性(图8)。

注：B—β-谷甾醇；S—豆甾醇；D—二氢山萘酚；Q—槲皮素。

表5 部分候选lncRNA对细胞色素P450基因表达的调控作用

2.6 qPCR验证

2.6.1lncRNA验证 为了进一步验证抗黄萎病相关lncRNA的表达，选取分别与激酶、转录因子和细胞色素P450基因表达相关的 lnc_011764、lnc_010753和lnc_013707进行qPCR定量分析。结果显示，黄萎病菌胁迫下lnc_011764在2、6、12和24 hpi显著上调表达，48 hpi显著下调表达；lnc_010753在2 hpi显著上调表达，48 hpi显著下调表达；lnc_013707在6 hpi显著下调表达，12和24 hpi显著上调表达(图9)。这些结果进一步揭示了lnc_011764、lnc_010753和lnc_013707参与棉花对黄萎病菌的防御反应。此外，qPCR检测结果与高通量测序结果基本一致(图5和9)，表明基于高通量测序的表达量分析结果质量较高。

2.6.2靶mRNA验证 根据lncRNA-mRNA靶关系预测结果，挑选lnc_011764、lnc_010753及lnc_013707的靶基因Gh_D05G2983、Gh_D05G1893、Gh_D13G0477和Gh_D13G0071进行qPCR分析。结果显示，Gh_D05G2983在2、6和12 hpi显著上调表达，48 hpi显著下调表达；Gh_D05G1893在2、6、12和24 hpi显著上调表达，48 hpi显著下调表达；Gh_D13G0477在2、6和12 hpi显著上调表达，48 hpi显著下调表达；Gh_D13G0071在12 hpi显著上调表达(图10)。对比图10和图9发现，Gh_D13G0071与其靶调控子lnc_013707的诱导表达规律不一致，Gh_D05G1893、Gh_D05G2983分别与lnc_011764及Gh_D13G0477与lnc_010753表达规律基本一致。这些结果为进一步研究lnc_011764与lnc_010753的抗病机制奠定了理论基础。

3 讨论

长链非编码RNA(lncRNA)作为一类重要的调控分子，通过多种机制广泛参与植物的生长发育和抗逆反应等生物过程[17, 20]。近些年，棉花中也报道了大量与纤维发育、生物胁迫和非生物胁迫相关lncRNA[2, 8, 21-24]。棉花黄萎病菌是一种广泛存在、极具破坏性的土传真菌病害，对棉纤维的产量和品质造成严重影响。黄萎病菌首先攻击寄主的根组织，因此前期研究主要集中在棉花根组织lncRNA对黄萎病菌的抗性。然而不少研究表明，茎组织在寄主拮抗黄萎病菌扩展过程中起重要作用，因此研究茎组织lncRNA的抗病功能亦十分重要。本研究首次表征了茎组织lncRNA在抗黄萎病过程中的重要作用。

黄萎病菌对植物的侵染包括根际附着、定植、孢子萌发、新菌丝延伸及向上部茎、叶组织扩展等过程，植物同时也进行一系列复杂有序的调控机制进行免疫反应[8-10, 25-26]。本研究显示，茎组织在抵抗黄萎病菌侵染过程中具有明显的时序性(图2)。与本实验室前期根组织抗病相关lncRNA进行比较，发现有474条(48.82%)lncRNA只在茎组织被诱导表达，表明棉花lncRNA在抵抗黄萎病菌侵染时具有时空特异性。因此，组织特异性诱导表达lncRNA和非特异性诱导表达lncRNA的抗病功能和机制可待进一步深入研究。

近些年研究鉴定了大量植物lncRNA，但具体功能大多还不十分清楚。本研究利用lncRNA靶基因进行GO富集分析，以帮助了解lncRNA潜在的抗病功能。并揭示了lncRNA在胞内、胞外和细胞膜多种细胞组分，通过多种机制响应病原菌诱导的各类刺激及激素信号，参与棉花对黄萎病菌的免疫反应。同时响应细菌、冷等生物、非生物胁迫的防御反应也被显著富集(未列出)，表明部分抗黄萎病相关lncRNA可能具有多抗性。另外，茎组织抗病相关lncRNA的数量显著少于根组织[8]，推测经过根组织对黄萎病菌的识别与互作，茎组织对病菌响应的目标性更强，可减少不必要的能量消耗，这更有利于寄主的存活。

本研究发现，lnc_004943等5条lncRNA对激酶CDPK4、CRK10、CRK15和LecRK_IX.1的编码基因具有潜在的调控作用；lnc_004653等5条lncRNA潜在调控转录因子WRKY6、WRKY33、ERF1A和NAC072的编码基因；lnc_007238等6条lncRNA通过反式作用正调控CYP82A3、CYP75B1和CYP710A1编码基因的表达(表2～4)。利用UniProt数据库对这些靶基因进行了功能预测，结果表明这些靶基因具有潜在抗病功能。据报道马铃薯CDPK4和CDPK5通过RBOH调节活性氧(reactive oxygen species，ROS)簇参与植物免疫[27]。LecRK-IX.1可诱导H2O2积累和细胞死亡，介导拟南芥对疫霉(Phytophthora)的抗性[28]。CRK5可正调控PR1表达并激活超敏反应引起的细胞死亡，提高拟南芥对丁香假单胞菌(Pseudomonassyringae)的抗性[29]；CRK4和CRK6通过调节ROS积累和气孔关闭，增强拟南芥对病原菌的免疫作用[30]；而CRK10、CRK15与CRK4、CRK5、CRK6属于同一CRK亚家族[18]。研究表明，拟南芥WRKY33通过激活类渗调蛋白(OLP)和PR3基因表达介导拟南芥抵御灰霉病菌(Botrytiscinerea)和黑斑病菌(Alternariabrassicicola)[31]。拟南芥WRKY6通过NPR1正调控病原菌防御基因PR1[32]。拟南芥乙烯响应转录因子ERF1通过碱性几丁质酶(CHIT-B/PR3)响应丁香假单胞菌[33]。拟南芥NAC072可能通过调控EDR1响应缺水胁迫[34]。据报道，大豆CYP82A3通过正调控PR3和PR4基因表达，提高转基因烟草对寄生疫霉(P.parasitica)的抗性[19]。拟南芥CYP710A1可催化β-谷甾醇合成豆甾醇，进而提高寄主对丁香假单胞菌的抗性[35-37]。CYP75B1可催化二氢山萘酚合成槲皮素[38]，而槲皮素具有真菌抗性[39]。结合这些激酶、转录因子和细胞色素P450家族基因的已有报道，本研究构建了lnc_004943、lnc_004653、lnc_007238等10条棉花茎组织lncRNA可能的抗病机制简图(图7～8)。采用qPCR技术分别对lnc_010753、lnc_013707和lnc_011764的靶基因进行了定量分析，结果进一步揭示了lnc_011764对棉花CYP82A3(Gh_D05G1893、Gh_D05G2983)及lnc_010753对棉花LecRK-IX.1(Gh_D13G0477)具有潜在的调控作用；而lnc_013707与预测靶基因Gh_D13G0071的诱导表达规律不相同，可能是因为lnc_013707不是Gh_D13G0071的主要调控因子，是否存在其他主效调控因子有待进一步研究。本研究为深度解析这些lncRNA的抗病机理提供了重要理论依据，将通过过表达和沉默实验进一步验证lncRNA功能。