土壤盐渍化是一种全球范围内的生态环境问题.在全世界范围内约有7%的土地和20%的耕地受到盐渍化的威胁^[1].盐胁迫可以影响植物的种子萌发、生长、分化以及发育的各个阶段^[2].提高植物内在耐盐基因的表达量，是盐渍化土地资源开发利用的有效方式之一.

棉花是我国尤其是新疆重要的经济作物之一，在我国棉花供应中发挥着重要的作用.但在棉花生产过程中，经常遭受来自于土壤的盐碱胁迫作用，土壤盐渍化一直是制约新疆的棉花质量和产量的重要环境限制因子.因此，通过挖掘和筛选棉花耐盐相关基因，阐明它们在应答盐胁迫中的分子调控机制，并利用相关耐盐基因资源进行耐盐棉花品种改良，为培育出对盐胁迫耐受性更强的棉花新品种具有十分重要的意义.

目前，利用转录组技术进行盐胁迫相关基因筛选和挖掘已经在多种植物中得到广泛的应用^[3]，LIU等^[4]对大豆叶片和根在盐胁迫下转录组进行分析，发现使碳代谢和氮代谢相关基因的表达发生了显著变化，表明碳代谢和氮代谢参与了大豆对盐胁迫的应答过程. PRIYANKA等^[5]对在150 mmol/L NaCl处理条件下的葡萄叶片进行转录组分析，发现其碳水化合物代谢、能量代谢、氨基酸代谢、脂类代谢途径和次生代谢产物的合成等代谢通路参与葡萄应答盐胁迫的过程.陆地棉全基因组测序的完成也极大地推进了棉花功能基因的筛选^[6].目前为止，有关棉花的抗逆基因研究取得了较大的进展，尤其是棉花耐盐相关的分子生物学研究取得了很多成果，如LONG等^[7]对亚洲棉、雷蒙德氏棉和陆地棉的NHX基因进行了全面、系统的比较研究，分别从这3种棉中鉴定了12、12和23个类NHX蛋白；基因结构分析表明，NHX基因内含子较多，可形成选择性剪接，有利于植株适应土壤高盐环境，其中GhNHX1定位于棉株的液泡系统，受盐胁迫的强烈诱导. GhNHX1的沉默增强了棉花幼苗对高盐浓度的敏感性.在转基因棉花中过表达K2-NhaD，能提高棉花植株的耐盐、耐旱性，增加胁迫相关基因的表达，改善SOS途径等代谢途径^[8].有关新疆主栽棉花品种的耐盐基因筛选研究鲜有报道，本文以新疆主栽棉花品种为研究对象，以转录组学方法为基础，筛选新疆棉花中的耐盐相关基因，为更深入地了解棉花的耐盐分子机制及选育棉花耐盐新品种提供理论依据.

1.   材料与方法

1.1   材料

以南疆地区主栽4个棉花品种(‘新陆中42号’、‘新陆中69号’、‘新陆中72号’及‘国审棉206-5’)为供试材料.

选取籽粒饱满、大小一致的棉花种子，70%的酒精消毒1 min后，用30%的双氧水浸泡1 h，再用无菌水冲洗3~5遍，无菌环境下放至种子萌发露白，在超净工作台上把露白的棉花种皮去掉，分别放入MS固体培养基和含0.50%NaCl的MS固体培养基上，对照为不含NaCl的培养基.每个培养基放3~4粒种子. 28 ℃培养，16 h光照/ 8 h黑暗，培养6 d.

1.2   棉花种子萌发期抗盐性分析

将棉花种子分别在不同盐质量分数(0.25%、0.50%、0.75%、1.00%、1.25%、1.50% NaCl)的培养基中培养，拍照并统计其萌发率.

1.3   总RNA提取

根据不同棉花品种的抗盐性，对棉花培养6 d的‘新陆中69号’棉花根系，用天根公司植物多酚多糖RNA提取试剂盒提取总RNA，电泳检测合格后，由上海派森诺公司进行转录组测序分析.

