澳大利亚东南部公牛海带（Durvillaea amatheiae）对未来海洋变暖的适应能力受损：基因组学证据与保护启示

本综述系统分析了澳大利亚东南部特有公牛海带（Durvillaea amatheiae）在快速海洋变暖背景下的遗传适应能力。研究通过中性（FST、HE）和适应性（RDA、LFMM2）基因组分析，揭示了该物种因洋流隔离和有限扩散力导致的种群近交与遗传多样性丧失。温度相关候选位点（SNP）沿热梯度呈现选择特征，尤其在温暖种群中等位基因固定化表明其适应未来变暖的能力严重受损。研究为气候脆弱性物种的保护策略提供了关键分子生态学依据。

　　

引言

海洋变暖正驱动许多海洋物种向极地迁移，导致局部海洋生态系统紊乱和区域生物多样性变化。在暖温带地区，海带等关键基础物种的丧失已在多个区域发生，有时导致海带森林生态系统崩溃。全球海带森林通过提供资源和生态系统服务为全球经济贡献约5000亿美元，而澳大利亚的大南部礁每年为澳大利亚经济贡献约100亿澳元。特别是，澳大利亚东部的海洋变暖速度比全球平均水平快两到四倍，气候模型预测，在更严重的未来气候情景下，到2100年一些海带物种将局部灭绝。因此，这些基础大型藻类物种的丧失对依赖健康海带生态系统的社区和工业以及以它们为食物和庇护所的物种具有重大影响。

澳大利亚东南部特有的公牛海带Durvillaea amatheiae预计在温和的代表性浓度路径（RCP）情景下将经历大幅范围收缩，但在剧烈变化下，预计到2100年在RCP 6.0下将在澳大利亚大陆灭绝。该物种沿澳大利亚大陆海岸从新南威尔士州（NSW）的Aragunnu到维多利亚州的威尔逊岬角呈斑块状分布，并向南进入东部巴斯海峡岛屿和塔斯马尼亚东部、南部和西部，在那里可能与D. potatorum同域分布。在新南威尔士州，D. amatheiae分布分散，位于低潮间带和浅潮下带的波浪暴露岬角边缘。其北部分布界限据称在近几十年已向南退缩约20公里，从Bermagui附近转移到Aragunnu。此外，新西兰的同属Durvillaea spp.因2017/2018年的海洋热浪而经历了局部灭绝，突出了该属的热敏感性。这种快速发生的区域灭绝揭示了新南威尔士州D. amatheiae面临类似命运的脆弱性。新南威尔士州南部海洋变暖加剧对该物种在该地区剩余的分散种群构成了相当大的风险。因此，迫切需要确定这一关键物种的适应能力以及未来在澳大利亚东部范围收缩的风险。

遗传适应能力在物种适应和生存于变化条件的能力中扮演重要角色。探索适应能力的一种方法是确定遗传多样性和结构水平，包括现存种群间的基因流动水平。遗传多样性本身，尤其是适应性多样性，对于赋予种群选择能力和适应不断变化的海洋条件非常重要。基因型与环境分析可用于识别与温度等环境压力相关的个体候选位点。通过检查这些遗传变异如何沿环境梯度变化，我们可以识别种群可能难以适应的关键阈值。当环境条件变化速度快于种群通过自然选择进化的速度时，种群崩溃可能发生，特别是在遗传隔离、无法从其他地方获得有益基因的种群中。在本研究中，我们检查了遗传多样性和结构，并评估了可能表明Durvillaea amatheiae大陆种群未来热适应能力受损的变异模式。鉴于这些种群之间的碎片化以及可能导致遗传隔离的生活史（基因流非常低）和近亲繁殖，我们检验了这些种群是遗传隔离的假设，并表明鉴于它们作为物种向赤道范围的边缘位置，一些种群可能处于其热适应能力的上限。

材料与方法

研究区域

2022年1月，从澳大利亚新南威尔士州南海岸和维多利亚州东部约230公里段的六个地点（图1，表1）随机采集了Durvillaea amatheiae的叶片样本（每个种群25-30个）。个体在低潮时从潮间带边缘采样，组织取自叶片远端边缘。注意分开叶片以避免两次采样同一个个体，样本来自相距1至约20米的植株。采样点覆盖了D. amatheiae大陆范围的大部分（图1），包括现今的北部分布界限Aragunnu, NSW，向南和向西至维多利亚州的Cape Conran。在该物种范围的这一部分，分布非常斑块化，种群间存在巨大的空间间隙；实际上，采样地点代表了新南威尔士州和远东部维多利亚州大部分已知、可到达的种群。该研究区域代表了预计在未来几十年内将发生变暖驱动范围退缩的部分海岸。

DNA提取与基因分型

