微生物群中单个细菌细胞的高通量鉴定和定量

编译：微科盟牛魔王，编辑：微科盟小编、江舜尧。

微科盟原创微文，欢迎转发转载，转载须注明来源《微生态》公众号。

导读  

细菌微生物群作为一个由许多个体生物(即细胞)组成的群落发挥作用。为了充分了解细菌微生物群的功能，必须对单个细胞进行鉴定；然而，目前的技术很难做到这一点。在本文中，我们开发了一种方法，即条形码鉴定和量化细菌(BarBIQ)，该方法基于具有单碱基精度的细胞条形码16S rRNA序列，将单个细菌细胞划分为分类群，本文命名为基于细胞的可操作分类单元(cOTU)，并以高通量方式量化微生物群中每个cOTU的细胞数。我们使用BarBIQ对小鼠盲肠微生物群进行测量，总共有3.4×105个细菌细胞包含810个cOTU。有趣的是，我们发现盲肠微生物群的位置依赖性整体差异取决于膳食维生素A的缺乏，并且近端位置的cOTUs比远端位置的cOTUs差异更丰富。重要的是，传统的16S rRNA基因扩增子测序方法无法清楚地显示这些位置差异，这些方法量化的是16S rRNA基因，而不是细胞。因此，BarBIQ能够通过识别和量化单个组成细菌来表征微生物群，这是微生物群研究的基础。

论文ID

原名：High-throughput identification and quantification of single bacterial cells in the microbiota

译名：微生物群中单个细菌细胞的高通量鉴定和定量

期刊：Nature Communications

IF：14.919

发表时间：2022.2

通讯作者：Katsuyuki Shiroguchi

通讯作者单位：日本RIKEN生物系统动力学研究中心(BDR)

DOI号：10.1038/s41467-022-28426-1

实验设计

结果

1 基于细胞的高通量方法——BarBIQ

在BarBIQ(图1b)中，我们首先在缓冲液中制备了细菌样本(模拟群落和小鼠盲肠内容物；见下文)，并通过涡旋打破了细菌团块(补充图1和补充注释3)。接下来，我们将细菌样本与包括细胞条形码、引物和DNA扩增试剂的溶液混合，然后使用Bio-Rad Droplet DigitalTM PCR(ddPCR)系统将它们包裹在直径约120微米的液滴中。我们根据泊松分布调整了条形码和细菌的浓度及其比例。对于条形码，约93%的条形码占用的液滴每个都有一个条形码；其余7%的液滴中的多个条形码不会对测量产生显著影响，因为细菌的数量平均变化系数为1.07(=0.93×1+0.07×2)。此外，这一因素对所有细菌都是通用的，因此它们的相对比例不会受到影响。我们发现细菌不影响条形码(83碱基单链DNA)的封装(补充图3a)，这应该遵循泊松分布。对于细菌来说，约90%的细菌液滴每个都有一个细胞；其余10%的液滴中的多个细胞大多是不同的细菌(补充注1)，可以根据它们的16S rRNA序列区分它们(补充注2)。通过对单个细菌的显微成像，我们确认盲肠样本细菌成功地被包裹到遵循泊松分布的液滴中(图1d和补充图2)。为了随后的测序，通过加热使细菌细胞裂解，扩增Barcode和16S rRNA基因(V3-V4区，约450个碱基)，连接测序接头，并且在一步中完成Barcode的扩增和其与16S rRNA基因间的连接(图1c和补充图4)。这一过程在很大程度上避免了不同细菌的嵌合体的产生，这些嵌合体主要通过16S rRNA基因扩增子测序中的批量PCR产生。我们发现，通过显微成像检测到细菌的液滴比例接近于盲肠样本扩增16S rRNA基因的液滴比例(图1e)。我们还通过调整细胞裂解的初始加热时间，证实了包括细胞裂解过程在内的16S rRNA基因的扩增是稳健的(补充图3b)。扩增后，我们破坏了液滴，纯化了文库(连接的扩增片段)，并使用高通量MiSeq测序仪对单个扩增分子的细胞条形码和16S rRNA序列进行了测序。我们分析了条形码的每种序列类型(即细胞)的测序分子(即Reads)，鉴定了每个细胞的16S rRNA序列(称为Bar序列)，并计算了每种16S rRNA序列类型的细胞数量(方法部分，补充图5和补充注释2)。在一次MiSeq运行中，确定了超过105个细胞。值得注意的是，这种分析也适用于其基因组中具有多种16S rRNA序列类型的细菌，因为相同的细胞条形码与来自同一细胞的多个扩增的16S rRNA序列相连，并且根据共存的量化和统计模型，这些相同的细菌16S rRNA序列与偶然共存的16S rRNA序列不同，后者来自液滴中的不同细菌(补充说明2，步骤15)。我们通过使用ddPCR测量的同一样本每单位重量或体积的总细胞数对测序确定的细胞数进行归一化，最终获得了样本中每种16S rRNA序列类型(称为cOTU，如下所述)每单位重量或体积的绝对细胞数。

