基于16S rRNA 測序分析阻塞性睡眠呼吸暫?；颊吣c道靶標菌群的變化

朱繼偉，盧曼路，焦倩倩，孫運良，劉 璐，丁紅紅，于 燕，潘 磊

濱州醫學院附屬醫院呼吸與危重癥醫學科，山東 濱州256603

阻塞性睡眠呼吸暫停（OSA）是一種常見的睡眠呼吸障礙類型，主要表現為睡眠時上呼吸道反復塌陷、氣流停止，導致間歇性缺氧和睡眠碎片化［1,2］。研究表明，全球目前有近10億人口患有OSA，其中我國患病人口最多［3］。OSA會導致氧化應激、炎癥反應、代謝紊亂等多系統功能障礙［4,5］，然而其機制目前尚不清楚。目前，臨床針對OSA患者治療方式大多數采用手術［6］或持續氣道正壓通氣治療［7］，患者診療費用較高，依從性較差［8］。

近年來，隨著高通量測序技術的發展，腸道菌群與OSA發生發展的關系也在逐漸被認識，OSA的發生發展與腸道菌群變化密切相關［9,10］。既往動物實驗研究發現間歇性低氧小鼠模型腸道菌群出現明顯變化，腸道菌群參與介導物質代謝和慢性炎癥［11-13］。臨床研究發現間歇性缺氧會改變腸道微生物組的生物節律［14］，OSA患者腸道菌群中厚壁菌門豐度升高，擬桿菌門豐度降低，炎癥相關細菌豐度明顯升高［15-17］。Baldanzi 等［18］研究發現OSA患者腸道微生物物種多樣性降低與呼吸暫停、缺氧密切相關。Ayyaswamy等［19］發現OSA誘導腸道微生態失調，導致出現神經免疫反應和腸道炎癥。有研究表明OSA患者腸道菌群變化誘導腸道通透性增加，高血壓合并OSA患者mTOR信號通路富集存在差異，與其病理生理變化密切相關［20］。目前涉及OSA患者的研究多為闡述OSA患者腸道菌群物種差異及與炎癥關系，關于人腸道靶標菌屬及其介導的相關差異代謝通路的研究尚未見報道。

本研究基于16S rRNA 高通量測序結果分析OSA患者及健康人群腸道菌群物種差異，基于多維度分析模型尋找驗證OSA患者腸道靶標菌屬，結合差異代謝通路分析標志菌群或腸道靶菌參與的代謝通路，以期OSA患者干預和治療提供策略。

1 資料和方法

1.1 研究對象

隨機選取2022年1～12月就診于本院呼吸睡眠醫學中心的39例OSA患者作為OSA組，同期檢查無OSA的20例健康志愿者作為對照組。本研究通過濱州醫學院附屬醫院倫理委員會批準（20220128-76），所有研究對象均知情同意。

1.2 納入、排除標準

OSA 組納入標準：經整夜多導睡眠監測初診為OSA的患者，根據美國睡眠醫學會睡眠及相關事件判讀診斷標準，計算呼吸暫停低通氣指數（AHI）為呼吸暫停和低通氣的總發次數除以全部睡眠時間的總持續時間，AHI≥5次/h可診斷為OSA［21］；未進行持續氣道正壓通氣治療。排除標準：年齡＜18歲；服用影響睡眠的藥物，如苯二氮?類藥物、巴比妥類藥物或鎮靜劑；接受相關治療，包括但不限于手術、呼吸機治療；合并胃腸道疾病，如潰瘍性結腸炎和腸易激綜合征等；入組前2月內服用腸道微生態制劑、抗生素或免疫抑制劑的患者；存在腫瘤、糖尿病、肝腎功能不全、風濕性免疫疾病或其他可能影響腸道菌群失調的疾病者。

