小编读完这篇文章,最大的感受就是:同一批 RNA-seq 数据,既能分析宿主基因表达,又能重建微生物组组成,还能挖掘宿主-菌群互作! 这对于那些已经做了大量 RNA-seq 但没有配套 16S 或宏基因组数据的研究者来说,简直是一个"存量挖掘"的福音!当然,不少数据可能没有挖掘价值的,比如微生物含量过少的组织等。

背景 为了在单一样本中同时描述宿主基因表达与微生物组组成,研究者传统上需要并行开展两套独立的实验与计算流程(mRNA测序,以及16S rRNA基因测序或宏基因组测序)。
结果 研究首先利用定量微生物标准品(模拟群落)和宏基因组分析对poly(A)富集mRNA测序在重建微生物组组成方面的可靠性进行了验证,随后对30例健康外阴样本的完整队列进行了分析。关键之处在于,整个队列的分析仅依赖mRNA测序,无需并行DNA宏基因组测序。这一统一的方法使我们不仅能够分析外阴细胞的转录组,还能够解析微生物群落的组成与动态,包括微生物基因表达特征。对同一批样本进行这三个层次的分析(宿主mRNA、单个细菌物种、细菌基因通路),进一步实现了宿主-微生物分子互作的基因水平探索。利用这一统一框架,揭示了外阴微生物组显著的异质性和高度的个体间差异,识别出与阴道中描述的社区状态类型(CST)相镜像的群落结构。重要的是,发现不同的微生物组成与特定的宿主转录程序相关联:Lactobacillus crispatus(卷曲乳杆菌)与上皮分化和屏障完整性正相关;而富含Gardnerella vaginalis(阴道加德纳菌)或其他与菌群失调相关分类群的群落,则表现出与炎症相关的转录特征。有趣的是,Lactobacillus gasseri(格氏乳杆菌)此前被认为保护效果较弱,在本研究中对外阴细胞表现出中间效应。
结论 本研究不仅为这一研究不足的解剖部位提供了新的生物学见解,还引入了一种具有广泛适用性的策略,对该领域具有重要影响。公共数据库中已存有数以万计的人类RNA-seq数据集,我们的方法能够从现有转录组数据中回溯提取微生物组信息和宿主-微生物互作信号,无需额外测序或专门的微生物组实验方案。这为重新审视跨组织、跨疾病及低生物量环境的存档RNA-seq研究提供了一种强大且具成本效益的途径,揭示此前无法获取的宿主-微生物组互作层次,并最大化已发表数据的科学价值。