BarBIQ与传统方法之间的一个本质区别是用于定义微生物群组成的单位和用该单位进行定量(图1a)。在传统方法中，常用的单位之一是操作分类单元(OTU)，它表示一组相似的16S rRNA 序列，这些序列是通过基于从大量样本中检测到的序列的身份进行聚类而获得的。传统方法中另一个广泛使用的单位是扩增子序列变异(ASV)，它用于检测样本中的独特16S rRNA序列。OTU和ASV并不总是代表细菌细胞，因为如上所述，有些细菌在其基因组中具有多种16S rRNA序列类型。然而，BarBIQ从每个带条形码的细菌细胞中识别出16S rRNA序列(即Bar序列)。为了将我们的基于细胞的方法与使用OTU的传统方法区分开来，我们将来自同一细胞的16S rRNA序列(即Bar序列)命名为“基于细胞的操作分类单元(cOTU)”。此外，BarBIQ量化每个cOTU的细胞数量，而传统方法测量扩增的16S rRNA 基因分子的数量(16S rRNA基因丰度)。

图1 BarBIQ及其质量控制。a. BarBIQ的主要概念及其与传统16S rRNA基因-扩增子测序方法的比较。b. BarBIQ的示意图。样品悬浮在溶液中后，进行涡旋以打破细菌团。细胞条形码，含有随机碱基和扩增引物位点的DNA分子；引物，用于扩增16S rRNA基因和细胞条形码的DNA引物，用于连接两种扩增产物，并用于连接测序接头；试剂，用于DNA扩增的试剂。文库生成、纯化和测序的原理图详见补充图4，数据处理的细节见补充图5。c. 液滴中的文库生成原理图。滴液中细菌基因组中的细胞条形码和16S rRNA基因(V3-V4区，约450个碱基)最初由含有测序接头和接头序列的引物扩增。随后，扩增的条形码通过接头序列与扩增的16S rRNA基因连接。d. 通过显微镜成像(补充图2)观察到的液滴中细菌数量的分布(条形)与根据泊松分布计算的理论分布(点)之间的比较。e. 通过ddPCR扩增细菌中16S rRNA基因的液滴比例(补充说明1)与显微镜成像观察到细菌的液滴比例(补充图2)；两个实验中的液滴都是用同一个粪便样本产生的。用相同的粪便样本和相同的稀释系数产生。数据以平均值±SD表示(扩增的n = 4，成像的n = 3)。 P 值通过Kruskal-Wallis秩和检验计算。(d)和(e)的源数据作为源数据文件提供。 

2 使用人类肠道细菌菌株检测BarBIQ的功效 

