棉花|二倍体和异源多倍体棉花染色质结构和重复基因表达的进化动态

最新推荐文章于 2026-08-19 17:21:09 发布

原创最新推荐文章于 2026-08-19 17:21:09 发布 · 468 阅读

9 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

收录于

文献学习

Evolutionary Dynamics of Chromatin Structure and Duplicate Gene Expression in Diploid and Allopolyploid Cotton

二倍体和异源多倍体棉花染色质结构和重复基因表达的进化动态

材料	基因组类型	倍性	说明
G. hirsutum cv. Acala Maxxa (AD₁)	AADD	异源四倍体 (2n=4X=52)	天然异源多倍体（陆地棉）
G. arboreum accession A2-101 (A₂)	AA	二倍体 (2n=2X=26)	A基因组二倍体祖先
G. raimondii (D₅)	DD	二倍体 (2n=2X=26)	D基因组二倍体祖先
F₁ 杂交种 (A₂ × D₅)	AD	二倍体 (2n=2X=26)	人工合成的A-D基因组杂交种（无染色体加倍）

############################################################################################################
1.DNS-seq
DNS-seq利用了DNase I酶的生物学特性。DNase I是一种核酸内切酶，它对DNA的切割具有选择性：

开放区域：在染色质上形成比较松散、没有核小体结合的区域（即调控因子结合的活性区），DNase I极易切割，所谓的“DNase I超敏感位点”（DHS）。

核小体包裹区域：在DNA缠绕在核小体上的区域，由于蛋白质的物理保护，DNase I的切割效率较低，但仍能以极小的概率在核小体内部进行切割。

############################################################################################################

摘要

多倍体化是植物物种形成和适应的重要机制，然而对于重复基因调控的机械机理理解仍不清晰。染色质结构的动态变化被认为主导了基因的调控过程。本研究利用 DNS-seq 技术，以异源四倍体陆地棉（Gossypium hirsutum，AADD, 2n = 4X = 52）及其两个二倍体亲本（AA 或 DD 基因组）和人工合成的二倍体杂种（AD）为对象，鉴定了全基因组范围内的核小体组织和染色质可及性。研究发现，在二倍体中，较大的 A 基因组表现出更宽的平均核小体间距，而这种基因组间的间距差异在异源多倍体中消失，但在杂种中依然存在。此外，全基因组多倍化还表现出启动子区域可及性的全基因组增长，以及亚基因组间顺式作用元件（cis-regulatory motifs）的同步化。研究推断顺式作用控制对染色质动态具有显著影响，并通过启动子区域转座元件（TE）的移除得到了证实。通过将可及性与基因表达模式相结合，我们发现杂交阶段与随后的全基因组加倍阶段具有不同的调控效应，包括同源偏向表达（HEB）和表达水平优势（ELD）的细微建立过程。组蛋白基因表达与核小体组织通过染色质可及性进行协同。本研究展示了追踪高分辨率染色质结构动态的能力，揭示了它们在多倍体顺式调控景观演变和重复基因表达中的作用，为深入理解亚基因组不对称性与优势现象的调控联系提供了新见解。

介绍

多倍体化是真核生物中一种普遍存在的生物学现象，在生物组织的各个层级都具有重要意义（Fox 等，2020）。全基因组加倍（WGD）在蕨类和开花植物中尤为盛行（Jiao 等，2011；Ruprecht 等，2017；一千种植物转录组计划，2019），由多倍体化导致的全基因组加倍对植物的生理、生态和进化产生了深远影响（Stebbins，1940；Levin，1983；Ramsey 和 Schemske，2002；Leitch 和 Leitch，2008；Wendel，2015；Soltis 和 Soltis，2016；Van de Peer 等，2017, 2021；Wendel 等，2018；Heslop-Harrison 等，2022）。

多倍体化可能与生态范围的扩大有关（Arrigo 等，2016；Coughlan 等，2017；Baniaga 等，2020；Wang 等，2021；Parshuram 等，2022；Zhao 等，2022；Elliott 等，2023；Mata 等，2023），并能增强植物对生物和非生物胁迫的耐受性（综述见 Van de Peer 等，2021），引起生理变化（Mishra，1997；Sugiyama，2005；Otto，2007；Knight 和 Beaulieu，2008；Coate 等，2012；Orr-Weaver，2015），以及改变生物合成途径（Combes 等，2022）。

这些变化可能赋予植物具有经济或生态价值的重要性状（Heslop-Harrison 等，2022）。不出所料，许多重要的作物物种都是相对年轻的多倍体（Olsen 和 Wendel，2013；Renny-Byfield 和 Wendel，2014；Zhang 等，2019；Heslop-Harrison 等，2022）。

由多倍体化引起的全基因组含量的增加通常与核型特征的变化相关，例如细胞大小、细胞核体积以及细胞周期时间（Wendel 等，2018；Doyle 和 Coate，2019）。这些基因组面积的变化还可能改变表观遗传动力学、表达、蛋白质组以及分子网络。

一个被广泛证实的效应是，在异源多倍体化过程中的基因组融合与加倍过程中，调整组发生深刻的重构（Grover 等，2012；Hu 和 Wendel，2019；Visger 等，2019；Shan 等，2020；Giraud这种全基因组范围的重构涵盖了多种现象，包括：在基因组水平上同源基因的不等量表达（称为“同源偏向表达”，HEB）（Flagel 等，2008；Grover 等，2012）；在基因组水平上的表达优势（基因组优势）（Schnable）等，2011）；不同组织或条件下同源偏向的不一致性（表达亚功能化和新功能化）（Adams 等，2003），甚至在单细胞水平上也是如此（Zhang 等，2023）；明显的重复基因表达反式控制（表达水平优势，ELD）（Rapp 等，2009；Grover 等）等，2012；Yoo 和 Wendel，2014；Yoo 等人。 2014）；以及改变的共表达基因网络（Gallagher 等，2016；Hu 等，2016）。

