干旱胁迫下野生大豆ABC转运蛋白转录组测序分析

张小芳 乔亚科 王冰冰 徐 燕 张 锴 李桂兰

(河北科技师范学院农学与生物科技学院,河北秦皇岛 066600)

野生大豆(Glycine soja)是栽培大豆的近缘野生种[1]。研究表明,野生大豆中携带着特有的新基因[2],而栽培大豆在人工驯化过程丢失了许多与环境适应相关的重要基因,包括一些与抗逆相关的基因[3]。因此对野生大豆抗逆基因的挖掘具有重要意义。

三磷酸腺苷结合盒(ATP-binding cassette, ABC)跨膜转运蛋白是植物细胞膜蛋白的重要组成成分,参与植物激素运输、气孔调节、次生代谢物的运输以及响应环境胁迫等过程[4]。最初人们对ABC 转运蛋白家族的研究多集中于多药耐药性研究[5-6],随着对ABC转运蛋白家族的进一步研究发现,植物中任何特定的生理过程都至少需要1 个ABC 转运蛋白[7],因此对植物ABC 转运蛋白的研究具有重要意义。

近年来,已有关于植物ABC 转运蛋白家族的结构和转运机制的研究报道[8-10],除对拟南芥[11]、水稻[12]的ABC 转运蛋白家族进行整体研究外,在铁皮石斛[13]、葡萄[14]和黄瓜[15]等中也开展了部分基因的研究,鉴定的ABC 转运蛋白主要在ABCB、ABCC 和ABCG 家族[16-19]中。目前大豆的基因组已经完成测序,但关于豆科ABC 转运蛋白各家族相关基因的研究尚鲜见报道[20-21]。

随着高通量转录组测序技术的发展,利用生物信息学对转录组测序的深入挖掘已成为研究重点[22-25],目前,已有多种具有特定功能的大豆或野生大豆基因被挖掘与利用[26-29]。本研究以抗旱性较强的野生大豆永46 为试验材料,对其转录组测序结果进行功能注释后,以ABC 转运蛋白为关键词进行筛选,对获得的ABC 转运蛋白进行多角度分析,探讨ABC 转运蛋白家族与干旱之间的关系,以期为ABC 转运蛋白家族基因的挖掘和利用提供一定的理论依据。

1 材料与方法

1.1 试验材料

以河北科技师范学院野生大豆种质资源课题组筛选鉴定的抗旱性强的野生大豆永46 为试验材料[30],采自河北省秦皇岛市昌黎县,并在河北科技师范学院植物细胞工程实验室保存。

1.2 转录组测序

将野生大豆永46 的种子播种于营养钵中,出苗后保留长势一致的幼苗10 株,设3 次生物学重复。定期浇定量霍格兰营养液,在苗龄达30 d 时,加入20%PEG 6000 进行干旱胁迫处理,分别在处理前(0 h)及处理后6、12、24 和48 h 取同一个叶位的复叶,按照RNA 提取试剂盒[TaKaRa(日本)]提供的方法提取野生大豆叶片总RNA。取2 μL 所得RNA 用1%琼脂糖凝胶电泳检测。检测合格后送至深圳华大基因科技有限公司进行测序,具体步骤参考张小芳等[31]的方法。

1.3 功能注释及ABC 转运蛋白关键词搜索

BLAST 比对到公共数据库,包括非冗余数据库(non-redundant, NR)、基因本体论数据库(Gene Ontology,GO)、京都基因与基因组百科全书数据库(Kyoto Encyclopedia of Genes and Gnomes,KEGG),之后以ABC 转运蛋白为关键词进行筛选,最终获得ABC转运蛋白家族相关信息。

1.4 数据分析

根据获得的ABC 转运蛋白家族信息,利用Excel 2007、ClustalX 和MEGA5.1 分析获得的ABC 转运蛋白基因簇的关系。

2 结果与分析

2.1 ABC 转运蛋白相关基因簇信息分析

