比较转录组分析显示二倍体D-基因组棉中的进化分歧和共用网络的冷和盐应激反应

背景

棉花野生种对非生物胁迫具有良好的抗逆性。二倍体d -基因组棉表现出丰富的表型多样性，可能是生产最大纺织天然纤维的异源四倍体棉的供体。

结果

通过mapping RNA-seq Illumina reads，在所有样本中共表达了41,053个基因G. Thurberi.(D₁），g . klotzschianum(D_{3 k}）, g . raimondii就(D₅),g . trilobum(D₈)到参考基因组。冷胁迫下差异表达基因(DEGs)数量显著高于盐胁迫下。4种植物在盐胁迫和冷胁迫下有34.1%的DEGs重叠。值得注意的是，通过基因共表达分析，我们发现了一个潜在的共享网络(包括16个基因的冷盐反应)。通过De novo组装，在4个二倍体物种中共鉴定出47180 ~ 55548个独特基因。阳性选择基因(PSGs)分别为163、344、330和161thurberi，g . klotzschianum，G. Raimondii.和g . trilobum通过进化分析，在冷或盐胁迫下，分别在相应物种中均为9.5-17％的四种物种。更重要的是，大多数PSG都是与刺激响应相关的最佳关系。g . klotzschianum显示这在冷、盐胁迫下具有良好的耐受性。有趣的是，我们发现一种ralf样蛋白编码基因不仅是PSGsg . klotzschianum但也属于潜在的共享网络。

结论

我们的研究提供了新的证据，即自然选择进化的基因表达变化是与进化期间与环境适应有关的形态变异的必要司机。此外，在寒冷和盐压力下存在共享调节网络，例如CA²⁺信号转导和氧化还原机制。本研究建立了四种二倍体d基因组棉花基因表达改变的转录组选择机制，为多非生物抗性棉花育种提供了可用的基因资源。

棉 （Gossypium)提供最天然的纤维，以制造纺织品[1，是全球重要的天然纤维作物。目前,g .分子和g .取得（同种异体一体化棉花）被广泛种植，并通过长期历史驯化。遗传证据表明，在大约1-2百万年前的二倍体A和D-基因组种杂交事件中形成了同种异体一体化棉花[2]．二倍体d基因组棉花包含至少13个品种。在这些物种中,G. Thurberi.，克劳茨夏努，雷蒙地弓形虫，和g . trilobum，分布在美洲的四个不同纬度地区[3.，形态性状差异明显。通过系统发育方法揭示了这四个物种的分子进化过程和系统发育[3.]．例如,G. Thurberi.和g . trilobum尽管自然分布的纬度有明显的不同，但显示出系统发育的密切关系。此外,G. Thurberi.D基因组的一个成员能忍受低温[4]．

目前，研究人员承认这一点g . arboreum和G. Raimondii.分别为优质异源四倍体棉花A和D基因组的假定供体[5]．A-基因组二倍体棉包含两种物种，g . arboreum和g . herbaceum，分布在非洲南部和亚洲。d -基因组二倍体棉花约有13种，主要分布在墨西哥，范围扩展到秘鲁、加拉帕戈斯群岛和亚利桑那州南部[6]．D-和a -基因组棉花基因组和形态特征的差异是由于地理隔离和划分造成的。约1-2百万年前，a -基因组二倍体与d -基因组二倍体棉花杂交，随后发生多倍体事件，异源四倍体棉花出现[7]．这就要求A和D基因组必须在物理上接近[8但是跨越太平洋联系是不可思议的。因此，现代异源四倍体棉花的起源和进化存在许多谜团，尽管目前有很多假说或理论来解释AD棉基因组的起源和进化。

低温和盐胁迫是严重限制世界棉花生产的重要环境因素[9]．植物对非生物胁迫的适应依赖于各种分子网络的激活和触发，这些分子网络涉及胁迫感知、信号转导、代谢物和特定胁迫相关基因的诱导[10]．应激诱导因子可以同时或依次发生，导致渗透胁迫、水分缺乏、离子失衡、过氧化损伤，最终导致生长抑制[11]．钙作为第二信使在非生物应激反应中起着重要作用[12，13]．盐和冷胁迫均能提高植物胞质游离钙的浓度。研究表明，OSCA1(降低高渗透压诱导的钙增加1)是渗透胁迫的假定传感器，涉及冷和盐胁迫反应[14，15，16]．MAPK(丝裂原活化蛋白激酶)级联，由钙等第二信使刺激，参与非生物应激信号转导[16]．此外，SnRK2(蔗糖非发酵相关蛋白激酶2)家族也参与了盐、渗透和干旱胁迫下的信号转导。植物中MAPK和SnRK2可被冷和盐快速激活[15，16]．

越来越多的棉花基因组测序和重测序、mRNA测序和表型评估[15[为研究棉花潜在的生物机制提供重要资源。通常用于通过基因表达改变构建调节模型的比较转录组分析[17]．系统转录组分析为研究驯化过程中的基因进化和表达变化提供了一种新的策略[14]．例如，在野生番茄亲属中鉴定了数百种已经进化的新蛋白序列或响应于自然选择改变的表达水平的候选基因，通过番茄素分析，表明人为和自然选择对番茄及其野生亲属的转录组表达改变在进化和驯化中发挥着重要作用[18]．加权相关网络分析(weighted correlation network analysis, WGCNA)是基因共表达网络(gene co-expression network, GCN)分析的一个R包，可作为数据挖掘工具或基因筛选(排名)方法，发现高度相关基因的聚类(模块)[19]．广泛用于发现生物医学科学中的集线基因[20.]．

在我们的研究中，为了揭示低温和盐胁迫下的遗传和表达多样性，我们对4个二倍体d -基因组物种进行了转录组测序，包括G. Thurberi.，G. klotzschianum, G. raimondii和g . trilobum．4个物种的系统发育关系与之前的研究一致[3.]．通过比较转录组学鉴定了6个物种特异性谱和4个物种特异性模块。在各基因的表达中，冷胁迫比盐胁迫下有更多的基因表达差异。基因进化分析在不同的基因组或亚基因组中发现了数百个PSGs(野生物种:G. Thurberi.，克劳茨夏努，雷蒙地弓形虫，和g . trilobum；栽培物种:g . arboreum，g .分子一个_T.和D_T.，g .取得一个_T.和D_T.）.我们还发现一个模块与WGCNA的盐和冷应力负相关。g . klotzschianum在这种物种中显示出在冷和盐胁迫下的抗性和盐胁迫下的良好抗性，并且在该物种中鉴定了171个常见的含量。总之，基因表达变化是与进化期间与环境适应有关的形态变化的必要驱动因素，并且存在参与冷和盐应激反应的共用网络，例如信号转导和氧化还原过程。通过组合蛋白质编码序列和基因表达多样性，我们的作品在演化期间对表达分歧，保护和对环境适应的反应提供了深入了解。

二倍体D基因组物种的表型多样性和RNA-seq数据

在转录组水平上观察野生二倍体d -基因组棉花的进化分化和保守G. Thurberi.(D₁），g . klotzschianum(D_{3 k}）, g . raimondii就(D₅),g . trilobum(D₈)，具有丰富的形态特征，以供进一步分析。这些物种在花的颜色和叶的形状上表现出巨大的表型差异，尽管它们在系统发育上彼此最接近(图。1a).此外，4个品种表现出不同的耐冷和耐盐性:g . klotzschianum和G. Thurberi.表现出良好的抗寒性和耐盐性G. Raimondii.(无花果。1b）。奇怪的是，g . trilobum表示最低的阻力，尽管它与G. Thurberi.．为了综合评价，我们对d基因组二倍体棉花进行了rna测序。四个物种,G. Thurberi.，克劳茨夏努，雷蒙地弓形虫，和g . trilobum，为方便起见，简称为GD1、GD3、GD5、GD8。在2个胁迫处理(冷T12和盐胁迫S12)后的不同时间间隔(分别为C0、C6和C12三叶期后0 h、6 h和12 h)，每个物种采集10个叶片样品，每个样品重复2次。共获得40个数据集，原始数据273.01 GB。在不同的数据集中，得到的GC含量率为44.41-46.64%。经过质量控制，共获得273.01 Gb未使用适配器、ploy-N和低质量读取的干净数据。每个数据集至少有1710万次干净读取，超过基准Q30的89.03%。之后，rna测序的干净序列被固定在参考基因组上，G. Raimondii.．mapping clean reads的范围为81.37 ~ 91.37%，其中只有81.37 ~ 91.37%的clean reads被特异地映射到参考基因组(表S)1）.根据作图结果，每个文库共鉴定出7374个SNPs(单核苷酸多态性)(表S2）.定量分析了41,053个基因的表达水平，其中3548个是新的转录基因。除GD1S12样品外，同一样品区域的R1和R2库相关性较高(R²> 0.7）;建议其他样本的数据集是可靠的(图。1c). GD1S12R1和GD1S12R2相关性较低的原因可能是发育和环境变异对基因表达的影响，虽然我们试图尽量降低这一影响。

四种二倍体D基因组物种的遗传差异

转录序列多态性分析有助于理解四种d -基因组二倍体棉花在转录组水平上的进化多样性和保守性。利用RNA-Seq数据进行序列多态性发现，在不同文库中鉴定了7374-404,737个SNPs(表S2）.超过59%的SNPs具有过渡类型。通过转录组测序获得的数据集显示，基因导入区染色体上有较多的SNPs分布。值得注意的是，在GD1、GD3和GD8的基因间区域发现了约9%的SNPs，这代表了其他d -基因组物种的一个新的基因区域。GD5图书馆(G. Raimondii.)的snp数量比其他三个物种少。一般来说，GD5样本的转录组应该与参考基因组具有相同的基因组。虽然GD5样本和参考基因组的转录组非常相似(大多数基因每Kb < 1 SNP)，但GD5和参考基因组之间存在一些SNPs，表明相同供体材料的植物之间存在转录组差异。在四种二倍体d基因组棉花中，预计有3,2651个基因通过单核苷酸多态性获得高功能效应。Go本体论(Go ontology, Go)富集分析表明，大量基因高度参与了氧化还原、应激反应、蛋白质修饰和蛋白质磷酸化等多种过程。研究表明，非生物胁迫在驱动4种棉花转录多样性方面发挥了重要作用(表S3.）.相反，也值得注意的是，8402个基因没有snp功能效应，包括许多管家基因，如(微管蛋白、核糖体蛋白、糖酵解酶编码基因)。共有8402个基因参与了包括光合电子传递链在内的多种生物途径。基于转录组数据，我们使用Neighbor-joining方法构建了4个物种的系统发育树(图。2a).并通过主成分聚类观察四种植物之间的关系(图4)。2b）。适度数量的SNP在四种物种中将它们分开（每kB大部分基因的5个SNP）。相同物种的数据集几乎彼此重叠。G. Thurberi.显示出与g . trilobum已经被证明。

转录组de novo集会

基于使用Trinity软件的RNA-SEQ数据进行转录组De Novo组件。范围从47,180到55,548个独特的基因，在四个二倍体物种中鉴定出来。大多数unigenes（> 40%)长度范围为200 ~ 500 bp。为了开始我们的进化分析，我们在7个物种中鉴定了严格同源的unigenes，包括3个栽培的(A₂-Genome：g . arboreum；广告₁-Genome：g .分子；广告₂-Genome：g .取得)和四种野生物种。为了进行严格的分析，将异源四倍体物种的基因组分离为A和D亚基因组，以鉴定同源转录本。共对47,119对同源基因对进行了表征，其中5312对单拷贝转录本对采用最大似然法构建了9个基因组的系统发育树。使用这些unigenes进行的系统发育分析显示了一个非常相似的拓扑结构(图。2c)基于转录组SNP数据的树(图。2a)，表明四个d基因组二倍体物种的系统发育牢固。与之前的研究相同，A亚基因组g .分子和g .取得从起源g . arboreum，在系统发育分析中形成单系。g .分子- d subgenome,g .取得d亚基因组与GD5形成单系。这些结果证实了异源四倍体AD的D亚基因组₁和广告₂基因组有相同的捐赠者，来自G. Raimondii.．

估计进化速率并识别积极选择基因(PSG)

采用非同义(dN)和同义(dS)估计各物种的进化速率和正选择基因(PSGs)。dN/dS值为> 1，这表示积极选择，而等于或大于1(> 1)表示净化或中性选择。因此，我们评估了9个基因组中已鉴定的同源单拷贝基因的dN、dS和dN/dS。基于上述构建的系统发育树，使用自由比率模型(模型= 1)评估各同源单基因对在不同分支中的dN/dS，该模型允许每个分支有单独的dN/dS比率。结果表明，野生二倍体棉花(GD1、GD3、GD5、GD5_ref和GD8)的dN/dS比栽培棉花(GAD1、GAD2和GA2)高。在自然选择压力下，野生物种在适应性选择的基础上表现出快速的进化比率(图2)。2d).通过PAML软件分别鉴定GD1、GD3、GD5、GD8的psg分别为163、344、330、161(表S)4）.二倍体棉花品种的PSGs富集了与蛋白质修饰、蛋白质泛素化、RNA加工和ncRNA加工相关的GO术语，涉及环境适应(表S5）.更多，我们还发现Ga2，GAD1_A，GAD1_D，GAD2_A和GAD2_D中的180,77,51,103和70psgs，这些PSG富集MRNA代谢过程（表S6）.

基因表达的分化与保守

我们在至少一个样本中检测了总共41,053个转录本的表达水平。40个数据集的整体表达水平分布相似。利用24组未经胁迫处理的叶片数据，对4种植物的基因表达差异和保守性进行了检测。以GD5C0R1和GD5C0R2作为对照组，鉴定不同的表达转录本。4个物种的差异基因表达模式被描述(图。3.一种）。基于基因表达模式，共聚集了29,512次的八种型材使用K-Means聚类方法。具有大量参数的配置文件，显示四种类似的表达模式（轮廓1,2和3）。曲线1和2显示出四种物种叶片的表达相当大的差异。在型材1中，DEG在下调表达水平随时间，并富集在光合生物中的碳代谢和碳固定。相反，概况2的含量上调并富含植物昼夜节律和核糖体生物发生在真核生物中。在轮廓3中，DEGs是显着富集的不饱和脂肪酸，DNA复制和光合作用 - 天线蛋白方法的生物合成（图。3.b).这些途径在叶片的生长发育过程中是必不可少的，而这些途径的转录本在棉花d -基因组物种分化过程中是保守的。Profile 4-9被观察到一个相当大的物种特异性表达模式。在基因型4和5中，3662和1048个DEGs显示了gd5特异性的表达模式。此外，在profile 6中，1516个DEGs显示了gd1特异性的表达模式。图7和图8显示了GD3-GD5和GD1-GD8的特殊表达模式。考虑到GD1和GD8的关系密切，GD1和GD8包含更保守的生物学途径是完全合理的。

基因共表达网络的进化保守与分化

虽然基因表达模式的分析为基因表达的分化和保守提供了思路，但由高连接基因构建的共表达基因网络更与生长发育的重要生物学过程和非生物胁迫反应的复杂调控过程相关。为掌握基因共表达网络的保守性和分歧性，采用加权基因共表达网络分析(WGCNA)方法，对40组叶片数据进行共表达网络构建。为了提高WGCNA分析的准确性，每千碱基外显子模型中每百万reads的片段小于一个的基因(FPKM <1）被移除。结合连接强度(软阈值功率:5，R²> 0.90）在12,110个基因中（图S1)，构建了四种物种间共表达网络拓扑的全局视图。同样显示模块的基因显示了更高的拓扑重叠(图S2）.最后，总共使用33个模块来研究基因共表达网络的进化保守和分化，基因共表达网络被定义为高度相互关联的基因簇(图)。4a).在这些模块中，基因数量在31个(暗橙色)到3128个(绿松石色)之间，各模块之间的相关系数较高。在每个模块上评估每个基因的连通性。高度连结基因(特制> 0.95)被鉴定为枢纽基因。最终共检测到425个hub基因(在模块内范围从1到293)。

在每个模块中，所有基因的表达水平通过热图显示出来，并通过特征基因值(模块表达谱的第一个主成分)进行总结。为确定的重要模块定义了七个样本条件(图。4b). Share.T.S组用于识别冷和盐胁迫反应的保守共享网络。此外，与冷和盐胁迫反应相关的模块分别被冷组和盐组识别。GD1、GD3、GD5和GD8组用于基因组特异性模块的鉴定。通过特征基因和样本条件之间的关联分析(7组:Share.T。通过Pearson相关系数分析，发现29个模块具有显著的基因组特异性和/或非生物胁迫调节的共表达模式(ANOVA，P < 0.05）.高度共表达基因的四个主要模块的相关性最强(Pearson相关性)r > 0.9）分别有四个基因组。最大的模块（绿松石），含有3128个高度连接的GD5基因，响应刺激和免疫系统过程而富集（表S.7）.第二个模块(Blue, 1193个基因)显示了GD3特异性。另外两个模块(黄色有420个基因，紫色有442个基因)包含的高连接基因较少，分别与GD1和GD8相关。有趣的是，四个主要模块的氧化石墨烯富集分析结果相似。近一半的氧化石墨烯富集基因与每个模块对刺激的反应和免疫系统过程有关，揭示了非生物和生物胁迫在这四个物种之间的转录变异中发挥了主要作用。我们注意到7个模块的特征与相应的组相关，其中5个模块是基因组特异性的。此外，大多数模块与四个基因组组的相关性最强。这些结果表明，转录变异主要与基因组差异相关。我们观察到一些模块在不同的组之间有重叠。重叠模块有助于理解两个或多个基因组之间的相似生物学特性。 For example, the dark green module was significantly showed a positive correlation with GD1 and GD3, but a negative correlation with GD5 and GD8, and enriched sesquiterpenoid and triterpenoid biosynthesis and flavonoid biosynthesis which were known related with abiotic and biotic stress resistance.

特别值得关注的是，由于四个二倍体物种呈现出不同的冷和盐胁迫性，与冷和盐胁迫相关的共表达网络。一个ssociation analyses between co-expression modules and abiotic stress (Cold and salt stress) revealed that 17 modules correlated with abiotic stress (Shar.T.S: 7 modules; Cold: 8 modules; Salt: 8 modules), thus representing suites of interconnected genes underlying the biological process of the abiotic stress response (Table S7）.值得注意的是，在四种物种中，GD3和Salt/Cold组有更多重叠模块。考虑到GD3比其他3种植物表现出更好的耐冷和耐盐性，表明GD3在非生物胁迫适应中进化出更完整的机制。相反，GD8几乎没有一个模块与非生物胁迫组存在重叠，可以预见GD8在冷和盐胁迫下表现出显著的敏感性。

有趣的是，SkyBlue3模块与份额显着与份额显着负相关。 - 冷和盐组。中的一个中心基因，快速碱化因子样(RALF -like)。来自拟南芥的同源基因可能参与调控植物的逆境反应和或适应、生长和发育。该模块中有15个基因与ralf样蛋白编码基因相互连接。除了5个非特征基因外，15个相互关联的基因中有10个与GO术语(对刺激的反应)有关，包括Gorai.005G234900，Gorai.007G094200，gorai.n023400.，和Gorai.011G238800．拟南芥的同源基因已被注释为4个重要的转录因子:MYB44、TCP9、TCP12和GATA8，主要参与非生物胁迫反应。另外，两个木葡聚糖内转葡萄糖基酶/水解酶蛋白22 (XTH22)编码基因，Gorai.003G052400和gorai.009g006400.被观察到参与碳水化合物运输和代谢。同源物的Gorai.005G094300在拟南芥中，EXORDIUM蛋白编码基因是叶片细胞扩张所必需的，并可能参与协调油菜素内酯(BR)对环境或发育信号的响应的信号传导过程(图。5）.15个基因在冷盐胁迫处理后大部分表达量下降。

盐碱应力条件下寒冷盐胁迫下患者胁迫的特征及对关键基因的验证

错误发现率(FDR)≤0.001和日志₂利率≥2(8组:处理组/对照组)对4种植物进行DEGs鉴定。8组处理间deg的数量从459个到3372个不等(图S3.）.共检出基因7515个，占检出基因总数的20.04%。有趣的是，在盐胁迫组中，4个物种的DEGs数量总是小于冷胁迫组(图2)。6a)，表明冷胁迫反应比盐胁迫反应有差异，且同一种非生物胁迫中有更多的独特的DEGs(6405个，占所有DEGs的85.2%)再次证明了这一点。但是，有1109个基因(16.1%的冷基因;在两种不同的非生物胁迫中发现了63.3%的Salt DEGs)，表明在冷胁迫和盐胁迫反应中存在潜在的共享调控途径。与预期一致，1109个基因在参与多种非生物胁迫反应的氧化还原途径中富集。

因此，我们对GD3中的DEGs进行了研究，发现GD3中的DEGs在冷、盐胁迫下表现出了较高的耐受性，最终鉴定出了2759个GD3中的冷、冷组DEGs，其中冷组DEGs 2300个，盐组DEGs 630个。在这些基因中，我们发现171个在盐和冷胁迫下共享的deg，其中80个下调，75个上调(图1)。6c).有趣的是，共有102个基因参与了对刺激的响应，其中32个基因与对冷或盐胁迫的响应有关，包括13个转录因子编码基因(NAC, ERF, MYB, G2, HD-ZIP)被认为与对非生物胁迫的响应有关。因为这些基因在拟南芥中有一些同源基因与非生物胁迫反应有关。例如,Gorai.002G073700编码拟南芥NaC72的同源物，它与抗旱响应粘合顺式脱水应激早期反应因子1 [21.]．同时,这两个Gorai.001G239000和Gorai.006G017400在拟南芥中编码参与植物防御反应的myb样蛋白同源物[22.]．显著的30℃和3℃的盐温均为PSGs。此外，skyblue3模块中有3个冷DEGs和12个盐DEGs。虽然deg更少(冷deg: 6.1%;在所有的盐类基因中，有部分基因与PSGs或skyblue3模块重叠(图3)。6b).这些结果再次表明适应性进化驱动转录多样性，并且参与冷和盐胁迫反应的共享调控网络。我们还发现Gorai.006G147500skyblue3模块的基因被积极选择并与中心类ralf蛋白编码基因。这进一步证实了前面分析中发现的潜在调控网络。通过反相遗传学对关键基因的验证表明，vigs植物对盐碱的适应能力受到了很大的损害。在盐胁迫下，vigs植株表现出萎蔫和叶缘逐渐干燥的现象。7A)，结果与先前的发现一致，即盐胁迫确实会对叶子产生负面影响，进而导致地上生物量的显著减少[23.). .此外，它们都表现出显著的低生物量(图。7b)，此外，各组织离子含量的测定，基因敲除的植株Na浓度显著降低⁺/ K⁺这表明两种离子之间的平衡受到盐碱毒性的影响(图。7c)。

在之前的研究中，一些研究集中在d基因组二倍体棉花的关系上，如d基因组二倍体棉花的关系g . trilobum接近G. Thurberi，相反g . klotzschianum和G. Raimondii.［3.]．通过将我们的转录组数据与参考基因组进行比对，检测到7374 - 4047.37亿个SNPs (G. Raimondii.变得更)。然后，利用这些SNPs构建系统发育树来研究4个物种之间的关系。4个物种之间的关系与前人研究一致，进一步验证了基于RNA-Seq de novo assembly的可行基因组序列分析。异源四倍体棉花有7种。基于异源四倍体和二倍体物种间染色体上45S和5S rDNA的相似分布，一些研究者认为不同的异源四倍体棉花物种有不同的供体物种，如G. Thurberi.，G. klotzschianum, G. raimondii和g . trilobum［24.]．但是,G. Raimondii.显示大多数同步块与D-亚基g .分子和g .取得通过比较基因组分析，提出G. Raimondii.可能的供体物种是g .分子和g .取得．因此，一些研究人员认为同种异体一体化棉可能是单巧妙的。在我们的研究中，5312个单拷贝基因用于构建系统发生树并估计进化率。一般形成的A-基因组和亚亚胍，以及单层形成的D-基因组和D-亚胍基。我们的结果支持G. Raimondii.是D亚基因组的供体种g .分子和g .取得．这与异源四倍体棉花的单系进化理论相一致[7]．

阳性选择在植物对生物和非生物胁迫的演化和适应中起重要作用，因为基因表达和通过阳性选择的调节变化已经假定为适应性进化率的关键决定因素[25.]．我们的结果表现出野生棉具有更高的DN / DS比（基因组加速的进化）而不是栽培棉花，表明野生棉可能会发生适应性的进化，使他们能够应对其极度广泛的可怕条件和环境。我们的研究发现在进化期间发现了一百个基因，这些基因也富集了与非生物或生物应激反应相关的术语。并确认阳性选择驱动野生物种的环境自适应演化。K-Means野生物种的聚类分析发现，富含物种 - 特殊型材与非生物应激反应相关的术语，表明四种野生物种的表达分析还证明了非生物胁迫在四种物种之间推动转录分类的转化多样性[26.]．

我们观察了四种植物的形态差异。四种不同的生态系统导致了形态上的差异。以往的遗传证据表明，基因表达的改变是在进化过程中驱动表型多样性的关键[27.]．我们的结果同意了这一结论，鉴定了五种特异性基因表达簇（简介4-8），大多数基因是与环境适应相关的富集术语的富集。同样，大多数野生物种与环境适应有关。超过9.6％的PSG与相应的DEG集合的基因重叠。因此，我们推测进化选择也可以驱使基因表达改变以适应环境。

形态学分析发现g . klotzschianum耐冷和盐胁迫。g . klotzschianum，这自然范围是加拉帕戈斯群岛，它可以自我调整以适应环境。因此，我们研究了g . klotzschianum在积极的选择。psg的数量很低g . klotzschianum(344个，占全部基因组的0.9%)。然而，33个PSGs在盐胁迫和冷胁迫下表现出差异表达，表明自然选择导致的基因表达改变可能在提高环境适应能力方面发挥了重要作用。

结合了前面的上下文的分析g . klotzschianum在寒冷和盐胁迫下表现出惊人的耐受性，我们推测了在自适应演进过程中形成了涉及冷和盐应激反应的共享网络g . klotzschianum．33个(~ 10.6%)PSGs在冷、盐胁迫下有不同表达。特别是在WGCNA中，我们发现了一个与盐、寒冷和共享t.s.条件负相关的模块。74%的基因富集了与刺激相关的GO项，36.3%的基因富集了与信号传导相关的GO项。这个结果进一步证实了我们的推测。同时，我们发现应激信号过程是盐和冷应激反应的一个初步共享调控网络。在WGCNA中发现了skyblue3中的一个中心基因。并鉴定出15个与中心基因显著相关的基因。有趣的是,Gorai.006G147500在15个高连接基因中发现。skyblue3模块中的16个基因在低温和盐胁迫下具有潜在的调控网络。研究者假定植物细胞必须能够感知各种环境信号。28.]，并且在以前的研究中也发现了一些假定的传感器，如OSCA(高渗透诱导钙增加减少1)[26.， G蛋白[26.]，冷却[29.]．盐和冷应激可能导致植物中的细胞溶质游离钙浓度的增加[30.]．的相同器官GORAI.010G168400.在拟南芥中,编码AtOSCA3.1蛋白，在高渗门控非选择性阳离子通道中发挥作用²⁺离子(31.]．和GORAI.010G168400.，在盐和冷胁迫下表达量高于对照组，是介导棉花冷和盐胁迫的潜在传感器，钙离子作为响应非生物胁迫的第二信使[32.]．Gorai.009G294400它编码cbl相互作用蛋白激酶18²⁺信号转导，在冷胁迫和盐胁迫下也有不同的表达。核心应激信号通路涉及与酵母SNF1和哺乳动物AMPK相关的蛋白激酶[33.]．SNF1/ ampk相关激酶介导各种非生物胁迫信号。在171种温度的盐和冷组中，gorai.005g081200.和Gorai.002G103300其中在激酶结构域中与SNF1 / AMPK分享同源性。

氧化还原过程也参与盐和冷应激反应，通过保持活性氧物种的内稳态[34.]．我们发现氧化还原酶活性相关基因如Gorai.004G093200，Gorai.013G025400，Gorai.013G059500和Gorai.013G176900在冷盐胁迫下表达上调，与氧化还原过程相关。PSGs和天空蓝模块基因也产生了许多重要的候选基因，如Gorai.004G227900和Gorai.002G226600，这与保持活性氧的内稳态有关。对这些基因在冷盐胁迫下的研究有助于我们了解DEGs如何调节活性氧的稳态和信号转导，以及它们在冷盐胁迫下的作用，从而帮助我们提高棉花产量。

对四种二倍体d基因组棉花进行转录组比较分析，揭示其序列和基因表达的差异。棉花野生品种和栽培品种的基因进化分析鉴定了参与棉花驯化和进化的正向选择基因，并估算了其进化速率。在本研究中，我们发现进化选择可以驱动基因表达改变以适应环境，基因表达变化是四个d -基因组物种分化过程中的主要进化事件。表达模式分析发现，6个谱图显示了不同的物种特异性特征，这些谱图的基因参与了对刺激的响应。因此，基因表达变异是进化过程中与环境适应相关的形态变异的重要驱动因素。与盐胁迫下相比，冷胁迫下发现了更多的基因表达变化，表明冷胁迫导致了更多的基因表达变化。G。显示klotzschianum表现出较好的抗寒性和耐盐性。与其他3种植物相比，在植物中检测到更多的PSGsg . klotzschianum， 9.6%的PSGs在冷和盐胁迫下有差异表达。在g . klotzschianum结果表明，27.1%盐胁迫下的DEGs与冷胁迫下的DEGs重叠，且skyblue3模块与冷、盐和shared . t.s条件组显著负相关。因此，在冷盐胁迫反应中存在信号转导和氧化还原过程等共享网络。在此基础上，我们推测出了一套参与冷盐胁迫响应的候选基因，为棉花抗多非生物育种提供了遗传资源。

植物生长和样品收集

本研究利用d -基因组棉的4个二倍体野生种G. Thurberi.(D₁），g . klotzschianum(D_{3 k}）, g . raimondii就(D₅),g . trilobum(D₈）.二倍体D-Genome棉有13-14种，它们最初从墨西哥西南到亚利桑那州分发，秘鲁和加拉帕戈斯群岛的额外分散物种分布。这里，G. Thurberi.，G. klotzschianum, G. raimondii和g . trilobum由中国农业科学院棉花研究所(ICR-CAAS)管理，现保存于海南省三亚市国家野生棉花苗圃。野生棉花种子来源于中国农业科学院棉花研究所管理的野生棉花苗圃。种子萌发温度为28℃，光照16 h /暗8 h，光照强度为150 μmol m^{- 2}年代^{- 1}在含水量为15%的砂土中。发芽3天后，将合适的植株放在土壤中盆栽，并置于相同条件下的生长室中。将幼苗分为2个单叶期和1个心形叶期(时间点:0 h)。第一组置于28°C, 16 h光照/8 h暗循环，光照强度150 μmol m的正常条件下^{- 2}年代^{- 1}和含水量15%),第二组与300毫米的氯化钠溶液(28°C, 16 h光/ 8 h暗周期和150μ摩尔强度m - 2 s - 1和15%含水量)和第三组在低温(4°C, 16 h光/ 8 h黑暗周期和150μ摩尔的光强度^{- 2}年代^{- 1}15%的水分含量)。分别在0 h、6 h和12 h采集幼苗叶片。第2组和第3组研究冷盐胁迫对叶片生长的调控作用，并于12 h采集叶片。这个实验重复了两次。

RNA提取，文库构建，RNA测序

样品中的RNA通过使用Trlzol试剂（Life Technologies，California，USA）在三个生物学重复作为教学手册中提取。通过使用Agilent 2100 BioAnalyzer（Agilent Technologies，Inc.，Santa Clara，Ca，Ca，USA）来验证RNA质量和浓度。纳米二玻璃2000分光光度计进一步检查RNA质量，并在琼脂糖凝胶上电泳。在260/280比例为1.8-280的比例为1.8-230，260/230比率≥2.0的RNA用于进一步分析。然后通过使用M-MLV转录酶试剂盒转录RNA，从大连，中国大连获得Takara Biotechnology。此外，通过Nebnext Poly（A）mRNA磁隔离模块（NEB，E7490）检索的MRNA。施用转录的CDNA以使用肌炎（NEB，E7530）和illumina（NEB，E7500）的NEBNEXT Ultra RNA文库预备试剂盒构建cDNA文库。简而言之，将富集的mRNA分离成具有约200nt的RNA，其用于合成第一链cDNA，然后是第二个cDNA。适配器连接程序用于分析双链CDNA，并通过Agencourt Ampure XP珠（Beckman Coulter，Inc。）检索符合片段，并通过聚合酶链反应（PCR）扩增来富集。通过使用Illumina Hiseq TM 2500测序平台在流动细胞上进行构建的cDNA文库的测序。 Beijing Biomarker Technologies (http://www.biomarker.com.cn)提供实验程序和商业执行。

RNA-seq数据处理

通过内部Perl脚本过程处理fastq格式的原始数据，在该脚本过程中，通过使用ploy-N、适配器和低质量读取来检索干净数据。此外，对检索到的clean reads进行gc含量、Q20、Q30和序列重复水平的评价。所有的下游分析都是基于高质量的干净数据。然后利用Hisat2工具将干净的reads映射到参考基因组序列，并对完全符合的reads进行进一步分析。

基因表达分析

对不同基因表达谱质量的评估，表达水平由每千碱基的转录本每百万片段(FPKM)分析[35.]．使用DESeq2评估两个条件/组的差异表达分析，其中P根据Benjamini和Hochberg的方法来调整错误发现的基因P-值小于0.01则称为差异表达基因(DEGs) [36.]．

SNP呼叫

Picard - 工具v1.41和samtools v0.1.18用于对重复的读取进行排序和删除并整合每个样本的BAM对齐。使用SAMTOOLS软件进行单一核苷酸多态性（SNP）呼叫。使用GATK规定的标准滤波方法和其他变量过滤原始的VCFF。标准和变量被编程为ClusterWindowsize：10;Mq0 > = 4和(mq0 / (1.0* dp)) > 0.1；Quent <10；战< 30.0或QD < 5.0或HRun > 5，仅使用距离为> 5的snp。

转录组组装

所有库中的左文件和右文件被合并成一个大的左文件和右文件。fq文件。转录组组装建立在左侧。fq和正确的。使用Trinity设置参数，min_kmer_cov默认设置为2，所有其他设置为默认[37.]．表达分析、SNP调用和转录组组装使用BMKCloud (www.biocloud.net）在这项研究中。

同源鉴定、系统发育分析、进化速率估计和正选择基因(PSG)鉴定

两个异源四倍体棉花(g .分子和g .取得)和两种二倍体棉花(g . arboreum和G. Raimondii.)从科通根数据库(https://www.cottongen.org）.异源四倍体物种的基因组包含两个亚基因组(A-和d -亚基因组)。为了保证同源鉴定的准确性和可靠性，分离了两个异源四倍体物种的亚基因组。结合我们的转录组组装结果(G. Thurberi.，g . klotzschianum和g . trilobum)，共9个基因组进行系统发育分析，a -基因组分支为外群(g . arboreum，g .取得subgenome和g .分子subgenome)。Orthofinder [38.]用于将基因与S设置钻石和设置为默认的所有其他参数进行群集。单复制的正交基因对具有来自每个基因组和亚基组的一种拷贝的单拷贝用于系统发育分析。单拷贝基因家族的蛋白质序列由肌肉排列（V3.8.1551）[39.]．使用Gblocks获得排列良好的蛋白质序列[38.，40]．基于4D位点从编码序列(CDS)比对，利用RA x ML构建系统发生树[41.]．此外，利用来自4个物种(G. Thurberi.，g . klotzschianum，g . raimondii就和g . trilobum）使用BMK云的SNPS（www.biocloud.net）.使用PAML中的codeml程序估计了九个侏儒的每个谱系的进化速率[42.[带有免费配给模式的包装(模型= 1）.我们从Codeml的结果中获取dN、dS和dN/dS。我们用dS = 0过滤基因。如果观测dN / dS > 1为PSGs。

基因聚类和可视化

为了确定每个物种中基因随时间的表达模式，使用K-means聚类可视化基因表达模式(log)₂-转换的FPKM值)使用BMKCloud (www.biocloud.net）.

加权基因共表达网络分析（WGCNA）

WGCNA [19]软件包鉴定基因共表达网络，研究性状相关模块。按照Xu等人的描述使用WGCNA [43.]．参数为相关矩阵的软阈值，强调的是基因之间的相关性[44.]．选取值5作为软阈值功率的决定因素和自由拓扑分析的评价尺度。动态树木切割算法和模块与特征的关联被Xu等人评价。[43.]．

关键基因功能分析的验证

病毒诱导候选基因的基因沉默（Vigs），RALF和FLA，g .分子玛丽-加兰特85号比赛采用规定的程序进行[43.，45.]．首先，我们构建了TRV: RALF、TRV: FLA和TRV: PDS的VIGS矢量，并将其引入农杆菌肿瘤术应变GV4104。将农杆菌培养物浸透到10日龄土壤育苗的两个膨大子叶中g .分子比赛Marie Galante 85（Mar85）。棉幼苗在温室中种植，温度设定在26°C，16小时光/ 8小时暗循环。每个构建体都接种了总共24个幼苗。14天后的土壤杆菌接种后，然后将Vigs-植物和非Vigs植物暴露于盐碱处理6天。然后收集样品以用于形态学和生理分析。

统计分析

本实验产生的所有数据均采用R (v3.5.0)软件进行分析。

数据沉积

本研究工作产生并报告的RNA-seq数据已保存在国家生物技术信息中心(NCBI)，登录号为PRJNA554555 (https://www.ncbi.nlm.nih.gov/sra/PRJNA554555)。

在目前的研究中生成和/或分析的数据集可在NCBI知识库中获得，https://www.ncbi.nlm.nih.gov/bioproject/PRJNA554555。

WGCNA:

加权基因共表达网络分析

NCBI:

国家生物技术信息中心

PSG：

积极选择基因

度:

差异表达基因

SNPS：

单核苷酸多态性

我们非常感谢整个研究团队在本次研究工作中给予我们的支持。

该研究工作得到了来自中国国家自然科学基金的资金支持;授予31621005,31530053和2016YFD0100203的国家重点研究和开发计划。新疆人才培养计划，培训青年优秀S＆T天赋2017Q050项目。资金编号为研究计划提供了财务支持，但没有涉及工作设计，数据收集，分析和稿件的准备。

作者指出

1. 徐彦超、马旺加、金鼎沙等人对这部作品贡献良多。

从属关系

1. 中国农业科学院棉花研究所，棉花生物学国家重点实验室，安阳455000

   徐彦超，马旺加，金鼎沙，蔡晓燕，侯玉清，王坤波，刘芳，周忠礼

2. 华中农业大学植物科学技术学院，湖北武汉40070

   验钞徐

3. 生物、物理、数学和精算科学学院(SBPMAS)，主校区，Jaramogi Oginga Odinga科技大学(JOOUST)，邮政信箱210-40601,Bondo，肯尼亚

   Richard Odongo Magwanga和Stephen Gaya Agong

4. 新疆省农业科学院经济作物研究所，新疆乌鲁木齐

   郑Juyun

5. 郑州大学农业科学学院，河南郑州450001

   方刘

贡献

YX, FL, KW, ZZ和ROM设计并构思了这项研究。YX、YCX、DSJ和ROM进行实验。YX, XC, ZZ, YH, ROM和FL贡献了材料/分析工具。YX和ROM进行生物信息学分析;FL、KW、ZZ、YC指导研究工作;YX、ROM、FL、ZJ、SGA、KW、ZZ对数据进行解释和修改;FL、KW、ROM和DJ对手稿进行了修改。所有作者都同意投稿。

相应的作者

对应到kunbeb王或方刘或Zhongli周．

伦理批准和同意参与

没有伦理或同意参与这项研究。这项研究工作是根据中国棉花研究所(CRI)的广泛授权进行的，CRI是负责发展、开展研究和批准中国所有棉花育种工作的国有研究机构。

同意出版物

不适用。

相互竞争的利益

作者声明没有任何形式的利益竞争。

出版商的注意

施普林格《自然》杂志对已出版的地图和机构附属机构的管辖权要求保持中立。

开放访问本文根据创意公约归因于4.0国际许可证，这允许在任何中或格式中使用，共享，适应，分发和复制，只要您向原始作者和来源提供适当的信贷，提供了一个链接到Creative Commons许可证，并指出是否进行了更改。除非信用额度另有说明，否则本文中的图像或其他第三方材料包含在文章的创造性公共许可证中，除非信用额度另有说明。如果物品不包含在物品的创造性的公共许可证中，法定规定不允许您的预期用途或超过允许使用，您需要直接从版权所有者获得许可。要查看本许可证的副本，请访问http://creativecommons.org/licenses/by/4.0/．创作共用及公共领域专用豁免书(http://creativecommons.org/publicdomain/zero/1.0/)适用于本文提供的数据，除非在数据的信贷额度中另有说明。

再版和权限