虽然这些研究揭示了多倍体转录组的进化动力学，但这些现象的机械强度（底层）仍然模糊不清，这一限制机制限制了我们对重复基因表达进化的理解，其次也限制了我们对伴随多倍体化产生的进化创新起源的认识。

染色质结构的研究已成为连接基因组进化与转录组进化之间鸿沟的一个新兴领域，为深入了解基因表达调控的动态变化提供了见解。染色质结构的景观图谱反映了精细调控基因表达的多个复杂调控层级（Talbert 等，2019；Ahmad 等，2022）。

核小体是染色质的基本结构单位，由包裹在核心组蛋白八聚体上的 147 个碱基对的 DNA 组成（Luger 等，1997）。核小体促进了基因组 DNA 向染色质的压缩，在控制 DNA 可及性方面发挥着至关重要的作用，进而影响基因转录、DNA 复制、修复和重组等过程（Kornberg，1974；Andrews 和 Luger，2011）。在转录激活过程中，核小体可以发生移动以暴露或遮蔽顺式调控 DNA 位点，或者在启动子区域发生瞬时失稳（被称为“脆弱”核小体）（Zlatanova 等，2008；Mieczkowski 等，2016；Klemm 等，2019）。

因此，核小体充当了染色质可及性的调节器，这种可及性本质上体现了组蛋白和 DNA 的无数种表观遗传修饰，这些修饰共同控制着基因的表达（Schmitz 等，2011；Jordan 和 Schmitz，2016；Kawakatsu 等，2016；Niederhuth 等，2016；Hofmeister 等，2017；Jackson，2017；Song 等，2017；Springer 和 Schmitz，2017；Giles 和 Taberlay，2019；Klein 和 Hainer，2020）。理解决定核小体特性的因素及其对染色质可及性和基因活性的影响，是当前生物学面临的一项核心挑战。

在过去的十年中，高通量技术已被广泛评估植物中，用于在全基因组测量上较弱的核小体激发度（核小体占据）和染色质可及性图谱（Tsompana and Buck 2014；刘等2015;张等2015;冯等2017;张和江2018;巴尔迪等2020;加利等2020；约旦等2020；赵等2020;巴比尔等2021）。

这些方法包括微球菌核酸酶扩增（MNase-seq）、DNase I超敏感位点扩增（DNase-seq）以及转座酶可及染色质鉴定酶（ATAC-seq），它们均基于对染色质的物理可及性。通过片段化（fragmentation）、转座切割加标签（tagmentation）或去除（elimination）等手段，利用核酸酶的模式切割来区分DNA的断裂与受核小体保护或受转录因子（TF）保护的区域。

自20世纪70年代以来，DNase I超敏感位点（DHSs）一直被认为是真核生物基因组中激活调节区域的标志（Weintraub 和 Groudine 1976;高通量 DHS 绘图已为多种植物物种的顺式调控元件（CREs）和调控因子结合位点（TFBSs）提供了全基因组层面的意见（Zhang 等 2012; Wu 等 1979a, 1979b）。江2015;沙利文等2015;邱等2016;赵等2018;韩等2020、2022）。

ATAC-seq作为DNase-seq一种更高效的替代方案，能够实现快速且低输入量的染色质可及性分析（Lu 等 2017），甚至可以达到单细胞水平（Dorrity 等 2021）。这些技术及其变体，为理解植物物种中的顺式景观和基因调控网络提供了高度视角（Lu 等 2019;利玛窦等2019;雷诺索等2022）。

另一方面，MNase-seq在历史上常用于分析核小体占据度（核小体占据度），并在拟南芥（Chodavarapu 等，2010；Li 等，2014；Liu 等，2015）和水稻（Wu 等，2014；Zhang）该技术的最新应用利用了两种微球菌核酸酶（MNase）消化条件——紫色（轻）和重度（重），能够同时提供核小体定位数据和染色质可及性/脊柱图谱（Vera 等，2014；Rodgers-Melnick 等，2016）。

也就是说，通过对核小体扭转度进行变异酶探针（DNS）分析，可以鉴定出不同水平的染色质可及性；这种方法最初是在玉米中基于DNA微阵列建立的（Vera等，2014），附加采用了高磁场幅度进行全基因组分析（Rodgers-Melnick等，2016）。与通过DNase-seq和ATAC-seq鉴定的DHS类似，来自差异性MNase-seq（DNS-seq）的MNase敏感足迹（MSFs）在基因的5’和3’边界富集，并且与基因表达、DNA低甲基化、保存非编码序列以及已知的转录因子结合位点（TFBSs）呈正相关。

在玉米中，MNase超敏感区域仅占基因组的1%，但它们与约40%表型性状变异的基因型变异相关，这与编码区（约48%）的影响相当（Rodgers-Melnick等，2016）。另外，通过MNase分析获得的顺式光照景观已与组织拓扑拓扑和环境响应建立了联系，突显了它们在模型模型变异中的作用（Pass）等，2017；Parvathaneni 等，2020）。一种基于休眠MNase消化产生的相关DNA短片段的分析方法MOA-seq最近被开发出来，用于异构小颗粒图谱，从而绘制出开放染色质区域内顺式元件上可能的艾滋病调控情况（Savadel 等，2021；Liang 等，2022）。

