2022《Nature Genetics》棉属植物演化过程中的基

导读

棉花是重要的经济作物之一，是植物纤维、油料和蛋白的重要来源。棉属（Gossypium spp.）已经发现有52个种，其基因组包括A、B、C、D、E、F、G、K共8个二倍体类型（45个种）和1个AD异源四倍体类型（7个种），广泛分布于美洲、亚洲、非洲、澳洲。

不同棉种植株表型差异明显，包括多年生高大树状、小灌木或匍匐生长；大部分纤维粗短，紧紧附着在种子上，只有 A 基因组二倍体和AD基因组类型的四倍体被驯化为一年生栽培种，其余棉种皆为野生种。

由于长期的驯化和育种选择，栽培种棉花的遗传多态性逐渐降低，亟需从野生棉中鉴定未被利用的优异基因，拓宽育种可用的基因资源。厘清棉属基因组结构演化与纤维性状形成的关系，对于另辟蹊径提升棉花纤维品质至关重要。

表型多样性和演化创新可追溯到基因组序列的变异和调控网络的重整。本研究基于十个具代表性的二倍体棉属物种基因组构建了棉属泛基因组并记录基因组演化历史和谱系特异性转座子扩增对不同基因组的影响。泛3D基因组揭示了转座子驱动的基因组大小改变与高阶染色质结构重组和染色质互作组重整之间的演化联系。在此，通过将染色质结构变化与棉花纤维表型差异联系起来，并对1,005个纤维发育时期转录组进行测序从而鉴定到解码纤维长度遗传基础的调控变异。本研究结果展示了如何将泛基因组、泛3D基因组和遗传调控数据作为描述可纺棉纤维演化基础的资源并提出了关于基因组结构组织和调控演化的新见解，以及使用基于调控组的方法为棉花改良提供有价值的信息。

论文ID
原名：Genomic innovation and regulatory rewiring during evolution of the cotton genus Gossypium
译名：棉属植物演化过程中的基因组创新和调控重整
期刊：Nature Genetics
IF：41.307
发表时间：2022年12月
通讯作者：张献龙和Jonathan F. Wendel
通讯作者单位：华中农业大学和爱荷华州立大学
DOI号：10.1038/s41588-022-01237-2
实验设计

结果

1. 组装注释 7 个 2 倍体棉花基因组

为构建棉属泛基因组图谱，本研究选择了10个绵种进行基因组分析，包括来自两个物种的三个A基因组类型绵种（A_1a、A₁和A₂）以及代表棉属所有其他二倍体基因组类型的七个棉种，即 B₁、C₁、D₅、E₁、F₁、G₁ 和 K₂ 基因组（图1a）。该抽样方案包含三个已发表基因组，并涵盖所有主要地理区域和二倍体基因组类型。通过整合Nanopore长读长（126–161×）、Illumina短读长（52–79×）和高通量染色体构象捕获（Hi-C）数据，对本研究新测序的基因组进行组装，即G.herbaceum（代表南非两个基因组类型：野生型‘A_1a’和驯化型‘A1’），G.anomalum（非洲佛得角群岛‘B1’），G.sturtianum（澳大利亚‘C1’），G.stockii（非洲-亚洲‘E1’），G.longicalyx（东非‘F1’）和G.bickii（澳大利亚中部‘G1’）。这些基因组组装得到的Contig N50大小在7.7-28.8 Mb之间，基因组大小为750 Mb到2444 Mb不等，其中98.4%–99.9%的Contig锚定在 13 条假染色体上。每个基因组的长末端重复（LTR）组装指数范围为11.3至14.4，表明重复序列区域的高度完整性。与最近发布的基因组组装结果相比，本研究的组装结果在连续性和Contig数量方面都有很大改进。

这些组装的基因组拥有40,321-41,949个蛋白编码基因和5,992-16,606个非编码RNAs。基因组中的重复序列比例从69.5%到86.6%不等。使用着丝粒相关LTR序列和基于CENH3染色质免疫沉淀法预测每个基因组中的着丝粒区域，然后进行（ChIP-seq）测序。

所用棉花材料包括7个代表性二倍体野生棉种和3个A基因组类型的棉花。

2. Gossypium（棉属）谱系的演化和历史种群规模