對照組納入標準：年齡＞18歲；無睡眠呼吸暫停、失眠及其他睡眠障礙疾病病史。排除標準：入組前2 月內有服用微生態制劑、抗生素或免疫抑制劑者；入組前1月內有較大飲食結構或生活環境改變；存在消化系統疾病、腫瘤、糖尿病、肝腎功能不全、風濕性免疫疾病或其他可能影響腸道菌群失調的疾病者。

1.3 研究方法

1.3.1 樣本采集與DNA提取 采集受試者清晨糞便樣本0.5 g，儲存于無菌無酶的糞便采樣管中，避免尿液、月經等污染，標記好信息后，迅速轉移至-80 ℃冰箱保存。所有糞便標本的凍融參照Wang等［22］的方法，采用Omega Mag-Bind Soil DNA Kit（200）（Omega）提取糞便樣本中總微生物DNA，采用0.8%瓊脂糖凝膠電泳分析樣本DNA 相對分子質量和純度，利用NanoDrop NC2000（Thermo Fisher Scientific）對DNA進行定量。

1.3.2 文庫構建及基因測序 以微生物核糖體RNA等能夠反映菌群組成和多樣性的目標序列為靶點，應用16S rRNA V3-V4 區域特異引物進行聚合酶鏈式反應（PCR）擴增。PCR擴增選用細菌16S rDNA V3-V4區特異性引物，上游引物338F（5'-barcode +ACTCCTACGGGAGGCAGCA-3'），下游引物806R（5'-GGACTACHVGGGTWTCTAAT-3'）。前引物中的barcode 是一個7 個堿基的寡核苷酸序列，用來區分同一文庫中的不同樣品。擴增結果進行2%瓊脂糖凝膠電泳，然后用凝膠回收試劑盒（AXYGEN）回收目的片段。采用Illumina 公司的TruSeq Nano DNA LT Library Prep Kit建立測序文庫。在Illumina MiSeq平臺上利用MiSeq Reagent Kit V3（600cycles）對16S rRNAV3-V4區進行2×250 bp的雙端測序。

1.4 生物信息學分析

1.4.1 腸道菌群測序質量評估 從每個樣本中隨機抽取一定數量的序列數據，繪制稀疏曲線，將各分組的分類操作單元（OUT）按其豐度從大到小沿橫坐標依次排列后，將豐度值經Log10對數轉換后的值作為縱坐標，在R軟件中編寫腳本繪制各樣本或各分組的豐度等級曲線，反映測序數據量的合理性及豐富度，當曲線趨向平緩時，測序質量良好，測序數據量充足，可以用于后續分析。

1.4.2 物種組成分析 統計抽平后樣本擴增子序列變異（ASV），獲得每個樣本的微生物群落的具體組成，對兩組在門、屬水平的物種進行可視化，使用circos圖用來展示樣本與物種的豐度關聯關系。

1.4.3 Alpha多樣性分析 使用QIIME 2進行Alpha多樣性分析，以Shannon 指數和Simpson 指數表征多樣性，以Chao1指數和Observed species指數表征菌群豐富度，以Pielou指數表征均勻度。

1.4.4 Beta多樣性分析 通過主坐標分析(PCoA)、非度量多維尺度分析(NMDS)非約束排序手段對多維的微生物數據進行降維，展示數據變化的主要趨勢?；趙eighted UniFrac距離矩陣分析，通過二維排序圖展示微生物群落的組成差異。

1.4.5 線性判別分析及影響因子（LEfSe）分析 LEfSe分析［23］對所有分類水平進行差異分析，識別組內組間物種豐度差異，來檢測不同組間的差異豐富分類群，區分不同生物學特征的物種，尋找分組之間穩健的差異標志物種。

