-
马铃薯(Solanum tuberosum L.)是茄科、茄属双子叶植物，是中国乃至全球次于小麦、水稻和玉米的第4大重要粮食作物，是支撑中国农村经济发展的重要组成部分[1].马铃薯高产、耐瘠薄、分布广泛、生育期短、适应性强，富含碳水化合物、蛋白质、膳食纤维、矿物质及维生素等营养物质，而且粮用、菜用、饲用兼优[2]. 2015年起国家农业部开始大力推广马铃薯主粮化工程[3].怀玉山高山马铃薯又称麻籽洋芋，2013年4月核准为国家地理标志农产品，主要种植于江西省上饶市怀玉山山区；其块茎呈卵圆形或长圆形，皮糙色黄肉白，淀粉含量10.3%~14.4%，蛋白质含量2.56%~2.60%，干物质含量16.8%~18.0%，维生素C 0.131~0.133 mg/g，未检测出还原糖，富含膳食纤维；有和中养胃、健脾利湿、降糖降脂、美容养颜、防治胃癌等功效[4].农业生产实践表明，怀玉山高山马铃薯不能在平原地区种植，否则将会出现形态变大、品质下降，且这些生产农艺形状移回高海拔地区也不可逆[5]，究其原因，尚无相关报道.
转录组学已逐渐成为基因功能研究的重要手段之一[6].近年来，基于高通量方法的转录组测序正成为植物生理生态研究领域的有力工具[7].转录组测序是对某一物种的mRNA进行高通量测序，而转录组测序的结果反映了特定条件和特定时间点的所有转录本序列信息考[8]，进而得到基因功能注释、蛋白质编码区序列、基因的表达量、代谢途径等大量信息，为进一步研究提供基础数据和重要参考[9].该方法不仅能够用于有参考基因组序列的物种研究，也能用于无参考基因组序列的物种研究，已经被广泛应用于各个研究领域[10-11].
本研究基于边合成边测序(Sequencing by synthesis，SBS)技术，利用Illumina Hiseq 4 000高通量测序平台对怀玉山高山马铃薯正常种及其平原种植退化种块茎进行转录组测序分析、基因表达量分析和差异表达分析、新基因预测与功能注释分析，以及差异表达基因功能注释分析，对怀玉山高山马铃薯平原种植出现形态变大、品质下降的机理进行研究，以期解析怀玉山高山马铃薯平原种植退化种的内在机理，并为其适应高海拔生境相关基因的挖掘与应用提供参考.
HTML
-
怀玉山高山马铃薯正常种(麻籽洋芋，编号CK)及其平原种植退化种(编号BZ)采后未萌芽的块茎，由江西省上饶市红日农业开发有限公司提供.
-
参照说明书用Trizol试剂提取怀玉山高山马铃薯正常种(CK)及其平原种植退化种(BZ)采后未萌芽块茎总RNA.所提RNA的质量由Nanodrop 2000和琼脂糖凝胶电泳检测，Agilent 2100测定RIN值.若RNA总量5 μg，浓度≥200 ng/μL，OD260/280介于1.8和2.2之间即为可用.
-
RNA测序送由上海美吉生物科技股份有限公司完成.利用带有Oligo(dT)的磁珠，将3′端带有PolyA尾结构的mRNA从总RNA中分离出来.然后利用TruSeqTM RNA sample prep kit试剂盒构建cDNA文库.利用Illumina hiseq 4000进行测序和分析(3次技术重复).使用Seq Prep工具(https://github.com/jstjohn/seqprep)和Sickle工具(https://github.com/najoshi/sjckle)将原始的FASTQ文件中包含一些带接头、低质量、N率较高序列及长度过短序列的reads进行过滤、修剪或删除，得到高质量的测序数据(clean reads)，以保证后续信息分析的质量.
-
使用Trinity软件(https://github.com/trinityrnaseq/trinityrnaseq/wiki)对clean reads进行de novo拼接，形成一个不能在两端扩展的Unigene序列.利用Trinity软件组装得到转录本，并与参考蛋白质数据库NR(Non-redundant protein，https://www.ncbi.nlm.nih.gov)进行BLAST比对，保留最优比对结果，挑选一条最长的代表序列作为Unigene.使用软件Trans Decode对组装结果进行CDS预测.再将序列注释到Swiss-Protein(https://www.expasy.ch/sprot)，KEGG(KyotoEncyclopedia of Genes and Genomes，https://www.genome.jp/kegg)和GO(GeneOntology Consortium，https://geneontology.org)数据库中获得基因功能信息.对所有单基因进行GO功能分类，在宏观层面上了解基因功能的分布.基因的表达水平由转录本的丰度表示，利用RPKM(每百万个reads映射到外显子区的每千个碱基上的reads数)对Unigene的read count值进行标准化. P值来识别差异表达的基因，在多次检验和分析中，用错误发现率来识别P值的阈值.根据错误发现率(false discovery rate，FDR≤0.001)和表达量倍数变化(|foldchange|>1.5)的阈值筛选出2个样本间差异表达的基因(Degs).对显著性差异表达基因进行KEGG Pathway显著性富集分析.
-
根据转录组试验结果，采用Primer BLAST在线工具对相关差异表达基因(GSC0003DMG400011811，PGSC0003DMG402005881，PGSC0003DMG400021727，PGSC0003DMG400000926，PGSC0003DMG400022022，PGSC0003DMG400013684，PGSC0003DMG400027631，PGSC0003DMG400004595，PGSC0003DMG400006448，PGSC0003DMG400003563)进行引物设计(表 1)，每个样品以GAPDH(GenBank登录号为NM_001288415.1)为内参，3个技术重复，采用1.2.1的方法提取2种材料的总RNA，按照PrimeScriptTM RT reagent Kit with gDNA Eraser(Perfect Real Time)试剂盒说明书合成cDNA的第1条链，在荧光定量PCR仪中进行扩增反应，每个反应重复3次.数据采用SPSS 19.0版本进行统计分析.
1.1. 试验材料
1.2. 试验方法
1.2.1. RNA提取
1.2.2. Illumina Hiseq测序
1.2.3. 数据分析
1.2.4. 相关差异表达基因的qRT-PCR
-
经过测序质量控制，共得到14.83 Gb Clean Data.其中BZ组Clean Data中pair-end Reads总数为25 461 192，总碱基数为7 689 279 984，G和C 2种碱基占总碱基的百分比为44.46%，质量值大于或等于30的碱基所占的百分比为93.76%；CK组Clean Data中pair-end Reads总数为23 643 517，总碱基数为7 140 342 134，G和C 2种碱基占总碱基的百分比为47.31%，质量值大于或等于30的碱基所占的百分比为93.53%. 2组样品Q30碱基百分比均不小于93.53%.结果表明，此次转录组测序数据质量较好，满足转录组分析的基本要求.
-
比对效率指Mapped Reads占Clean Reads的百分比，是转录组数据利用率的最直接体现.经统计，BZ组Total Reads(即Clean Reads数目，按单端计)为50 922 384，比对到参考基因组上的Reads数(Mapped Reads)为37 005 516，在Clean Reads中占的百分比为72.67%；比对到参考基因组唯一位置的Reads数(Uniq Mapped Reads)为35 067 236，在Clean Reads中占的百分比为68.86%；比对到参考基因组多处位置的Reads数(Multiple Mapped Reads)为1 938 280，在Clean Reads中占的百分比为3.81%；比对到参考基因组正链的Reads数为17 764 983，在Clean Reads中占的百分比为34.89%；比对到参考基因组负链的Reads数为18 197 367，在Clean Reads中占的百分比为35.74%. CK组Total Reads(即Clean Reads数目，按单端计)为47 287 034，比对到参考基因组上的Reads数(Mapped Reads)为17 179 940，在Clean Reads中占的百分比为36.33%；比对到参考基因组唯一位置的Reads数(Uniq Mapped Reads)为16 150 402，在Clean Reads中占的百分比为34.15%；比对到参考基因组多处位置的Reads数(Multiple Mapped Reads)为1 029 538，在Clean Reads中占的百分比为2.18%；比对到参考基因组正链的Reads数为8 093 153，在Clean Reads中占的百分比为17.11%；比对到参考基因组负链的Reads数为8 399 367，在Clean Reads中占的百分比为17.76%.比对效率除了受测序质量影响外，还与指定的参考基因组组装的优劣、参考基因组与测序样品的生物学分类关系远近(亚种)有关.通过本次试验的比对效率可知，所选参考基因组组装能满足后续信息分析的需求.
-
从图 1可知，CK组和BZ组的密度分布范围在0~0.5. CK组和BZ组表达量FPKM的对数值基本为0~2.
-
基因差异表达火山图(图 2)中的每一个点表示一个基因，横坐标表示某一个基因在2个样品中表达量差异倍数的对数值；纵坐标表示基因表达量变化的统计学显著性的负对数值.横坐标绝对值越大，说明表达量在2个样品间的表达量倍数差异越大；纵坐标值越大，表明差异表达越显著，筛选得到的差异表达基因越可靠.差异表达基因MA图(图 2)中每一个点代表一个基因.横坐标为A值：log2(FPKM)，即2个样品中表达量均值的对数值；纵坐标为M值：log2(FC)，即2个样品间基因表达量差异倍数的对数值，用于衡量表达量差异的大小.
对筛选出的差异表达基因做层次聚类分析，将具有相同或相似表达模式的基因进行聚类，差异表达基因聚类结果如图 3.经分析，BZ组和CK组差异表达基因数目为3 604个；其中，上调基因的数目为1 331个，下调基因的数目2 273个. BZ组和CK组差异最显著的10个基因分别为PGSC0003DMG400000028，PGSC0003DMG400000046，PGSC0003DMG400000048，PGSC0003DMG400000049，PGSC0003DMG400000064，PGSC0003DMG400000066，PGSC0003DMG400000069，PGSC0003DMG400000070，PGSC0003DMG400000075，PGSC0003DMG400000079.其中相比CK组，BZ组双功能L-3-氰基丙氨酸合酶/半胱氨酸合酶2、类CCR4-非转录复合物亚单位1等2个基因上调，LOC101261174、类甘油-3-磷酸2-O-酰基转移酶4、类WRKY转录因子22、乙烯反应性晚期胚胎发生丰富蛋白、植物细胞壁蛋白SlTFR88、非特征蛋白LOC101265683、非特征蛋白LOC101245497、类BAG家族分子伴侣调节因子等8个基因下调.
-
图 4展示的是在差异表达基因背景和全部基因背景下GO各二级功能的基因富集情况，体现2个背景下各二级功能的地位，具有明显比例差异的二级功能说明差异表达基因与全部基因的富集趋势不同.通过富集，镉离子响应、根毛伸长、对缺水的反应、真菌防御反应、盐胁迫反应、脱落酸反应、单元部分、质膜、细胞内膜结合细胞器、叶绿体、高尔基体、光系统Ⅱ、光系统Ⅰ、蛋白质结合、电子载体活性、血红素结合、叶绿素结合、谷胱甘肽结合等GO term富集显著(表 2).这些功能类别可能在怀玉山高山马铃薯的平原种植出现形态变大、品质下降中起着重要作用.
-
差异表达基因COG分类统计结果如图 5所示，BZ组和CK组差异表达基因General function prediction only(一般功能预测)，Transcription(转录)，Signal transduction mechanisms(信号转导机制)，Replication，recombination and repair(复制、重组和修复)和Secondary metabolites biosynthesis，transport and catabolism(次级代谢产物的生物合成、运输和分解代谢)等5大功能的基因数较多.
-
对BZ组和CK组差异表达基因KEGG的注释结果按照KEGG中通路类型进行分类，分类图如图 6(左)所示.结果表明，BZ组和CK组差异表达基因注释到植物激素信号转导(Plant hormone signal transduction)、苯丙酸生物合成(Phenylpropanoid biosynthesis)、植物病原相互作用(Plant-pathogen interaction)、碳代谢(Carbon metabolism)、氨基酸生物合成(Biosynthesis of amino acids)的基因个数及其个数占被注释上的基因总数的比例分别为48(7.92%)，47(7.76%)，45(7.43%)，32(5.38%)和31(5.12%).
差异表达基因KEGG通路富集分析结果见图 6(右)，呈现了显著性Q值最小的前20个通路.其中，差异表达基因富集q value≤0.25的KEGG通路主要有苯丙氨酸代谢(Phenylalanine metabolism)、植物病原相互作用(Plant-pathogen interaction)、苯丙素生物合成(Phenylpropanoid biosynthesis)、谷胱甘肽代谢(Glutathione metabolism)、半乳糖代谢(Galactose metabolism)、光合作用(Photosynthesis)、氮代谢(Nitrogen metabolism).这些途径可能对怀玉山高山马铃薯的平原种植出现形态变大、品质下降起着重要调节作用.
-
从RNA-seq数据中挑选10个在KEGG通路注释数量较多的差异表达基因PGSC0003DMG400011811(ferredoxin-NADP，reductase铁氧还蛋白-NADP还原酶)，PGSC0003DMG402005881(ferredoxin-3，铁氧还蛋白-3)，PGSC0003DMG400021727(photosystem Ⅱ oxygen-evolving complex protein 3，光系统Ⅱ析氧复合蛋白3)，PGSC0003DMG400000926(xygen-evolving enhancer protein 2，产氧增强蛋白2)，PGSC0003DMG400022022(photosystem I reaction center subunit IV B，光系统Ⅰ反应中心亚基Ⅳ B)，PGSC0003DMG400013684(trans-cinnamate 4-monooxygenase-like，反式肉桂酸4-单加氧酶样)，PGSC0003DMG400027631(gibberellin 2-oxidase，赤霉素2-氧化酶)，PGSC0003DMG400004595(ABC transporter B family member 11-like，类ABC转椅子B家族成员11)，PGSC0003DMG400006448(anthocyanin O-methyl transferase，花青素O-甲基转移酶)，PGSC0003DMG400003563(flavanone 3 beta-hydroxylase，黄烷酮3β羟化酶)，并利用qRT-PCR技术对其表达量进行验证(表 3)，结果显示，10个基因表达量的变化趋势均与RNA-seq结果一致，说明本次转录组测序数据可信度高. qRT-PCR也进一步验证了与怀玉山高山马铃薯平原平原种植退化种有关的PGSC0003DMG400000926，PGSC0003DMG400003563，PGSC0003DMG400006448，PGSC0003DMG400022022，PGSC0003DMG400021727基因在怀玉山高山马铃薯平原种植退化种中下调表达，而PGSC0003DMG400004595，PGSC0003DMG400011811，PGSC0003DMG400027631，PGSC0003DMG402005881，PGSC0003DMG400013684基因在怀玉山高山马铃薯正常种中上调表达.
2.1. 测序数据产出统计
2.2. 转录组数据与参考基因组序列的比对效率
2.3. 基因表达量分析
2.4. 基因差异表达分析
2.5. 差异表达基因GO分类
2.6. 差异表达基因COG分类
2.7. 差异表达基因KEGG注释
2.8. 相关差异表达基因的qRT-PCR分析验证
-
怀玉山高山马铃薯种植于平原地区即低海拔区域容易产生平原种植退化种，其主要原因是怀玉山高山马铃薯的高海拔生境发生了改变.怀玉山高山马铃薯在高山地区的长期栽培已经产生一整套适于高海拔生境的生理机制，将其移种到低海拔区域，必然打破其原有的生理适应机制，产生一整套适于低海拔生境的生理机制.高山地区对植物生理机制产生影响的因素主要有低温、紫外线、水分等，这些因素的胁迫对高山地区植物的生理机制将会产生一个综合性影响.张飞等[10]对水分胁迫下高粱差异基因KEGG分析结果显示，苯丙烷类生物合成、苯丙氨酸代谢、类黄酮生物合成3个途径中差异基因表达较多.白子彧等[12]的研究也指出水分胁迫与苯丙氨酸代谢密切相关.此外，有学者指出水分胁迫将造成光合系统发生变化，影响光合物质生产[13].王晓娇等[14]发现“克新1号”马铃薯在萌芽出苗期遇到水分胁迫时，显著富集的有谷胱甘肽代谢、植物激素信号转导、苯丙素生物合成、苯丙氨酸代谢和植物-病原菌相互作用等途径.本研究与上述研究结果基本一致.在本研究中，怀玉山高山马铃薯平原种植退化种差异表达基因富集的KEGG通路主要有苯丙氨酸代谢、植物-病原菌相互作用、苯丙素生物合成、谷胱甘肽代谢、半乳糖代谢、光合作用、氮代谢.这些途径对怀玉山高山马铃薯的平原种植退化起着重要调节作用.苯丙氨酸代谢途径是植物重要的次生代谢途径，主要包括苯丙烷代谢和异黄酮合成代谢[15-16].植物的抗病与病原菌的致病是生物界一种十分复杂的相互作用系统.当植物组织受到病原真菌侵染后，寄主体内的防御体系会发生一系列的改变，具有明显的时空性，且受到环境条件、代谢产物与基因表达的调控[17].怀玉山高山马铃薯移种到低海拔区域，受低温、紫外线、水分等胁迫减少，病原菌侵染加强，植物与病原菌相互作用可能也是怀玉山高山马铃薯平原种植退化的主要原因.苯丙素是具有C6-C3芳香核骨架的一类天然化合物的统称，包括简单苯丙素类、木质素和木脂素类、香豆素类及黄酮类等许多天然芳香族化合物[18-19].其中香豆素类具有多种生理活性[20]，如抗菌、抗病毒等；木质素作为植物细胞壁的主要组成部分，与纤维素、半纤维素、果胶等物质交叉联接，其合成能够对植物细胞起重要的保护作用，增强植物细胞抵抗病害和逆境的能力[21].由于苯丙素生物合成和苯丙氨酸代谢均属于次生代谢产物生物合成，这些次生代谢物质对植物抵抗逆境胁迫具有重要作用，说明怀玉山高山马铃薯移种到低海拔区域对马铃薯的生长发育造成一定程度的伤害.谷胱甘肽是一种含有巯基的三肽化合物，在植物体内是一种尤为重要的抗氧化物质[22]，在减缓细胞氧化损伤上具有重要的作用，可降低怀玉山高山马铃薯移种到低海拔区域对马铃薯造成的毒害作用[23].维生素C含量的多少可以作为衡量园艺产品品质的重要指标之一，在植物中维生素C的合成中，L-半乳糖内酯脱氢酶作为关键酶可直接氧化半乳糖内酯形成L-抗坏血酸(VC)[24].抗坏血酸是植物生长发育所必需的微量营养素，在生命过程中发挥重要生理功能，特别是在抵抗逆境胁迫反应中具有重要作用[25].研究表明，低温胁迫下秋茄叶片上述抗氧化酶活性和抗氧化剂含量表现为抗坏血酸/氧化型抗坏血酸、还原型谷胱甘肽/氧化型谷胱甘肽比值增加，从而抑制活性氧(ROS)大量产生，减少膜系统伤害[26].怀玉山高山马铃薯移种到低海拔区域，降低了抗坏血酸-还原型谷胱甘肽循环中抗氧化酶活性，降低清除ROS能力，增强了叶片PSII光抑制，从而降低了植株光合作用，这可能也是怀玉山高山马铃薯移种到低海拔区域出现形态变大、品质下降的原因之一.氮在植物生长中起着重要的作用，氮素作为植物蛋白质、核酸、酶及叶绿素的组成成分，又参与多种内源激素及其前体的合成[27].研究表明，在逆境条件下，氮的有效性对于植物的光合能力和植株生长至关重要，植株中硝态氮含量呈降低的趋势，但铵态氮含量呈增加的趋势，氮代谢相关酶活性均有不同程度的降低，植物的硝态氮更多的分配到了根部，以增强植物的抗逆性[28].怀玉山高山马铃薯移种到低海拔区域出现形态变大、品质下降等现象，可能是降低了氮代谢相关酶活性，促使碳流由光合碳代谢转向氮代谢，破坏了怀玉山高山马铃薯碳氮代谢的正常进行，减弱了怀玉山高山马铃薯抵御逆境胁迫的能力.
本试验研究了怀玉山高山马铃薯正常种及其平原种植退化种块茎涉及苯丙氨酸代谢、植物-病原菌相互作用、苯丙素生物合成、谷胱甘肽代谢、半乳糖代谢、光合作用、氮代谢等相关基因的差异表达，仅是关于怀玉山高山马铃薯正常种及其平原种植退化种块茎转录组的变化的研究；下一步将对怀玉山高山马铃薯正常种及其平原种植退化种块茎试管苗转录组变化进行系统研究，以获得对怀玉山高山马铃薯适应高海拔生境的全面认识.
DownLoad: