湿实验室方案至关重要：DNA提取与文库构建方法的选择显著影响古口腔微生物组的回收效果

本研究系统比较了两种古DNA（aDNA）提取方法（QG法与PB法）和两种文库构建方法（DSL与SSL）对考古牙结石中古口腔微生物组回收的影响。结果表明，实验方案的选择会显著影响aDNA的片段长度、GC含量、克隆性、内源性含量、脱氨损伤模式及微生物组成，且其效果受样本保存状态的调控。研究强调，在古宏基因组学研究中需审慎选择湿实验室方案以提高结果的可比性和可重复性。

　　

1 引言

考古牙结石中的古DNA（aDNA）分析为了解古代健康、人口统计和生活方式提供了丰富信息。然而，古代宏基因组学的工作流程仍在不断演变，引发了关于结果可重复性的担忧。尽管有来自现代群体的证据表明DNA提取和文库制备方法会影响微生物组的回收，但很少有系统研究探讨这些方法如何影响古代口腔微生物组的恢复，这使我们对湿实验室方案如何影响牙结石中aDNA回收的理解存在空白。

2 材料与方法

2.1 出口许可与样本来源

本研究经阿德莱德大学人类研究伦理委员会（H-2012-108）评估并批准。分析的古代牙结石样本来自匈牙利（n=3）和尼日尔（n=3）。所有三个尼日尔样本均从撒哈拉沙漠中的Gobero考古遗址回收，该遗址在距今9500年和4500年前被密集占据。这些个体估计生活在大约公元前5200年至公元前2200年之间，跨越了该地区的中全新世至晚全新世时期。一个匈牙利样本采集自Balatonszárszó–Kis-erdei-d?l?的个体，其年代可追溯至新石器时代早期（约公元前5500–5000年），而另外两个牙结石样本采集自Alsonyek的个体，其年代在公元前5000年至4500年之间。

2.2 样本收集程序

使用无菌牙科刮匙从牙齿样本表面取下龈上牙结石沉积物。为避免釉质损伤，施加与牙齿表面平行的轻柔压力。收集的碎片随后储存在无菌密封自封袋中，运输至阿德莱德大学澳大利亚古代DNA中心（ACAD）的古代DNA设施。这些样本被选中进行比较是因为，尽管两组样本的年代相似（约7000年前），但它们的aDNA保存状态预期会有所不同。这一预期源于先前的发现，即从匈牙利样本的相应牙齿中回收到了真实的古代人类DNA，但未从尼日尔的样本中回收到了。

2.3 设施协议

所有样本处理和实验室程序均在阿德莱德大学的专门aDNA设施（ACAD）中进行。该设施配备正压空气，每天用3%次氯酸钠（NaClO）进行消毒，并每晚暴露于紫外线进行消毒。本研究中的所有实验均在隔离的静止空气房间内经过紫外线处理的静止空气罩中进行。所有人员通过专用的单通道房间进入设施，并穿着一次性全身工作服、手套和面罩，遵循严格的设施工作流程。

2.4 样本去污染与消化

为去污染，我们将每个样本分成两个视觉上相等的部分，并对其进行以下处理：将样本置于无菌塑料培养皿中，每面暴露于紫外线辐射15分钟。接下来，将每个样本浸入5%漂白剂中，然后浸入纯水中各5分钟，随后浸入90%乙醇中3分钟以去除任何残留的漂白剂。样本在无菌容器中风干5分钟，然后转移到无菌塑料袋中。在袋内，用钢锤压碎牙结石。然后剪掉塑料袋的一角，将样本倒入无菌的2 mL螺旋盖管中。

2.5 QG方法

使用QG方法制备的粉碎牙结石样本遵循Brotherton等人先前描述的方案，该方案改编自Rohland和Hofreiter的方法。简而言之，通过将粉末化样本添加到含有1.8 mL 0.5 M乙二胺四乙酸（EDTA）、100 μL 10%十二烷基硫酸钠（SDS）和20 μL 20 mg/mL蛋白酶K的无菌管中，在55°C下旋转20小时进行脱钙。消化一天后，将释放的DNA与3 mL改良QG缓冲液（Qiagen）中的硅胶悬浮液结合，离心沉淀，并用80%乙醇洗涤两次。接下来，将清洁干燥的硅胶重新悬浮在100 μL Tris-HCl缓冲液中两次以洗脱DNA。洗脱的DNA分装并在-20°C下储存直至扩增。值得注意的是，用于制备试剂的所有无DNA和RNA认证水（Invitrogen超纯蒸馏水）都是新鲜开封的或在分装后冷冻以防止微生物生长和污染。