总之，利用MNase作为染色质结构所表现出的特性，已被证明在表征染色质景观、核小体定位、核小体稳定性以及功能性顺式调节元件（CREs）的鉴定方面具有极高的信息价值。

棉属（Gossypium）是研究多倍体进化基因组学的公认模型。目前已知有 50 个多个物种（Wendel and Grover 2015;胡等2021; Viot 和 Wendel 2023），并且不断有新的棉花物种多样性（Stewart 等 2015；系统发育分析（Wendel and Cronn 2003; Gallagher 等 2017）。温德尔等2010; Chen 等 2017）和基因组序列数据（Huang 等 2021）表明，该属起源于约 500 万至 1000 万年前（mya）。

异源多倍体棉花（AD基因组）起源于更新世（更新世），当时A基因组的祖先跨洋传播至新大陆，并与当地的D基因组二倍体发生了杂交。异源多倍体由此产生了现在的七个物种，其中包括具有重要商业价值的移植棉（G. hirsutum）和海岛棉（G. barbadense），两者均在过去价值7000全年被驯化（Wendel and Grover 2015）。

与D基因组祖先亲缘关系最近的最初物种是雷蒙德氏棉（G. raimondii），而统一的A基因组物种——亚洲棉（G. arboreum）和草棉（G. herbaceum），在模拟杂交过程中的母本（种子亲本）方面都是优秀的模型（Wendel 等 1989）。棉属这种清晰的进化历史已成为研究异源多倍体化的卓越模型。

以往的研究强调了棉属（Gossypium）中重复基因表达进化的多个方面，包括“同源偏向表达”（HEB）——即两个同源基因中其中一个的表达水平明显另一个；以及“表达水平优势”（ELD）——这是一种神秘的现象，即两个同源基因的总表达量在统计上与仅其中一个亲本的表达水平异无（Rapp 等，2009；Grover 等，2012；Hu等，2013, 2014, 2015；Yo 等，2013；Gallagher 等，2020）。

中部还对异源多倍体棉花表达的顺式（顺式）和反式（反式调控）进行了研究，发现在约5,000至8,000年的驯化过程中，反式调节变异发生了优先积累（Bao等，2019）。棉花中的这些质及其他调控变化存在与染色景观的各个方面相关联或因果关系，包括DNA甲基化（宋）等，2017）、组蛋白修饰（Zheng 等，2016）、染色质可及性（Han 等，2022）以及三维组基因拓扑结构（Wang 等，2018）。但截至目前，染色质结构的分子机制及其重复对基因表达的影响很大程度上仍是未知的。

在本研究中，我们应用 DNS-seq 技术，对异源四倍体陆地棉（G. hirsutum）及其模型二倍体祖先，以及模拟 100 万至 200 万年前的杂交自然过程的人工合成二倍体 F1杂种，进行了全基因组染色质可及性和核小体组织的全面分析。除了重复表征基因组与叠加而来的染色质结构动态变化外，我们还考察了基因组的表达模式，以揭示异源多倍体棉花中染色质重综上所述，我们的详细研究展示了染色质结构和顺式调控景观的演化动力学，强调了基因融合组与加倍如何这些景观，并阐明了它们在重复基因表达演化中的调控改变作用。

结果

1.植物材料

本研究使用了四种棉属（Gossypium）基因型，包括：一种天然异源四倍体（AD 基因组）——本土棉（G. hirsutum）品种 Acala Maxxa (AD1)；以及其模型二倍体祖先（A 基因组和 D 基因组）——即亚洲棉（G. arboreum）材料 A2-101 (A2) 和雷蒙氏棉德（G. raimondii ））（D5）。

A和D这两个二倍体基因组群最后一次拥有共同祖先是在500万至1000万年前（Wendel 和 Albert 1992），此后它们发生了杂交，导致基因组大小（GS）外侧了2倍。因此，研究中还包括了相应的种间二倍体F1杂种（A2 × D5），旨在研究两个杂交的基因组融合后的直接后果（在不存在基因组加倍以及多倍化后的演化时间的情况下）。

多重基因型种植4到5株，生长于美国爱荷华州立大学（Ames, IA, USA）Bessey Hall温室，环境控制为短日照条件（10小时天线，下午5点至次日早晨7点为黑暗；夜间/日间温度为22/28°C）。在开花枝条上采集集群组织，采集时间为下午5点，随后立即在液氮中速冻并储存于-80°C。

2.数据处理

在使用CutAdapt（Martin 2011）进行质量过滤并引引序列后，将来自不同棉属物种的双端读取段（paired-endreads）比对到从CottonGen数据库（Yu等2014）下载的相应参考基因组上。这些基因组包括：棉花品种TM1 UTX v2.1（Chen等2020）、亚洲棉品种SXY1 WHU更新版v1.0（Huang 等 2020）以及雷蒙德氏棉 JGI v2.0（Paterson 等 2012）。F1 杂志的读段则比对到亚洲棉和雷蒙德氏棉的合并参考基因组上。

使用 Bowtie2 (v2.5.1) 进行比对，参数设置为“no-mixed”（禁止单端比对）、“no-discordant”（禁止不一致比对）、“no-unal”（不记录未记录未比对读段）和“dovetail”（允许交尾比对）。并保留比对质量评分（MAPQ）≥20 的结果用于后续分析。