使用Qiagen DNEasy Plant Kit从叶片子样本中提取基因组DNA。处理遵循制造商针对富含酚类化合物的样品的方案，没有变化。样品在最后一步使用50 μL提供的洗脱缓冲液洗脱，并使用NanoDrop、Qubit测定和凝胶电泳进行质量检查。

对于单核苷酸多态性（SNP）基因分型，将每个样本的20 μL提取DNA送至Diversity Arrays Technology Pty Ltd (DArT; Canberra, Australia)。DArT组织提供了一个全基因组分析流程，不需要参考基因组。使用高通量DArTseq技术对Durvillaea amatheiae DNA进行基因分型。采用基于PstI的复杂性降低方法用于富集单拷贝序列的基因组表示。该方法涉及用切割酶PstI消化DNA样本，与一组次级频繁切割限制性内切酶配对，与位点特异性适配器连接，并扩增连接适配器的片段。限制酶对消化后，连接与PstI悬垂兼容的寡核苷酸适配器，并按照标准方案扩增连接适配器的片段。为了开发SNP，使用两个PstI兼容适配器对应两个不同的限制酶悬垂来优化DArTseq技术。按照所述程序生成基因组表示。使用HiSeq2000实施下一代测序技术以检测SNP标记。使用DarTsoft14和DArTdb分析序列数据。

SNP调用与质量控制

最初使用DArTseq?平台对161个Durvillaea amatheiae个体样本进行了基因分型，产生了总共60,547个SNP位点，平均读取深度为14.76，缺失数据为29.23%。原始数据通过DArT的几个专有过滤器以生成最终用户的高质量数据集：对条形码区域应用严格阈值（75%的条形码Phred ≥ 30；50%的整个读取Phred ≥ 10），这确保了可靠的样本分配。然后将序列截短至69 bp，使用3-bp汉明距离聚类，错误校正（Phred < 20 校正为Phred > 25的标签），并通过DArT GS14流程进行SNP调用，使用技术重复和孟德尔验证来确保标记质量和可重复性。作者收到后，应用了额外的质量控制过滤器以提高SNP质量，同时优化可用于群体基因组分析的个体和位点数量。基于DArTseq?流程的描述性统计的阈值使用R包dartRverse v.2.7.2应用于原始的DArT提供的SNP数据矩阵。质量控制过滤器部署如下：对于读取深度，为了保留最小深度为5到最大深度为30读取的位点（图S1），我们使用函数gl.filter.rdepth移除覆盖不足（可能导致基因分型错误）和覆盖过高（可能指示重复区域或比对伪影）的位点。然后使用函数gl.filter.reproducibility，阈值0.99，仅保留在重复样本中至少99%时间被相同调用的SNP——这是基因分型一致性和可靠性的指标。基于调用率过滤，使用函数gl.filter.callrate移除缺失数据超过80%的个体（阈值=0.2；图S2），并排除调用率低于0.7的位点，因此位点需要在至少70%的个体中被基因分型才能保留（图S3）。调用率过滤减少了低质量样本和稀疏基因分型标记的影响，这些可能会偏倚群体分析。然后应用次要等位基因频率（MAF）阈值0.05，使用函数gl.filter.maf移除稀有变异。稀有SNP通常对群体水平分析信息量较少，并且可能更容易出现基因分型错误。为了最小化冗余或非独立位点，使用函数gl.filter.secondaries过滤次级SNP，仅保留每个69 bp序列标签的一个SNP，基于进行1000次模拟后选择信息量最大的变异（图S4）。过滤次级有助于确保SNP不是物理连锁的，并降低违反独立性假设的风险。最后，使用函数gl.filter.allna移除具有100%缺失基因型（所有位点NA）的个体，因为它们没有为下游分析提供有用信息，并且可能扭曲摘要统计。过滤后，数据集1有16,237个SNP，跨越154个个体。

第二个数据集，数据集2，是通过进一步过滤偏离哈迪-温伯格平衡（HWE）在α=0.05的SNP生成的。这产生了一个包含8323个推定的中性、非适应性SNP位点的数据集，保留在154个个体中用于下游群体遗传分析。由于Durvillaea没有注释的参考基因组，排除了同义替换或非编码区域的识别，我们使用HWE偏差作为代理来识别可能受选择或其他非中性过程影响的位点。尽管我们承认基于基因组注释的替代方法在这些资源可用时更可取，但我们的方法与非模型系统中的既定实践一致。

