基于16S rRNA 高通量測序技術(shù)的腸道菌群多樣性推斷死亡時間

曹潔，李文晉，王一飛，安國帥，盧曉軍，2，杜秋香，李晉，孫俊紅

1.山西醫(yī)科大學(xué)法醫(yī)學(xué)院，山西 太原 030001；2.包頭市公安局刑事偵查支隊，內(nèi)蒙古 包頭 014030；3.山西醫(yī)科大學(xué)第二醫(yī)院，山西 太原 030001

死亡時間又稱死后間隔時間（postmortem interval，PMI），通常指死亡發(fā)生到尸體被發(fā)現(xiàn)的間隔時間[1]。死亡時間推斷是法醫(yī)學(xué)鑒定中需要解決的重要問題之一，對確定案發(fā)時間、認(rèn)定和排除嫌疑人、劃定偵查范圍具有重要意義。傳統(tǒng)的PMI 推斷主要依據(jù)死后尸體現(xiàn)象變化規(guī)律粗略判斷。隨著分子生物學(xué)技術(shù)的發(fā)展，實時熒光定量聚合酶鏈反應(yīng)（real-time fluorescent quantitative polymerase chain reaction，qPCR）技術(shù)[2]、雙向電泳（two-dime nsional electrophoresis，2-DE）和高效液相色譜-質(zhì)譜（highperformance liquid chromatography-mass spectrometry，HPLC-MS）技術(shù)[3]、傅里葉變換衰減全反射紅外光譜（attenuated total reflectance-Fourier transform infrared，ATR-FTIR）技術(shù)[4]等越來越多的新型技術(shù)被應(yīng)用到PMI 推斷研究中。大量相關(guān)研究為尋找便捷快速的技術(shù)方法、敏感的指標(biāo)以及建立合適的數(shù)學(xué)模型以縮小PMI 窗口提供了實驗基礎(chǔ)和新的思路。

微生物在周圍環(huán)境和個體中無處不在，生態(tài)環(huán)境中的微生物在地球化學(xué)循環(huán)中發(fā)揮著重要的生態(tài)作用，人體內(nèi)的微生物更是廣泛參與免疫、消化、代謝等生理過程[5]。腸道微生物在食物分解、藥物降解、脂肪積累、維生素和氨基酸合成中扮演著重要的角色[6]。腸道微生物不僅參與體內(nèi)的穩(wěn)態(tài)代謝，而且影響尸體分解的發(fā)展和進(jìn)程，在尸體腐爛過程中發(fā)揮著重要作用[7]。機體死亡后，微生物既幫助分解機體組織，又通過分解代謝完成物質(zhì)的自然循環(huán)[5]。

既往對腸道微生物的研究主要基于培養(yǎng)法，而能夠培養(yǎng)的微生物僅占微生物種類總數(shù)的1%[8]。隨著高通量測序技術(shù)的成熟與發(fā)展，細(xì)菌群落結(jié)構(gòu)分析被廣泛應(yīng)用。不同環(huán)境中的微生物多樣性被逐一解讀，大樣本、高通量的測序可對微生物變化進(jìn)行動態(tài)的、系統(tǒng)的分析。有研究[9]表明，尸體降解過程中腸道厚壁菌門（Firmicutes）和擬桿菌門（Bacteroidetes）的數(shù)量減少，而變形菌門（Proteobacteria）的數(shù)量增加，提示腸道微生物群落變化具有較強的時間變化規(guī)律，而腸道微生物在尸體腐敗破壞之前，處于相對封閉的環(huán)境，其菌群的變化主要與PMI 相關(guān)，有望成為精確推斷PMI 的有效方法。

本研究收集了大鼠死后30 d 內(nèi)不同時間點盲腸內(nèi)容物，采用16S rRNA 高通量測序技術(shù)探究大鼠腸道菌群與PMI 之間的關(guān)系，從多方面分析比較各時間點間差異，建立偏最小二乘（partial least squares，PLS）回歸模型，尋找時序性變化特征，從而更準(zhǔn)確地推斷PMI。

1 材料與方法

1.1 大鼠死亡模型制備及檢材提取

健康成年雄性Sprague-Dawley大鼠84只，10～12周齡，體質(zhì)量250～300 g，由山西醫(yī)科大學(xué)實驗動物中心提供。溫室飼養(yǎng)2 d 后，采用斷頸方式處死大鼠，死后存放于16 ℃人工氣候箱。于死后0、1、2、3、5、7、9、12、15、18、21、24、27 和30 d（每個時間點6 只大鼠）共14 個時間點分別取大鼠盲腸內(nèi)容物置于干凈的微量離心管內(nèi)，立即放入液氮中冷凍，隨即將樣本保存至-80 ℃冰箱待檢。