系统发育：以 Gossypioides kirkii 作为外群，使用11,630个单拷贝直系同源基因进行棉属系统发育树的重建。结果表明，Gossypioides kirkii 与棉属物种大约在8.5-10 百万年前（MYA）发生分歧（图1a）。棉属 D 基因组类型的物种类群与该属其余物种为姐妹群，这与此前的相关研究结果相一致。而与先前结果矛盾的是，本研究显示F₁和B₁与A基因组类型（的物种类群）亲缘关系更近，G₁、K₂和C₁共同属于另一个包含G₁、K₂、C₁、D₅和E₁谱系的类群（图1a）。这两个主要的棉属类群之间的分歧发生在大约5百万年前，随后出现了更多近期的分歧，最近的一次是A₁和A₂的分岐，发生在它们独立驯化之前，距今大约0.57-0.85百万年。

祖先核型和基因组重排：每个谱系的祖先核型是使用10个棉属物种和 Gossypioides kirkii 之间28,556-31,720个直系同源基因从Gossypioides kirkii的12条染色体中推断出来的（图1a）。正如G.kirkii中已知的降序异常倍体所预期的那样，所有棉属基因组中的2、4和6号染色体在着丝粒区域附近出现重排，尽管从系统发育上讲，这些被认为属于祖先配置。相反，与Theobroma cacao中频繁发生的染色体重排事件相比，棉属内的重排是有限的，包括从祖先到澳大利亚进化支（C、G和K）中的单个相互易位、F谱系中的非相互易位（染色体5至9）以及A₂谱系中1号染色体和2号染色体之间的相互易位（图1a）。

图1 | 二倍体棉属基因组的演化和历史种群规模统计。