基于对的段覆盖度，利用deepTools (v2.5.2)（Ramírez 等 2014）的plotCorrelation命令plotPCA来评估生物学重复之间的一致性，以及不同的读取MNase实验的恐情况；使用computeMatrix和plotHeatmap来可视化感兴趣的基因组区域（如中继起始位点TSS和中继终止位点TTS）的信号富集情况。

读段覆盖度数据使用UCSC基因组生物信息学工具集（https://github.com/ucscGenomeBrowser/kent）中的“bedGraphToBigWig”代码转换为bigWig文件，并在Broad研究所开发的集成基因组浏览器（IGV）上进行可视化（Robinson 等 2011）。

3.核小体鉴定、分类与预测

在重度 MNase 消化实验中，将过滤后的 MNase-seq 比对数据导入 R/Bioconductor 框架（版本 3.5.0），并使用nucleR架构（Flores 和 Orozco 2011）进行分析。长度在 260 bp 以下的双端读取段被裁剪至以 DNA 片段中心为原点的 50 bp计算全基因组每百万映射读段数（RPM）的覆盖度，并利用每个样本的总比对数进行归一化处理。使用以下nucleR参数进行噪声过滤和峰检测（Peak Calling）：pcKeepComp = 0.02，峰宽 = 147 bp，峰检测阈值 = 35%，最小重叠度 = 50 bp。

如果识别出的峰宽超过150 bp，则认为该峰包含两个以上重叠的核小体中心（dyads）。在峰宽低于150 bp的非重叠核小体位点中：

定位良好的核小体（Well-positioned, W）：定义为峰高得分为0.6且峰得分为0.4。

模糊定位的核小体（Fuzzy, F）：其余不满足上述条件的核小体均归为此类。

**核小体重复长度（Nucleosome Coverage，NC）**定义为被核小体重复长度（Nucleosome Repeat Length，NRL）**定义为包裹在组蛋白八聚体上的DNA长度连接DNA（linker DNA）的长度，即相邻核小体中心之间的距离；该指标使用NucTools脚本nucleosome_repeat_length.pl和plotNRL.R（Vainshtein等2017）进行说明。此外，使用R包NuPoP（Xi等2010）根据基因组DNA序列进行核小体定位预测，该工具利用四阶或一阶马尔可夫链（HMM）对隐连接DNA长度进行显式建模。NuPoP输出基于Viterbi算法预测的最佳核小体位置图谱，并据以计算预测的NC和NRL值

4.通过 DNS-seq 异构开放染色质量区域图谱

MNase敏感足迹（MNase Sensitive Footprints, MSFs）生物学重复和技术重复之间高度的可恢复性（皮尔逊相关系数） $r > 0.9$ ），将每个基因组及每次 MNase 消化实验的比对结果进行合并，生成计算每个基因组的 DNS 概况图。参照先前建立的差异 MNase-seq 数据处理流程（Turpin 等，2018），执行了以下连续步骤：归一化：将磁场（轻）和重度（重）MNase消化样本的映射读取段覆盖度归一化为RPM。计算DNS评分：将消化样本的读取段覆盖度下降重度消化样本的覆盖度，所得差值即为DNS评分（Differential Nuclease Sensitivity Score）。轨迹生成文件：制作使用基因组浏览器（Genome Browser）读取的数据轨道（Data Tracks）。峰检测（Peak Calling）：利用基因组件分段算法iSeg (v1.3.4)（Girimurugan 等，2018）识别正向峰（即 MNase 敏感位点）和负向峰即（ MNase 耐受位点）。

为了实现物种间及（亚）基因组间的比较，在运行 iSeg 之前增加了一个**分补归一化（Quantile Normalization）**步骤，对二倍体基因组（A2 和 D5）以及杂种和四倍体棉花的亚基因组（At 和 Dt）的全基因组 DNS 评分进行标准化。在鉴定 **MNase 敏感足迹（MSFs）**和MNase 耐受足迹（MRFs）（分别由正向和负向 DNS 峰代表，术语沿用自 Vera 等，2014）时，测试了一系列生物学截断值（Biological Cutoff, BC）的严谨度。最终采用优化的严谨度参数BC = 6.0（除补充材料在线文本 1）生成了最终的 MSF 列表。

5.亚核小体颗粒占用度（Subnucleosomal Particle Occupancy）

正如先前研究所报道的（Grossman 等 2018；Savadel 等 2021），来自突变 MNase 消化的短序列片段（0 到 130 bp）也可用于直接分析参与调控的亚核小体尺寸颗粒的处理情况。

使用 awk 和 BEDTools (v2.27.1)（Quinlan 2014），紧急通过提取了来自睡眠消化的每个短片段比对（0 到 130 bp）的几何中心，将其与步长为 5 bp、窗口大小为 21 bp 的滑动基因组窗口进行交集缺血。片段中心的平滑图谱被短一化为RPM，作为全基因组的亚核小体颗粒引发度（SPO）评分。

与 DNS 的相对评分不同，对各基因组间的 SPO 评分进行补充归一化（Quantile Normalization）会导致明显的信号损失。因此，对每个基因组生成的 BedGraph 文件分别使用了优化的严谨度参数进行 iSeg (v1.3.4) 分析（见补充材料在线文本 1）。得到最终的碎片列表代表了通过 SPO 鉴定的开放染色质域（ACRs）。