1.4   测序文库构建及转录组测序数据分析

选择4个样本进行实验，包含2个重复对照(CK1、CK2)和2个处理(NaCl1、NaCl2)，将待测序列样品用Illumina Hiseq测序平台上机测序之前把每个样本的mRNA打断成300 bp左右大小的片段，然后反转录成cDNA，最后每个片段加上一个标签序列，就构建成样本文库，文库名按照CK1、CK2、NaCl1和NaCl2的顺序分别命名为LRA01810、LRA01811、LRA01812和LRA01813.待序列下机之后就可以按照文库的标签序列区分不同的样本列.

对下机数据进行CASAVA碱基识别分析，获得测序数据序列信息，利用fastQC软件对原始数据进行质量分析^[9]，而后利用NGSQC Toolkit去除低质量碱基后得到Clean Reads^[10].使用Hisat2软件将Clean Reads比对到陆地棉基因组上，而后再使用Stringtie软件对转录本进行组装^[11].使用Htseq软件对得到的转录本进行表达定量分析. Q20和Q30分别代表被99%和99.9%测准确的碱基数. N代表测序仪不能确切地测定出某一个碱基时的百分比(%)，测定不出的碱基就会标注为N.

1.5   差异表达基因筛选及功能富集分析

采用R语言DESeq程序包对表达基因进行差异分析，筛选差异表达基因条件为：表达差异倍数|log 2(Fold Change)|>1，显著性P < 0.05^[9].对得到的差异基因进行GO注释和KEGG通路分析.

1.6   基因表达检测

利用上述方法提取‘新陆中69号’棉花根系RNA，5X All-In- One RT MasterMix反转录试剂盒(购买自Applied Biological Materials (abm) Inc.)进行cDNA合成.根据已知的基因序列，设计RealTime PCR引物(表 1)，按照TaKaRa公司的SYBR Premix Ex Taq^TM II说明书进行Real Time PCR. PCR反应条件为：95 ℃ 30 s，1个循环；95 ℃ 5 s，60 ℃ 34 s，共39个循环.利用2^-ΔΔCT方法分析基因的相对表达量.

表  1  检测基因表达的引物序列

Table  1.  The primer of gene expression detected

基因编号	正向序列(5′-3′)	反向序列(5′-3′)	PCR产物长度/bp
UBQ7	AGGCATTCCACCTGACCAA	TCGCACCTCCCCTGAGAC	131
08054	GAGGTGAATAAAGGGAGGATGG	GCCTGTCGGGCTTGTGAT	166
03482	TTTTGTGACGGAGCACGATT	ACCATCTACTACCTGCCCACTG	141
30608	CATGTCTCCTCGGATTGCC	TGGTATGCCGAATTGTTGTTAA	126
18245	GCGTTGCTGTGGCTTTATCT	TCGCTTACCGTCACTTCTGG	193
33843	CCTGTTTACGGCTGTGCG	TGGAGGAGACGACGGTGC	168
44291	GGTGCCTGAAAACATTGTGAA	TCTTCCAGTCTCCAACCTTCC	149
05721	GAAACAATGGAAGAATGGGACA	GAGCATCATAGCAATGGTAGCC	160
49287	AGGTGGCTCCCAAGGTTCT	TTGAGTGCATCTCCAATCAAAA	148
55203	GCAAGAGGATGTGGAGGGC	TTTGGGGACATTGTGATGGA	118
39457	TAGGTGGAGGAATTGAAAGGG	GGATCAGACCAAGCAGAGGG	159
54195	GGGAACGAAACCTCACCAAG	GGAAGCAACGGCGACAAT	134
56209	TTCCCACTGTTGTCATTATCGG	TTGTGCTGCCTCACTTGCTT	152
17346	AAGATGTTCGATCAAAAGACGG	CAAATGCCCAGCAGAGTGC	161
50458	AACAACCCTTCTACGCCACG	AATGCCCGCAATGCTGATA	160
51234	GAGTTCGTGGATACGCCATTAC	CATTTGGGTTTGTTGATTGAGC	156
34390	GAGCCAGGGAAGGTGTTAGC	GGATCTCCCCTGTTTGATGC	115
47521	TGCCCAACAACCTTTACTACCC	CCATCCTCAGACCAGCTCCTAT	181
55296	ATCGTTGTGGCTCCTACCG	GCTCCCATCTACTGTTTCTTCG	138
59144	ATTTCCTAGAGTGCAGCCACAT	ACAAGCCCGTTCCCCAGT	198
42639	CCCAATGGAAAGGGTGCA	TCAGGGTCCCCAGGAGAAG	164