a，二倍体棉种的系统发育地位。使用 Gossypioides kirkii 作为外群进行每个棉属基因组的祖先染色体重建推断。根据棉属谱系之间的Ks峰进行分歧时间推断，同义替换率（r）的扩张度为3.48×10⁻⁹。每个物种中总基因家族的扩张和收缩分别以蓝色和红色标示。如图所示，比例尺表示的是1cm或5mm尺度。
b，不同棉属植物的地理分布。蓝色箭头表示多个棉属物种之间的基因流，箭头大小表示D值大小。根据历史记录，红色箭头表示从A₁物种到A₂基因组物种的两次人类介导的传播。该地图经R包ggplot2“maps”和公共领域的自然地球数据集(http://www.naturalearthdata.com)修改。
c，来自25个物种70个棉种的主成分分析。虚线框表示四个组：A基因组物种（A，第I组）；B，F基因组物种（BF，第II组）；C、K、G-基因组物种（CKG，第III组）；D基因组物种（D、E第IV组）。
d，D 值用于检测棉属谱系中古基因流的相对渗入程度。较小的点代表13条染色体的 D 值，较大的红点代表平均值。
e，除K₂外，二倍体棉种的历史种群规模。使用PSMC模型估计有效种群大小（Ne）（突变率 = 3.48×10⁻⁹）。

群体基因组和基因流：为探索棉属的历史种群规模，本研究整合了涵盖 25 个棉属物种的 70 个棉种的深度重测序数据（平均测序深度 23×）（图1b）。使用20,084,560个单核苷酸多态性位点（SNPs）进行的群体基因组分析表明，70个棉种被聚类为四个具有可变核苷酸多样性的类群（图1c）。通过使用四分类树进行D统计，在棉属物种中发现了广泛的基因流，通常主要涉及B、F和E基因组（图1d）。

历史种群规模：通过估算动态有效种群大小（Ne）对历史种群规模的分析表明，每个棉属物种都表现出显著的种群规模降低（开始于>100,000年前，部分与末次盛冰期相吻合），随后在约20,000年前恢复（图1e）。

通过不同棉种基因组的比较，阐明了棉属的演化历史：所有二倍体棉种在5百万年前（MYA）开始发生物种分化，A基因组与B1、F1基因组棉花分化时间更近，G1、K2和C1基因组棉花与A、B1、F1基因组棉花在4.3 MYA分化，E1和D5基因组棉花形成了相对独立的进化分支。通过与棉花的近缘物种叉柱棉（Gossypioides kirkii）进行比较分析，发现棉属与叉柱棉存在2–4次染色体重排事件，重排位置发生在着丝粒附近。

同时，发现A、B、F等基因组间存在广泛的基因流，并发现所有棉种在末次盛冰期（Last Glacial Maximum）发生了种群规模的缩小事件。

此外，研究发现转座子元件在棉花与叉柱棉分化后，发生了大规模的扩增，导致了基因组大小的不均等扩增，揭示了为什么棉属各基因组呈现至少三倍大小的变化（图1）。

3 Gossypium 的泛基因组

SV：为构建棉属泛基因组，本研究使用A₂基因组作为参考，鉴定到二倍体基因组之间的大量结构变异（SVs）。该分析中，在其他 9 个基因组中的每一个中都识别出150-746 Mb的缺失、2-68 Mb的重复、114-605 Mb的倒位和2-53 Mb的易位（图2a）。在2到9个基因组中共有的倒位变异（>1kb）共有14,827个。基于共线性比对、Hi-C热图、断点长reads比对和聚合酶链式反应（PCR）扩增结果，共在6条染色体上发现共有的倒位和特异性倒位存在（图2b）。在13号染色体和5号染色体之间发现了一个大的染色体间易位，这一变异似乎只出现在G₁、K₂和C₁基因组中，因此推测可能首先出现在G₁、K₂和C_1d的祖先基因组中（图2c）。还检测到将A₂与所有其他基因组区分开来的1号染色体和2号染色体之间的易位，这一结果与近期比较基因组学研究中报道的倒位一致。

图2 | 二倍体棉属的泛基因组和转座子扩增

a，A₂和其他九个棉属二倍体基因组之间的基因组比对，以识别共线性（SYN）块和SVs，包括缺失（DEL）、重复（DUP）、倒位（INV）和易位（TRA）。
b，基于A₂参考基因组的棉属基因组中的大型染色体（Chr）倒位。共线区块上的粉红色和黄色框分别代表共享和可变的倒位。Ma，百万年前；Ka，千年前。
c，棉属谱系间的基因组共线性。共线性区块以灰色连接。A₂和K₂中1号染色体/2号染色体和5号染色体/13号染色体之间的两个大的相互易位以绿色标示。一些大的倒位以红色标示。
d，10个亚基因组和12个二倍体基因组中棉属泛基因家族的存在/缺失变异（PAV）信息。热图显示了22个（子）基因组中的核心基因家族、20或21个（子）基因组中的软核基因家族、2-19个（子）基因组中的可变基因家族和特有基因（包括仅存在于一个（子）基因组中的基因家族）在不同的（亚）基因组中。
e，基于22个基因组和8个模拟基因组的成对基因家族比较估计的泛基因和核心基因家族的数量。每个黑点对应于由特定基因组组合估计的泛基因或核心基因家族大小。
f，每个基因组中的Gypsy、LARD、DIRS、Copia、TIR和LINE大小统计。
g，10个基因组中棉属逆转录转座子gypsy-like元件（Gorge3）的系统发育关系。以Gossypioideskirkii（G.kirkii）作为外群。
h，棉属谱系中特异性flLTRs和共有flLTRs的插入时间与染色体分布之间的关系。染色体位置被标准化为相应染色体长度的百分比。

核心和非核心基因集：结合部分已发表的基因组，本研究构建了一个包含42,394个基因家族和12,007个特异性基因家族的棉属泛基因组基因集（图2d）。共包括22个（亚）基因组中的17,079个（40.3%）核心基因家族、20个或21个（亚）基因组中的7,832个（27%）软核基因家族、2-19个（亚）基因组中的16,803个可变家族和680个基因组特异性家族（图2d）。每个基因组有20,661-21,652个核心基因、7,423-9,520个软核基因、5,590-12,346个可变基因和225-3,170个特异性基因。用迭代采样对泛基因大小进行建模，结果表明核心基因家族的数量随基因组数量（n）的增加而减少，但泛基因家族的数量在n = 15时达到平台期（图2e）。

转座子：记录不同种类转座子对棉属基因组大小演化的贡献（图2f）。共鉴定出3,773–32,684个全长LTR（flLTR）。在棉属与Gossypioides kirkii（8-9MYA）发生分歧后，在棉属中检测到几次flLTR的爆发，包括在~1MYA发生在三个A基因组祖先基因组中的一次fLTR爆发和~4MYA发生在C₁、G₁和K₂共同祖先基因组中的一次fLTR爆发。与G₁、K₂和C₁谱系相比，棉属逆转录转座子gypsy-like元件（Gorge3）的不同进化支在A基因组谱系中得到扩增（图2g）。还观察到Ty3-Gypsy、Copia和DIRS元件的谱系特异性或共有的扩增模式。研究发现特异性flLTR往往位于染色体的中间位置（图2h），这表明染色体臂中基因富集区域中的基因间部分在不同谱系分化后相对保守。共有的或谱系特异性的flLTR对基因类别和表达有不同的影响。

4. 泛基因组和泛 3D基因组的联合

为了研究不同棉种演化过程中基因组结构变异对同源基因转录调控分歧的影响，该研究构建了棉属的泛三维基因组（pan-3D genome）图谱。发现不同基因组结构变异（插入或缺失）导致拓扑相关结构域（TAD）发生重组（新的TAD形成或TAD融合），棉种特异的TAD边界形成更多地发生在不活跃的染色质区域，与棉种特异的转座子插入相关。

为探索转座子驱动的基因组演化对3D染色质结构组织的影响，本研究对包含所有核型的8个物种进行了Hi-C、ATAC-seq 和 RNA-Seq测序，以构建高分辨率的全基因组染色质互作图谱（A₂、B₁、C₁、D₅、E₁、F₁、G₁和K₂）。

AB：这些数据能够识别基因组区域中36%的具有A/B（即打开/关闭染色质）区室状态变化的单拷贝直系同源基因对，其中49%的变化是物种特异性的（图3a)。