对转录组中39 183 条基因簇以ABC 转运蛋白为关键词在NR、GO、KEGG 3 个数据库中进行筛选,结果共有337 条基因簇被注释,这些基因簇的表达量FPKM 值介于0～180.57 之间,序列拼接长度介于165～6 827 bp 之间。这些被注释的基因簇中,87.83%在3 个数据库中均有注释,3.86%只注释到了1 个数据库中。有299 条被GO 数据库注释、321 条被KEGG 数据库注释、146 条被NR 数据库注释。GO 本体中注释的细胞组分数据库有13 个,分子功能数据库有34 个,生物过程数据库有68 个;NR 数据库中注释的ABC 转运蛋白家族有8 个,包括ABCA、ABCB、ABCC、ABCD、ABCE、ABCF、ABCG 和ABCI(表1)。KEGG 注释的有28 个代谢途径分支,排名前10 的为K05658(226 个)、K00924(102 个)、K12843(81 个)、K05681(80 个)、K04733(54 个)、K05666(40 个)、K02065(20 个)、K13430(15 个)、K05643(13 个)、K05663(13 个)。

表1 ABC 转运蛋白家族相关基因簇统计Table 1 Statistical of Unigene about ABC transporter family

2.2 ABC 转运蛋白相关差异表达基因分析

对获取的原始基因簇信息进行差异表达基因(differentially expressed genes,DEGs)筛选,筛选阈值为P-value ≤0.05, | log2Fold Changes | ＞1,其中Fold Changes 为处理材料与未处理材料中基因表达量的比值,对于差异表达基因,其log2Fold Change＞0 时,则认为该差异表达基因是上调的, 反之, 若log2Fold Changes＜0,则认为该差异表达基因是下调的。筛选获得的DEGs 共182 个。干旱胁迫处理12 h 时涉及到的差异表达基因最丰富。在干旱胁迫处理12 h 时,差异表达基因为118 个,其中上调基因93 个,下调基因25个。干旱胁迫处理24 h 时,差异表达基因涉及范围变小,与未经干旱胁迫(0 h)相比,仅有36 个表现为差异表达,与干旱胁迫处理12 h 相比,上调基因减少,下调基因增多,说明部分基因在干旱胁迫处理24 h 时表达量上调,之后恢复原始表达量或低于原始表达量。干旱胁迫处理48 h 时,其变化趋势与干旱胁迫处理24 h相同,涉及的差异表达基因数无明显变化,较干旱胁迫处理24 h 表现为差异表达的基因仅有14 个(图1)。

图1 不同胁迫处理时期的差异表达基因统计分析Fig.1 Statistical analysis of DGEs under different stress treatment periods

2.3 ABC 转运蛋白相关差异表达基因功能分析

2.3.1 GO 功能显著性富集分析 GO 共有3 个本体(ontology), 分别描述基因的细胞组分(cellular component)、分子功能(molecular function)、生物过程(biological process)。在被注释的13 个细胞组分数据库、34 个分子功能数据库和68 个生物过程数据库中分别有11 个、28 个和50 个数据库注释到差异表达基因。由图2 可知,多数基因有分子功能和生物学过程注释。分子功能主要为 ATP 酶活性(ATPase activity)、蛋白激酶活性(protein kinase activity)、转运活性(transferase activity)、氧化活性(oxidoreductase activity)等;生物学过程主要为嘌呤核苷酸分解过程(purine ribonucleotide catabolic process)、磷酸化合物代谢过程( phosphate-containing compound metabolic process)、蛋白质修饰(cellular protein modification process)、囊泡运输(vesicle-mediated transport)、羧酸转运(carboxylic acid transmembrane transporter activity)、激素应答(response to hormone stimulus)等。

图2 3 个GO 本体中数量排名前10 的GO 分类统计情况Fig.2 The number of top 10 in 3 gene ontology classification

由表2 可知,6 个细胞组分数据库、4 个分子功能数据库和17 个生物过程数据库中注释的基因均表现为差异表达,在细胞组分本体中,主要定位在细胞膜系统中,与ABC 转运蛋白作用场所相符;在分子功能本体中,功能分类为转运活性或转移酶活性,与其作用功能相一致;在生物学过程本体中,注释的生物学过程多样,充分体现了ABC 转运蛋白的功能多样性。

2.3.2 KEGG 功能显著性富集分析 在KEGG 数据库中注释的28 个途径分支中,有9 个途径分支被注释到了ABC 转运蛋白途径(map02010)。这9 个途径分支分别是 K05641 ( ABCA1)、K05643(ABCA3)、K05658 ( ABCB1)、K05663 ( ATM)、K05656 ( ABCB9)、K05657 ( ABCB10)、K05666(ABCC2)、K05677(ABCD3)、K05681(ABCG2), 共包括ABCA、ABCB、ABCC、ABCD 和ABCG 5 个亚家族。K05658 途径在各个时期均出现差异表达基因,K05681 途径在各个时期均有上调基因分布,K05666 途径在各个时期均有下调基因分布,K05663 途径仅在干旱胁迫处理48 h/6 h 有下调基因分布,且下调基因仅有1 个。在各时期中,干旱胁迫处理6 h 和12 h 的差异表达基因最多,之后基因表达量变化范围变小,直至恢复到原始水平(表3)。