2.6 PB方法

使用PB方法制备的样本采用与QG方法相同的方案进行去污染和消化。消化后，PB方法的步骤遵循Moore和Weyrich描述的方案，该方案改编自Dabney等人。简而言之，PB结合缓冲液由12.2 mL Qiagen PB缓冲液、7 μL Tween-20和378 μL NaOAc（3 M）组成。该方案所需的额外试剂包括100 μL硅胶溶液、200 μL TLE缓冲液和1.8 mL分子级水中的80%乙醇。将每个压碎的牙结石样本在静止空气罩中加入到15 mL锥形管中。接下来，加入12.6 mL PB结合缓冲液和100 μL硅胶。所有样本在微量离心机中以19,500 g离心3分钟。接下来，将样本的上清液转移至无菌15 mL管中。将管置于旋转混合器上，在室温下混合1小时。此步骤之后，以4400 g离心5分钟。弃去上清液，向每个沉淀中加入900 μL 80%乙醇，使用长程移液器吸头通过上下吹打进行重悬。将重悬的溶液转移至DNA LoBind管中。所有样本随后在微量离心机中以14,000 rpm离心1分钟。向每个干燥的沉淀中加入100 μL TLE，并使用涡旋振荡器重悬。可以将30 μL DNA提取物移入管中。

2.7 DSL方法

使用DSL方法构建的DNA文库按照Meyer和Kircher描述的方案制备，并进行了微小修改。大约20 μL来自PB和QG方法的DNA提取物用T4多核苷酸激酶和T4 DNA聚合酶（New England Biolabs）在25°C下处理15分钟。使用MinElute反应纯化试剂盒（Qiagen）清洁反应。使用T4 DNA连接酶（Fermentas）在22°C下将带有5 bp独特条形码的截短Illumina接头序列连接到双链DNA分子上60分钟。使用Qiagen MiniElute纯化试剂盒去除多余的连接酶；使用Bst DNA聚合酶（New England Biolabs）在37°C下填充接头序列30分钟，然后在80°C下使聚合酶变性10分钟。将所得反应物作为模板，在五个独立的PCR反应（12 μL无DNA dH₂O，2.5 μL 10×缓冲液，2.5 μL 25 mM MgCl₂，0.25 μL HiFi DNA聚合酶和5 μL Bst反应混合物）中使用以下条件进行：94°C 12分钟；94°C 30秒，60°C 30秒，72°C 45秒，共13个循环；最后72°C 10分钟。将PCR反应合并，并使用Ampure PCR纯化（Agencourt）进行清洁。然后使用GAII索引引物在相同条件下重新扩增文库，以包含每个样本独特的单个P7（3'）索引序列，重新合并，并再次使用Ampure进行清洁。

2.8 SSL方法

单链DNA文库按照Gansauge等人的方案制备，并进行了微小修改。首先，使用FastAP热敏碱性磷酸酶（ThermoFisher）对来自PB和QG方法的每个样本的20 μL DNA提取物进行去磷酸化处理，以去除5'和3'磷酸基团。每个反应在37°C下孵育30分钟，然后在75°C下热灭活10分钟。去磷酸化后，通过加热至95°C持续2分钟使DNA变性，并立即在冰上速冷以保存单链片段。使用T4 DNA连接酶（ThermoFisher）将生物素化的接头连接到单链DNA的3'端。将连接产物固定在Dynabeads C1链霉亲和素包被的磁珠（Invitrogen）上，并彻底洗涤以去除未结合的DNA。然后退火与生物素化接头互补的引物，并使用Klenow片段（3'→5'）（ThermoFisher）合成第二链。接下来，使用T4 DNA连接酶将平末端双链接头连接到新合成的链上，随后进行第二次珠固定和洗涤步骤。接头连接后，将文库热变性以从珠上释放文库分子，并保留上清液用于扩增。