TAD：结果表明在8个物种中共鉴定到 15,039（1,038–2,763）个拓扑关联域（TADs），包括17,859（1,288–3,191）个TADs边界，关联域长度在300到3,000 Kb之间。两个物种的保守TADs比率与相对基因组大小呈负相关（P = 0.0075）。该信息指向一个包含6,155个边界（核心：129，可变：4,388和特异性：1,638）的泛TADs边界数据集，表明特异性边界更有可能位于染色体中间，并包含更多物种特异性基因。TADs边界的保守程度（代表2到8个物种的保守性）与TAD分离分数、基因数量和染色质可及性呈正相关。与A2基因组中的TADs相比，在其他基因组中发现了242个neo-TAD（新TAD的形成）和193个TAD融合事件（图3e）。

图3 | 棉属中TADs的演化基因组学。

a，八个基因组中直系同源基因的A/B区室状态变化。
b，TAD边界在不同物种中的累积贡献。
c，Circos图显示了A₂上不同保守程度的TAD边界情况。从外到内的轨迹分别代表特定物种的边界、两到八个物种的共有边界和转座因子(TE)密度。中间曲线表示A₂上TAD边界的密度分布，具有不同的保守程度。d中显示的区域由红色箭头指示。
d，2号染色体上共线区块的代表性Hi-C热图。数字和箭头分别代表保守的和谱系特异性的TAD边界。蓝色和绿色方块分别代表正链和负链的基因。
e，基因组大小比率与两个事件的数量之间的相关性。
f，通过双侧Fisher精确检验测试的不同类型TAD边界中稳定和差异表达基因（DEG）的比例。
g，TAD结构变化示意图。该比率是不同分类中的SV长度除以总SV的长度。
h，在不同的比较组中，不同分类的数量（图3g）。
i，在g中显示的不同分类中缺失（蓝色）和插入（卡其色）的实际观测结果（箭头）和预期结果的分布（直方图）。双侧Wilcoxon秩和检验，P < 2.2×10⁻¹⁶。
j，基于ATAC信号、基因数（GN）和基因表达（GE）的TAD边界分类。
k，三种类型特定边界的比例。
l，不同保护程度的TAD边界周围的平均TE覆盖度分布图。
m，相对于保守边界，特定边界处不同类型flLTR的富集。
n，两组TE在不同物种中的扩增程度。虚线表示TE与基因组长度成比例的扩增程度。
o，A₂基因组中共有的、特定的和随机的TAD边界上年轻和古老flLTR的覆盖率。中心线，中位数；框限制，上四分位数和下四分位数；误差棒，1.5×四分位数间距；异常值，灰点。双侧Wilcoxon秩和检验，****P*< 2.2×10⁻¹⁶。

TAD与SVs：发现大型或小型基因组具有更多的neo-TAD或TAD融合事件，这与基因组之间发生的大量SVs一致（图3e）。这些事件中涉及的基因和位于保守边界内的基因相比更容易出现差异表达（图3f）。由于SVs对TAD重组的影响主要包括neo-TAD、TAD融合和大小变化，因此根据SVs在TAD内部或边界的位置将TAD分为六类（图3g），这些类别与A₂中7.5%–16.8%的TADs和35.1%–55.5%的边界有关（图3h）。结果表明7.1%的基因组插入和6.1%的缺失发生在TAD内部，9.6%的插入和9.1%的缺失发生在TAD边界（图3g）。其中，0.7%的插入和0.4%的缺失伴随着neo-TAD和TAD融合事件，分别显著高于随机改组边界中的比率（图3g，i；P < 2.2×10⁻¹⁶）。