本研究证明了BarBIQ对一个模拟群落的功效，该群落包含10种已知丰度的可培养的人类肠道细菌菌株(方法部分，补充表1)；这些菌株已被用作模型群落，代表了健康人体肠道微生物群中4种最重要的细菌门类(放线菌门(1株)、拟杆菌门(3株)、厚壁菌门(4株)和变形菌门(2株))。包括革兰氏阳性(5株)和革兰氏阴性(5株)细菌。我们从模拟群落中发现了16个条形码序列(图2a 和补充数据1)。这些确定的条形码序列在三个技术重复中完全一致(补充数据1)。所有16个Bar序列与我们从10种培养菌株中通过Sanger测序鉴定的16S rRNA序列(San序列)之一相同(图2a，b，方法部分，补充数据2)。另一方面，还有另外13个San序列具有上述16个序列中的一个或两个替换(San-Bar匹配序列)(图2b和补充数据2)。我们发现这13个San序列可以通过Sanger测序扩增中的PCR错误和我们使用的聚合酶的错误率得到合理解释(补充图6a)，这表明San序列和Bar序列是完全一致的。然后，我们根据细胞条形码(方法部分，补充说明2)从16个Bar序列中鉴定出10个cOTUs，每个对应于10个菌株中的一个，这也表明两对Bar序列在一个碱基上各不相同并且来自同一个细菌细胞(图2a)。因此，我们得出结论，BarBIQ在模拟群落中识别出正确的菌株数量(cOTU)，并且对16S rRNA序列的鉴定具有单碱基准确性和分辨率。 相比之下，我们通过使用传统方法，基于ASV的DADA219 的分析和使用基于OTU的Mothur的分析(方法部分)分析了相同的模拟群落，并确定了28个ASVs和12个OTUs(补充数据3和4)。ASVs和OTUs的数量大于模拟群落中的菌株数量，表明至少有一些ASV和OTU不代表菌株。此外，我们还比较了ASVs和OTUs的代表性序列(OTU-RepSeqs)与10个菌株的San序列和模拟群落的Bar序列的相似性。对于ASVs，虽然15个ASV与San-Bar匹配序列之一相同，但其他13个ASV与任何San 序列都不相同(图2b)。我们发现13个不同的ASV的丰度变化大于估计的采样噪声，而15个相同的ASV很好地遵循了估计的采样噪声(补充图 6b)，并且13个不同的ASV的丰度与15个相同的ASV的丰度有很大不相同(补充图6b和补充数据5)，表明13个不同的ASV不是来自10个菌株的正确16S rRNA序列。对于OTU来说，12个OTU-RepSeqs中只有2个与来自不同菌株的San-Bar匹配序列之一相同，其他10个OTU-RepSeqs与任何San序列都不相同(图2b和补充图6c)，这意味着基于OTU的分析没有从10株模拟群落中的至少8株中检测到正确的16S rRNA序列。重要的是，与任何Bar序列不相同的13个San序列不匹配任何ASV或OTU-RepSeqs(图2b)，这支持了我们对这13个San序列包含扩增错误的解释。因此，与我们在此使用的传统方法相比，BarBIQ对模拟群落的16S rRNA序列识别准确性最高。 通过BarBIQ，我们随后测量了模拟群落中每个cOTU( [C] BarBIQ )的每单位体积的细胞丰度(细胞数量)(补充数据6)，并将 [C] BarBIQ 与每单位体积的细胞丰度进行了比较，后者通过每个菌株的显微镜成像( [C] Microscope )测量的丰度的稀释因子计算得出(图2c，方法部分，补充图7和补充说明4)。首先， [C] BarBIQ  的3个技术重复之间的可重复性很高(标准偏差/平均值：0.01~0.25，N=3；图2c)。其次， [C] BarBIQ 和 [C] Microscope 之间的Pearson's r为0.98， [C] BarBIQ /[C] Microscope 平均比率为0.92(95%置信区间：0.68-1.25；图2c)。因此，BarBIQ以高重现性准确测量了模拟群落中每种细菌菌株的细胞丰度。 然后，我们将15个San序列相同的ASV( [C] ASV )的16S rRNA基因丰度与 [C] Microscope (具有多个16S rRNA序列的菌株与多个相同的ASV进行比较)进行了比例归一化。正如预期的那样， 成比例的 [C] ASV 与 成比例 [C] Microscope 没有很好的相关性(图2d)。为了理解这种差异，我们从 [C] ASV ( [C] ASV-estimated )估计了细胞丰度；10株中有6株的基因组中有16S rRNA基因拷贝数，这些基因组在rrnDB数据库中注册，我们总结了与同一菌株的San序列相同的多个ASV的丰度，并使用其16S rRNA基因拷贝数计算了由 [C] ASV 估计的比例 [C] ASV-estimated (补充图8)。成比例的 [C] ASV-estimated 与成比例的 [C] Microscope (0.97)之间的Pearson相关系数相较于成比例的 [C] ASV 与成比例的 [C] Microscope (0.91)之间的Pearson相关系数有所增加，这与传统方法测量基因拷贝数而非细胞数的设计原则一致。