1.4.6 關聯網絡分析 采用ASV在本項目所有樣本中的豐度數據，過濾掉序列總數少于10、出現樣本數少于5 的ASV，采用SparCC算法，構建相關性矩陣，使用隨機矩陣理論確定相關性數值的過濾閾值。使用subgraph函數過濾共現網絡計算各樣本子網絡的圖形水平拓撲學指數。根據Erdos Renyi模型以及當前網絡的節點數目和邊數目，構建隨機網絡。分別計算當前網絡、隨機網絡中每個節點的相連節點數（度）。根據當前共現網絡的模塊化切割結果，計算當前網絡中每個節點的Zi、Pi score值。根據Zi、Pi score值確定每個節點在關聯網絡中的角色［24］。

1.4.7 隨機森林分析 使用未抽平的ASV，調用q2-sample-classifier中的“Classify_samples_ncv”函數進行隨機森林分析以及巢式分層交叉檢驗［25］。進行10倍交叉檢驗，確定各物種在組間鑒別的重要性。

1.4.8 代謝通路預測分析 構建微生物基因組16S rRNA進化樹，推測特征序列的最近序列物種，結合各樣本特征序列豐度，計算基因家族拷貝數。將基因家族映射到MetaCyc數據庫中，獲得各樣本中代謝通路的豐度數據。使用歸一化的pathway豐度表，計算第二層次通路的平均豐度。運用MetagenomeSeq分析調用fitFeatureModel 函數使用zero-inflated log-normal model 對每個pathway的分布進行擬合，并使用該模型的擬合結果判別差異的顯著性。依據選擇的通路抓取分層的樣本代謝通路豐度表中的相應數據，繪制物種組成柱狀圖。

1.5 統計學分析

采用SPSS26.0軟件對數據進行統計學分析。用QIIME 2 和R（3.2.0）軟件進行序列數據分析。符合正態分布的計量資料采用均數±標準差表示，兩組間比較采用獨立樣本t檢驗；非正態分布資料應用非參數檢驗分析，兩組間比較采用χ2檢驗。以Ρ＜0.05為差異具有統計學意義。

2 結果

2.1 臨床資料比較

本研究最終獲得59例有效樣本，其中OSA組樣本39例，對照組樣本20例，兩組人群性別、年齡差異均無統計學意義（Ρ＞0.05），與對照組相比，OSA組BMI、AHI指數明顯升高，夜間平均血氧降低（Ρ＜0.01，表1）。

表1 OSA組和健康對照組一般資料Tab.1 General data of OSAgroup and healthy control group(Mean±SD)

2.2 腸道菌群測序質量評估

兩組樣本的稀釋曲線及豐度等級曲線均趨向平坦，說明物種豐富度足夠，測序樣本量及數據量漸進合理（圖1）。腸道菌群測序質量良好，測序數據量充足，可以用于后續分析。

圖1 稀釋曲線及豐度等級曲線Fig.1 Rarefaction curve and RankAbundance curves.OUT:Operational taxonomic unit.

2.3 物種組成分析

對受試者樣本進行16S rRNA 基因測序，得到有效序列總數為12 615 425。在門水平上，測序結果顯示兩組菌群主要是由厚壁菌門（Firmicutes）、擬桿菌門（Bacteroidetes）、變形菌門（Ρroteobacteria）、放線菌門（Actinobacteria）4類菌群組成，與對照組相比，OSA組Firmicutes、Ρroteobacteria、Actinobacteria豐度比例升高，而Bacteroidetes豐度比例降低（圖2A）。在屬水平上，樣本菌群測序結果相對豐度前10的菌群依次為擬桿菌屬（Bacteroides）、糞桿菌屬（Faecalibacterium）、羅斯氏菌屬（Roseburia）、雙歧桿菌屬（Bifidobacterium）、志賀氏菌屬（Shigella）、布勞特氏菌屬（Blautia）、巨單胞菌屬（Megamonas）、瘤胃球菌屬（Ruminococcus）、普雷沃氏菌屬（Ρrevotella）、戴阿里斯特桿菌屬（Dialister），與對照組相比，OSA組Faecalibacterium、Roseburia豐度明顯降低，Megamonas豐度明顯升高（圖2B）。