本實驗由山西醫(yī)科大學(xué)科學(xué)研究倫理審查委員會批準(zhǔn)（審批號2018LL351）。

1.2 DNA 提取及濃度、純度檢測

取200 mg 盲腸內(nèi)容物于新微量離心管中，使用QIAamp PowerFecal Pro DNA 試劑盒（德國Qiagen 公司）提取總DNA，具體操作步驟參照試劑盒說明書。使用Infinite?200 PRO 酶標(biāo)儀（瑞士TECAN 公司）檢測總DNA 純度及濃度，光密度D260/D280在1.8～2.0 的DNA 樣本用于后續(xù)實驗。

1.3 文庫構(gòu)建及高通量測序

DNA 質(zhì)檢達(dá)要求后，以此為模板，采用帶barcode的引物對16S rRNA 的V3～V4區(qū)進(jìn)行第一輪PCR 擴增，擴增引物為341F（5'-CCTACGGGNGGCWGCAG-3'）和805R（5'-GACTACHVGGGTATCTAATCC-3'），反應(yīng)體系30 μL：15 μL 即用型的常規(guī)PCR 預(yù)混合溶液（上海翊圣生物科技有限公司），上下游引物各1 μL，10～20 ng DNA 模板，9～12 μL ddH2O。隨后引入二代測序接頭引物試劑盒（Illumina 平臺Index Primer13-24，北京百奧萊博科技有限公司）進(jìn)行第二輪PCR 擴增，反應(yīng)體系30 μL。用20 g/L 瓊脂糖凝膠電泳檢測文庫大小，檢測合格后，為了得到均勻的長簇效果和高質(zhì)量的測序數(shù)據(jù)，使用Qubit?3.0 熒光定量儀（美國Thermo Fisher Scientific 公司）進(jìn)行文庫濃度測定。通過Illumina Miseq PE300 測序平臺（美國Illumina 公司）對DNA 片段進(jìn)行雙端高通量測序，獲得原始數(shù)據(jù)。

1.4 數(shù)據(jù)預(yù)處理

得到原始數(shù)據(jù)后，用Cutadapt v1.2.1（https://cutadapt.readthedocs.io/en/stable/）軟件去除接頭序列，將成對的短序列用PEAR v0.9.6（https://cme.h-its.org/exelixis/web/software/pear/）軟件依據(jù)重疊關(guān)系進(jìn)行數(shù)據(jù)拼接。采用PRINSEQ v0.20.4（http://prinseq.sourceforge.net/）軟件切除barcode 序列和引物序列，并采用滑窗法檢查各個序列的質(zhì)量值，以去除含N 部分序列、短序列以及低復(fù)雜度序列，長度閾值為200 bp。使用USEARCH v8.1.1831 軟件（www.drive5.com/usearch）去除預(yù)處理后序列中非擴增區(qū)域序列，并用UCHIME v4.2.40 軟件（http://drive5.com/usearch/manual/uchime_algo.html）檢測嵌合體并予以刪除，隨后去除嵌合體的序列與核糖體數(shù)據(jù)庫項目（Ribosomal Database Project，RDP；http://rdp.cme.msu.edu/misc/resources.jsp）代表性序列進(jìn)行比對，剔除相似度低于80%的序列，得到合格序列。USEARCH v8.1.1831軟件將上述數(shù)據(jù)進(jìn)行97% 相似性運算分類單元（operational taxonomic unit，OTU）聚 類，采用RDP Classifier v2.12 軟件（https://sourceforge.net/projects/rdp-classifier/files/）將各樣本序列進(jìn)行物種分類，選擇豐度最高的序列作為OUT 代表性序列。

1.5 統(tǒng)計分析

在門、屬水平對微生物群落構(gòu)成及分布豐度進(jìn)行分析。用R v3.2（https://www.r-project.org/）對物種分類學(xué)統(tǒng)計結(jié)果做柱狀圖，云工具（https://www.lc-bio.cn，杭州聯(lián)川生物技術(shù)股份有限公司）繪制?；鶊D觀察微生物群落構(gòu)成?；谇笆鯫TU 聚類結(jié)果進(jìn)行α多樣性分析[10]，使用Mothur v1.30.1 軟件（https://mothur.org/）計算各個樣本α多樣性指數(shù)值（ACE 指數(shù)和Shannon 指數(shù)），檢驗水準(zhǔn)α=0.05。線性判別分析效應(yīng)大小（linear discriminant analysis effect size，LEfSe）是一種線性判別分析（linear discriminant analysis，LDA）方法[11]，結(jié)合Wilcoxon 和Kruskal-Wallis 非參數(shù)檢驗可以篩選出最能解釋組間差異的微生物菌群，以及其對組間差異的影響程度，本研究使用LEfSe v1.1.0 軟件（http://huttenhower.sph.harvard.edu/lefse/）選擇各時間點LDA 分值＞2 的差異細(xì)菌群落。利用SIMCA?14.1（Q845）軟件（德國Sartorius 公司）建立PLS 回歸模型[12]，得到?jīng)Q定系數(shù)（R2）和交叉驗證均方根誤差，分析PMI 與腸道菌群之間的線性關(guān)系。