使用数据集2按照在R studio中部署的dartR.popgen和dartR.spatial分析框架分析了群体间的中性遗传分化。使用gl.dist.pop函数处理群体距离矩阵（欧几里得距离），并使用gl.plot.heatmap可视化。使用dartR函数gl.pcoa进行主成分分析，并使用gl.pcoa.plot函数可视化前两个主成分。使用R中的poppr包生成几个群体水平遗传多样性的估计，包括观察到的（H_O）和预期的（H_E）杂合度以及等位基因丰富度（A_R）。对于A_R，等位基因计数通过使用R包hierfstat中的allelic.richness命令通过最小基因分型个体数进行稀化。使用poppr中的private_alleles函数计算每个群体的私有等位基因。使用R包adegenet中的Hs.test函数估计地点间遗传多样性测量的统计差异。使用hierfstat中的basic.stats函数计算偏离随机交配的F_IS（近交系数）。

随后通过计算全局群体分化（F_ST）及其95%置信限，以及使用dartR中的gl.fst.pop函数通过置换（999次）确定显著性的群体 pairwise F_ST测量来进行群体遗传结构测试。使用Benjamini-Hochberg错误发现率（FDR）程序校正多重测试效应。使用gl.amova函数进行分子变异分析（AMOVA），使用一个在样本地点之间划分变异的模型，基于999次置换的随机化检验显著性（表S1）。使用函数gl.ibd进行距离隔离（IBD）的Mantel检验（图S5）。最后，在R包LEA中实现了稀疏非负矩阵分解（SNMF）算法。该算法估计每个样本的遗传祖先成分。对于本研究，对每个K值（1到6）执行了10次运行，α=100。选择最佳推定祖先群体数量的指导是交叉熵准则（其中，对于K，交叉熵曲线图形成一个“膝”或“平台”），并使用barplots函数可视化最佳运行的结果。

环境数据

鉴于驱动澳大利亚东南部大型藻类分布的最强环境因素是水温，我们将用于基因型-环境关联（GEAs）的变量限制为水温。这些变量源自涵盖2001-2020年期间（每日时间分辨率）的海面温度栅格数据集，由Nimbs, Champion等人开发和使用。为初步评估GEAs，为三个遥感温度变量提取了点数据（最接近采集点纬度和经度的栅格单元格）。选择的变量是最大海面温度（MaxSST；°C）、最小海面温度（MinSST；°C）和平均海面温度（MeanSST；°C）。第四个变量，海面温度范围（SSTRange），通过从每个采样点的MaxSST减去MinSST计算。

该环境数据集首先使用R包LEA中的函数pairs.panels通过成对散点图分析以降低维度/共线性（图S6）。进行此操作是因为一些GEAs在使用高度相关的预测因子时会遇到问题（回归参数方差膨胀）。通常，r > 0.7被视为移除相关预测因子的良好阈值。使用R包vegan中的vif.cca函数通过分析方差膨胀因子（VIF）进一步检查共线性。通过共线性和方差膨胀检验的变量减少导致保留两个温度相关环境变量用于进一步分析：MaxSST和SSTRange（表1）。

全基因组范围的适应性变异特征

进行基因组扫描以识别可能处于选择压力下的候选SNP，使用数据集1。此步骤使用了两种不同的模型：冗余分析（RDA）和潜在因子混合模型（特别是LFMM2）。冗余分析用于执行约束排序，模拟遗传变异和环境预测变量之间的线性关系。这种方法能够识别与环境因素同时表现出关联的协同等位基因频率，从而捕获局部适应的多基因特征。为了检测可能处于选择下的位点，检查了RDA加载分布的尾部——代表极端关联——使用偏离平均值2.5个标准差的截止值（相当于双尾p值0.0005）。作为补充，使用LFMM2来测试遗传变异和环境预测因子之间的线性关联，同时考虑可能混淆的群体结构。LFMM2的算法将未观察到的潜在因子纳入最小二乘回归框架，这有助于减少由于隐性群体结构导致的假阳性。使用R包LEA（版本3.10.2）中的snmf函数实现的SNMF来估计潜在群体结构。这生成了个体祖先系数，用于推断祖先群体数量（K）并控制LFMM2分析中的结构。通过生成一个熵准则来评估统计模型对数据的拟合度，使用交叉验证技术为1到13个祖先群体（K）确定祖先系数。从100次重复中，选择具有最低交叉熵得分的K值。接下来，使用最优因子K=4来告知LFMM，以识别等位基因频率是否与任何环境变量相关。通过使用dartR包中的gl.impute函数使用最近邻选项插补缺失基因型数据，增加了关联的统计功效。随后，使用函数lfmm_ridge计算使用岭惩罚的正则化最小二乘估计。使用通过基因组膨胀因子校准的统计检验（函数lfmm_test）评估每个SNP频率与每个环境变量之间的个体关联。使用Benjamini-Hochberg算法应用多重比较校正，FDR阈值为5%。当在两个方法中都识别出时，保留候选SNP位点用于下游分析。