圖2 OSA組與對照組物種組成分布Fig.2 Composition and distribution of TOP10 species in terms of relative abundance in OSA group and control group.A:Phylum;B:Genus.

2.4 腸道微生物Alpha多樣性分析

Alpha多樣性分析結果顯示，OSA組物種多樣性Shannon 和Simpson 指數低于健康對照組（Ρ＜0.01）；OSA組物種豐富度Chao1 和Observed species低于健康對照組，Observed species兩組間差異具有統計學意義（Ρ＜0.05）。OSA組菌群均勻度低于對照組（Ρ＜0.05，表2）。

表2 OSA組和對照組的Alpha多樣性分析Tab.2 Alpha diversity analysis in OSAgroup and control group(Mean±SD)

2.5 腸道微生物Beta多樣性分析

PCoA分析結果顯示，OSA組與對照組相距較遠，兩樣本差異較大（Ρ＜0.05，圖3A）；NMDS分析應力值Stress=0.118（圖3B），當應力值＜0.2時差異分析結果可靠，具有統計學意義（Ρ＜0.05）。

圖3 OSA組和對照組的Beta多樣性分析Fig.3 Beta diversity analysis in OSA group and control group.A:Principal co-ordinate analysis;B:Nonmetric multidimensional scaling(NMDS).

2.6 物種差異與標志物種分析

2.6.1 組間差異分析 Anosim分析顯示，兩組腸道菌群組間差異大于組內差異，差異有統計學意義（R=0.268，Ρ=0.001，圖4）。

圖4 OSA組和對照組的Anosim分析Fig.4 Anosim analysis in OSAgroup and control group.

2.6.2 物種差異分析 在門水平上，Firmicutes、Ρroteobacteria、梭桿菌門（Fusobacteria）、TM7、酸桿菌門差異有統計學意義（Ρ＜0.05，圖5A）。在屬水平上，OSA組和健康對照組在26個菌屬中的差異有統計學意義（Ρ＜0.05），豐度較高的前7位菌屬中，OSA組假單胞菌屬（Ρseudomonas）、Megamonas、Ruminococcus、毛螺菌科（屬）（Lachnospiraceae）、多爾氏菌屬、丁酸球菌屬高于對照組（Ρ＜0.05，圖5B）。

圖5 組間物種差異分析Fig.5 Species difference analysis at the genus level.A:Phylum;B:Genus.

2.6.3 LEfSe分析 LEfSe分析柱狀圖及進化分支圖結果顯示，兩組腸道菌群中23種菌群豐度比較，差異有統計學意義（Ρ＜0.05）。其中OSA 組韋榮氏球菌科（Alicyclobacillaceae）、Megamonas、假單胞菌科（Ρseudomonadaceae）、假單胞菌目（Ρseudomonadales）、Ρseudomonas、腸球菌科（Enterococcaceae）、腸球菌屬（Enterococcus）、梭桿菌綱（Fusobacteriia）、梭桿菌目（Fusobacteriales）、梭桿菌門（Fusobacteria）、梭桿菌屬（Fusobacterium）等17 種腸道菌群豐度高于對照組（Ρ＜0.05）；對照組棲糞桿菌屬（Faecalibacterium）、桿孢囊菌屬（Virgisporangium）、RF3、1_68、脂環酸芽孢桿菌屬（Alicyclobacillus）、脂環酸芽孢桿菌科（Alicyclobacillaceae）腸道菌群豐度高于OSA 組（Ρ＜0.05，圖6）。

圖6 兩組人群腸道菌群LEfSe 分析柱狀圖及進化分支圖Fig.6 LEfSe analysis histogram and evolutionary branching diagram of gut microbiota in the two groups.LDA:Linear discriminant analysis.