图2 BarBIQ对模拟群落的功效以及与传统方法的比较。a. 通过Sanger测序和BarBIQ鉴定的16S rRNA序列的比较。Edit distance、Levenshtein distance，定义为替换、插入和缺失的最小数目；San序列，Sanger测序鉴定的16S rRNA序列；ATCC/JCM/DSM <number>,菌株ID；A、B或C，每个菌株的San序列；Bar序列-MK-<number>，BarBIQ-鉴定序列(Bar序列)；cOTU-MK-<number>，基于细胞的可操作分类单元(cOTUs)；红星，有一个碱基差异的Bar序列。b. 从10个菌株中每个菌株鉴定的San序列、Bar序列、扩增子序列变异(ASVs)和从模拟群体中鉴定的操作分类单位的代表序列(OTURepSeqs)之间的比较(文氏图)。圆圈，每种方法的总序列；圆圈中的数字，由给定方法检测到的独特或相同序列的数量；括号中的数字，由给定方法检测到的序列总数。c. 使用BarBIQ([C]BarBIQ)(补充数据6)和显微镜成像([C]Microscope)(补充表1)比较模拟群落中10株的单位体积的绝对细胞丰度。数据以平均值±SD表示([C]BarBIQ为n=3，[C]Microscope为n=5)。蓝色细线，考虑到[C]BarBIQ和[C]Microscope的标准误差，用固定斜率为1的对数比例拟合线(截距：-0.035)，表示平均比率，表明[C]BarBIQ/[C]Microscope的平均比为0.92；灰色粗线，拟合线的95%置信区间，表明[C]BarBIQ/[C]Microscope的95%置信区间为0.68~1.25。d. 由Conv_ASV(比例[C]ASV)测量的模拟群落中15种ASV的比例丰度(补充数据5)与显微镜成像测量的比例[C]Microscope进行比较。数据以平均值±SD表示(比例[C]ASV为n=2，比例[C]Microscope为n=5)。显示了具有常见检测序列的菌株。与多个相同的ASV比较的菌株用颜色表示。通过与c相同的拟合(截距：-0.28)，成比例的[C]ASV和成比例的[C]Microscope的平均比率为0.52，其95%置信区间为0.28-0.9。

3 小鼠盲肠微生物群的测量条件

作为BarBIQ在生物样品中的应用，我们主要根据饮食调查了C57BL/6J雄性小鼠盲肠中两个位置(远端和近端，图1b)

的微生物群落。分析了在以下三种条件下饲养的小鼠(方法部分)：(i)来自日本CLEA 的三只6周龄的小鼠(CEa、CEb和CEc)，由营养均衡的饮食CE-2维持他们整个生命周期(CE2-营养组)；(ii)从日本SLC公司购买的四只8周龄小鼠(VSa、VSb、VSc和VSd)，并通过喂食成分明确的营养均衡饮食再维持5周(补充表2)(VA-充足组) (图3a)；(iii)四只小鼠(VDa、VDb、VDc和VDd)与(ii)中的小鼠经历相同，只是在最后 3周的饮食(补充表2)中不包含维生素A(VA-缺乏组)(图3a)。研究膳食维生素A对肠道微生物群的影响很重要，因为细菌对维生素A的吸收和储存至关重要，而维生素A缺乏会影响人体的视力、生长和免疫功能，从而导致公共健康问题。对于BarBIQ测量，我们首先通过过滤从每只小鼠的每个位置的盲肠内容物中分离细胞和细胞外DNA(ecDNA)，然后测量过滤器残留物(cellsample)和流通物(ecDNA-sample)(方法部分，补充图9和补充说明5)，因为正如最近报道的那样，细胞外细菌DNA可能会影响肠道微生物群的定量。我们注意到，即使需要纯化和/或浓缩，这种过滤程序也可用于不同类型的样品，例如口腔和皮肤样品。我们在小鼠CEa的每个位置对细胞和ecDNA样本进行了3次重复，并对所有其他位置进行了一次测量；每次测量需要<1 mg含量。 

4 准确鉴定小鼠盲肠cOTU的16S rRNA序列 