预扩增的文库使用含有样本特异性条形码的P5和P7引物进行索引PCR。进行13个循环的PCR，根据定量确定适合每个样本的循环数以实现平衡扩增。

2.9 文库合并与测序

所有DSL和SSL索引文库均使用AMPure XP磁珠（Beckman Coulter）进行纯化，使用Agilent TapeStation进行定量，并以等摩尔浓度合并。最终库使用qPCR（Applied Biosystems）进行定量。文库在Illumina NextSeq平台（Illumina，USA）上使用2×150 bp配置进行测序。

2.10 生物信息学处理与分析

使用Illumina bcl2fastq（v1.8.4）软件将测序数据转换为FASTQ文件格式。然后使用AdapterRemoval v2对原始FASTQ文件进行解复用、修剪和合并，基于独特的P5/P7条形码接头（–minlength 25, –minquality 25, –trimns, –trimqualities, –collapse），产生分析就绪的读段。使用Seqkit（v2.6.1）计算每个样本的读段数量、平均读段长度和GC含量。使用内部脚本计算克隆性（即独特读段的百分比）。我们使用R（v.4.1.1）中的aov()函数进行双向ANOVA来测试这些评估的统计显著性。这些分析的脚本可以在GitHub页面上找到。

使用MEGAN比对工具（MALTn；v.0.3.8）中的核苷酸比对选项对分析就绪的读段进行 taxonomic 分箱。将读段与一个内部RefSeq（'RefSeq'）数据库进行比对，该数据库先前已发布，包含47,696个古菌和细菌基因组组装（ scaffold、染色体和完整水平）。使用MEGAN6（v 6.11.1）程序中包含的blast2rma脚本将生成的基于比对的blast-text文件转换为RMA文件，并使用以下最后共同祖先（LCA）参数：Weighted-LCA = 80%，最小比特值= 42，最小E值= 0.01，最小支持百分比= 0.1。将生成的RMA6文件导入MEGAN6。我们使用MEGAN6 CE也导出了物种水平的BIOM表格，使用taxonNameToCount汇总选项，并排除了未分配的读段。

2.11 保存状态评估

为了评估样本中内源性含量的保存情况，我们使用了SourceTracker2（v2.0.1）。该分析的比较数据来自研究现代牙结石、现代牙菌斑、土壤和皮肤的研究，以及本研究中的EBC。选择比较源的原因如下：它们使用Illumina鸟枪法测序策略生成，并且至少产生1000个映射到我们NCBI-RefSeq数据库的读段。将源的物种水平表格与本研究的样本合并。我们使用以下参数运行SourceTracker2：alpha2 = 1.0，sink_rarefaction_depth = 5573，source_rarefaction_depth = 1489。sink的稀疏深度设置为5573个读段，基于样本18393的测序深度，该样本在至少有1000个读段的牙结石样本中读段数最低。源样本的稀疏深度设置为1489个读段，基于样本19028，它代表了至少有1000个读段的参考源数据集中的最低测序深度。计数少于1000的样本被排除在此分析之外。使用R（v.4.1.1）可视化此分析的结果。

使用MapDamage（v2.0）通过将分析就绪的读段映射到Anaerolineaceae bacterium oral taxon 439 (ASM171754v1)、Methanobrevibacter oralis DSM 7256 (ASM163927v1) 和 Olsenella sp. taxon 807 (strain F0089) (ASM118951v2) 来评估单个微生物基因组的死后损伤。选择这些微生物是因为它们在尼日尔和匈牙利数据集的许多样本中含量丰富，并且在其他古代牙结石研究中有充分记载。对于单链文库，我们在mapDamage中加入了单链选项，以确保正确检测读段两端的C→T替换。

除了mapDamage，我们还使用了ChangePoint。与MapDamage不同，ChangePoint不依赖于将读段映射到参考基因组，而是使用似然比检验来评估DNA片段末端是否在5'端富集胸腺嘧啶（T），在3'端富集腺嘌呤（A）。因为它以无参考方式操作，ChangePoint假设当在古代宏基因组数据集中观察到这种核苷酸模式时，反映了一个以真实aDNA为主的群落，脱氨信号分散在多个分类群中。如果样本的调整后p值小于0.05，则认为它们代表了真实的古代宏基因组。为了优化ChangePoint的性能，我们将每个样本二次采样至100,000个读段。