2 結(jié)果

2.1 樣本測序結(jié)果

本研究對大鼠死后14 個時間點共84 個樣本進(jìn)行測序，獲得原始數(shù)據(jù)8 282 115 條，經(jīng)數(shù)據(jù)預(yù)處理獲得有效數(shù)據(jù)為6 591 968 條，每條序列平均長度為（416.11±4.09）bp。根據(jù)相似性將上述序列聚集成為23 759 個OTU，采用RDP Classifier v2.12 軟件對OTU代表序列進(jìn)行物種分類，共得到28個門、53個綱、87個目、163個科、347個屬，用于后續(xù)不同水平的菌落組成分析。

2.2 門、屬水平微生物群落構(gòu)成及分布豐度分析

在門、屬水平分別選擇豐度最高的前20 個（＞1%）物種做柱狀圖進(jìn)行分析（圖1A～B）。門水平（圖1A），厚壁菌門在0 d 平均相對豐度最高，達(dá)到近80%，隨后整體呈下降趨勢，但在各時間點相對豐度始終最高；擬桿菌門相對豐度在死后0～7 d 上升，9 d 后下降，變形菌門在死后12、18、21 d 相對豐度高于其他時間點。厚壁菌門、擬桿菌門、變形菌門在整個樣本中的相對豐度最高，其余前10 大菌門分別為放線菌門（Actinobacteria）、疣微菌門（Verrucomicrobia）、糖念珠菌門（Candidatus Saccharibacteria）、螺旋體門（Spirochaetes）、廣古菌門（Euryarchaeota）、浮霉菌門（Planctomycetes）、酸桿菌門（Acidobacteria），在分解過程中所有時間點相對豐度都比較低。屬水平（圖1B），乳酸菌（Lactobacillus）在0、3、5、7、9、12、15、21、24、27、30 d 相對豐度最高，其余前10 大菌屬分別為羅姆布茨菌屬（Romboutsia）、腸球菌屬（E.enterococcus）、擬桿菌（Bacteroides）、雙歧桿菌（Bifidobacterium）、巴內(nèi)氏菌屬（Barnesiella）、阿克曼氏菌屬（Akkermansia）、梭狀芽孢桿菌XlVa（ClostridiumXlVa）、梭狀芽孢桿菌（Clostridium sensu stricto）、假單胞菌（Pseudomonas），其相對豐度都比較低。

圖1 不同水平微生物群落相對豐度Fig.1 Relative abundance of microbial communities at different levels

?；鶊D（圖1C）中，不同顏色彩色條帶代表不同菌屬，可以直觀地看到相對豐度前10 的菌屬都屬于相對豐度前5 的菌門，左邊與各時間點相連的條帶端點寬度表示菌屬在該樣本中的豐度，不僅反映了每個時間點的優(yōu)勢菌群組成比例，同時也反映了各優(yōu)勢菌群在不同時間點之間的分布比例。

2.3 α 多樣性分析

由表1 可知，除死后5 天微生物群落總數(shù)升高（P＜0.05）外，其他PMI 點在數(shù)量上差異無統(tǒng)計學(xué)意義；與0 d 相比，死后9、12、15、24、27、30 d 的Shannon指數(shù)降低（P＜0.05）。

表1 腸道微生物菌群α 多樣性分析Tab.1 Alpha diversity analysis of intestinal microbial communities（n=6，xˉ ±s）

2.4 細(xì)菌群落結(jié)構(gòu)差異統(tǒng)計分析

圖2 為LDA 分值分布柱狀圖，對14 個時間點的細(xì)菌群落進(jìn)行比較，在13 個時間點（除第24 天）共篩選出119 個具有差異的細(xì)菌群落。

圖2 線性判別分析效應(yīng)Fig.2 LDA effect

2.5 PLS 回歸模型建立