2.3.3 NR 功能聚类分析 对符合筛选条件的182 个差异表达基因在NR 数据库注释功能,获得差异表达基因的蛋白信息(有14 个没有注释,共统计168 个)。用邻接法对差异表达基因的氨基酸序列进行进化分析。分析结果表明,可将这168 个差异表达基因分为三大类25 小类。第一大类有3 小类,包括43 个差异表达基因,第二大类有6 小类,包括43 个差异表达基因,第三大类有16 小类,包括82 个差异表达基因。第一大类主要是ABC 转运蛋白家族的ABCE(1 个)、ABCF(5 个)、ABCG 家族(23 个)以及多药耐药蛋白(pleiotropic drug resistance protein);第二大类主要是ABC 转运蛋白家族的ABCB(21 个)与ABCC(17 个)家族;第三大类涉及基因最为广泛,除ABCA 家族(4个)、ABCD 家族(2 个)外,还包括锌指结构域蛋白(zinc finger CCCH domain-containing protein)、U-box 结构域蛋白(U-box domain-containing protein)、L 型外源凝集素受体蛋白(L-type lectin-domain containing receptor kinase)、富亮氨酸重复受体蛋白(leucine-rich repeat receptor-like serine/threonine-protein kinase)以及其他不同的蛋白激酶家族成员等。在进化分析时,NR注释功能为ABCI 家族的蛋白序列共有3 个(ABCI10、ABCI11、ABCI19),这3 个序列分别聚类在了第二大类、第一大类和第三大类中(图3)。

表2 不同GO 下的差异表达基因注释统计Table 2 Statistical of DEGs in different GO

对注释功能为ABC 转运蛋白家族的基因进行统计分析。除ABCC 家族外,各家族中涉及的差异表达基因多为上调基因;ABCA 家族共涉及4 个基因,3 个为上调基因,1 个为下调基因;ABCB 家族共涉及21 个基因,其中glysoja_016804 表达量差异倍数达-10.75,glysoja_ 023579 表达量差异倍数达10.54, 为注释到的基因中表达量变化最大的2 个基因;ABCC 家族共涉及17 个基因,表达量变化明显,最高差异表达倍数为9.32;ABCD 家族共涉及2 个基因,表达量变化幅度较小;ABCE 家族仅涉及1 个基因,且该基因为上调基因;ABCF 家族共涉及5 个基因,均表现为上调表达;ABCG 家族共涉及23 个基因,上调表达较多;ABCI家族共涉及3 个基因(表4)。

表3 差异表达基因的KEGG 通路富集分析Table 3 Statistical of DEGs in different pathway

表4 ABC 转运蛋白家族的差异表达基因分析Table 4 Statistical of DEGs in ABC transporter family

图3 差异表达基因的亲缘关系分析Fig.3 Genetic relationship analysis of DEGs

2.4 基因表达量分析

对同时注释到GO 和KEGG 数据库,且注释基因均差异表达的基因在干旱胁迫处理12 h 前的表达量变化进行汇总分析。由表5 可知,不同基因的原始表达量存在差异,glysoja_008864 和glysoja_010956 在干旱胁迫处理0 h 的FPKM 值均在100 以上,而glysoja_017978 和glysoja_047568 在干旱胁迫处理0 h 的FPKM 值均为0;多数基因在干旱胁迫处理12 h 表达量变化最明显,glysoja_005213 和glysoja_010956 则在6 h 的表达量变化最明显,而glysoja_008864、glysoja_010547 和glysoja_024284 在干旱胁迫处理12 h 前的表达量变化无明显差异;glysoja_010547 和glysoja_013434 的注释均为APK1A,其原始表达量差距不大,glysoja_002097 和glysoja_010956 的注释均为ABCC 家族,但原始表达量差异较大。

3 讨论