6.DNase-seq

本研究使用了之前报道的棉花幼叶公共数据（王等，2017，2018；韩等，2022），这些数据已从NCBI数据库下载（参见补充材料在线表S1）。

7.数据处理

对原始的ATAC-seq和DNase-seq读段进行接头去处和质量，并使用“Trim Galore”(v0.4.5) (Krueger 2012)进行过滤。然后，使用Bowtie2 (v2.3.4) (Langmead and Salzberg 2012)将干净的读段比到相应的参考基因组上，比对参数设置为“–no-mixed” --no-discordant --no-unal --dovetail”。针对ATAC-seq，测试了不同的峰检测（Peak Calling）方法（具体补充材料在线文本2），而DNase-seq则采用了MACS2方法。

8.HOMER 与 MACS2 峰检测（Peak Calling）

使用 Picard (v2.17.0) 默认参数（http://broadinstitute.github.io/picard/）去除重复读段（Duplicate Reads）。仅保留比对评分质量（MAPQ）至少为 20 的唯一比对读段用于峰检测。利用Phantompeakqualtools (v1.14)（Landt 等 2012）计算链交叉相关性（strand cross-correlation），并使用deepTools (v2.5.2)（Ramírez 等 2016）计算生物学重复之间的一致性。本研究采用了两种峰检测算法：HOMER (v4.10)（Heinz 等 2010）：运行其findpeaks工具，模式设为“region”，并设置峰之间的最小距离为150 bp。MACS2 (v2.1.1)（Zhang 等 2008）：运行其callpeak命令，使用参数“-f BAMPE”以仅分析正确配对的比对序列，并使用默认设置和假发现率（FDR）< 0.05 过滤推定的峰值。然后deepTools计算出的比对高度高度（皮尔逊相关系数） $r = 0.99$ ，斯皮曼相关系数 $r = 0.77$ ），使用BEDTools（v2.27.1）（Quinlan 2014）将每个工具在不同重复样本中检测到的峰进行合并。其次，使用BEDTools对HOMER和MACS2识别的峰区域取交集，仅保留被两种工具同时鉴定的区域作为后续分析的ATAC ACR（开放染色质区域）。

┌─────────────────────────────────────────────────────────────────────────────┐
│                         ATAC-seq Raw Data (BAM)                             │
└─────────────────────────────────┬───────────────────────────────────────────┘
                                  │
                                  ▼
┌─────────────────────────────────────────────────────────────────────────────┐
│                    Pre-processing & Quality Control                         │
│  ┌──────────────────────┐    ┌──────────────────────┐    ┌──────────────┐  │
│  │   Picard v2.17.0     │───▶│   MAPQ ≥ 20 Filter   │───▶│ Unique Reads │  │
│  │ (Remove Duplicates)  │    │ (Default Parameters) │    │              │  │
│  └──────────────────────┘    └──────────────────────┘    └──────┬───────┘  │
└─────────────────────────────────────────────────────────────────┼──────────┘
                                                                  │
                    ┌─────────────────────────────────────────────┼─────────┐
                    │                                             │         │
                    ▼                                             ▼         │
┌─────────────────────────────────────┐     ┌─────────────────────────────┐ │
│  Phantompeakqualtools v1.14         │     │    deepTools v2.5.2         │ │
│  (Strand Cross-correlation)         │     │  (Replicate Concordance)    │ │
│                                     │     │                             │ │
│                                     │     │  Pearson r = 0.99           │ │
│                                     │     │  Spearman r = 0.77          │ │
└─────────────────────────────────────┘     └─────────────────────────────┘ │
                                                                            │
                                                                            ▼
                           ┌────────────────────────────────────────────────┐
                           │           Peak Calling Algorithms              │
                           │           (Parallel Processing)                │
                           └──────────────┬─────────────────┬───────────────┘
                                          │                 │
                          ┌───────────────▼──┐   ┌─────────▼────────────────┐
                          │   HOMER v4.10    │   │     MACS2 v2.1.1         │
                          │   findPeaks      │   │     callpeak             │
                          │                  │   │                          │
                          │  - mode: region  │   │  -f BAMPE                │
                          │  - minDist: 150bp│   │  - FDR < 0.05            │
                          └────────┬─────────┘   └──────────┬───────────────┘
                                   │                        │
                                   │                        │
                                   ▼                        ▼
                          ┌─────────────────┐      ┌──────────────────┐
                          │ HOMER Peaks     │      │ MACS2 Peaks      │
                          │ (Replicates)    │      │ (Replicates)     │
                          └────────┬────────┘      └────────┬─────────┘
                                   │                        │
                                   ▼                        ▼
                          ┌─────────────────┐      ┌──────────────────┐
                          │  BEDTools       │      │   BEDTools       │
                          │  v2.27.1        │      │   v2.27.1        │
                          │  (Merge reps)   │      │   (Merge reps)   │
                          └────────┬────────┘      └────────┬─────────┘
                                   │                        │
                                   └──────────┬─────────────┘
                                              │
                                              ▼
                               ┌──────────────────────────────┐
                               │        BEDTools              │
                               │      Intersection            │
                               │  (HOMER ∩ MACS2)             │
                               │                              │
                               │   High Confidence Peaks      │
                               └──────────────┬───────────────┘
                                              │
                                              ▼
                               ┌──────────────────────────────┐
                               │     Final ATAC ACRs          │
                               │  (Open Chromatin Regions)    │
                               │                              │
                               │   ✓ High Reproducibility     │
                               │   ✓ Dual Algorithm Support   │
                               │   ✓ Stringent Filtering      │
                               └──────────────────────────────┘