下载: 导出CSV 

| 显示表格

2.   结果与分析

2.1   棉花品种间的耐盐性比较及耐盐品种筛选

通过对4种棉花种子在不同盐质量分数条件下培养，发现‘新陆中69号’棉花种子萌发最多(图 1A).随着处理的盐质量分数增加，4种棉花种子的萌发率都下降，其中‘新陆中69号’棉花种子的萌发率最高(图 1B)，表明其具有较强的耐盐性.

图  1  棉花品种间耐盐性比较

Figure  1.  The comparison of salt tolerance among cotton varieties

下载: 全尺寸图片 幻灯片

2.2   盐胁迫处理的棉花根系测序质量分析

选择经0.50% NaCl处理的新陆中69号品种根系提取RNA，进行转录组分析. Illumina Hiseq平台测序结果表明：CK1、CK2、NaCl1、NaCl2获得的原始Reads分别是49 841 082、41 348 682、42 165 466、45 793 106条，去除接头序列和低质量碱基后，获得的Reads分别是49 665 682、41 187 572、41 939 240和45 577 200条，且Q20和Q30分别达到97%和94%以上，说明转录组测序质量较高(表 2).

表  2  数据过滤统计

Table  2.  The data filtering statistics

样本	Reads/条	碱基数/bp	Q30/bp	N/%	Q20占碱基数的 百分数/%	Q30占碱基数的 百分数/%	Clean Reads/条	Clean Reads 占Reads的百分数/%
CK1	49 841 082	7 526 003 382	7 151 088 931	0.001 900	98.01	95.01	49 665 682	99.64
CK2	41 348 682	6 243 650 982	5 938 337 538	0.001 053	98.03	95.11	41 187 572	99.61
NaCl1	42 165 466	6 366 985 366	5 992 756 836	0.002 118	97.61	94.12	41 939 240	99.46
NaCl2	45 793 106	6 914 759 006	6 525 598 019	0.001 877	97.74	94.37	45 577 200	99.52

下载: 导出CSV 

| 显示表格

2.3   盐胁迫下棉花根系差异表达基因分析

与对照材料相比(图 2)，经过0.50% NaCl处理6 d后的棉花植株中有1 557个基因表达量发生显著变化.其中，在盐胁迫条件下，有891个基因表达量增加，而666个基因表达量出现下调.

图  2  棉花根系差异表达基因的火山图

Figure  2.  The volcano plot of differentially expressed genes in cotton roots

下载: 全尺寸图片 幻灯片

2.4   差异表达基因GO和KEGG注释分析

对得到的棉花根系差异表达基因进行GO注释分析(图 3A)，结果表明：在细胞组分分类中，差异基因共注释到65个条目，其中注释到最多的是细胞膜组分，共有173个基因；在分子功能分类中，差异基因共注释到261个条目，其中共有162个基因功能注释为氧化还原酶类，这说明氧化还原酶类基因在对于棉花应答盐胁迫是十分重要的；在生物学过程分类中，注释到最多的是代谢进程，共有440个基因，当植物受到盐胁迫后，会通过改变相应的代谢进程对盐胁迫做出应答，从而提高自身对盐胁迫的抵御能力.因此，参与植物各种代谢进程的基因在植物应答盐胁迫过程中有着重要的意义.

图  3  棉花根系差异表达基因的GO和KEGG注释分析

Figure  3.  The GO and KEGG pathway annotation analysis of differentially expressed genes in cotton roots