亮点
微生物组组成分析传统上采用两种主要策略:扩增系统发育标记基因(通称为宏条形码,例如原核生物的16S rRNA基因)和全基因组鸟枪法宏基因组测序(WMS)。宏条形码依赖使用简并寡核苷酸进行靶向PCR,从基因组DNA中扩增核糖体DNA位点的高变区;而宏基因组测序则不聚焦于单一区域,而是靶向样本中存在的所有DNA。
应用最为广泛的方法是16S rRNA基因分析,尽管快速且成本低廉,但存在若干局限:分类分辨率低(通常仅限于属级)、PCR引物偏差、非靶标扩增,以及无法表征真菌和病毒(真菌依赖18S/ITS标记,病毒则完全缺乏核糖体基因)。此外,不同16S DNA模板的扩增效率因所选变区(如V1、V2等)、基因拷贝数变异及PCR扩增偏差而存在差异。引物与保守区之间的细微错配会引入扩增偏差,导致某些区段被优先扩增而其他区段受到抑制,造成特定分类群的系统性高估或低估,甚至完全丢失。
相比之下,WMS能够捕获样本中所有DNA,包括细菌、病毒和真核生物序列,提供更全面、偏差更少的信息,可将分类分辨率提高至种乃至株级,同时减少与PCR相关的偏差。然而,WMS成本更高,且并不总是适用于低生物量样本,因为它需要相对大量的高质量微生物DNA,理想情况下宿主DNA污染要尽可能少。
这两种策略可在优化的实验设计中结合使用,充分发挥各自优势,以应对不同的研究问题,并兼顾样本量、纵向采样和稀有分类群检测等因素。与此同时,鉴于每种技术的局限性,还需要额外的多组学数据才能全面理解特定背景下微生物的功能效应。事实上,16S rRNA基因测序和宏基因组测序都只描述群落的分类组成,因而无法衡量其在特定条件下的功能活性。
为克服这一局限,近年来宏转录组学方法(分析给定多物种生态位中所有RNA的综合)对微生物群落的功能活性提供了更深入的认识。正如基因组学扩展到宏基因组学以分析混合微生物基因组,转录组学同样扩展到宏转录组学,以表征整个微生物组的群落水平基因表达模式。
宏转录组学利用微生物群落的全部RNA,直接读取功能活性,并捕获群落对特定条件或扰动的响应。在人类微生物组研究中,这种基因表达分析对于理解单个分类群如何影响健康与疾病至关重要,因为即使是同一物种,在不同转录状态下也可能对宿主组织产生截然不同的生物学效应。
这种组学方法减少了16S rRNA基因测序固有的PCR偏差,既能最小化扩增错误(尤其对于低生物量样本),又能提高物种和株级检测的准确性。然而,它也面临重大挑战:细菌rRNA相对于mRNA丰度极高、mRNA半衰期极短,以及大多数原核生物转录本缺乏poly(A)尾。此外,人类样本中微生物与宿主细胞的比例通常较低,需要对样本采集和处理进行仔细优化,以获得足够的微生物数据。
尽管上述方法在科学研究中已被广泛采用,但每种方法都依赖专门的提取和文库制备方案。特别是为了成功从不同微生物物种中提取基因组DNA,通常需要采用珠磨法进行强力机械裂解,继以酶解裂解,这与RNA提取及后续分析不兼容,因为RNA的稳定性较低。总体而言,多项研究表明,不同的基因组提取方法会对观测到的群落组成产生显著影响,尤其对于难以裂解的生物体,如革兰氏阳性菌和真菌。
对于宏基因组分析,文库由片段化DNA制备,经末端修复、接头连接后进行测序。而对于宏转录组工作流程,从微生物组样本中制备RNA文库是关键步骤,因为在保存原核mRNA的同时,还需耗尽大量丰富的核糖体RNA。为减少测序数据中的宿主污染和rRNA含量,通常采用杂交捕获技术。
然而,针对多个研究问题的研究往往需要在独立样本上并行分析,这本身构成一种局限。事实上,即使样本采自同一受试者的同一时间点,它们在数量和组成上也可能存在差异,从而引入偏差。因此,从单一样本中获取多层次信息是确保一致性和数据可比性的最佳策略。
分析表明,一套统一的方案即可在物种水平评估特定样本的微生物组成,同时保留宿主细胞和微生物表达谱的信息。为验证这一假设,我们以外阴生态系统为模型系统,因为目前对这一器官的了解非常有限。
实验室专注于健康与疾病状态下的外阴生态系统研究。在微生物组研究背景下,女性盆腔疾病常与微生物失衡(即菌群失调)相关联。在女性生殖道中,外阴和阴道微生物组是两个既独特又相互关联的生态系统,对女性健康发挥着至关重要的作用。
尽管影响阴道和外阴的疾病负担大体相近,但人类阴道微生物组的研究远比外阴微生物组深入得多。已知阴道是一个复杂的动态生态系统,主要由乳杆菌属物种主导,与其他身体部位的微生物组相比,其α多样性相对较低。这些乳杆菌通过产生乳酸和杀菌化合物发挥保护作用。一项针对394名育龄健康女性的大型研究将人类阴道微生物组分为五种社区状态类型(CST):CST-I由L. crispatus主导,CST-II由L. gasseri主导,CST-III由L. iners主导,CST-V由L. jensenii主导;CST-IV则由α多样性高、厌氧菌比例高的多种非乳杆菌分类群主导,包括Gardnerella、Atopobium、Prevotella、Dialister、Aerococcus、Megasphaera、Peptoniphilus、Sneathia、Eggerthella、Finegoldia和Mobiluncus。
由乳杆菌主导的稳定酸性阴道环境与抵御细菌性阴道病、性传播感染及宫颈癌等疾病相关。在乳杆菌物种中,L. gasseri与健康相关,但似乎处于保护谱系的中间位置。关于HPV和早产的流行病学研究表明,L. gasseri主导的群落提供的保护程度适中,其风险介于L. crispatus和L. iners之间。在这一框架下,L. gasseri可被视为有益物种,有助于维持有利的阴道环境,但其保护效果属于中等水平。相比之下,菌群失调通常以乳杆菌减少和厌氧菌扩增为特征,使生态系统向CST-IV偏移。
尽管以乳杆菌为主,阴道微生物组在个体之间及同一个体随时间的变化中都存在相当大的差异,这造就了其特有的低α多样性和高β多样性。影响这种变异的因素包括内在因素(如种族、遗传背景、年龄、激素状态和个体宿主因素)以及外在因素(如性行为、卫生习惯、吸烟、饮食、压力和肥胖)。类似地,外阴微生物组包含多种细菌属,包括Lactobacillus、Corynebacterium、Staphylococcus、Prevotella和Actinobacteria等。
这一格局提示受邻近阴道群落的影响,同时还有来源于皮肤和粪便共生菌的额外多样性。此外,现有外阴微生物组研究大多依赖16S rRNA基因测序,很少采用宏基因组方法。迄今为止,仅有一项已发表的研究将鸟枪法宏基因组测序应用于外阴硬化性苔藓(LS)和高级别鳞状上皮内病变(HSIL)背景下的外阴微生物组分析。该研究发现,LS病变以Prevotella和Bacteroidales减少为特征,而HSIL病变与健康对照相比Fusobacteria增加、Actinobacteria减少。此外,LS和HSIL外阴病变中的乳头瘤病毒科(尤其是α乳头瘤病毒)含量均高于健康皮肤。
鉴于现有微生物组分析策略的优势与局限,我们着力通过实施优化的宏转录组方法来填补外阴生态系统功能表征方面的现有空缺。具体而言,我们开发了一套统一的提取和文库制备方案,能够从单一样本中同时分析外阴微生物组和宿主上皮基因表达。这一整合策略使我们能够重建外阴微环境的微生物组成,同时捕获在同一生物学背景下发生的宿主转录应答。宿主细胞与微生物组通过持续而复杂的互作共存并相互影响。这种互作可显著影响宿主和微生物细胞的基因表达,这一现象在外阴组织中尤为缺乏深入研究。
通过将微生物信息与宿主表达数据相结合,我们的方法提供了一个全面框架,用于研究特定细菌群落如何影响外阴上皮活性,从而为研究这一尚未充分认识的解剖部位的宿主-微生物互作提供了更高效、生物学上更连贯且更具成本效益的方法。
为验证RNA-seq数据重建微生物组成的能力,我们首先采用了一种商业化的微生物模拟群落(ZymoBIOMICS微生物群落标准品;Zymo Research),其中包含3种易裂解的革兰氏阴性菌、5种难裂解的革兰氏阳性菌和2种难裂解的酵母菌,其丰度根据细胞数量和基因组DNA含量确定。