根據(jù)全部時間點實際值與預(yù)測值建立的PLS 模型（圖3A），其R2為0.795，交叉驗證均方根誤差為6.57 d。根據(jù)9 d 前樣本和12 d 后樣本分別建立PLS回歸模型（圖3B～C），兩者R2、交叉驗證均方根誤差分別為0.767、1.96 d 和0.445、5.37 d。

圖3 PLS 回歸分析Fig.3 PLS regression analysis

3 討論

隨著16S rRNA 高通量測序技術(shù)的成熟與發(fā)展，腸道菌群的研究在疾病診斷、臨床治療等[11，13-14]方面均有應(yīng)用。近年來，有學(xué)者將死后微生物的時序性變化用于PMI 推斷的研究。METCALF 等[5，15]研究了小鼠和人尸體皮膚、腹腔、周圍土壤在不同季節(jié)下微生物隨PMI 變化的關(guān)系，建立回歸模型推斷PMI。LIU 等[16]評估了死后15 d 內(nèi)微生物群落結(jié)構(gòu)，基于隨機森林（random forest，RF）、支持向量機（support vector machine，SVM）、人工神經(jīng)網(wǎng)絡(luò)（artificial neural network，ANN）等算法建立模型推斷PMI。本研究通過分析大鼠死后30 d 內(nèi)不同時間點腸道菌群測序結(jié)果，探索其與PMI 之間的關(guān)系。

桑基圖是樣本與物種的共現(xiàn)性關(guān)系圖，描述樣本與物種間對應(yīng)關(guān)系的可視化圖，從圖中可見相對豐度最高的前10 個菌屬都屬于豐度最高的前5 個菌門，反映每個樣本的優(yōu)勢物種組成比例，同時也反映各優(yōu)勢物種在不同樣本之間的分布比例。本研究結(jié)果顯示，在門水平厚壁菌門、擬桿菌門、變形菌門相對豐度最高，且各時間點α多樣性指數(shù)表明，隨時間推移，微生物群落總數(shù)未發(fā)生明顯變化，而菌群多樣性卻在多個死后時間點明顯升高（圖2）。上述針對各采樣時間點微生物群落在門、屬水平發(fā)生的變化，與LI 等[17]的研究結(jié)果基本一致，但本研究探索了死后30 d 內(nèi)更多時間點的腸道菌群變化規(guī)律。

LEfSe 分析可得到LDA 分值分布柱狀圖和組間差異具有統(tǒng)計學(xué)意義的差異物種在不同組別的豐度。本研究發(fā)現(xiàn)，在14 個時間點中除第24 天外其余13 個時間點共篩選出119 個具有差異的細(xì)菌群落，第24 天沒有差異的細(xì)菌群落的原因可能為隨著時間的推移，菌群多樣性下降，沒有LDA 分值＞2 的差異菌群。

如何通過合理的建模方式提高PMI 推斷的準(zhǔn)確性是一個亟待解決的問題，PLS 回歸是目前最常用的建模方法[12，18-19]。為明確OTU 水平與PMI 的相關(guān)性，本研究建立了PLS 回歸的PMI 推斷模型，結(jié)果顯示，R2為0.795，但交叉驗證均方根誤差為6.57 d。若將間隔時間較短的9 d 前和間隔時間較長的12 d 后分別做PLS 回歸，其R2、交叉驗證均方根誤差分別為0.767、1.96 d 和0.445、5.37 d，說明通過分段建模的方式可以分別提高9 d 前和12 d 后的模型推斷準(zhǔn)確性，提示在PMI 推斷中應(yīng)避免采用單一統(tǒng)計方法進(jìn)行PMI 推斷，運用多段建模的方式可提高其準(zhǔn)確性。

腸道菌群多樣性容易受個體飲食習(xí)慣、臨時性用藥及周圍環(huán)境變化、死后腐敗等多種因素的影響[7]，如何有效避免上述因素對與PMI 具有高度相關(guān)的菌群變化的影響是腸道菌群用于PMI 推斷的關(guān)鍵。本研究在實驗設(shè)計中并未考慮動物模型的致死方式（如麻醉藥）及周圍環(huán)境對腸道菌群多樣性的影響，因此在今后的研究工作中應(yīng)考慮尸體致死方式、周圍環(huán)境（溫度、濕度等）、土壤類型、土壤質(zhì)地、昆蟲等[9，20-21]影響因素，同時繼續(xù)延長死后時間，增加時間點，提高實驗數(shù)據(jù)準(zhǔn)確度，并盡可能減少誤差。

綜上所述，利用16S rRNA 高通量測序技術(shù)檢測的腸道微生物在死亡30 d 內(nèi)呈現(xiàn)一定時序性變化，采用分段建模的方式可以提高死亡9 d內(nèi)和死后12～30 d的準(zhǔn)確性。