使用R中的gradientForest包进行梯度森林（GF）分析，以模拟沿环境梯度的等位基因频率转换。这种机器学习方法构建回归树组合以识别环境变量和等位基因频率之间的非线性关联，同时排序每个位点的变量重要性。累积重要性度量捕获了每个环境预测因子在SNP集中解释了多少基因组变异。然后使用GF算法描述环境梯度和遗传变量之间的关联。GF使用基于回归树的方法来拟合基因组数据和环境变量之间的响应模型。使用最初通过RDA和LFMM2识别的候选SNP作为响应变量，在预测变量上模拟适应性遗传变异的转换。机器学习算法沿每个环境梯度在多个分割值处划分等位基因频率，并计算每个分割的等位基因频率变化。分割重要性沿环境梯度累积求和，并在等位基因间聚合以构建非线性转换函数，以识别受预测变量显著影响的位点。该分析针对环境预测变量运行超过500棵回归树，所有其他参数设置为默认设置。为了评估假阳性的可能性，我们作为对照进行了一次GF分析，使用一组139个SNP，这些SNP随机选自既未在LFMM2也未在RDA中识别为异常值的位点池。这些SNP与候选位点频率匹配，并在可能的情况下从相同染色体采样以控制基因组背景。

结果

中性群体结构

Durvillaea amatheiae表现出非常强的遗传结构，在主成分分析图中聚类为六个良好分离的组（图2b）。分子变异分析表明群体间存在高度显著的基因组变异（F_ST=0.1925, p<0.001, df_groups=5, df_error=148, 1000次置换）（表S1），但IBD分析（Mantel检验）返回了不显著的结果（p=0.264, y=0.22 + 2.1e?07x, r²=0.0131；图S5）。所有六个群体都表现出低等位基因丰富度和近交（正F_IS值；表1），这也表现为杂合子缺陷（H_E > H_O）和整体杂合度匮乏（低H_O值）。私有等位基因计数（A_P）在两个范围边缘群体Aragunnu（882）和Tathra（728）最低，其次是Merimbula（1900）和Cape Conran（2252），最高计数在Mallacoota（4330）和Green Cape（2875）。

pairwise F_ST分析（图2a）和显著性检验（使用置换检验（1000次置换的p值，df=n?2））揭示了每个群体对之间存在高度显著差异（Bonferroni校正后p<0.05），其中最大的分化出现在Mallacoota与两个范围边缘群体之间：Aragunnu和Tathra（F_ST=0.48和0.43）。第二大的分化出现在Mallacoota和Merimbula之间（F_ST=0.42）。较低的分化出现在Tathra和Merimbula之间（F_ST=0.20），Merimbula和Green Cape之间（F_ST=0.16），以及Aragunnu和Tathra之间（F_ST=0.10），最小的分化出现在Cape Conran和Green Cape之间（F_ST=0.07）。后两个地点共享了相当大比例的遗传物质（与其他群体对相比），而Mallacoota似乎高度隔离（图2a）。

祖先分析显示最可能的祖先群体数量是四个（交叉熵膝或平台=4, K=4；见图S7）。观察到四个祖先群体以不同比例出现在每个群体中（图2c）。群体1（黄色）构成了Aragunnu和Tathra的大部分祖先比例，但后者来自其他群体的混合增加。Green Cape和Cape Conran都表现出类似的混合比例，祖先以群体3（红色）为主。Mallacoota和Merimbula分别由单独的祖先群体主导，分别是群体2（蓝色）和群体4（紫色）。

温度相关适应的特征

冗余分析和LFMM2独立地识别了与两个环境预测因子都表现出显著GEAs的SNP。在RDA中，置换方差分析证实约束排序模型整体高度显著（p=0.001, F=6.049, df_model=2, df_residual=151, 999次置换；表S2）。此外，两个约束典型轴RDA1和RDA2都单独显著（RDA1: F=7.885, p=0.003, df=1, df_residual=151; RDA2: F=4.214, p=0.003, df=1, df_residual=151–999次置换；表S3）。特征值分解表明RDA1和RDA2分别解释了65.17%和34.83%的约束遗传变异（图S8）。使用偏离SNP加载平均值±2.5个标准差的截止阈值，RDA识别了332个候选SNP。LFMM2分析也揭示了SNP变异与两个环境预测因子之间的显著关联。基于校正模型，LFMM2识别了607个SNP为显著异常值。通过直方图检查每个环境变量的p值分布，以确认零假设下的均匀性和备择假设下低p值的富集（图S9和S10）。总共有139个SNP在RDA和LFMM2方法的结果中重叠。这些位点——在个体SNP水平定义，不是基因——被保留用于使用GF建模进行进一步分析（表2）。为了评估方法间SNP重叠是否大于随机机会预期，我们使用base R中的phyper函数进行了超几何检验。检验产生了高度显著的结果（p<0.000），表明观察到的方