9.Genrich峰检测（Peak Calling）

比对后续的后续处理以及多样本重复的集合峰检测，均通过Genrich (v0.6.1)（https://github.com/jsh58/Genrich）的一个命令统一完成。该工具由哈佛大学 FAS 信息学小组开发并经过了广泛测试。来自两个重复样本的比对文件由Genrich进行集合分析，并启用了以下选项：消除PCR重复 (-r)：排除增加产生的干扰。处理单端比对(-x)：通过延伸至平均片段长度来保留非成对的比对序列。修复问题区域 (-E blacklist.bed)：过滤掉已知的基因组黑名单区域（如高重复序列或人工造成的信号富集区）。峰检测阈值：设置最大 $q$ 值为0.05 (-q 0.05)，且最小曲线下面积（AUC）为20.0 (-a 20.0)。Genrich输出的文件为ENCODE slimPeak格式，其中列出了每个识别出的峰的基因坐标组、峰顶（Peak Summit）位置以及各项统计数据。

########################################## ACR 特征分析 ##################################################

10基因组注释

如前所述，本研究识别了多种来源的开放染色质区域（ACRs），包括MSF（MNase敏感足迹）、SPO区域（亚核小体颗粒对象区）以及ATAC-seq峰。此外，还应用了一个额外的过滤步骤，剔除雷蒙德氏棉（G. raimondii）中的基因组黑色名单区域（详见补充材料在线文本3）。

根据与最近基因的距离，这些ACR被分为以下三类：

基因内ACR（gACRs）：与基因区域重叠。

近端ACR（pACRs）：位于基因2 kb范围内。

最后ACR（dACRs）：基因距离超过2 kb。

为了比较 ACR 与非开放基因区域组的GC 结构，使用 BEDTools 的shuffle命令生成了对照区域：末端对照（排除基因及其 2 kb 侧翼区域）和基因内/近端对照（包含基因及其 2 kb 侧翼区域）；并nuc命令使用计算了每个 ACR 及置换对照区域的 GC 。

利用 R 封装ChIPseeker (v1.18.0)（Yu 等 2015），将 gACRs 和 pACRs 合并，并进一步细分为以下子类别：启动子（<1 kb, 1-2 kb, 2-3 kb）、外显子、内含子、下游区域（<1 kb, 1-2 kb, 2-3 kb）以及基因间隔区（距离 TSS 上游 >3） kb 且距离 TTS >3 kb）。

11.与转座元件（TE）的相关性分析

使用 EDTA (v1.9.5)（Ou 等 2019）流程对所有参考基因组进行了全基因组转座元件（TE）的注释。当 ACR 的坐标与 TE 区间发生重叠时，计算各 TE 超家族中 ACR 的比例。为了代表背景噪音，利用 BEDTools 的 shuffle 命令模拟了随机对照区域（这些区域在数量、区间宽度以及远端与基因内/近端区域的组成比例上均与真实的 ACR 保持一致）。通过置换检验（Permutation tests, $n = 1000$ ），评估了每个 TE 超家族中 ACR 的富集程度是否显著高于对照区域的零假设分布。**富集得分（Enrichment scores）**计算为：观测到的 ACR 比例与基于置换检验得出的平均比例之比，并进行 $log_2$ 转换。

12.差异可及性分析

本研究参照已建立的差异可及性（DA）流程（Reske 等，2020），利用 R 软件包 csaw (v1.16.1)（Lun 和 Smyth 2016），检测了由杂交和异源多倍体化引起的染色质可及性差异。

为了在不同棉种之间进行直接比较，所有 MNase-seq 数据均比对至相同的参考基因组，即 AD1 参考基因组或 A2 与 D5 的合并参考基因组；通过检查两种参考基因组得出的 DA 结果以减轻比对偏好性。

将经过质量过滤的比对读段对计入**滑动窗口（sliding windows）或特定的峰集（peak set）**中，以量化全基因组范围内的 MNase 信号，随后采用 TMM 或 Loess 方法进行归一化；研究中评估了多种分析方法，以确定最合适的 DA 工作流程（详见补充材料在线文本 4）。

所得的计数矩阵随后被导入 edgeR (v3.24.3) 统计框架，通过**经验贝叶斯（Empirical Bayes）估算离散度，并利用拟似然广义线性模型（Quasi-likelihood GLM）**拟合进行假设检验。具体实验设计如下：

二倍体中的轻度 vs. 重度：确定基础开放区域。

F1 杂种中的轻度 vs. 重度：确定杂种中的开放区域。

AD1（四倍体）中的轻度 vs. 重度：确定多倍体中的开放区域。

F1 (轻-重) vs. 二倍体 (轻-重)：代表杂交效应（Hybridization effect）。

AD1 (轻-重) vs. F1 (轻-重)：代表多倍体化效应（Polyploidization effect）。

13.基序发现与富集分析

本研究使用了 MEME Suite (v5.4.1)（Bailey 等 2015）并在默认设置下进行分析。利用 FIMO（Grant 等 2011）在 1 kb 启动子区域内扫描已知基序（Motif）的出现情况；同时，结合 XSTREME（Grant 和 Bailey 2021）和 AME（McLeay 和 Bailey 2010）进行基序的从头发现（de novo discovery）与富集分析。