2.12 微生物组成分析

样本18393_Niger_PB_SSL因其分配计数低（400）而从所有组成分析中移除。然后我们对剩余的23个样本进行下游分析。虽然合并尼日尔和匈牙利数据集可以增加统计效力，但这种方法可能会掩盖我们旨在评估的样本内变异性，因为区域水平的微生物组成差异将主导信号。通过按国家分层分析，我们能够更好地分离和解释DNA提取和文库制备方案对关键变量的影响，例如保存状态、DNA片段长度回收、GC含量、克隆性、内源性含量、DNA脱氨以及微生物多样性和组成，所有这些已被证明是评估aDNA和古代微生物组回收的关键指标。

2.13 分类学分析与污染物过滤

为了评估湿实验室方案对微生物组成的影响，我们使用未过滤的数据集进行了分类学分析。使用未过滤的数据集来评估某些DNA提取或文库制备方法是否更容易回收已知的污染物分类群，这可能表明方案特异性的污染易感性偏差。同时，我们使用decontam包（v.1.26.0）根据频率模式识别和去除污染物物种。我们测试了三个decontam污染分数阈值——0.10、0.25和0.50。0.50阈值在先前的研究中经过验证，识别出两个口腔分类群——Ottowia sp. oral taxon 894和Actinomyces israelii——作为污染物，尽管该数据集中的大多数物种是非口腔物种。基于此以及先前的文献，我们继续使用以0.50阈值过滤的数据集。此外，我们移除了在整个数据集中未在至少两个样本中发现且计数少于10的特征，因为这些可能是虚假比对。经过此过滤后，样本18398_Niger_PB_DSL、18398_Niger_QG_DSL和18400_Niger_PB_DS的计数少于1000。因此，它们被排除在过滤数据集的所有分析之外。这个精炼的数据集使我们能够确定在排除可能的污染物后，观察到的微生物多样性和组成模式是否仍然成立。通过比较未过滤和过滤数据集的结果，我们能够评估方案相关趋势的稳健性，并了解 taxonomic 谱如何随着污染物过滤而变化，以及整个宏基因组（内源性和外源性特征）如何受这些方法的影响。

使用Kruskal-Wallis检验评估alpha多样性（观察到的物种）的差异。该检验应用于整体组和成对比较。对于beta多样性，使用QIIME2（v.2024.10）中的DEICODE插件将物种水平表格进行中心对数转换，并将其转换为Aitchison距离矩阵。我们使用QIIME2中的adonis函数执行了置换多元方差分析（PERMANOVA），使用了999次置换，应用以下模型：Aitchison距离矩阵 ~ “ExtractionType” × “LibraryMethod”。我们进一步使用基于物种水平Aitchison距离的主坐标分析（PCoA）探索了beta多样性。所有检验的统计显著性设定为< 0.05。

2.14 微生物丰度与实验室方法之间的关联

我们采用微生物组多变量线性模型关联（MaAsLin2）（v.1.20.0）对从MEGAN导出的物种水平表格进行差异丰度分析。使用了默认的MaAsLin2参数（分类学特征在至少10%的所有样本中普遍存在，最小相对丰度百分比0.01），但校正后的max_significance阈值设置为0.05。元数据中的“ExtractionType”和“LibraryMethod”列作为模型中的“fixed_effects”输入。

3 结果

3.1 片段长度回收

古代DNA数据集的特点是高度碎片化的序列，导致读段长度短。为了评估实验室方案如何影响片段长度回收，我们首先评估了DNA提取和文库制备方法的效果。在尼日尔数据集中，DNA提取和文库制备方法均对片段长度回收没有显著影响。对于DNA提取比较，使用QG方法制备的样本回收的片段（67.60 bp，标准差[SD] = 25.471）比PB方法（61.42 bp，SD = 18.060）长，但差异不显著（F = 0.199, p = 0.667）。对于文库制备的比较，DSL方法回收的片段（69.25 bp，SD = 30.211）比SSL方法（59.77 bp，SD = 5.587）长，但差异不显著（F = 0.469, p = 0.513）。当比较DNA提取和文库制备方法的组合时，两者之间的相互作用不影响片段长度回收（F = 0.008, p = 0.930）。