2.6.4 關聯網絡分析 構建相關性矩陣確定閾值為0.565（圖7A），拓撲指數分析顯示，腸道菌群度分布結果符合無標度網絡，說明影響宿主穩態的主要因素為差異性標志菌群（圖7B）。ZiPi圖結果顯示，在門水平上厚壁菌門作為主要影響宿主穩態的關鍵菌群（圖7C），在屬水平上影響宿主穩態的關鍵物種為棲糞桿菌屬、志賀菌屬、假單胞菌屬（圖7D）。

圖7 關聯網絡分析Fig.7 Analysis of the association networks.A: Correlation matrix threshold;B: Topological exponential analysis;C:ZiPi;D:Key species that affect the host homeostasis.

2.6.5 差異標記物種分析 隨機森林模型分析結果顯示Pseudomonas在鑒別兩組中的重要性最高（圖8A）。在屬水平上對LEfSe分析差異物種進行受試者工作特征（ROC）曲線分析，根據曲線下面積評估差異性菌群診斷OSA的敏感性。其中Pseudomonas AUC=0.798[95%CI（0.674-0.921），Ρ=0.0003]、Megamonas AUC=0.678[95%CI（0.537-0.818），Ρ=0.029]，可作為重要鑒別意義的生物標記物（圖8B）。

圖8 差異標記物種分析Fig.8 Analysis of the differentially labeled species.A:Random forest model;B:ROC analysis.

2.7 差異生物功能代謝通路功能預測分析

基于MetaCyc數據庫，通過PICRUSt2預測分析顯示，在宿主體內腸道菌群穩態起主要作用的腸道菌群的功能是生物合成功能，起主要作用的為氨基酸生物合成、核苷酸生物合成等。其他預測功能包括降解、利用、同化、解毒、前體代謝產物和能量的產生、聚糖途徑、大分子修飾、代謝簇（圖9A）。MetagenomeSeq分析結果顯示，在芳香族生物胺降解、酮葡萄糖酸代謝（Ρ＜0.001），葡萄糖降解、沒食子酸酯降解、L-精氨酸和L-鳥氨酸降解超途徑、L-酪氨酸降解I、腺苷鈷胺生物合成II、L-亮氨酸降解I、甲苯降解IV（Ρ＜0.05，圖9B）。物種組成結果顯示Pseudomonas可能是影響物質通路代謝的關鍵菌屬（圖9C、D）。

3 討論

本研究基于16S rRNA高通量測序，分析OSA患者與健康對照組腸道菌群物種差異，發現兩組人群腸道菌群在組成結構和物種多樣性方面存在顯著差異。OSA組潛在致病菌屬（如假單胞菌屬、巨單胞菌屬）豐度增加，多維度分析模型及差異代謝通路分析發現假單胞菌屬可能通過參與芳香族生物胺降解和酮葡萄糖酸代謝通路影響OSA的發生發展，可望作為OSA患者腸道菌群的潛在靶標。

基于Shannon稀釋曲線及豐度等級曲線分析結果表明本研究測序數據量完整，測序質量良好，研究結果可靠。α多樣性分析結果顯示OSA組物種多樣性、豐富度、均勻度均低于健康對照組，β多樣性分析和Anosim分析結果表明兩組間菌群結構存在差異且組間差異大于組內差異。物種組成差異分析結果顯示，OSA組厚壁菌門、梭桿菌門、假單胞菌屬、瘤胃腸球菌屬等豐度明顯升高，產短鏈脂肪酸豐度明顯降低。這與Collado等［26］研究結果相似，研究發現厚壁菌門/擬桿菌門比例明顯升高，OSA患者腸道微生物群物種豐富度明顯低于對照組。有研究同樣證明OSA兒童患者腸道菌群組成發生明顯改變，產短鏈脂肪酸菌屬豐度降低［27］。結合既往研究，表明OSA患者腸道菌群穩態發生改變，菌群微生態紊亂失衡。