对于所有细胞样本，我们总共计数了3.4×10 5 个细菌细胞并鉴定出810个含有954条Bar序列的cOTUs(补充数据1)。此外，我们还从ecDNA样本中鉴定了独特的50个Bar序列(补充数据1)；我们没有将这些Bar序列定义为cOTUs，因为ecDNA不代表细胞。重要的是，383个已识别的Bar序列(1004个中的38%(954+50))未在三个广泛使用的公共数据库(GreenGenes46、Ribosomal Database Project47和Silva48)中注册，并且所有Bar序列与这三个数据库中最接近的16S rRNA序列具有大于93%的同一性(图3b，方法部分，补充数据1)。为了确认数据库未注册的Bar序列的准确性，我们从小鼠VDd(VDdprox)近端位置的细胞样本中随机选择单个细菌细胞，对其16S rRNA基因的PCR扩增子进行Sanger测序，并确定了34个独特的San序列(方法部分，补充数据2)。34个San序列中有17个(6个未在数据库中注册)与从同一样本VDdprox中鉴定的Bar序列之一相同(图3c)；其余17个San序列在Bar序列中具有1个、2个或3个替换或1个缺失，这可以通过Sanger测序的模拟PCR错误和我们使用的聚合酶的错误率合理解释(补充图10a)。这些结果表明所识别的Bar序列是准确的，上述模拟群落实验也证明了这一点。 相比之下，我们还分别使用基于ASV和OTU的分析测量了来自VA-充足组和VA-缺乏组(VA组)的16个细胞样本(4只小鼠×2个位置×2个饮食条件)(补充数据3、4、5和7)。为了研究此处基于ASV和基于OTU的分析的序列识别的准确性，我们再次使用上述细胞样本VDdprox的结果。对于ASV，16个San序列与其中一个ASV相同，其他18个不相同的San序列通过Sanger测序的模拟PCR错误得到合理解释(图3c和补充图10b)。然而，对于所有16个细胞样本，4%的ASV(21个ASV)与数据库中最接近的注册16S rRNA序列具有较低(<70%)的同一性(7%ASV具有<93%同一性)(图3b)，而数据库中几乎所有注册的16S rRNA序列对之间的同一性>70％(补充图10c)，表明这些低同一性ASV不是16S rRNA序列，可能是PCR错误的序列，例如嵌合体。对于OTU，虽然13个San序列与其中一个OTU-RepSeqs相同，但其他21个不相同的San序列与OTU-RepSeqs有很大差异，这不能用全球Sanger测序的PCR错误来解释(图3c和补充图10d)，表明基于OTU分析的16S rRNA序列识别准确度不高。重要的是，与任何Bar序列不同的17个San序列不匹配任何ASV或OTU-RepSeqs(图3c)，这支持了我们对这17个San序列包含扩增错误的解释。总的来说，对于模拟群落和小鼠盲肠样本，BarBIQ显示出最高的16rRNA序列识别准确度，这使BarBIQ能够识别未知的16S rRNA序列而无需任何预先确定的序列。

图3 膳食维生素A的实验设计和小鼠粪便微生物群的16S rRNA序列鉴定。a. 膳食维生素A实验设计示意图(方法)。b. Bar序列和ASVs的序列同一性概况；每个Bar序列或ASV与其在三个公共数据库 (GreenGenes、Ribosomal Database Project和Silva)中最接近的16S rRNA序列之间的同一性。原始数据见补充数据1和3。c. 小鼠VDd近端位置细胞样本中San序列、Bar序列、ASVs和OTU-RepSeqs之间的比较(维恩图)。

5 小鼠盲肠中cOTUs的高度可重复的细胞丰度量化

对于每个细胞样本中每个鉴定的cOTU，我们通过BarBIQ(补充数据8)量化每单位重量的绝对细胞丰度，并使用比率中位数方法通过样本之间的归一化计算相对细胞丰度(方法部分)。我们发现每单位重量的总绝对细胞丰度在小鼠之间存在差异(最大倍数变化：2.6(远端)和3.2(近端))，并且远端位置的丰度始终高于近端位置(1.1-3.4倍)(补充图11)。对于每个cOTU，我们通过对同一样品的技术重复，证实了细胞相对丰度和绝对丰度的结果都具有高度可重复性(图4a和补充图12a、b)，并且量化的噪声主要来自模拟采样(补充图12c-e)。 为了研究BarBIQ与基于ASV和OTU分析的量化差异，我们将cOTUs的细胞丰度比例(除以总细胞数)与ASV和OTUs的16S rRNA基因丰度比例与VA组16个细胞样本中的每一个中通常检测到的16S rRNA序列进行比较；将包含多个Bar序列的cOTU与每个相应的ASV或OTU进行比较。正如预期的那样，差异很大(图4b、c和补充图13a、b)：在所有样本中，在cOTU和ASV之间，25%的比较对差异大于2倍，2%的差异大于10倍；在cOTU和OTU之间，26%的差异大于2倍，5%的差异大于10倍。此外，对于常见的16S rRNA序列，ASV和OTU的16S rRNA基因丰度也不同(12%大于2倍和2%大于10倍)(补充图13a-c)。为了理解这种差异，我们考虑了模拟群落中描述的基因组中16S rRNA基因的拷贝数(补充图14)，这再次与这些方法的设计原则一致，其中BarBIQ测量细胞丰度以及基于ASV和OTU的分析测量16S rRNA基因拷贝。