TAD与转座子：为进一步探索转座子对TAD组织的影响，根据其局部的染色质状态和转录活性将所有TAD边界分为三种类型（图3j）。结果观察到代表相对不活跃边界的第3型物种特异性边界的比例高于分别具有高度和适度染色质可及性的第1型和第2型边界（图3k）。随着边界保守性的增加，转座子的覆盖率显著下降（图3l）。与保守边界相比，四类flLTR（Gypsy、DIRS、SINE和LARD；即第1组）在物种特异性边界中发生富集，而Copia、LINE和TIR（第2组）显示出相反的模式（图3m）。因此，第1组中flLTR诱导区域基因组大小的增加大于整体基因组扩展的预期，而对于第2组，观察到相反情况，表明快速扩增的转座子被富集在特定的TAD边界处（图3n）。具体而言，A2中的物种特异性边界包括较年轻的flLTR，而相对保守的边界具有更古老的flLTR（图3o）。

5. 染色质 loop 调控的互作组的重整

loop与互作：使用Hi-C矩阵，以5kb的分辨率水平在这8个基因组的染色体内共鉴定到1,012,351个（每个基因组74,379-226,135个）染色质loops。将染色质loops的锚定位置与ATAC-seq峰相结合，在8个物种中共推断出225,786个假定的长距离转录顺式调控元件（long-range transcriptional cis-regulatory elements）。与基因或顺式调控元件（CREs）相关的loop中共有348,966（36,879-69,431）个基因-基因（G-G）、103,954（8,280-23,497）个基因-CRE（G-E）和13,654（690-3,731）个CRE-CRE（E-E）loops；其余loops与基因或CREs无关。

为识别两个物种之间的保守loops，根据序列保守性以及是否发生染色质相互作用，可将两loops间关系归为四类情况（图4a和方法）。第I类主要为E-E相互作用，第II类中三种可能的相互作用间无较大差别，第III类和第IV类主要为G-G相互作用（图4b）。

• 对于第 I 类，83.1% loops锚定位置包含的PAVs显著高于其他类别（总体54.5%，卡方检验，P < 2.2×10⁻¹⁶）。
• 对于第 III 类，其他基因组中两个锚点之间跨越的距离大于D₅中的距离，甚至比不同基因组大小的差异程度更大，而第IV类则相反（图4c）。
• 第 II 类的ATAC信号强度明显低于第IV类（图4d）。
  图4 | 染色质loop介导的调控网络重整

a，示例展示了两个loop之间的四类情况。RC代表Hi-C reads数量。FDR，错误发现率。
b，折线图展示了两个基因组之间三种loop对中上述四类的比例。橙色、红色和蓝色圆点分别代表基因-基因（G-G）、基因-CRE（G-E）和CRE-CRE（E-E）相互作用。
c，散点图表示III型（红色）和IV型（蓝色）loop在两个物种间的距离。虚线表示log₁₀（另一个基因组的大小）除以log₁₀（D5基因组的大小）。
d，II型（蓝色）、IV型（红色）和随机（黄色）loop锚点的ATAC信号强度（log₂（RPKM（ATAC）+0.1））。双侧Kruskal–Wallis检验，****P* < 2.2×10⁻¹⁶。中心线，中位数；框限制，上四分位数和下四分位数；误差棒，1.5×四分位数间距；异常值，灰点。
e，八个基因组中核心和特异性CREs的数量。
f，具有高、中、低保守性的octad基因的环数（LN）和保守性得分（CS）的标准化卡方值。保守分数表示每个octad基因的平均保守物种数。
g，AUX/IAA相关基因家族的染色质loop介导网络。中间是A₂基因组中的环状网络，虚线圈出的是最大的网络。红线代表特异性互作，橙色线代表保守互作。红色、蓝色和灰色点分别代表直系同源基因、假定的直系同源CREs和缺失元件。其他基因组中最大网络的情况展示在外侧，其中实线和虚线分别代表存在和缺失的互作。外侧网络由f中的红色箭头指示。

8个物种染色质loops的比较表明，随着保守程度的增加，G-G loops的比例增加，转座子在锚点上的覆盖率降低。两个实例证明了直系同源基因之间距离的变化往往伴随着loops的变化和基因的差异表达。对8个基因组中CREs进行比较的结果表明，共有165个保守核心CREs、29,418个可变CREs和136,920个物种特异性CREs（图4e）。与物种特异性CREs相比，保守和可变CREs倾向于调节更多高表达水平的基因，并且涉及的loops距离更远。此外，53.9%的物种特异性CREs与转座子重叠，比例显著高于变量（14.1%）和保守CREs（9.41%；卡方检验，P < 2.2×10⁻¹⁶）。正如预期的那样，更高比例的特异性CREs与PAVs重叠，暗示转座子扩增和PAVs在CREs演化中的作用。根据涉及的loops和保守基因的数量变化（以卡方值表示），在23,693个octad基因（8个基因组中直系同源基因一一对应）中分别有7,126个、10,433个和6,134个octad基因参与高、中、低保守的染色质互作网络（图4f）。与AUX/IAA基因家族相关的调控网络如Fig.4g所示，它详细展示了相对保守的调控网络中ARF3基因（图4g）。这一网络中82%（23/28）的loops和所有（17个）CREs在至少两个物种中是保守的。

