《mSphere》：Co-occurrence network analysis reveals novel associations between the neonatal airway microbiome and bronchopulmonary dysplasia risk: an observational, population-based study
编辑推荐：
本研究通过16S rRNA测序和微生物共现网络分析，首次揭示了新生儿气道微生物生态系统的网络结构（包括网络密度和关键类群）与支气管肺发育不良（BPD）严重程度显著相关。研究发现，随着BPD严重程度增加，微生物网络复杂性显著降低，而较高的网络密度（OR=0.12）是BPD的独立保护因素。该研究为BPD的早期风险分层和微生物生态干预提供了新视角。
ABSTRACT
本研究旨在评估早产儿出生时呼吸道微生物与后续支气管肺发育不良（BPD）发生及严重程度的关联。这项前瞻性队列研究纳入了98名早产儿（胎龄<32周，出生体重<2,000克），在出生后2小时内通过气管插管收集气管吸出物样本，采用16S rRNA测序技术分析气道微生物组，并运用成分稳健的方法进行共现网络分析。在98名早产儿中，BPD发生率为68.4%，其中Ⅰ级31例、Ⅱ级20例、Ⅲ级16例。BPD患儿的气道微生物群呈现明显的严重程度分级模式：Escherichia-Shigella和Streptococcus在Ⅰ级显著富集，而Ⅲ级中Chryseobacterium显著增加，同时Streptococcus明显减少。微生物共现网络分析得出三个关键发现：（1）网络复杂性随BPD严重程度急剧下降，Ⅲ级最为稀疏；（2）不同组别存在独特的关键类群：非BPD组为Acinetobacter和Fusobacterium；Ⅰ级为Brevundimonas和Fusobacterium；Ⅱ级和Ⅲ级为Fusobacterium和Acinetobacter；（3）经多变量模型校正关键临床混杂因素后，出生时较高的微生物网络密度与BPD风险显著降低独立相关（OR=0.12）。新生儿出生时气道微生物组的生态结构（由网络复杂性和关键类群定义）与BPD严重程度相关，这凸显了微生物网络稳定性作为新因素以及生态相互作用作为未来研究靶点的重要性。
INTRODUCTION
支气管肺发育不良（BPD）是新生儿慢性肺疾病，也是早产儿最常见并发症。尽管新生儿护理取得进展，BPD发病率仍居高不下，约影响56.2%的胎龄（GA）小于28周和29.2%的GA小于32周的婴儿。BPD对早产儿心肺功能和神经发育有不良影响，加重公共卫生负担。BPD发病机制涉及遗传易感性、宫内和产后炎症、感染及围产期因素。系统综述提示气道微生物组失调可能与BPD进展和严重程度相关。严重BPD早产儿气管支气管吸出液中革兰氏阴性菌检出率显著较高，且气道革兰氏阴性菌定植与严重BPD患儿死亡或出院时需氧风险增加以及长期不良结局风险增加相关。微生物群落中微生物相互作用促进群落稳定性，但微生物如何相互作用影响BPD发生尚不清楚。共现网络分析利用图论说明群落间潜在相互作用，是研究群落内不同微生物群之间复杂相互作用的强大工具。本研究旨在描述早产儿出生时下呼吸道（LRT）微生物组特征，并应用共现网络分析探讨其与BPD严重程度的关系，假设特定微生物定植模式（尤其是关键类群）与BPD发生和严重程度相关。
MATERIALS AND METHODS
研究设计与参与者
这项前瞻性队列研究于2022年2月至2024年3月在厦门市妇幼保健院（中国福建厦门）进行。
纳入与排除标准
纳入标准包括：（i）GA在24+0/7至31+6/7周之间，出生体重400–2,000克；（ii）出生后2小时内需要气管插管；（iii）单次呼吸道分泌物样本量≥0.5 mL；（iv）获得婴儿法定监护人的知情同意和书面授权。排除标准包括：（i）严重先天性呼吸道结构异常；（ii）非呼吸系统疾病死亡且校正GA低于36周；（iii）样本收集前接受全身抗菌治疗；（iv）样本收集前气管内给药；（v）样本收集前开始肠内喂养；（vi）羊水胎粪污染或血性羊水。随访至出院或校正GA达36至40周。
样本量考量
所有参与者均在出生后2小时内在NICU接受气管插管。根据2021年中国新生儿协作网（CHNN）报告，中国大陆极早产儿BPD发病率约为29.2%，本探索性研究计划纳入100例以确保足够大的队列进行初始微生物生态学分析。
标本收集与处理
LRT分泌物标本在出生后2小时内无菌条件下通过气管插管封闭抽吸收集，每样本最小体积0.5 mL以保证足够微生物生物量用于下游分析。所有收集样本立即储存于-80°C直至DNA提取。
临床结局
主要结局是参与者BPD严重程度。本研究遵循美国尤尼斯·肯尼迪·施莱弗国家儿童健康与人类发展研究所制定的BPD诊断和严重程度分级标准。BPD定义为持续肺实质疾病（经影像学检查诊断），在妊娠36周时需要任何形式的呼吸支持（包括有创或无创通气或补充氧气）。BPD根据不同吸入氧浓度（FiO2）范围分为Ⅰ、Ⅱ、Ⅲ级，还包括因持续肺实质疾病和呼吸衰竭导致的早期死亡。
提取与测序
使用FastDNA Spin Kit for Soil从103个临床样本和一个试剂盒对照中提取基因组DNA，其中98个样本产生高质量基因组DNA。为监控潜在污染，在DNA提取和PCR扩增过程中包含无菌水作为阴性对照。通过琼脂糖凝胶电泳确认阴性对照中无扩增，支持临床样本微生物谱的特异性。对16S核糖体RNA基因的V3和V4区域（338F: 5′-ACTCCTACGGGAGGCAGCAG-3′, 806R: 5′-GGACTACHVGGGTWTCTAAT-3′）在Illumina MiSeq平台上进行测序。PCR条件：95°C 3分钟，随后30个循环（95°C 30秒、55°C 30秒、72°C 45秒），最后72°C延伸10分钟。所有测序文库使用NEXTFLEX Rapid DNA-Seq Kit统一制备，配对末端测序在单个Illumina MiSeq PE300平台上使用相同试剂盒类型进行，以进一步减少技术变异。
生物信息学处理
序列数据进一步预处理在QIIME2流程上进行。基于测序质量进行配对末端读长质量控制和过滤，并根据配对末端读长间重叠关系进行拼接，获得质量控制拼接后的优化数据。然后使用序列去噪方法（DADA2）处理优化数据，获得扩增子序列变异（ASV）代表性序列和丰度信息。为减轻潜在污染物和测序伪影影响，采用保守过滤方法：去除总频率为1的所有特征（单例）以及分类为线粒体或叶绿体的序列。在此过滤步骤后更新分类分配。ASV针对Silva参考数据库（版本138.2）进行聚类，并使用RDP分类器（版本2.11）通过朴素贝叶斯算法分配分类身份，该算法在SILVA v138.2参考序列上训练。
统计分析
评估BPD严重程度与气道细菌微生物组关键特征之间的关联：群落组成、α多样性、β多样性和细菌共现网络。使用Mothur软件（版本1.48.0）计算α多样性指数（包括Chao1物种丰富度估计、Shannon指数、Simpson指数和Good's覆盖度指数）。使用Mothur中Bray-Curtis相异矩阵通过相似性分析（ANOSIM）评估不同BPD分级间的β多样性。ANOSIM P值经Bonferroni校正进行多重比较校正。通过R统计软件（版本4.4.2）中vegan包（版本2.6.0）实现的capscale函数，使用约束主坐标分析（CPCoA）进行Bray-Curtis距离排序。分析期间统计控制潜在临床混杂因素。首先使用PERMDISP检验评估组间离散均匀性，确保后续PERMANOVA结果有效性。使用随机森林插补方法（通过R中missForest包（版本1.5）实现）插补缺失数据。使用SPSS（版本26.0）进行统计分析。连续变量经Shapiro-Wilk检验正态性和Levene检验方差齐性后，使用Kruskal-Wallis H检验评估；分类变量通过Pearson卡方检验或Fisher精确检验比较。属水平差异丰度使用ANCOM-BC（偏倚校正微生物组组成分析）评估，能够获得跨BPD组各属的偏倚校正效应大小估计和置信区间。P值小于0.05认为有统计学意义。
网络分析中仅保留相对丰度总和大于0.01%且在5%样本（5/98）中存在的特征，以最小化潜在测序伪影和低丰度污染物影响，并增强后续成分相关性分析的数值稳定性。为减轻共现网络推断中因组大小不等产生的偏差并确保计算稳定性，对组间样本量进行标准化。每组随机下采样至N=16，匹配最小组大小。为考虑微生物组数据的成分性质，使用稀疏成分数据相关性（SparCC）算法（v0.1.0）构建微生物共现网络，这是在Python中实现的成分稳健相关性方法。SparCC相关性从标准化相对丰度数据计算，进行20次迭代以确保收敛。通过100次自助重采样评估相关显著性和边稳定性，仅保留显著边（P< 0.05）用于后续网络构建。网络构建采用|r| >0.3阈值定义显著微生物关联。该阈值基于微生物生态学中SparCC的既定实践及其在成分数据中捕获强生物学合理相互作用的适用性选择。通过多阈值敏感性分析（使用|r| > 0.2和|r| > 0.4阈值）验证结果稳健性。为评估方法选择影响，系统比较SparCC推断的网络与相对丰度上Spearman相关性生成的网络，使用Jaccard指数评估相似性。使用igraph分析网络拓扑属性。网络中ASV代表节点，节点大小表示相对丰度，ASV间相关性显示为边。本研究使用参数包括节点数（ASV数量）、边数（所有节点间连接数）、平均度（每个节点平均连接数）、平均路径长度（节点间最短路径平均步数）、直径（任意两节点间最长最短路径）、边密度（实际边数与最大可能边数之比）、聚类系数（网络中节点与其他节点关联程度）、模块性、介数中心性（节点作为其他两节点间最短路径桥梁的次数）、接近中心性（节点到所有其他节点距离的倒数）和特征向量中心性（中心性集中于少数主导节点的程度）。此外，计算网络中单个节点的度、介数中心性、接近中心性和同配性。
通过度、介数和接近中心性评估节点级中心性。微生物生态学广泛研究表明，共现网络中的关键物种往往表现出独特的拓扑模式，具体表现为高度中心性、低介数中心性、高接近中心性和高传递性。因此，根据其特定生态学解释为每个中心性度量分配差异权重以精炼识别。该复合关键得分前10%的类群被指定为关键类群。
基于识别的关键类群，进一步构建有序逻辑回归模型，将关键类群相对丰度与BPD严重程度关联。纳入基于临床重要性和差异变量选择的11个潜在混杂因素：性别、GA、BW、分娩方式、抗生素类型、产前皮质类固醇（ACS）、绒毛膜羊膜炎、胎粪污染羊水（MSF）、宫外生长受限（EUGR）、宫内窘迫、住院时间。对丰度数据进行对数转换并标准化为均值为0、标准差为1的变量，同时对数值变量进行标准化。由于样本量限制和多重共线性问题，采用前向逐步变量选择。从仅包含ASV的基线模型开始，按临床重要性顺序依次添加协变量。每步检查模型收敛性和数值稳定性，排除导致Hessian矩阵奇异的变量，最终保留保持模型稳定性的协变量组合。
使用Gephi可视化共现网络。利用网络相异系数评估不同BPD严重程度间整体网络相异性。进一步使用NetCoMi评估网络拓扑特征参数差异。Jaccard指数用于评估度和中心性相似性，范围0至1，值越高表示相似性越大。对每组进行1,000次随机排列，统计显著性阈值设为P<0.05。使用ComplexHeatmap包可视化各组网络拓扑特征差异。
使用R MASS包中polr函数开发多变量有序逻辑回归模型，采用前向逐步选择方法评估网络拓扑指标与BPD严重程度之间的关联。所有连续网络拓扑指标（包括边数、直径、平均路径长度、平均最近邻度、介数中心性、密度、度中心性、度同配性和传递性）使用Z分数标准化以消除尺度效应。模型纳入以下协变量作为潜在混杂因素：性别、GA、BW、分娩方式、抗生素类型、ACS、MSF、EUGR、宫内窘迫、住院时间。
RESULTS
临床特征
研究人群男性多于女性，非BPD组GA为30.57周（标准差1.405），Ⅰ级29.13周（1.700），Ⅱ级29.48周（1.758），Ⅲ级27.78周（2.425）。BPD患病率为68.4%（67/98），其中Ⅰ级31例（31.6%）、Ⅱ级20例（20.4%）、Ⅲ级16例（16.3%）。不同BPD严重程度婴儿的GA、BW和分娩方式存在显著差异。所有婴儿因肺炎接受抗生素单药或联合治疗，单药治疗仅使用β-内酰胺类，联合治疗指β-内酰胺类与其他抗生素（如糖肽类和大环内酯类）联用。大多数参与者（91.8%）因呼吸窘迫接受肺表面活性物质治疗。基于治疗类型，住院时间和抗生素治疗类型存在差异。不同结局婴儿在症状性动脉导管未闭（hsPDA）和EUGR频率上存在显著差异。
微生物组谱
98个样本经质量控制过滤和去除潜在人DNA污染后，获得6,182,749条测序读长用于16S rRNA分析，每个样本读长范围22,026至72,630。通过de novo过程推断出2,894个ASV，在22,000稀释深度后区分差异仅一个核苷酸的序列变异。但出生时收集样本的α和β多样性无显著差异。Chao1丰富度Kruskal-Wallis检验得R2=4.22（P=0.24），表明组间差异极小。同样，Shannon多样性R2=4.11（P=0.25），Simpson多样性R2=5.49（P=0.14），覆盖指数R2=4.71（P=0.19）。表明α多样性变异中不到5%可归因于组间差异。组间均值差异的95%置信区间均包含零，进一步支持出生时样本多样性无有意义变异。PERMDISP检验表明组内变异性相似（F=0.93，P=0.43），使得PERMANOVA检验适用。PERMANOVA分析显示BPD组对微生物组成的独特贡献无显著差异（F=0.89，P=0.64）。进一步检查个体协变量显示抗生素类型是微生物变异主要驱动因素，解释19.3%的变异（R2=0.19）。
所有样本中鉴定出的33个细菌门中，变形菌门相对丰度最高（均值87.9%，中位数98.7% [四分位距93.6–99.3]），其次是厚壁菌门（8.2%，0.4% [0.1–4.1]）、拟杆菌门（1.9%，0.3% [0.2–0.6]）和放线菌门（0.7%，0.2% [0.1–0.3]）。按BPD严重程度分组的四组间这些门相对丰度无显著差异。所有样本检测到的566个属中，Acinetobacter最丰富（41.5%，47.9% [34.9–54.0]），其次是Pseudomonas（19.3%，21.8% [17.1–24.8]）、Sphingomonas（7.6%，6.6% [3.5–8.8]）、Pelomonas（4.0%，1.6% [1.1–3.9]）和Stenotrophomonas（3.5%，3.4% [2.7–4.0]）。随BPD严重程度增加，Escherichia-Shigella相对丰度呈增加趋势，而Bacteroides呈下降趋势。Staphylococcus仅在非BPD组检测到，而Enterococcus在除Ⅲ级外的BPD组中相对丰度增加。值得注意的是，Aquabacterium和Devosia相对丰度在不同BPD严重程度间保持稳定，无统计学显著波动。ANCOM-BC分析鉴定出三个属在BPD严重程度组间存在差异丰度。Escherichia-Shigella在Ⅰ级组与非BPD组相比丰度显著增加（lfc=2.09，q=0.004）。Streptococcus在Ⅰ级（lfc=1.78，q=0.048）和Ⅱ级（lfc=2.09，q=0.028）组均显著增加，但在Ⅲ级组与Ⅰ级（lfc=-2.30，q=0.006）和Ⅱ级（lfc=-2.61，q=0.003）组相比显著减少。Chryseobacterium在Ⅲ级组与Ⅱ级组相比显著增加（lfc=1.83，q=0.011）。
微生物共现网络
进行网络分析以探索ASV水平下气道微生物组的微生物共现模式。鉴于样本量，以下网络比较作为与BPD严重程度相关微生物群落结构的探索性分析呈现。敏感性分析显示，健康对照组中共现网络的全局拓扑结构和关键类群识别对相关性阈值选择高度稳健。相比之下，BPD组网络显示出更大变异性，强调了疾病的生态影响。基于这些发现，所有下游分析采用中等阈值（|r| > 0.3，FDR < 0.05）。
SparCC网络识别出保守数量的强相关性（组间边数范围35至72），而Spearman方法检测到大量关联（边数范围2,044至7,989）。量化两种方法识别边相似性的Jaccard指数极小（0.003至0.023）。最终稳健的SparCC网络可视化显示。随后发现微生物共现网络复杂性在BPDⅢ级组急剧下降。确定|r| >0.3相关性阈值产生最稳健生态网络，因其最佳平衡了避免假阳性关联与保留有意义的生态信息。使用|r| >0.3阈值，识别非BPD组关键类群为Acinetobacter和Fusobacterium，Ⅰ级为Brevundimonas和Fusobacterium，Ⅱ级和Ⅲ级为Fusobacterium和Acinetobacter。对于四个关键物种，模型排除了引入时引起数值问题和收敛失败的抗生素类型变量。因此，所有模型包含其余10个协变量。关键物种丰度对BPD水平无显著影响：ASV 1303（OR=1.15，P=0.52）、ASV 708（OR=0.85，P=0.46）、ASV 1702（OR=0.98，P=0.94）和ASV 2215（OR=0.89，P=0.56）。
网络相异系数表明BPDⅢ级组微生物相互作用网络与其他组显著不同