图4 BarBIQ的高度可重复性定量以及与传统方法的比较。a. 成对技术重复间的cOTU相对丰度的比较(点)(其他对见补充图12)。橙色线，基于泊松分布和相对丰度归一化的采样噪声理论置信区间(99.9%)(b和c的比例丰度)；青色线，2倍变化；品红色线，10倍变化。蓝点，显示出比采样噪声和2倍变化更大且小于10倍变化差异的cOTU；括号中的比例和数字分别是相应颜色的比例和点的数量；CEa，小鼠；dist和prox，位置；1和2，技术重复。b. 所有16个细胞样本的cOTUs和ASV间的丰度比例比较，格式与a相同。红点表示cOTUs和ASVs的差异大于采样噪声，10倍变化。c. 所有16个细胞样本的cOTUs和OTUs间的比例丰度比较，格式如a和b所示。源数据以Source data文件提供。
6 基于cOTU的α和β多样性

BarBIQ能够基于细胞(生物体)对细菌微生物群进行α和β多样性分析；大约半个世纪前提出的基于生物体定义α和β多样性对于理解细菌生态系统很重要。作为α多样性的一个例子，我们将基于cOTU的分类单元丰富度(cOTU richness)定义为从样本中检测到的一定总数的细胞中观察到的cOTU数量；类群丰富度是一种常见的生物多样性评估。我们发现CE2营养素组和VA组之间的cOTU丰富度差异很大(大约2倍)(图5a)，这可能是由于组间的饮食、维护设施或年龄不同造成的。另一方面，VA-充足组和VA-缺乏组之间以及每组内远端和近端位置之间的cOTU丰富度是一致的(图5a)。 作为β多样性的一个例子，我们根据检测到的cOTUs的相对细胞丰度(方法部分)计算了Bray-Curtis相异性(基于丰度的beta多样性)，它量化了每对样本的总体差异，并对CE2-营养组(补充图15a)和VA组(图5b)进行了主坐标分析(PCoA)。首先，对于CE2-营养组、VA-充足组和VA-缺乏组，远端和近端位置组没有分开。详细地说，来自远端或近端位置的不同小鼠的样本之间的差异(补充图15b中的洋红色符号)高于来自同一只小鼠的不同位置的样本(补充图15b中的蓝色符号)，表明小鼠总体差异大于位置差异。其次，VA-充足组和VA-缺乏组在远端和近端位置上被很好地分开(图5b)。有趣的是，近端位置的VA-充足组和VA-缺乏组之间的差异显著大于远端位置的差异(图5c)。基于cOTU的绝对细胞丰度计算的Bray-Curtis差异也发现了这种现象(补充图15c)。这种位置差异表明，VA缺乏饮食饲喂3周后，小鼠盲肠近端位置的微生物群总体上比远端位置的微生物群受到的影响更大。

图5   小鼠盲肠微生物群落cOTU丰富度与Bray-Curtis差异。a. 每个细胞样本的cOTU丰富度，通过使用R包Vegan中的函数rarefy对6608个细胞进行子采样确定。CE2，CE2-营养组；VA-suf，VA-充足组；VA-def，VA-缺乏组；dist和prox，位置。b,d和f. Bray-Curtis差异的主坐标分析(PCoA)，根据VA组每对细胞样本间cOTUs (b)、ASVs (d)、OTUs (f)的相对细胞丰度计算Bray-Curtis。标签，与(a)相同；灰线，来自同一只小鼠的联系；圆圈，每组的95%置信椭圆。c、e和g，分别是b、d和f中Bray-Curtis相异性的定量比较。远端，分别为VA-sufdist和VA-defdist的所有可能的配对；近端，分别为VA-sufprox和VA-defprox的所有可能的配对。a、c、e和g中的方框代表第25至75个百分位数(四分位数范围)，水平黑线表示中位数，须表示1.5倍的四分位数范围(CE2dist和CE2prox的n=3;VA-sufdist、VA-sufprox、VA-defdist和VA-defprox的n=4；远端和近端n=16)。 P 值通过Kruskal-Wallis秩和检验计算。 

7 盲肠位置之间cOTU的细胞丰度差异 

我们首先比较了每只小鼠远端和近端位置之间每个cOTU的相对细胞丰度(图6a中的一个示例；均在补充图16a中)。结果表明，差异丰富cOTU比例(差异比采样噪声大2倍)在同一小鼠的CE2-营养组(4–13%；22–75cOTU)、VA-充足组(7–16%；15–43个cOTU)以及VA-缺乏组(2–8%；5–18个cOTU)都显著大于技术重复之间的差异(0.2–0.9%；1–5个cOTU)，且VA-缺乏组的小鼠的差异丰富cOTUs比例最低(图6b)。我们还发现了绝对细胞丰度的类似结果(补充图16b，c)。这些结果表明，每只小鼠的两个位置之间的细菌组成略有不同，这与上面的Bray-Curtis差异一致(图5b和补充图15b)。然而，同一只小鼠中不同位置之间的大多数差异丰度cOTUs在小鼠之间并不一致(补充图17)，这表明应该研究每个cOTU以准确表征微生物群。