XSTREME 结合了 MEME 和 STREME 两种算法进行从头基序发现，随后利用 SEA（Bailey 和 Grant 2021）进行富集分析；而 AME 则用于识别与对照序列相比，在给定序列中相对富集的已知基序。分析中，各（亚）基因组的启动子区（<1 kb）ACR 序列及其对应的启动子全长序列分别被用作输入序列和对照序列。

已知功能基序参考了 JASPAR 2018 植物核心非冗余基序库以及来自 plantTFDB v5.0（Jin 等 2017）的拟南芥基序。为了对富集的基序进行聚类，使用了 RSAT 矩阵聚类工具（Castro-Mondragon 等 2017），参数设置如下：-hclust_method average -calc sum -metric_build_tree Ncor -lth w 5 -lth cor 0.6 -lth Ncor 0.4 -quick。最后，利用 R 软件包 pheatmap（Kolde 2019），基于欧几里得距离生成热图并进行层级聚类。

14.RNA-seq 分析

使用 Sigma Spectrum 植物总 RNA 试剂盒（货号 STRN50）进行总 RNA 提取，并使用 BioAnalyzer（安捷伦，帕洛阿尔托，加利福尼亚州）进行定量。mRNA 文库使用 Illumina TruSeq RNA 文库制备试剂盒（Illumina，圣迭戈，加利福尼亚州，美国）制备，并在三条 HiSeq 4000 车道上进行双端 150 循环测序。共生成了来自 A2、D5、F1 和 AD1 样本的 12 个文库，每个样本平均包含 1100 万个读段对（参见补充材料在线表 S1）。在使用 TrimGalore（Krueger 2012）进行质量过滤并去除接头序列后，使用 Kallisto（Bray 等 2016）将双端读段伪比对（pseudo-aligned）至参考转录组。在 R 环境 3.5.0 版本下，使用 DESeq2（Love 等 2014）进行差异基因表达分析，并要求假发现率（FDR） $α<0.05\alpha < 0.05$ 以鉴定显著变化。

为了优化推断重复基因表达模式的方法，我们测试了以下比对策略：

D5-ref（以 D5 为参考）：利用雷蒙德氏棉（D5）参考基因组（Paterson 等 2012）和先前生成的物种特异性 SNP 索引（Page 等 2013）构建参考转录本序列。

D5 样本：直接比对至 D5 转录本。

A2 样本：比对至“伪 A2（pseudo-A2）”转录本（通过将 D5 基因模型上的物种特异性 SNP 替换为 A2 特异性 SNP 生成）。

F1 杂种：比对至“伪 A2”与 D5 转录本的合并库。

AD1（四倍体）：比对至“伪 AD1-At”与“伪 AD1-Dt”转录本的合并库（同样利用 AD1 特异性 SNP 生成）。

AD1-ref（以 AD1 为参考）：将所有物种的读段分别比对至陆地棉（AD1）的转录本序列（Chen 等 2020）。

individual-ref（各自独立参考）： F1 的读段比对至亚洲棉（A2）与雷蒙德氏棉（D5）合并的转录本；而 A2、D5 和 AD1 的读段则分别比对至它们各自对应的参考基因组转录本。

使用 pSONIC 流程（Conover 等 2021）推断异源多倍体基因组内以及不同参考基因组之间的共线性直系同源/同源基因关系（即 A2、D5、F1:At、F1:Dt、AD1:At 和 AD1:Dt 之间的关系），并以此为基础对不同策略产生的读段计数（Read Counts）进行比较。

基于 At 和 Dt 同源基因的总表达量（两者之和），将 F1 杂种和 AD1 四倍体的基因表达水平与 A2 和 D5 亲本进行比较，并随后归类为以下几种模式（Rapp 等 2009）：
加性表达（Additivity）：杂种或异源多倍体中的总表达量在统计学上等于双亲表达量的中值（Mid-parent value）。

A 基因组表达水平优势（A-genome ELD）：总表达量在统计学上与 A2 亲本相当，但显著区别于 D5 亲本和中值表达水平。

D 基因组表达水平优势（D-genome ELD）：总表达量在统计学上与 D5 亲本相当，但显著区别于 A2 亲本和中值表达水平。

超越亲本上调（Transgressive up-regulation）：总表达量显著高于 A2 和 D5 两个亲本。

超越亲本下调（Transgressive down-regulation）：总表达量显著低于 A2 和 D5 两个亲本。

基于 At 和 Dt 同源基因各自的拆分表达量，通过评估 F1 和 AD1 中两个同源基因（At 和 Dt）之间的差异表达来分析同源基因表达偏好性（HEB）。

参照先前报道的方法（Bao 等 2019）对**顺式（cis-）和反式（trans-）**调控分化进行归类：
首先，通过 A2 和 D5 表达量的 $log_2$ 比值测量整体调控差异 $A = \log_2(A2/D5)]$ ；

其次，通过对应同源基因表达量的 $log_2$ 比值测量顺式效应 $B = \log_2(At/Dt)]$ ；

最后，通过 $A$ 减去 $B$ 得到反式效应。

根据 $A$ 、 $B$ 以及 $A - B$ 的统计显著性，将调控演化划分为六个类别（如原文图 6c 所示）。此外，参照 Hu 和 Wendel（2019）的方法（如图 6c 所示），确定了杂交（Hr）、**异源多倍体化（Pr）以及基因组加倍（Wr）**对演化的具体影响。

                            RNA-seq 数据分析流程