下载: 全尺寸图片 幻灯片

KEGG分析的结果显示：差异基因富集到20条代谢通路上，主要集中在木质素生物合成途径，该途径在植物体的机械支持、水分运输以及抵抗各种不良环境影响中起着很重要的作用；其中，有64个差异基因富集到苯丙烷合成代谢途径，是富集差异基因最多的代谢途径.此外，半乳糖代谢、胡萝卜素生物合成、淀粉和蔗糖代谢、鞘脂代谢、类黄酮生物合成、MAPK信号通路等代谢途径也富集到较多的差异基因(图 3B)，说明这些代谢途径可能在棉花应答盐胁迫过程中具有重要的作用.

2.5   棉花响应盐胁迫期间差异基因表达变化

利用差异倍数大于2、校正后的P小于0.01进行差异基因的筛选，选取对照和盐胁迫条件下基因表达差异性最大的10个差异基因(表 3)，分别在NCBI数据库中进行序列比对和基因功能注释分析，结果表明：编号为Cot-AD 05721、Cot-AD 54195和Cot-AD 47521的基因参与了ABA信号通路，其中，编号为Cot-AD 05721的基因属于富含半胱氨酸类受体激酶(CRK)基因家族，该基因家族编码的蛋白是植物体内普遍存在的一类蛋白激酶，是许多信号识别传递途径中的关键组分，作为识别信号的质膜受体，能够感受外界刺激，通过磷酸化作用参与胞内信号传递^[12]；Cot-AD 54195基因编码的蛋白主要参与调控环核苷酸离子通道(CNG)，在Ca²⁺响应以及抗盐胁迫信号传导中发挥重要作用^[13]；Cot-AD 47521基因编码的蛋白是受体蛋白激酶家族中LRR型类受体蛋白激酶，能够参与植物胁迫应答过程^[14]. Cot-AD 49287编码角鲨烯环氧酶(SQE)，可催化鲨烯转变为环氧化鲨烯.拟南芥SQE1可参与调节根和种子萌发^[15]. Cot-AD 56209可编码硫氧还蛋白. Cot-AD 17346编码的蛋白具有丙氨酰-tRNA合成酶活性. Cot-AD 50458编码SPX1-like蛋白，AtSPX1在拟南芥叶片衰老中发挥重要作用^[16]. Cot-AD 34390编码的是ABCG家族基因，可参与对ABA转运. Cot-AD 55296编码AMSH-like泛素硫酯酶，拟南芥AMSH3参与泛素化膜蛋白的降解过程^[17]. Cot-AD 59144编码的是细胞色素P450家族基因，是一种末端加氧酶，可参与生物体内甾醇类激素合成等过程.

表  3  盐胁迫下筛选的10个差异基因

Table  3.  The 10 DEGs screened under salt stress

编号	mRNA表达	染色体	起始位点	终止位点	长度/bp	预测蛋白	预测蛋白的功能
CotAD_05721	Up	At_chr9	84082252	84092166	2 067	CRK37; BRI1	Protein kinase, catalytic domain; Serine/threonine- / dual-specificity protein kinase, catalytic domain; Gnk2-homologous domain; Serine/threonine-protein kinase, active site; Protein kinase, ATP binding site; Tyrosine-protein kinase, catalytic domain
CotAD_49287	Up	At_chr8	68076229	68080579	1 605	SQE1;SQLE; ERG1	Aromatic-ring hydroxylase-like; Squalene epoxidase
CotAD_54195	Down	Dt_chr8	34057090	34062151	2 160	CNGC1; CNGC	IQ motif, EF-hand binding site; Cyclic nucleotide-binding domain; Ion transport
CotAD_56209	Down	At_chr5	41436225	41438394	1 713	At1g60420; CPT2	C1-like; Thioredoxin-like fold; Thioredoxin, conserved site
CotAD_17346	Down	scaffold568.1	799268	801407	804	alaS; AARS	Alanyl-tRNA synthetase, class IIc, N-terminal; Alanyl-tRNA synthetase, class IIc, core domain
CotAD_50458	Down	Dt_chr9	68556782	68559227	792	SPX1; PHO81	SPX, N-terminal
CotAD_34390	Down	Dt_chr7	35388212	35394633	2 022	ABCG11; ABCG2; CD338	ABC transporter-like; AAA+ ATPase domain; ABC-2 type transporter; ABC transporter, conserved site
CotAD_47521	Down	Dt_chr1	77812160	77813734	1 575	FER; IRAK4	Peroxidase, active site; Malectin-like carbohydrate-binding domain
CotAD_55296	Down	Dt_chr5	21522746	21527212	894	AMSH2; STAMBP; AMSH	JAB1/Mov34/MPN/PAD-1
CotAD_59144	Down	scaffold2496.1	211466	214066	1 536	CYP72A1; CYP735A	Cytochrome P450; Cytochrome P450, E-class, group I; Cytochrome P450, conserved site