图6 每只小鼠位置依赖性相对细胞丰度的比较。a. 小鼠VSa远端和近端cOTUs相对细胞丰度的比较，如图4a所示(补充图16a为其他小鼠)。b. 每只小鼠中位置相关差异丰度cOTUs的比例(差异大于采样噪声，是采样噪声的2倍)。技术上，所有对来自小鼠CEa远端和近端位置的三个技术重复；CE2，CE2-营养组小鼠；VA-suf，VA-充足组小鼠；VA-def，VA-缺乏组小鼠。方框表示第25至75个百分位(四分位数范围)，水平黑线表示中位数，须表示1.5倍的四分位数范围(技术重复n= 6，CE2组n= 3, VA-suf和VA-def组n = 4)。 P 值采用Kruskal Wallis秩和检验计算。 

8 有关饮食维生素A缺乏的cOTUs细胞丰度差异 

为了研究膳食维生素A缺乏对每个cOTUs的影响，我们使用远端和近端位置的cOTUs的相对细胞丰度，在VA充足组和VA缺乏组之间进行了丰度差异分析。在去除低丰度噪声cOTUs后，比较了远端位置的153个cOTUs和近端位置的150个cOTUs(方法部分，补充图18a、b)。我们发现VA-充足组和VA-缺乏组的来自5个属( Acetatifactor、Barnesiella、Lactobacillus、Marvinbryantia 和 Romboutsia )和2个门的8个cOTUs(由每个cOTU的Bar序列预测；方法部分)在近端位置存在显著差异(FDR<0.05和倍数变化>2)，而只有一个(cOTU-CM-2002)(上述8个中)在远端位置显示出显著变化(图7a，b和补充图18a，b)。这一结果与一项研究一致，该研究表明，当小鼠中定植人类肠道衍生细菌群落(92株)且被喂食维生素A缺乏饮食3周时，只有少数菌株显著增加或减少。VA组中所有16个细胞样本中8个cOTU的最大组成分别在0.4%(COTU-CM-0025)到6.5%(cOTU-CM-2002)的范围内(补充数据8)。在VA-缺乏组中，cOTU-CM-2002在近端位置表现出最大的增长(倍数变化：320-1100(标准误差的log2倍数变化：9.2 ± 0.9)和远端位置140-320(7.7 ±0.6))(图7b和补充图18a，b)。该cOTU来自 Barnesiella 属；然而，在 Barnesiella 属中分析的13个cOTU中有11个没有表现出显著差异。同样，在 Acetatifactor 、 Lactobacillus 和 Marvinbryantia 属中，一些cOTU表现出显著差异，而另一些则没有(图7a)。这些结果表明，维生素A缺乏在cOTU水平上塑造了微生物群，并且近端位置的cOTU水平差异大于远端位置。 我们再次在维生素A饮食实验中，对BarBIQ与传统方法进行比较，分别利用每个位置的ASVs和OTUs的相对16S rRNA基因丰度进行VA-充足组和VA-缺乏组之间的丰度差异分析；(近端位置：236个ASVs和138个OTUs；远端位置223个ASVs或131个OTUs)(方法部分；补充图18c-f)。ASV或OTU数量在VA-充足组和VA-缺乏组之间显著不同(FDR<0.05和倍数变化>2)，但在远端位置(6个ASVs或6个OTUs)和近端位置(8个ASVs或9OTUs)相似，这与cOTUs的结果不一致。即使我们改变了定义显著性的FDR和倍数变化的阈值，cOTUs和ASVs或OTUs之间的这种差异也是稳定的(图7c)。这些结果再次表明，细胞丰度与16S rRNA基因丰度之间的不一致性，对于理解这种生物学现象是不可忽视的。重要的是，通过基于ASV(ASV28和ASV31)和基于OTU的分析(OTU029和OTU032)将由BarBIQ(如上所述)确定的最大差异丰度cOTU(cOTU-CM2002)确定为两个独立的单元，这也表明BarBIQ基于细胞的鉴定与传统方法基于基因的鉴定的本质区别。此外，OTU032的OTU-RepSeq与ASV31有一个碱基差异，这是基于ASV分析与基于OTU分析不一致的一个示例。值得注意的是，BarBIQ检测到远端位置的6个差异丰富ASVs中的5个和近端位置的8个ASVs中的7个，而BarBIQ检测到远端位置的6个OTUs中的4个和近端位置的9个OTUs中的5个(补充数据1、3和4)。