本研究系统分析了野生大豆8 类ABC 转运蛋白基因对干旱胁迫的响应,结果表明,不同处理时期对应基因种类和数目不同,同一基因在不同处理时期表达量也不同。在氨基酸进化分析时,聚类结果并未被分为8 类,说明在进化上有一定的保守性,但功能发生了改变,说明其结构域发生变化,对基因产物的功能产生了一定的影响。通过对野生大豆干旱胁迫下的ABC转运蛋白家族基因的系统研究,为挖掘具有抗旱功能的ABC 转运蛋白的研究提供了新思路。

本研究在GO 生物学过程聚类分析中,部分GO条目未被注释到,差异表达基因在5 个以上的GO 条目增加了GO：0009154(嘌呤核苷酸分解代谢过程)、GO：0006796 (含磷酸盐化合物代谢过程)、GO：0006464(细胞蛋白质变性过程)、GO：0006468(蛋白质磷酸化)、GO：0006952(防御反应)和GO：0010038(金属离子响应);在KEGG 途径中,ABCA、ABCB、ABCC、ABCD 和ABCG 家族成员注释为差异表达,表明这5个亚家族与干旱相关,注释的GO 本体与朱璐等[32]注释结果略有不同。

表5 差异表达基因表达量分析Table 5 Analysis on the expression of DEGs

Byrne 等[33]发现黑麦草在亚硒酸钠的诱导下ABCA 家族的基因表达量显著增加,Maathuis 等[34]发现AtATH14 与AtATH15 在盐胁迫下表达下调;Wu等[35]发现AtABCB19 可以使种子萌发后下胚轴快速的伸长。Lee 等[36]发现AtABCB14 可通过将苹果酸盐从质外体运输到保卫细胞, 提高渗透压来调控气孔对CO2的响应, 从而调控气孔的开合;AtABC5 也可以通过调节气孔开闭来增强植物的抗逆性[37];拟南芥中AtABCG40、AtPDR12 基因的表达能够提高植物的抗旱能力[2];水稻OsABCG36 基因也参与了植物的抗旱反应[38];大豆GmPDR12 基因在水杨酸及功能类似物质的诱导下能够快速大量表达[7]。本研究中,注释到的ABCA 家族基因既有上调也有下调,注释基因有ABCA1、ABCA2 和ABCA7;ABCA1 主要参与胆固醇的逆转运过程[39],AtABCA2(别名AtATH1)与拟南芥的开花相关[40],ABCA7 与人体的阿尔茨海默病相关[41],但在植物中的功能尚不清楚;注释到的ABCB 家族基因上调和下调明显,注释基因有ABCB1、ABCB4、ABCB6、ABCB11、ABCB13、ABCB25、ABCB26 等,与徐杏等[42]的研究结果一致;注释到的ABCC 家族除已知的与干旱胁迫相关的ABCC5 外, 还包括ABCC8、ABCC10、ABCC13 和ABCC14 等,AtABCC10 证实与病原反应相关[13], AtABCC13 与营养生长相关[32];注释到的ABCG 家族包含的基因种类最多,关于ABCG11、ABCG14、ABCG22、ABCG32、ABCG36 基因功能已有研究报道[11,43],本研究中还注释到了ABCG3、ABCG7、ABCG15、ABCG21、ABCG41(PDR13)等,AtABCG3 参与PC-Cd 的运输[44];VvABCG7 与葡萄的胚珠发育有关[11],注释到的ABCD 家族仅包含2 个基因(ABCD1和ABCD2),研究表明ABCD1 可能对植物与病原菌、有害动植物以及共生体间相互作用有关键作用[45],ABCD2 的功能仍未知,尚需要进一步的研究证明。

4 结论

本研究结果表明,野生大豆ABC 转运蛋白基因在干旱胁迫12 h 差异表达数量最多,且在24 h 和48 h无明显差异。ABC 转运蛋白的不同家族基因对干旱胁迫的响应存在差异,ABCB、ABCC 和ABCG 家族的响应基因最多,ABCE 家族的响应基因最少。出现差异表达的基因也是ABCB、ABCC 和ABCG 家族最多,ABCE 家族最少。在GO 生物学过程中,嘌呤核苷酸分解代谢过程的注释基因最多。在注释到KEGG 途径的ABC 转运蛋白途径中,K05658(ABCB1),K05681(ABCG)和K05666(ABCC)的注释基因最多,K05663(ATM)最少。本研究结果为抗旱基因的挖掘提供了新的基因资源。