下载: 导出CSV 

| 显示表格

利用qRT-PCR检测棉花盐胁迫处理后上述基因的表达变化，结果表明：编号为Cot-AD 05721和Cot-AD 49287的基因上调表达；编号为Cot-AD 54195、Cot-AD 56209、Cot-AD 17346、Cot-AD 34390、Cot-AD 50458、Cot-AD 55296、Cot-AD 59144和Cot-AD 47521的基因下调表达. 10个基因的表达量变化趋势均与转录组结果一致，统计相关性系数R²为0.923 4，具有良好的线性关系，说明转录组测序的质量较高(图 4).

图  4  qRT-PCR检测在棉花盐胁迫时差异基因的表达

Figure  4.  The expression analysis of different genes detected by qRT-PCR

下载: 全尺寸图片 幻灯片

3.   讨论

当植物受到盐胁迫后，植物会通过调控耐盐相关基因表达来应对盐胁迫的威胁.植物的根部是决定整个植株生长发育的关键器官，由于根部直接接触土壤中的盐离子，也是接受盐胁迫信号的最初位点，因此对盐胁迫的反应也最为敏感、准确^{[2, 18]}.基于此，本文选择棉花根部作为实验取材部位，旨在通过转录组学方法从棉花根系中筛选出棉花耐盐相关基因.所备选的4个新疆主栽棉花品种中，‘新陆中69号’品种表现出较好的耐盐性，适合作为转录组筛选耐盐基因的实验材料.通过将对照组和盐胁迫处理组的棉花根系转录组测序结果进行分析，结果从‘新陆中69号’中筛选到1 557个差异表达基因，其中有891个基因表达量上调，666个基因表达量下调.有64个差异基因富集到苯丙烷类合成代谢途径中，是富集差异基因最多的代谢途径，这说明在棉花抵御盐胁迫过程中，苯丙烷类合成代谢途径发挥了重要的调控功能.其他差异表达基因多集中在半乳糖代谢、胡萝卜素生物合成、淀粉和蔗糖代谢等途径.

通过qRT-PCR方法对棉花响应盐胁迫时差异性较大的10个基因表达分析的变化结果与转录组分析结果一致.其中，编号为Cot-AD 05721的基因属于CRK家族基因，可参与细胞生物和非生物应激反应.拟南芥CRK2基因能提高种子萌发期幼苗的耐盐性，并调节幼苗的根长. CRK2是盐诱导胼胝质沉积所必需的关键基因，通过胼胝质沉积对盐分增加的响应作用，证明CRK2基因在种子萌发过程中耐盐性的重要性^[19].本文中转录组结果显示：棉花在盐胁迫后编号为Cot-AD 05721的基因表达量上调，推测该基因在棉花耐盐胁迫应答分子调控过程中发挥一定的作用.编号为Cot-AD 54195的基因属于调控环核苷酸离子通道CNG家族基因，CNG家族在抗盐胁迫信号传导中发挥重要作用.编号为Cot-AD 47521的基因属于LRR型受体蛋白激酶家族成员，JUNGA等^[20]发现辣椒中LRR受体蛋白激酶基因CaLRR1的表达受高盐的胁迫诱导，说明LRR可作为信号识别受体在多种信号途径中起作用.