同时,我们还从商业化人类外阴癌上皮细胞系(VLS-VSCC)中提取了总RNA。从所有样本中,我们使用TRIzol和适用于真核细胞的标准方案提取总RNA,不进行机械或酶解的细菌细胞壁破碎,以模拟传统的哺乳动物基因表达分析。我们按照制造商说明,分别用poly(A)富集和总RNA试剂盒制备文库,分别对一个微生物RNA样本和三个含不同比例人源与微生物RNA混合物的样本(分别含25%和60%人源RNA)进行处理。
我们对样本进行了测序,并将原始读段比对至人类基因组以去除宿主污染。剩余未比对读段经自定义Kraken2/Bracken流程处理,以获得每个样本的分类丰度。如预期所示,5种革兰氏阳性菌株中有4种相对于其理论浓度被低估,这可能是因为我们未进行专门的细胞壁裂解。重要的是,10种物种相对丰度的检测结果在样本间保持一致,与人源和微生物读段的比例无关。
我们的结果与使用人类微生物组计划(HMP)方案获得的结果十分接近,与模拟群落制造商文档中的报告数据相符。尽管与真实丰度存在偏差(但这种偏差在各样本间是一致的),我们的初步结果表明,真核mRNA测序工作流程也可成功地应用于重建微环境生态位组成,如两种方法均能可靠地检测到优势物种所示。
在评估了该方法的可靠性之后,我们对采集自健康个体前庭的外阴拭子进行了RNA-seq分析。样本通过轻柔拭取小阴唇内侧(特别是外阴前庭区域)采集,同时注意避开阴道口。
我们选择外阴作为研究系统,是因为尽管阴道微生物组已被广泛研究,外阴微生物组的复杂性仍知之甚少,是一个新兴研究领域。共纳入30名年龄在23至44岁之间、无复发性盆腔疾病病史的健康女性。对于部分受试者,我们还提取并处理了DNA用于宏基因组分析,以比较微生物组分析结果。我们未进行16S rRNA宏条形码测序,原因如下:样本生物量低、该方法固有的PCR偏差,以及需要物种级微生物分辨率以与宏转录组数据进行比较。两种DNA和RNA文库的质量控制指标、基因组比对效率及人源唯一比对、微生物Kraken指定和未比对读段的分布见图S2A-D及补充材料。利用这一体内数据集,我们随后系统研究了在poly(A)选择文库中回收细菌mRNA的机制。我们对捕获的读段进行了计算扫描,寻找内部poly-A/T片段(8、10和12 nt),这些片段可能在RNA片段化之前促进oligo-dT引物错误引导。我们将观测频率与经验背景概率进行比较,如预期所示,人类宿主读段中A/T片段的频率较高。相反,在细菌组分中,8 nt A/T片段的观测频率(约2.5%)显著低于预期背景概率(约5.5–6%),反映了细菌编码序列中长A/T同聚物在进化上受到抑制。这表明,尽管内部位点的机会性错误引导确实存在,但它不能完全解释所回收微生物组的全部。细菌RNA更可能是通过非特异性物理携带(physical carryover)被保留下来的,这一捕获机制解释了我们定量标准品中poly(A)和rRNA耗竭文库之间观测到的组成一致性。

为初步验证我们方法在离体样本中的可靠性,对6个样本同时进行了DNA和RNA提取后测序。宏基因组(DNA)和宏转录组(RNA)数据集的比较分析证实了两种技术获得的分类学图谱之间高度一致(R≈0.8),尽管读段在基因组中的分布存在较大差异)。主坐标分析(PCoA)显示配对的DNA和RNA样本紧密聚集,两种数据类型之间的分离极小。层次聚类进一步支持了这种一致性,因为所有病例中配对的DNA和RNA样本均聚集在一起,反映了宏基因组和宏转录组图谱在群落水平的高度相似性。分类组成比较显示,Lactobacillus crispatus和Lactobacillus iners等优势物种在DNA和RNA数据集中均得到一致检测。

然而,在相对丰度方面存在一些差异,尤其是对于丰度较低的分类群(Prevotella bivia),以及G. vaginalis(其在RNA图谱中相对于DNA方法表现出更高的丰度)。我们还分析了每个基因/通路的RNA/DNA丰度一致性。我们发现,特定的两条通路中的基因在RNA分析中被过度表达(核糖体和代谢通路),另外两条通路在RNA分析中被低表达(核苷酸切除修复NER和碱基切除修复BER)。尽管每个样本由不同的物种主导(L. iners、L. crispatus和G. vaginalis),这一趋势在所有样本中总体上得到维持,提示这是不同物种间的共同特征。