对于来自匈牙利的样本，QG方法回收的片段（82.72 bp，SD = 30.625）比其PB对应方法（56.13 bp，SD = 8.871）长，并且这些差异是显著的（F = 53.92, p = 0.0000805）。而对于文库制备分析，DSL方法回收的片段（86.57 bp，SD = 26.303）显著长于SSL方法（52.28 bp，SD = 6.597）（F = 89.68, p = 0.0000127）。DNA提取和文库制备方法之间的相互作用也是显著的（F = 31.6, p = 0.000498）。总体而言，DNA提取和文库方法对这些保存较差的样本的片段回收没有显著影响，但对这些保存良好的样本有影响。

3.2 GC含量

我们探讨了DNA提取和文库方案是否影响样本的GC含量。对于尼日尔数据集，QG方法（平均值= 48.23%，SD = 7.110）和PB方法（平均值= 47.90%，SD = 6.483）产生相似的GC含量（F = 0.045, p = 0.838）。然而，DSL方法（平均值= 53.85%，SD = 1.466）和SSL方法（平均值= 42.28%，SD = 3.194）之间的GC含量存在显著差异（F = 52.303, p = 0.0000896）。DNA提取和文库制备方法之间的相互作用不显著（F = 0.003, p = 0.959）。

对于匈牙利数据集，使用QG方法提取的样本的平均GC含量为48.14%（SD = 5.363），使用PB方法的样本为47.95%（SD = 6.503），显示无显著差异（ANOVA: F = 0.009, p = 0.92）。相比之下，文库制备方法之间观察到显著差异：使用DSL方案制备的样本具有更高的GC含量（平均值= 52.70%，SD = 2.796） than SSL（43.39%，SD = 3.351）（ANOVA: F = 23.187, p = 0.001）。两种方法的相互作用在GC含量上没有差异（F = 0.479, p = 0.51）。这些发现表明，无论保存状态如何，文库制备方案的选择都会影响GC含量的回收，而DNA提取方法没有明显影响。

3.3 克隆性

我们探讨了实验室方案是否影响从牙结石样本中回收DNA的克隆性。对于来自尼日尔的样本，DNA提取方法和文库制备方法均未显著影响克隆序列的回收，但观察到轻微差异。QG方法回收的独特读段百分比（平均值= 12.71%，SD = 6.247）低于PB方法（平均值= 27.61%，SD = 31.22），但这些差异不显著（F = 1.121, p = 0.303）。虽然使用DSL方法制备的样本产生的独特读段百分比（平均值= 17.00%，SD = 5.229）低于SSL方法（平均值= 23.31%，SD = 33.103），但这些差异也不显著（F = 0.218, p = 0.653）。两种方案之间的相互作用没有显著的交互作用（F = 1.015, p = 0.343）。

在来自匈牙利的样本中观察到类似的趋势。QG（平均值= 78.12%，SD = 32.936）和PB方法（平均值= 94.22%，SD = 2.315）之间的独特读段百分比存在微小差异（F = 1.331, p = 0.282），DSL（平均值= 92.34%，SD = 5.662）和SSL方法（平均值= 80.00%，SD = 33.500）之间也是如此（F = 0.782, p = 0.402）。方法的组合对回收的独特读段百分比影响也很小（F = 0.540, p = 0.483）。这表明DNA提取和文库制备方法对古代宏基因组文库中观察到的克隆性影响很小。

3.4 微生物源追踪