為進一步確定影響OSA患者腸道菌群穩態的靶標菌屬，本研究運用多維度模型分析尋找驗證OSA患者腸道菌群生物標記物。LEfSe分析結果表明在屬水平上OSA組假單胞菌屬、巨單胞菌屬豐度顯著升高。關聯網絡分析結果顯示，影響宿主穩態的關鍵物種為棲糞桿菌屬、志賀菌屬、假單胞菌屬等。同時隨機森林模型分析機器學習算法及ROC曲線分析結果同樣顯示假單胞菌屬在組間差異菌屬中重要性和敏感度最高。這與Li等［28］研究結果相似，在OSA患者腸道菌群中同樣假單胞菌屬豐度升高，假單胞菌屬豐度和D-乳酸呈正相關，D-乳酸是組織缺氧的敏感指標，間接證明假單胞菌屬和組織缺氧密切相關。Ren等［29］研究發現假單胞菌屬與腸道屏障功能障礙和感染有關，可能導致敗血癥和多器官功能障礙綜合征。需要注意的是，有研究發現瘤胃球菌屬同樣與腸道屏障功能密切相關，本研究OSA組瘤胃球菌屬豐度明顯升高［30］。Wen等［31］研究發現假單胞菌屬感染會引起機體的Th17細胞反應和系統性炎癥反應。同樣有研究發現假單胞菌屬與炎癥性腸病發病密切相關，通過影響腸道菌群的穩態來導致腸道免疫系統出現異常反應［32］。有研究發現假單胞菌屬參與腸-肺軸免疫調節，影響呼吸系統疾病的發生發展［33］。結合上述研究說明假單胞菌群可能通過介導免疫炎癥反應，參與腸-肺軸免疫調節導致OSA患者的腸道屏障功能受損和代謝紊亂。

在差異生物功能代謝通路功能預測分析中，二級功能通路分析顯示在生物合成途徑的相對豐度最高，占比最多的是氨基酸生物合成。組間顯著差異通路分析結果顯示OSA組芳香族生物胺降解、酮葡萄糖酸代謝通路顯著上調。分層樣本代謝通路分析結果發現在這兩個通路中起主要作用的均為假單胞菌屬。研究發現芳香族生物胺是生物體內活性成分，在物質合成及生物膜穩定性方面起重要作用，具有抗氧化及調節腸道生理功能的作用［34-36］。因此，OSA患者芳香族生物胺降解通路豐度升高，芳香族生物胺降解增多，影響OSA患者物質合成、物質代謝及抗氧化能力。Sun等［37］研究發現假單胞菌屬與酮葡萄糖代謝密切相關，在缺氧條件下假單胞菌屬、巨單胞菌屬通過酮葡萄糖酸代謝通路氧化生成2-酮基-D-葡萄糖酸，利用2-酮基-D-葡萄糖酸對細胞生長及代謝產生負面影響。因此，推測假單胞菌屬可能通過上述代謝通路，介導OSA患者的病理生理學改變，有助于解釋腸道菌群與OSA患者的內在關系，有望為干預OSA提供新的策略。

綜上所述，本研究發現與健康人群相比，OSA患者腸道菌群的整體結構失衡，物種組成紊亂且有差異，考慮差異與OSA的病理生理如間歇性缺氧及睡眠碎片化等因素有關；假單胞菌屬可作為重要鑒別意義的生物標記物，作為腸道靶標菌群可能參與芳香族生物胺降解和酮葡萄糖酸代謝通路，進而影響OSA的發生發展。另外，本研究存在一定局限性：首先，兩組存在性別差距，OSA組男性患者較多，這也與OSA的發病人群有關［38］；其次，兩組年齡和BMI有一定差異，可能影響兩組腸道菌群物種組成，后續研究將進一步探討年齡和BMI對OSA腸道菌群的影響。最后，腸道菌群容易受到環境、遺傳、飲食行為習慣等因素的影響，本研究地域單一，且受限于方法學，本研究所取樣品可能無法反映腸道內厭氧菌狀態。然而有關OSA腸道菌群研究方興未艾，尤其專注腸道菌群代謝通路研究，將為OSA的早期干預和治療提供新思路。