拟南芥植物中角鲨烯环氧化酶(SQE)可以环氧化角鲨烯为氧化三萜类化合物，在植物根和下胚轴伸长过程中发挥重要作用^[15].本研究表明：编号为Cot-AD 49287的基因是一类编码SQE的基因，该基因表达量显著上调，表明其功能可能与棉花种子萌发期的根系伸长调控过程有密切关系.编号为Cot-AD 56209的基因属可编码硫氧还蛋白，调控细胞氧化还原稳态，而Cot-AD 17346编码的蛋白则具有丙氨酰-tRNA合成酶活性.目前，有关这2个基因在棉花根系发育和盐胁迫中的调控过程如何发挥作用，未见相关文献报道，作用机理尚不明确.

盐渍环境下，植物根系的盐分与营养元素之间产生协同或拮抗作用.在盐碱条件下，植物根系会降低对必需的大量营养元素磷的吸收能力，使磷的吸收和运输以及在体内的分配遭到破坏, 导致营养失衡并加速植物叶片的衰老过程^[21].拟南芥衰老叶片中存在显著的磷酸盐饥饿诱导基因AtSPX1. AtSPX1是拟南芥叶片衰老、Pi饥饿和SA信号通路之间的关键调控因子^[16].本文研究表明：Cot-AD 50458是一类编码SPX1-like蛋白的因子，该基因表达量在盐胁迫环境下表达量显著下调，推测可能是在盐胁迫的作用下抑制了棉花对磷元素的吸收功能，增加棉花叶片的衰老过程. ABCG家族基因，可参与植物对ABA的转运，其中ABCG25和ABCG40是高亲和脱落酸(ABA)转运蛋白. ABCG14与细胞分裂素生物合成高度共表达，是根至芽的主要细胞分裂素转运体.编号为Cot-AD 34390的基因属于ABCG家族基因，推测Cot-AD 34390可能参与了棉花在盐胁迫时期的ABA转运功能. Cot-AD 55296是一类可编码AMSH-like泛素硫酯酶，KATSIARIMPA等^[17]研究表明：拟南芥AMSH3基因调控泛素化膜蛋白的降解过程，泛素翻译后修饰在调节内吞降解过程中起关键作用；拟南芥AMSH3是一种去泛素酶，amsh3缺失突变体导致幼苗的致死率大幅提高，并在多种细胞内转运途径中都表现出缺陷.

编号为Cot-AD 59144的基因是细胞色素P450家族基因，是一种末端加氧酶，可参与生物体内甾醇类激素合成等过程. MAGWANGA等^[22]研究表明：棉花CYP450基因在响应植株干旱和盐胁迫等逆境环境刺激过程中，发挥着重要的调控作用.同样，在拟南和水稻中，细胞色素P450蛋白家族可通过赤霉素失活在植物发育过程中发挥着靶特异性调控作用^[23]. PtCYP714A3在形成层-韧皮部组织中广泛表达，且在杨树中受盐胁迫和渗透胁迫影响较大. PtCYP714A3异位表达导致水稻半矮化表型，分蘖促进，籽粒变小^[24].过量表达PtCYP714A3基因的转基因株系较野生型植株，GA积累水平更低，部分GA生物合成基因的表达受到明显抑制，表明PtCYP714A3在水稻幼苗对盐毒反应中起着重要的作用^[23].

本研究对盐胁迫下棉花根系进行转录组分析，通过与正常培养条件下比较，初步筛选出的差异表达基因，通过qRT-PCR方法，对部分差异显著的基因表达进一步验证.后续实验会克隆这些经过验证的完整cDNA序列，通过构建过量表达载体，将基因转化到拟南芥和棉花植株中进行基因功能的遗传学分析，其结果为进一步克隆棉花关键抗盐基因和认识其功能提供了依据，也为揭示棉花的抗盐机制奠定了基础.