发现不同棉种间基因与基因之间的染色质互作比非编码调控区与基因之间的互作更保守，发现结构变异、转座子插入和染色质开放程度变化驱动了非编码调控元件与同源基因之间调控网络的演化，为棉花中非编码调控区的结构和功能研究提供了丰富资源

6. FL 的调控变异和遗传基础

尽管棉属植物具有广泛的基因组多样性，但研究者们对与纤维质量相关基因型的功能性遗传变异知之甚少。本研究对G. arboreum的 216 个不同个体进行基因组重测序（平均测序深度为20.9×）共鉴定出142万个 SNPs（次要等位基因频率 > 0.05）。基于这些 SNPs，全基因组关联分析研究（GWAS）共鉴定到11个纤维长度（FL）的QTLs，发现这些 QTLs 有利等位基因的积累能够影响FL。此外，还鉴定到7个与纤维强度（FS）、2个纤维伸长率（FE）和6个纤维均匀性（FU）相关的 QTLs。在这26个与纤维质量相关的QTLs中，有18个此前未被报道。

在跨越初级细胞到次级细胞伸长的五个时间点（开花后（DPA）4天、8天、12天、16天和20天；图5a）对来自这些个体的1,005个纤维样本进行转录组测序。表达基因的主成分分析（每百万个mapped fragments中对比到每1 Kb外显子上的fragments数量 > 0.1）结果表明，属于每个发育时期的转录组大致分别聚类在一起。

图5 | eQTL定位和候选基因的发现

a，J85和J98中发育的纤维特征。比例尺为1cm尺度。
b，转录组的主成分分析。
c，在每个纤维采样阶段具有不同数量独立eQTLs的eGenes数量。
d，在不同数量的采样阶段中共有的顺式eQTL的比例，橙色表示具有一定效果的eQTLs，青色表示有效应的eQTLs（不考虑效果大小）。
e，在5个采样阶段间成对共有eQTL的效果大小情况。
f，具有Hi-C loop支持的eQTL和基因之间调控关系的示意图。
g，非编码loop锚点中基因间eQTLs的富集。箱线图展示真实eQTLs与随机变异集（n = 100）的比值。
h，三个发育阶段中Hi-C loop支持的和模拟的eQTL-eGene对比例。卡方检验，P < 2.2×10⁻¹⁶。
i，eQTL热点的基因组分布。外侧堆叠直方图表示2 Mb窗口中的热点数。在circos图内部展示的是一个由eQTL热点（Hot145）介导的调节网络。黑色比例尺标示的是10个热点尺度大小。
j，受Hot145调控的52个eGene表达水平的聚类。
k，具有QTL FL5和FL6中前导SNP不同基因型的样品的纤维长度。
l，在4DPA和8DPA时，基因Garb_05G040220的FL和eQTL区域关联图。用紫色菱形突出显示FL前导SNP（Chr5:90700977）。
m，按eVariant Chr5:90700977的基因型分类的两个棉花组中Garb_05G040220的表达水平。使用双侧Wilcoxon秩和检验进行统计检验，精确P值如图所示。
n，顺式eQTL介导的Garb_05G040220调控由Hi-C loop证据支持。连接eQTL和目标基因的loop以红色突出标示，其他loop以灰色标示。黄色条带代表eQTL和loop锚点的共同定位。对于g、k和m，中心线、中位数；框限制，上四分位数和下四分位数；误差棒，1.5×四分位数间距。

eQTL分析：为表征纤维发育中基因表达的遗传调控，本研究进行了基因表达QTL（eQTL）分析。结果表明，每个时期平均有6,150个基因被鉴定为顺式eGenes，其中17.8%具有两个或多个条件独立的eQTLs，81.3%的eGenes可以在至少两个时期中出现（图5c）。位于基因及其上下游2-kb区域的顺式eQTLs占所有独立顺式eQTLs的49.3%。平均鉴定出3,444个反式eGenes，其中30.7%具有两个或更多反式eQTLs。为比较不同发育阶段的顺式eQTLs交集情况，本研究使用了多变量自适应收缩（mash）方法，重点关注所有阶段每个基因的最强顺式eQTL。仅考虑效应方向时，72.7%的顺式eQTLs出现在所有时期的共有交集中；然而，当关注效果大小时，只有20.3%出现在所有时期的共有交集中（图5d）。不同发育阶段之间的成对eQTL分析显示12、16和20 DPA的交集程度较高（图5e）。顺式eQTL效应大小和顺式eGene表达水平的相关性分析表明，15.3%（1785/11670）的顺式eQTL具有跨阶段的显著相关性。