为了验证牙结石样本中的微生物含量，我们使用了SourceTracker2。SourceTracker2预测了在尼日尔和匈牙利牙结石样本中发现的微生物DNA的潜在来源。对于来自尼日尔的样本，DNA提取方法对口腔DNA（即与现代牙结石或牙菌斑源匹配的DNA序列）的回收影响很小，因为QG（平均值= 7.2%，SD = 0.002）和PB（平均值= 6.4%，SD = 0.002）方法的产量相似（双向ANOVA，F = 0.039, p = 0.850）。文库制备方法的情况则不同。SSL方法在回收口腔DNA（平均值= 14.8%，SD = 0.002）方面比DSL方法（平均值= 0.05%，SD = 0.0002）有效得多（双向ANOVA，F = 12.847, p = 0.009）。值得注意的是，PB+SSL方法对样本18393（12.0%）和18398（25.20%）产生的口腔DNA最多，而QG+SSL方法对样本18400（15.50%）产生的口腔DNA最多。然而，DNA提取和文库制备之间没有显著的交互作用（双向ANOVA，F = 0.373, p = 0.561）。此外，DSL方法倾向于回收更多归因于污染源的DNA，例如EBC、皮肤和土壤（平均值= 66.20%，SD = 25.70%），而它们的SSL对应方法则较少（平均值= 16.50%，SD = 34.20%）。

与尼日尔数据集类似，DNA提取方法的选择对内源性含量的回收没有显著影响（双向ANOVA，F = 0.027, p = 0.873）。然而，与尼日尔数据集不同，文库制备不显著（双向ANOVA，F = 0.334, p = 0.579），尽管DSL方法对匈牙利样本产生的内源性含量百分比更高（平均值= 46.70%，SD = 0.005） than SSL方法（平均值= 37.4%，SD = 0.015）。值得注意的是，对于匈牙利组，所有样本中归因于污染源的不到1%（DSL平均值= 0.04%，SD = 0.0002；SSL平均值= 0.05%，SD = 0.0003）。总之，这些结果表明，文库制备在从保存较差的样本中回收内源性含量方面起着重要作用，但对于保存良好的样本则不然，并且可以影响污染DNA的回收。

3.5 实验室方案影响样本损伤谱的回收

我们探讨了DNA提取和文库制备方法是否影响使用MapDamage2分析脱氨率的能力。我们将每个样本的分析就绪读段映射到三种微生物：Anaerolineaceae bacterium oral taxon 439 (ASM171754v1)、Methanobrevibacter oralis (strain DSM 7256) (ASM163927v1) 和 Olsenella sp. taxon 807 (strain F0089) (ASM118951v2)。选择这些微生物是因为它们在样本中相对丰度较高。对于尼日尔数据集，使用DSL方法制备的样本没有表现出古代样本预期的脱氨模式。这一结果可能是由于唯一映射读段数量不足。使用DSL方法制备的样本都没有产生超过两个映射到这三种微生物的唯一读段。使用SSL方法制备的样本产生了更多唯一映射读段，尽管基于提取方法的组间存在显著差异。使用QG提取和SSL方法制备的样本产生了更多唯一映射读段（平均值= 1415.44；SD = 1702.23） than 使用PB提取和SSL方法制备的样本（平均值= 116.56；SD = 113.86）（ANOVA, F = 5.217; p = 0.036）。这些结果表明，QG+SSL方法是唯一一种始终回收足够数量读段的方案——这里定义为至少100个唯一映射读段——以便对保存较差的样本进行脱氨模式的初步评估。虽然这个阈值提供了一些aDNA真实性的指示，但它仍然低于通常接受的1000个唯一映射读段的基准，以自信地确定一个微生物是否是真正的古代微生物。

DNA提取和文库制备方法似乎对保存良好的匈牙利样本影响较小。大多数使用DSL方法制备的匈牙利样本在5'和3'端分别表现出预期的A→T和G→C错误掺入。只有使用PB+DSL方法制备的样本18427产生的唯一映射到M. oralis基因组的读段少于100个。与尼日尔数据集的结果不同，DNA提取和文库制备方法在产生更多唯一映射读段方面都没有发挥显著作用。然而，SSL方法确实产生了显著更多的唯一映射到M. oralis基因组的读段（平均值= 6034.83） than DSL方法（平均值= 678.67）（ANOVA: F = 6.284; p = 0.03）。

为了评估匈牙利和尼日尔样本的整体损伤谱，我们使用了ChangePoint，一种无参考的损伤分析器。我们的ChangePoint结果表明，尼日尔的样本支持MapDamage结果，即样本中具有5'端富集T和3'端富集C的DNA片段数量不显著。另一方面，