综上所述,这些数据证实了宏基因组和宏转录组方法在确定微生物组成方面的高度一致性。
在评估了我们方法的可靠性之后,我们将其应用于从30名绝经前健康个体采集的30份外阴拭子的mRNA-seq数据。众所周知,外阴微生物组与阴道微生物组非常相似,但外阴微生物组通常因皮肤和粪便来源共生菌的存在而表现出更高的多样性。为进一步研究微生物组成的异同,我们对批次校正后的数据进行了NMDS分析,比较了来自已发表宏基因组数据集和我们自身样本的阴道、外阴和肛门微生物组。

分析揭示,每个身体部位的样本聚集在一起,每组的质心明显分离。这种明显的分离强调了每个解剖环境特有的微生物群落。值得注意的是,我们的宏转录组样本与已发表的外阴宏基因组数据紧密聚集,进一步验证了我们方法的可靠性。个体间差异可能受采样位置的影响(如外阴黏膜、皮肤皱褶、肛门或肛周区域),在我们的数据集中采样限于前庭区,而已发表数据集包含了更多样化的部位。
有趣的是,外阴微生物组的质心占据了肛门和阴道群组之间的中间位置,反映了其中央的解剖学位置。在检查三个环境中共同存在的物种时,我们发现在每个区室前100个表达物种中,外阴与肛门共享的物种多于与阴道共享的物种。值得注意的是,没有任何物种仅在肛门和阴道之间共享而外阴没有,进一步支持了外阴作为这两个在其他方面截然不同的环境之间微生物中介者的概念。相对丰度分析同样表明,外阴代表了一个中间环境,同时共享来自肛门和阴道的微生物特征。
这些发现进一步验证了我们宏转录组方法的稳健性,并强化了外阴微生物组是一个由阴道、皮肤和粪便来源的微生物群落共同塑造的异质性环境这一假设。我们着力寻找外阴和阴道微环境之间的共同特征。我们因此假设,在不同个体阴道微生物组中发现的社区状态类型(CST)也可以在外阴微环境中得到体现。
为验证这一假设,我们对未标准化数据进行了聚类分析,分析了阴道微环境中丰度最高的前100个物种,从而找到由最优势微生物驱动的群簇。分析揭示出五大微生物群落群组中的四个,其中三个由乳杆菌主导,一个由非乳杆菌属细菌组成,没有任何个体由L. jensenii主导。具体而言,我们发现16个样本由L. crispatus主导,9个样本由L. iners主导,3个样本由L. gasseri主导,3个样本由其他物种主导。这四种CST不仅与阴道研究中先前确立的CST相吻合,而且这些组内的样本分布也与阴道CST类型相似,在健康个体中CST-I和CST-III最为常见,相对于CST-II和CST-V更为丰富,再次强化了这些解剖部位之间的生物学关联。在尝试校正任何裂解和提取偏差时,我们应用了Rauer等人描述的基于形态学的计算校正(见方法部分)。
我们因此估算了不同模拟群落(图1中的poly(A) mRNA样本)的校正偏差。由于驱动外阴生态系统差异的关键物种不受提取偏差的严重影响(图S3D),且任何基于形态学的校正都会去除病毒和真菌,我们在下游分析中保留了原始丰度。从全局视角来看,NMDS分析也证实了CST分析识别的主要物种是样本分离的主要驱动因素(图S3E)。从热图来看,我们未观察到主要的共现细菌亚群,除了L. crispatus群组(粉色)中的一个共现物种群簇。有趣的是,该群簇主要由乳杆菌和芽孢杆菌(Bacillus)物种组成(包括B. subtilis、B. cereus和B. velezensis),这些物种已知在阴道微环境中具有保护作用。此外,观察不同属在各CST中的相对富集情况,我们未发现显著差异,唯一例外是上述Bacillus属在CST-I和CST-II中相对于CST-III和CST-IV更为富集。
我们还尝试识别其他相关变量(饮食、吸烟、年龄等)是否能解释队列中观测到的部分变异,但未发现任何显著性(数据未显示)。
对每位个体,我们评估了生物多样性(以每个样本的α多样性衡量),并在4个CST群簇之间进行比较。如预期所示,CST-IV因其定义上的"其他"物种富集,与L. iners和L. crispatus主导的CST相比,表现出较高的生物多样性趋势。
在优势属(定义为相对丰度>5%)中,我们观察到显著的多样性。厌氧和皮肤相关属(如Anaerococcus、Veillonella和Cutibacterium)较为普遍。此外,还检测到常与阴道或胃肠道微生物组相关的属,包括Fannyhessea、Gardnerella和Prevotella。有趣的是,Gardnerella表现出与优势属相互排斥的趋势,而Prevotella和Ureaplasma则在样本间呈现更为均匀的丰度分布。
在真菌群落中,只有Malassezia的相对丰度超过>5%的阈值,其次是Saccharomyces(丰度较低,<1%)。这些发现凸显了外阴微生物生态系统的丰富性和复杂性,涵盖来自多个微生物生态位的贡献,以及我们方法同时研究外阴微环境中细菌和真菌组分的能力。
为阐明外阴微生物组与宿主细胞过程之间的功能互作,我们通过计算单个人类基因与微生物之间的相关性构建了整合网络。这一方法尤为有利,因为宿主基因表达和微生物组数据均来自相同的RNA-seq样本。同步采样两个参数不仅减少了技术批次效应的影响,还增强了相关性的可靠性。