Hi-C和eQTL数据的整合分析：为研究远距离顺式调节的潜在3D空间基础（图5f），本研究生成了来自两种栽培类型的G. arboreum（A2基因组）【不同之处在于具有短（J85）或长（J98）纤维（图5a）】三个发育阶段（4 DPA、12 DPA 和 20 DPA）纤维的原Hi-C图（5 kb分辨率） 。通过对Hi-C和eQTL数据的整合分析，发现基因间区域的eQTLs在染色质loops的非编码锚点位置被富集（图5g），表明染色质loops与远距离遗传调控相结合。这一结论也被进一步证实：与模拟数据相比，实际eQTL-eGene对更可能与loops重叠。在五个时间点共鉴定到490个eQTL热点，调节4,832个eGenes。其中，一个调控热点与5号染色体上FL相关QTL（FL5）重叠（图5i）。据推测，该热点可调节52个基因，这些基因被富集到生物过程植物型细胞壁组织（GO：0009664）、蛋白质转运（GO：0015031）和KEGG途径核糖体（ko03010）（图5j）。

FL5的有利等位基因可使FL增加6.7%（图5k）。根据这52个基因的表达模式，棉种可分为FL显著差异的两组。通过整合全转录组关联研究（TWAS）和共定位分析，以优先考虑FL相关GWAS QTLs假定的结果基因。在TWAS分析中，共有191个基因在至少一个时间点达到FDR阈值0.05，其中12个通过Bonferroni校正显著性。通过评估GWAS和eQTL关联的共定位，9个基因被优先考虑与FL有关。此外，本研究注意到基因Garb_05G040220的eQTL信号在4 DPA和8 DPA时与FL QTL（FL6）共定位的概率很高（PP.H4 = 0.98和0.88）（图5l）。TWAS还在五个纤维发育阶段中的四个阶段鉴定到该基因（Bonferroni调整的P值在4 DPA、12 DPA、16 DPA和20 DPA时分别为0.004、0.019、0.006和0.046）。Garb_05G040220的表达水平与FL呈负相关（图5k，m），其在拟南芥中的直系同源基因（AT1G28570）编码了一种SGNH水解酶型酯酶超家族蛋白，表明它可能是一种负调节因子，用以响应纤维发育过程中的含氧化合物。特别是，eQTL和Garb_05G040220之间的调控关系得到染色质互作结果的支持，染色质互作与该基因在短纤维绵种J85中的高表达相吻合。

7. 棉花纤维的演化创新

为深入了解长纤维的基因组创新，本研究首先将所有二倍体物种分为三组，包括无纤维组（K₂）、稀疏短纤维组（C₁、G₁、E₁、D₅、F₁和B₁）和长纤维组（A₁和A₂）（图6a）。

泛基因分析表明，在无纤维物种、短纤维物种和长纤维物种中分别发现了1,077、1,697和458个基因。在pan-LTR分析中发现3,136、16,472和2,943个基因受到特定转座子插入的影响。A₁和A₂中存在而B₁、C₁、D₅、E₁、F₁、G₁和K₂中不存在的964个直系同源基因广泛参与生长素和糖代谢。同时，五个已报道与纤维发育相关的功能基因在不同的棉花属间存在PAV。

研究棉花纤维演化和遗传调控的框架，包括基因组多样性、表达差异和遗传调控差异。

对A₂和B₁物种纤维转录组的比较分析表明，6,947个直系同源基因在纤维发育过程中具有相似的表达模式，9,603个具有不同表达模式，1,499个具有相反的 表达模式。具有不同表达模式的基因被富集到微管运动活性、细胞壁、糖代谢和脂肪酸合成方面（图6b）。