图7 膳食维生素A缺乏导致cOTUs细胞丰度差异。a. 检测到的cOTU和分类(由RDP分类器预测，参见方法部分)。VA-充足组和VA-缺乏组在远端或近端至少有一个差异丰富cOTU的属被显示(其他在补充数据8中)。FDR<0.05和倍数变化> 2。b. 通过DESeq2估计远端(Dist)或近端(Prox)位置的差异丰富cOTUs(FDR <0.05，fold change >2)。COTU-CM-0025在远端位置丰度太低，所以没有进行差异分析(方法部分)。c. 近端(Nprox)和远端(Ndist)差异丰富cOTUs、ASV或OTUs数量之间的比值，作为FDR(左)和倍数变化(右)阈值的函数。FDR<0.05，fold change>2(灰色虚线)在其他分析中被认为有显著差异。(b)和(c)的原始数据以源数据文件的形式提供。 

9 小鼠盲肠中ecDNA的定量 

虽然我们分离ecDNA的主要目的是准确测量细胞丰度，但ecDNA的测量可能对某些研究有用，例如，通过监测死亡细胞中的ecDNA来测量特定细菌的药物杀伤效率。本研究利用BarBIQ对ecDNA样本中含有16S rRNA基因(s)的DNA片段进行了检测(方法部分，补充说明5)，检测到664个由细胞样本确定的cOTUs，并鉴定了50条独特的Bar序列。对于这些cOTU和Bar序列，我们量化了DNA片段的数量(补充数据8)。我们发现CE2-营养组中每个位置和小鼠的每单位重量的总DNA片段丰度(补充图19a)和总细胞丰度(ecDNA/细胞)的比率(补充图11)(0.01-0.04)远小于VA组(0.08-1.3)(补充图19b)。然后，我们比较了在所有细胞和ecDNA样本中通常检测到的5种cOTU的ecDNA/细胞比率。有趣的是，CE2-营养组中cOTU-CM-0074和cOTU-CM-2009的比例明显小于VA组，这与总ecDNA/细胞比例的趋势一致，而cOTU-CM-0823在CE2-营养组和VA组之间具有可比性(补充图19c)。显示cOTU差异的结果表明，ecDNA与细胞之间的比率的明显趋势，可能代表不同细菌中破碎和未破碎细胞之间的比率，可能与细菌-微环境相互作用有关。

讨论

在本研究中，BarBIQ是一种先进的高通量和精准方法，在单细胞水平上确定了微生物群。通过BarBIQ提供的基于cOTU的分析，我们发现膳食维生素A缺乏会影响小鼠盲肠近端和远端位置的微生物群，并且在喂食3周后对近端位置的影响更大。基于小鼠肠道的结构，这一观察到的现象表明盲肠内的内容物并未完全均质化，并且饮食对近端位置微生物群的影响在3周内部分转化为对远端位置微生物群的影响。重要的是，这种现象并没有被基于ASV或OTU的传统方法清楚地识别出来。因此，与传统的16S rRNA基因扩增子分析相比，BarBIQ能够准确可视化整体微生物群和单个细菌成员，从而有效地了解微生物群的功能。 BarBIQ将继续为未来的研究做出贡献，因为它可以提供包含cOTU、每个cOTU的量化以及任何微生物群落中新发现的 16S rRNA序列的数据库。重要的是，cOTU可能与具有全基因组序列和/或功能的生物直接相关。此外，基于cOTU的分析还可以通过对多种元组学、成像和/或使用一个共同单位(细胞的计算建模)的定量整合，从当前的关联研究到所需的机制研究，成功地促进各种微生物研究的进展。因此，BarBIQ可用于与微生物群相关的广泛研究领域，例如专注于肠道、口腔、皮肤、海洋、土壤、植物和其他陆地环境的研究，以阐明上述整体微生物群和单个细菌成员，这将为生物学现象提供新的见解，例如我们研究中所显示的例子。此外，通过设计靶向特异性引物，BarBIQ还可以进一步扩展到癌细胞、真菌和病毒的定量分析。事实上，高通量单细胞基因组测序，例如分析肿瘤中的细胞基因变异，一直很困难。因此，BarBIQ将有助于微生物群和高通量基因组测序相关的广泛研究。