微生物互作网络揭示,尽管乳杆菌属物种的平均相对丰度最高(节点大小),但与丰度较低的分类群相比,其中心性较低,后者被识别为网络的主要结构枢纽。这些分类群与G. vaginalis正向连接,形成一个整合良好的条件致病菌和菌群失调物种枢纽。虽然网络中的连接以正向为主,但少数负向相关特异性地将L. crispatus与L. iners和Gardnerella属联系起来,代表了相互排斥的拮抗关系。

以属特异性相关群簇为特征的微生物群落结构随后被映射到宿主转录组上,揭示了一个复杂的相互作用组,其中特定微生物充当宿主基因模块的关键调控因子。这些转录组模块与CST之间的关系通过模块-性状分析进一步得到验证,突出了CST组成与角化(green模块)和免疫激活(turquoise模块)等生物过程之间的显著关联。对功能景观的深入研究揭示,微生物影响高度专门化且具有两极分化特征。通过量化特定基因本体(GO)条目中互作的百分比,我们观察到,潜在的条件致病菌(即L. gasseri或G. vaginalis)对上皮完整性(角化和上皮细胞分化)呈负向调控,同时伴随免疫和炎症通路的激活。这种专门化调控通过其协同互作图谱得到印证,其中两种菌通过重叠基因靶向共同通路。特别是,G. vaginalis与一般免疫相关基因高度连接,而两种细菌都与先天免疫基因正向互作。