A₂和B₁中的直系同源表达模式映射到A₁和E₁。结果表明，在长纤维（A₁、A₂）和短纤维（B₁、E₁）物种中，1,460和40个表达模式差异的基因被富集到微管运动活性和细胞氧化还原稳态中。接下来，使用TWAS和A₂共定位分析确定的195个候选基因，将与纤维质量相关的基因组变异和调控关系制成表格。基于基因组变异数据，对其他二倍体棉花基因组中的285个eGene-eQTL关系进行了该分析。仅存在于A₂基因组中的21个eGene-eQTL代表A₂特异性调控关系，例如AGAL2、CESA7和ACT3。存在于A₁和A₂中但在其他基因组中不存在的64个eGene-eQTL代表了A基因组特异性调控关系。9个基因组中与PAV相关的178个eGene-eQTL代表了可变的调控关系，例如CSLG3、SGNH、MYB4和EXPA4。存在于所有10个基因组中的其他22个eGene-eQTL显示出保守的调控关系（图6c）。以FL的两个串联重复基因（Garb_05G040220、Garb_05G040230）为例，其表达与FL呈负相关（图5m）。这两个基因都存在于A₁、A_1a、B₁和F₁基因组中，但在D₅、E₁和G₁基因组中丢失，而Garb_05G040230在K₂和C₁基因组中丢失（图6d）。进一步发现，Garb_05G040220在B₁纤维发育过程中高表达，但在A₂中低表达（图6e），这与A₂中FL变异和基因表达变化之间的关系一致。这些结果表明该基因可能是与棉花纤维演化相关的重要候选基因。

图6 | 纤维演化过程中的基因组和调控创新

a，研究棉花纤维演化和遗传调控的框架，包括基因组多样性、表达差异和遗传调控差异。
b，在纤维发育过程中，A₂和B₁物种之间具有相似和不同表达模式的直系同源基因的GO富集。比例尺为1cm尺度。
c，通过TWAS鉴定eGenes和eQTL的PAV并在10个二倍体基因组中对这些PAV进行共同定位分析。
d，10个棉属基因组中Garb_05G040220和Garb_05G040230的微共线性分析。
e，纤维发育过程中A₂和B₁基因组中Garb_05G040220和Garb_05G040230的表达水平。

研究发现基因组结构变异对纤维品质性状有影响。该研究对216份亚洲棉材料组成的群体进行纤维发育5个时间点的1005份转录组进行研究，鉴定了引起材料间基因差异表达的遗传调控位点（eQTL），解析了纤维长度变化的遗传基础。通过整合泛基因组、泛三维基因组、群体遗传等数据，揭开了棉花纤维性状从“从无到有”演化的基因组基础，探讨了纤维伸长相关基因在不同棉种中的序列变异和表达模式变化，发现SGNH基因的拷贝数变异和调控演化与不同棉种纤维品质的差异相关（图3）。

讨论

来自比较基因组学的研究结果正以越来越快的速度出现，这得益于代表多个物种和/或同一物种的多个品系的基因组组装以及转录组学和其他组学资源。这项努力的一个基本认识是，基因组比以前想象的更具活力，基因内容、非基因部分和多层次的监管控制都在快速演化。棉属就是这些进步的例证，迄今为止已生成超过15个棉属基因组以及3D和表观基因组查询工具的快速发展。本研究利用这些发展来构建棉属的泛基因组，阐明在不同分类水平上保守的基因和基因组特征以及那些更不稳定或“可有可无”的特征。系统发育和分子演化分析结果增加了对棉属生物演化的理解，包括起源于它在全球范围内传播到多个大陆的干旱地区（主要是热带地区）的基因组多样性。

鉴于谱系特异性转座子更有可能与棉属中物种特异性无活性TADs的形成相关，但这些谱系特异性基因组结构是否通过创建物种屏障或环境适应与物种形成相关仍然未知。这一发现明确指出了转座子在塑造植物高级染色质结构方面的功能必要性，这在哺乳动物中是一致的。尽管转座子通过启动子或局部表观遗传沉默与转录调控建立了联系，但本研究提供了一个额外的维度来了解转座子相关远距离调控元件的演化及其对表达的影响。因此，就一个属的基本或最小基因组而言，应该包括基因间序列，例如支撑染色质结构或充当与编码基因相关的远距离调控元件的转座子。与之前主要关注基因组序列变异变化的泛基因组研究相比，本研究的泛3D基因组研究为了解转座子驱动的基因组扩展和染色质拓扑结构创新之间的关系提供了一个窗口。综上所述，本研究阐明了棉属植物的演化史，为植物三维基因组的演化研究提供了参考，为该领域的研究提供了丰富的可用资源。

这项研究构建了首个覆盖所有二倍体棉花基因组类型的泛基因组和泛三维基因组图谱，揭开了棉花种质资源多样性形成的基因组结构演化特征，推动了棉花纤维品质等性状演化的遗传机制研究。研究鉴定了野生棉中可供棉花育种利用的优质基因，将助力破解优质棉育种的瓶颈问题。研究也为其他植物种属水平的种质资源演化和性状形成机制研究提供了参考。

转自：
https://new.qq.com/rain/a/20221208A06Y1C00
https://mp.weixin.qq.com/s/WYepuZWCCrX4Wb81VyA6iw