================================================================================

   样本材料 (12个文库)
   ┌─────────┬─────────┬─────────┬─────────┐
   │   A2    │   D5    │   F1    │   AD1   │
   │(G.arboreum)│(G.raimondii)│(杂交种) │(G.hirsutum)│
   └────┬────┴────┬────┴────┬────┴────┬────┘
        │         │         │         │
        └─────────┴────┬────┴─────────┘
                         │
                         ▼
               ....................
              ┌─────────────────────┐
              │   质控 (TrimGalore) │
              │   过滤+去接头       │
              └──────────┬──────────┘
                         │
                         ▼
         ┌───────────────────────────────┐
         │   Kallisto 伪比对到参考转录组  │
         └───────────────┬───────────────┘
                         │
         ┌───────────────┼───────────────┐
         │               │               │
         ▼               ▼               ▼
   ┌──────────┐   ┌──────────┐   ┌──────────┐
   │ (i) D5-  │   │(ii) AD1- │   │(iii) Ind-│
   │  ref策略 │   │  ref策略 │   │ividual-  │
   │          │   │          │   │  ref策略 │
   └────┬─────┘   └────┬─────┘   └────┬─────┘
        │              │              │
        ▼              ▼              ▼
   D5→D5转录本    所有样本→AD1   A2→A2, D5→D5
   A2→pseudo-A2   转录序列       AD1→AD1
   F1→pseudo-A2              F1→A2+D5拼接
       +D5
   AD1→pseudo-At
       +pseudo-Dt
        │              │              │
        └──────────────┼──────────────┘
                       │
                       ▼
         ┌───────────────────────────────┐
         │   pSONIC pipeline 推断        │
         │   直系同源/同源异型体关系      │
         └───────────────┬───────────────┘
                         │
                         ▼
         ┌───────────────────────────────┐
         │      DESeq2 差异表达分析       │
         │        (FDR < 0.05)           │
         └───────────────┬───────────────┘
                         │
         ┌───────────────┴───────────────┐
         │                               │
         ▼                               ▼
   ┌──────────────┐              ┌──────────────────┐
   │  总表达量分析 │              │  同源异型体分离  │
   │  (At+Dt总和) │              │   表达量分析     │
   └──────┬───────┘              └────────┬─────────┘
          │                               │
          ▼                               ▼
   ┌──────────────┐              ┌──────────────────┐
   │  表达模式分类 │              │  HEB评估 (F1和   │
   │              │              │      AD1)        │
   │ • Additivity │              └────────┬─────────┘
   │ • A-ELD      │                       │
   │ • D-ELD      │                       ▼
   │ • 超亲上调   │              ┌──────────────────┐
   │ • 超亲下调   │              │ 顺反式调控分析:  │
   └──────────────┘              │                  │
                                 │ A = log2(A2/D5)  │
                                 │ B = log2(At/Dt)  │
                                 │ Trans = A - B    │
                                 └────────┬─────────┘
                                          │
                                          ▼
                                 ┌──────────────────┐
                                 │  6类调控演化分类  │
                                 └────────┬─────────┘
                                          │
                                          ▼
                                 ┌──────────────────┐
                                 │   演化影响评估    │
                                 │  ┌────────────┐ │
                                 │  │ Hr: 杂交    │ │
                                 │  │ Pr: 多倍化  │ │
                                 │  │ Wr: 基因组  │ │
                                 │  │     加倍    │ │
                                 │  └────────────┘ │
                                 └──────────────────┘

================================================================================
注: 
- A2 = G. arboreum (二倍体A基因组)
- D5 = G. raimondii (二倍体D基因组)  
- F1 = A2 × D5 杂交种 (二倍体)
- AD1 = G. hirsutum (异源四倍体)
- At/Dt = AD1中的A/D亚基因组
- ELD = Expression Level Dominance (表达水平显性)
- HEB = Homoeolog Expression Bias (同源异型体表达偏好)

15.组蛋白基因家族分析

从 HistoneDB 2.0（Draizen 等 2016）和 Probst 等（2020）研究中获取拟南芥（Arabidopsis thaliana）的组蛋白序列，并以此为查询序列，通过 BLASTP（设置 $e^{-5}$ 为截断值）在大隐棉花编码基因中进行搜索。利用 Seaview 第 5 版软件（Gouy 等 2021）的内置功能，使用 MUSCLE (v3.8.31)（Edgar 2004）进行多序列比对，并采用**邻接法（NJ）和最大似然法（ML）**进行系统发育分析。NJ 树：基于“泊松校正（Poisson correction）”模型构建，并进行 1,000 次重复的 Bootstrap 检验。ML 树：使用 PhyML (v3.0)（Guindon 等 2010）构建，采用默认的“LG”模型和 100 次非参数 Bootstrap 重复。对于每个组蛋白家族，在 MEGA11（Tamura 等 2021）中计算家族成员间的平均演化散度，即所有序列对之间平均每位点的氨基酸替换数（即总平均距离）。计算采用泊松校正模型，并针对每个序列对移除所有模糊位点（采用成对删除/pairwise deletion 选项）。

16.数据与代码可用性

本研究产生的数据已存入 NCBI 短读段存档库（SRA）：MNase-seq：项目编号 PRJNA529909ATAC-seq：项目编号 PRJNA1018916RNA-seq：项目编号 PRJNA529417所有使用的数据详情见补充材料在线表 S1。自定义脚本可在以下 GitHub 仓库获取：https://wendellab.github.io/cottonMNase-seq/。

标签

#网络