最后,我们识别出一组"桥接基因"(定义为与多个微生物分类群互作的宿主节点),它们充当外阴环境中的核心分子开关。这些桥接基因的共调控热图揭示了一种二元控制机制:尽管大多数宿主过程对微生物变异具有一定弹性,但炎症应答和上皮分化中的关键检查点(如NF-κB信号、IL-1β介导通路)受到共生菌和条件致病菌相反调控压力的共同影响。这一整合分析使我们能够赋予那些在更广泛数据集中可能被忽视的单个基因深刻的生物学价值。一个典型案例是ZDHHC5——一种棕榈酰基转移酶,是激活gasdermin D介导的焦亡(pyroptosis)所必需的,焦亡是一种细胞自主的抗病原体和炎症疾病防御机制。
值得注意的是,我们的分析将ZDHHC5识别为G. vaginalis、L. gasseri和L. crispatus之间的桥接基因。ZDHHC5与G. vaginalis和L. gasseri的正相关表明,这些微生物可能使外阴上皮进入焦亡信号的准备状态。相反,其与L. crispatus的负相关表明对ZDHHC5-GSDMD轴的保护性抑制,可能有助于维持屏障完整性并限制细胞因子释放。综上所述,这些数据提示外阴健康是通过微生物输入的精妙平衡来维持的,这些输入汇聚于特定宿主桥接基因上,调节组织稳态和免疫耐受。
随着实验策略和计算方法的协同进化,微生物组分析变得日益精细、灵敏和特异,能够对复杂微生物群落进行高分辨率表征。然而,大多数现有方案都有其固有的优势和局限,因此针对特定生物学问题审慎选择实验和分析流程至关重要。
我们实验室主要专注于人类外阴疾病中的RNA扰动研究,我们近来认识到我们的测序数据集中还包含了此前未经探索的微生物群落组成信息。
通过在体外和离体模型中验证基于DNA和RNA的分析之间的一致性,我们证明了常规poly(A)富集或核糖体RNA耗竭的mRNA-seq可以被重新用作一种统一的、类宏转录组的策略,用于识别微生物物种并估计其相对丰度。
从RNA和DNA分析推断的微生物丰度之间的轻微差异也可能反映了特定分类群内的转录激活。在活细胞中,调节蛋白质合成最快速的机制之一是通过核糖体RNA含量的变化,以适应不断变化的环境条件。核糖体基因以及代谢通路相关基因在RNA水平上的相对高丰度表明,它们对于细胞在动态外阴微环境中的存活至关重要。
我们的软裂解RNA测序方案的效率偏向于易裂解的革兰氏阴性菌,而革兰氏阳性菌通常被低估。重要的是,这种偏差在不同条件和样本间是一致的。尽管存在这一局限,一套方案实现了对外阴微生物组和宿主上皮转录组的同步分析,来自同一样品。为将软裂解方案推广至具有更复杂生态系统的其他身体部位,整合基于形态学的计算校正可能对于确保准确的绝对丰度估算至关重要,尤其是在难以裂解细菌比例较高的情况下。
另一个局限在于低生物量样本,这类样本本质上放大了环境和试剂污染的挑战。完全匹配的物理阴性对照是直接扣除背景的金标准,我们成功地将这一策略应用于我们的DNA文库。然而,这种方法在涉及复杂、多批次、多地点样本采集的临床研究中往往在操作上具有挑战性。在持续物理对照分散(如我们RNA拭子的情况)的情况下,使用强健的统计方法(如decontam R包)对于准确识别和去除污染物变得至关重要。根据每个数据集的操作限制来定制去污策略是确保最高准确性的关键步骤。
我们将方法应用于外阴微环境,这是一个新兴研究领域,目前只有少数近期研究开始探索其复杂性。尽管研究有限,现有证据表明外阴微生物组与阴道环境共享若干特征。然而,它通常表现出较高程度的皮肤和粪便来源共生菌定植,赋予外阴独特的微生物特征。本研究表明,外阴社区状态类型(CST)与规范的阴道CST相对应,同时占据肛门和阴道群落之间的中间生态位,支持外阴是一个独特但相互关联的微生物生态系统的概念。值得注意的是,CST-I群簇不仅富含Lactobacillus crispatus,还富含与阴道保护相关的芽孢杆菌属物种,提示CST-I通常具有的保护特征可能并非仅由乳杆菌驱动,而可能因其他保护性物种的同时存在而得到强化。有趣的是,外阴菌群失调与更高的生物多样性相关(图S3G、H),这与肠道中更高的微生物组异质性总体上与健康状态相关的规律恰好相反,提示外阴微环境存在更严格的控制机制。
我们的发现表明,不同的微生物群落与特定的宿主转录程序相关联。Lactobacillus crispatus主导的群落与上皮分化和屏障完整性相关基因正相关,而富含Lactobacillus gasseri的群落则与炎症转录特征正相关,与皮肤角化负相关。这些趋势支持了L. crispatus的保护作用和L. gasseri的中间角色,与文献数据一致,同时提供了宿主与微生物组之间更直接的分子联系。
我们方法的样本内框架整合了微生物丰度与宿主基因表达,揭示了特定的宿主-微生物互作模块,将乳杆菌主导的群落与上皮分化和屏障维持联系起来,将菌群失调相关分类群与炎症转录程序联系起来。
这些发现共同提供了一个概念和技术框架,用于从现有RNA-seq数据集中挖掘微生物组和宿主-微生物互作信号。尽管在poly(A)选择文库中回收细菌转录本看似违反直觉,我们的数据表明这可能是由罕见内部错误引导和非特异性物理携带的综合作用驱动的(图S1C),这通过poly(A)富集和rRNA耗竭文库之间一致的微生物图谱得到证实(图1C)。
因此,尽管存在这些技术特殊性,该方法确保核心生物学信号被准确捕获,并提供了一种可扩展、具成本效益的途径,用于重新审视存档的转录组队列,揭示微生物组在健康和疾病中的功能影响。
尽管对外阴微生物组功能活性的全面分析并非本研究的主要目标,我们采用的方法也可能有助于研究外阴微生物的动态基因表达谱,从而提供外阴环境中微生物活动更准确的快照。
总之,我们的工作表明,常规mRNA-seq数据集可被可靠地重新用于同步分析微生物组成和宿主上皮转录组,在同一样品中提供宿主-微生物互作的整合视角。这一整合策略:(1) 消除了制备独立RNA和DNA测序文库相关的批次效应;(2) 最大化单次实验的信息产出;(3) 降低实验成本;(4) 允许研究宿主基因表达与微生物物种及其功能通路之间的相互关系。将这一框架应用于外阴微环境,我们表明外阴社区状态类型反映了类似阴道的结构,同时定义了一个独特的生态位,特定的微生物配置与控制屏障完整性和炎症的宿主转录程序紧密关联。超越对外阴生态系统的表征,本研究建立了一种可扩展且具成本效益的概念和方法策略,用于挖掘现有RNA-seq队列中的微生物组成、活性和宿主-微生物互作信息,涵盖健康和疾病状态。
样品制备中,我们同时采用了VLS-VSCC人类外阴鳞状细胞癌细胞系与ZymoBIOMICS微生物群落标准品(Zymo Research;欧文,加州,美国)的共提取和提取后混合两种方式。共生成了四种用于RNA提取的输入样品:(1) 含60%人源细胞和40% ZymoBIOMICS标准品的共提取样品;(2) 仅含Zymo的阳性对照;(3) 含75% ZymoBIOMICS和25%人源RNA的提取后混合样品;(4) 条件3的技术重复。所有四种输入样品均采用基于TRIzol的方案提取总RNA,并用不含核酸酶的水洗脱。两个提取后混合重复样品通过手动混合人源RNA和Zymo洗脱液制备。
为评估细菌mRNA在poly(A)选择文库中的捕获机制,对比对到宿主和微生物的读段进行了计算扫描,寻找内部poly-A/T片段(8、10和12 nt)。为建立预期随机发生的基准,通过序列置换计算经验背景概率,在严格保留读段长度和GC含量的前提下,对细菌读段进行随机化处理。随后将实际细菌读段中A/T片段的观测频率与这一置换基准进行比较,以估计机会性内部错误引导与非特异性物理携带的相对贡献。
30名健康志愿者于2023年2月至2023年7月在盆腔健康中心(Centro Salute Pelvi)招募。纳入标准:无外阴疼痛或外阴症状的绝经前女性(年龄18至50岁)。排除标准包括:绝经后、妊娠、哺乳,或产后4个月内,以及在访视前1周内接受过系统性或阴道抗菌或益生菌治疗者。本研究经都灵大学伦理委员会审查批准(批准号:Prot. N. 0214340,2023年4月6日)。所有参与者提供了参与本研究的书面知情同意。
我们用温和的无菌干拭子拭取前庭区域采集外阴样品,随后将拭子置于含RNA/DNA保护试剂的1.5 ml离心管中以保存核酸。RNA提取通过加入1 ml QIazol/TRIzol裂解试剂(Qiagen)进行。在室温孵育3分钟后,加入200 μl氯仿,旋涡振荡15秒,室温孵育3分钟。随后以最高速度在4°C离心15分钟,将获得的水相转移至新管中,用500 μl异丙醇和1 μl糖原在室温沉淀。10分钟后,在最高速度下4°C离心15分钟。将沉淀用70% EtOH洗涤,重悬于DNase/RNase-free水中,在Nanodrop上定量。
外阴拭子用无菌棒状拭子采集后立即储存于-80°C。将拭子浸泡在1 ml 0.9% NaCl中,旋涡振荡5分钟以重悬样品。在室温下以13,000 rpm离心10分钟后,弃去上清,将沉淀重悬于80 μl TrisHCl中。随后加入20 μl溶菌酶(10 mg/ml),在37°C孵育30分钟使细胞壁降解。孵育后,使用TissueLyser II在30 Hz下以1 mm金属珠进行珠磨2分钟。随后,加入1 μl蛋白酶K(20 μg/μl)去除蛋白质,在56°C孵育30分钟。加入500 μl苯酚-氯仿-异戊醇和500 μl 0.9% NaCl提取DNA,室温孵育10分钟,以11,000 rpm离心5分钟。相分离后,向水相加入1 μl糖原、30 μl 3M醋酸钠和800 μl无水EtOH,在-20°C沉淀过夜。次日,以13,000 rpm在4°C离心30分钟纯化DNA,弃去上清,用500 μl 70% EtOH洗涤样品。以13,000 rpm在4°C再次离心5分钟后,弃去上清,将沉淀重悬于30 μl无RNase水中。随后在Qubit仪器上定量。
RNA-seq文库制备时,用Qubit和TapeStation(Agilent)检测起始RNA的数量和质量。图1中体外实验的文库由500 ng总RNA使用Watchmaker mRNA文库制备试剂盒和带Polaris®耗竭的Watchmaker RNA文库制备试剂盒按制造商说明制备。
图2、3和4中离体样本的文库专门使用Illumina Stranded mRNA制备试剂盒制备。
DNA-seq文库时,100 ng DNA在Pico Bioruptor(Biosense srl)上进行超声处理,3个循环(开机20秒,关机30秒)。文库使用Watchmaker DNA文库制备试剂盒按制造商说明制备。
文库在TapeStation(Agilent)上进行质量评估,在Illumina NextSeq1000系统上测序(配对末端60 bp读段)。
核心生物信息学分析使用自定义Snakemake流程执行。简要地说,流程按顺序执行以下规则:
i) 使用fastp工具(https://github.com/OpenGene/fastp)的snakemake包装器对FastQ数据进行预处理和质量控制。
ii) 为解决旧参考基因组中常见的高度重复区域、着丝粒间隙和结构变异等问题,我们采用了两步比对策略。RNA-seq配对末端读段分别比对至两个人类基因组:hg38(GRCh38)和T2T-CHM13(端粒到端粒)组装,使用STAR比对工具。这种双重比对方法确保即使低复杂度或非规范人类读段也能被识别。仅未能比对到任一人类组装的读段被保留用于微生物分析。
iii) 剩余读段使用Kraken2与预配置参考数据库Kraken2 PlusPF(标准版加上RefSeq原生生物、真菌和人类)进行分类。在该数据库中加入人类基因组充当次级"安全网":任何由于序列差异或文库伪迹而绕过初始比对过滤的人类读段,都会与内部人类k-mer库进行交叉比对,被标记为"人类"或"未分类"。Kraken2分类使用附加选项(--minimum-hit-groups 3)执行,以实现更严格的分类。通过利用Kraken2的最低共同祖先(LCA)算法,我们的流程通过将含有冲突k-mer特征的读段分配到高层分类学节点或将其丢弃为未分类,从而确保对嵌合伪迹具有鲁棒性,防止特定微生物丰度被错误地膨胀。
iv) 使用预配置数据库(Kraken2 PlusPF)通过Bracken处理Kraken2报告文件,参数:-r 50 -l S -t 5。Bracken通过分析特定分类学水平上的k-mer分布重新估计分类学丰度。输出包括详细的丰度报告和含有丰度数据的附加文件用于下游分析。
v) 为监测双重转录组学方法的性能和捕获效率,在流程的每个步骤计算读段分布统计。具体而言,追踪了原始读段总数、经fastp质量控制保留的读段(QC reads)、通过STAR成功比对到人类参考基因组的宿主来源序列的比例,以及剩余未比对部分。最后,通过统计Kraken2成功分类的读段数量化微生物含量。每个处理阶段的完整读段统计数据针对整个队列进行了制表,并在补充图S2A–D中进行了可视化展示。
为处理数据集中的潜在污染,使用R中的decontam包进行了去污染处理。污染物通过基于频率的方法识别,该方法评估分类群丰度与DNA浓度之间的关系。对于DNA浓度测量(Qubit定量)中含有非数字或"过低"值的样本,通过将"过低"替换为最低值0.01 ng/μl进行调整,以允许正常分析。RNA提取批次被指定为批次变量,以解释处理过程中的潜在差异。识别为污染物的物种被排除在数据集之外。最后,对数据集进行标准化以确保分类群的比例表示,仅将非污染物种保留用于下游分析。
为处理DNA样本中来自阴性对照的潜在污染,采用归一化程序旨在去除阴性对照样品中污染物引入的背景信号,确保只有生物学相关信号被保留用于下游分析。数据集根据两个不同的测序文库被分成两个子集。对于两个子集,每个样本列中的值通过减去相应阴性对照的值进行调整。具体而言,对于每个样本列,减去阴性对照列的值,任何由此产生的负值被设为零。这确保在归一化过程中不引入人为的负值。所有样本中总丰度为零的行被排除,确保只有具有有意义信号的物种保留下来。
数据集经过过滤,以仅保留在至少1个样本中相对丰度超过0.1%的物种。在显示两种技术比较的NMDS图中,不进行丰度过滤。
聚类在未标准化数据上进行,以相对丰度作为输入,有意地优先考虑优势细菌的贡献,并防止"双零效应"(稀有分类群之间共同缺失)干扰分析。使用欧氏距离计算不相似性。为确保识别的群簇不是该度量标准的伪迹,样本分组使用Bray-Curtis不相似度进行独立验证(图S3B、E)。随后使用Ward最小方差法(ward.D2)进行层次聚类,在每一步最小化群簇内的总方差。随后生成树状图可视化;为提高可解释性,树状图的分支经色彩编码以代表四个群簇。
为评估微生物组的致病潜力,我们建立了一个专门基于细菌性阴道病的参考菌群失调环境。我们查询了MicroPhenoDB数据库,为每个识别的分类群分配一个关联评分,反映其与BV的已知关联。随后对每个样本通过对对应分类群相对丰度加权的关联评分求和计算累积菌群失调评分。评分越高表示微生物图谱越接近BV相关群落。
非度量多维尺度(NMDS)用于探索不同研究中样本间的组成差异。使用sva包中的ComBat函数校正批次效应(本研究与其他研究之间)。计算Bray-Curtis不相似度以量化样本间的两两差异。NMDS排序以最多100次迭代执行以确保收敛。样本根据其NMDS1和NMDS2坐标绘制,点按采样位置色彩编码,按批次区分形状。代表每个位点95%置信区间的椭圆叠加在图上,以突出各组间的聚集模式。沿各轴的边际密度图提供了额外背景,显示NMDS1和NMDS2上的样本分布。
主坐标分析(PCoA)用于评估来自RNA或DNA测序样本之间的组成不相似性。使用Bray距离度量量化样本间的不相似性。为解决潜在的负特征值,在排序过程中应用了Cailliez校正。
首先确定过滤后数据集中的最小和最大测序深度。为标准化文库大小,使用vegan R包中的"rrarefy"函数将样本稀疏至最小测序深度。随后使用diversity函数对亚采样群落矩阵计算Shannon多样性指数。
为识别外阴微生物群落内的依赖关系,我们使用中心对数比率(CLR)转换构建相关网络,以减轻成分数据中固有的虚假相关。随后对转换后的值应用Spearman秩相关,以捕获单调依赖关系,确保分析对微生物组数据集中常见的非正态分布和离群值保持稳健。这种整合方法为检测独立于测序深度变化和非线性增长动态的生物学关联提供了数学上严格的框架。
相关计算前,我们保留在至少一个样本中相对丰度达到至少0.1%的微生物,共84种微生物。相应的相关网络加载到R中并进行过滤,仅保留绝对值至少为0.4且FDR小于5%的相关。随后使用Cytoscape(https://cytoscape.org/)可视化网络,并使用"microbetag" Cytoscape插件(https://hariszaf.github.io/microbetag/)注释。使用MCODE插件(https://mcode.readthedocs.io/en/latest/)进行网络聚类,参数:K-core=2;节点评分阈值=0.2;最大深度=100;度阈值=2。宿主基因-微生物组相关网络首先通过标准化原始计数构建,使用DESeq2包的normalization方法和varianceStabilizingTransformation函数进行标准化和转换。标准化数据随后使用WGCNA R包的corAndPvalue函数进行相关分析。保留绝对值至少为0.6且FDR小于5%的相关用于进一步分析。基因模块和模块-性状关系使用WGCNA R包获取。基因模块的功能富集使用enrichR(https://maayanlab.cloud/Enrichr/)进行。生成的网络随后导出并使用Cytoscape进行可视化。通过选择属于特定模块的所有基因并包括每个基因的第一个微生物邻居,检查基因模块间的互作。为评估微生物组对宿主的功能影响,与最丰富微生物分类群显著相关(Spearman |R|>0.5)的宿主基因使用基因本体生物过程(GO BP)条目进行注释。单个微生物与特定宿主上皮功能之间正向和负向相关的相对比例经计算并以发散条形图可视化。此外,进行了"桥接"分析,以识别与多个不同微生物分类群(n≥2)表现出显著共调控的宿主基因。这些多分类群互作基因的相关图谱使用层次聚类热图(Ward.D2方法)可视化,按其指定的生物功能分组。