多殺性巴氏桿菌脂多糖刺激羊支氣管上皮細(xì)胞后的轉(zhuǎn)錄組測序與分析

吳艷茹，陳巧玲，劉志勇，李崇瑞，黃惠嫻，翟 哲，王雪梅，滿初日嘎，杜 麗，王鳳陽

（海南大學(xué)動物科技學(xué)/海南省熱帶動物繁育與疫病研究重點實驗室/海口市動物基因工程重點實驗室，海南 ?？?570228）

多殺性巴氏桿菌（Pasteurella multocida，Pm）是一種革蘭氏陰性、兼性厭氧、不產(chǎn)生孢子的球狀桿菌，能夠感染多種家養(yǎng)和野生哺乳動物、鳥類和爬行動物以及人類[1]。Pm 可分為5 個莢膜血清組（A、B、D、E 和F）和16 個脂多糖（LPS）血清型（血清型1～16），禽霍亂主要由A：1、A：3 和A：4 血清型引起，牛出血性敗血癥主要由B：2 和E：2 血清型引起，豬的進(jìn)行性萎縮性鼻炎主要由A 和D 血清型引起，牛呼吸道疾病主要由A：3 血清型引起，兔巴氏桿菌病主要由B：2 和E：2 血清型引起，羊巴氏桿菌病主要由A 和D 血清型引起[2-5]。在山羊和綿羊中，Pm 是類似于溶血性曼氏桿菌的條件性致病菌，主要引起機體以肺部病變?yōu)橹鞯陌褪蠗U菌病[6]。對Pm感染后有臨床癥狀羊的檢測結(jié)果顯示，病原菌主要存在于其內(nèi)臟、淋巴結(jié)、血液及部分病變組織中，而對Pm 感染后無臨床癥狀羊的檢測結(jié)果顯示，病原菌主要存在于其黏膜、上呼吸道以及扁桃體中[7]。Pm嚴(yán)重危害羊的健康，養(yǎng)羊業(yè)也因此遭受重大經(jīng)濟損失[8]。盡管如此，Pm 的致病機制仍尚未明確。

脂多糖（Lipopolysaccharide，LPS）是Pm 的主要毒力因子，也是革蘭氏陰性菌細(xì)胞外膜的主要成分[9]，因此，LPS 在細(xì)菌與環(huán)境的相互作用中起至關(guān)重要的作用，并與膜屏障功能、細(xì)菌毒力和宿主免疫密切相關(guān)[10]。在病原體中，LPS 通過直接與宿主固有免疫系統(tǒng)（如抗菌肽和補體成分）的相互作用，在疾病的發(fā)展過程中發(fā)揮重要作用，如LPS 可以通過與Toll 樣受體4（TLR4，在某些細(xì)菌中為TLR2）的相互作用引起宿主強烈的固有免疫反應(yīng)，激活一系列宿主免疫細(xì)胞，從而引起細(xì)胞因子風(fēng)暴和宿主死亡[11]。疫苗是防控疾病最經(jīng)濟、有效的手段，LPS具有良好的免疫原性，其特異性抗體可介導(dǎo)針對同源菌株和LPS 結(jié)構(gòu)相似菌株的免疫保護(hù)，因而LPS成為Pm 疫苗的重要研究對象[12]。羊的巴氏桿菌病主要由莢膜血清D 型和A 型引起，目前對羊源Pm LPS研究較少，對其在Pm 侵入羊機體過程中發(fā)揮的具體作用以及引起的宿主免疫學(xué)變化還有待研究。

本實驗利用莢膜血清D 型和A 型Pm LPS（D-LPS和A-LPS）刺激羊支氣管上皮細(xì)胞6 h 后，通過轉(zhuǎn)錄組學(xué)測序與分析獲得的轉(zhuǎn)錄差異基因，并對轉(zhuǎn)錄差異基因進(jìn)行GO功能注釋和KEGG富集分析，為Pm LPS與宿主相互作用及其致病機制和疫苗的研發(fā)提供參考依據(jù)。

1 材料與方法

1.1 主要實驗材料D 型Pm HN01 株（cp037861）和A 型Pm HN02 株（cp037865.1）由本實驗室分離；LPS Extraction Kit 購自iNtRON 公司；羊支氣管上皮細(xì)胞SHP-iCell-a008、雙抗（青霉素/鏈霉素）和胰酶均購自賽百慷（上海）生物技術(shù)股份有限公司；RPMI 1640 液體培養(yǎng)基（含雙抗）購自博士德生物工程有限公司；胎牛血清（用于配制含10% FBS 的培養(yǎng)基）購自美國Gibco 公司；O111：B4 血清型大腸桿菌LPS凍干粉購自Sigma 公司；cDNA 文庫建庫試劑盒購自諾唯贊公司；熒光定量PCR 引物由武漢天一華煜基因科技有限公司合成。D-LPS 和A-LPS 由實驗室制備，并由銀染法鑒定[13]。

1.2 LPS 刺激羊支氣管上皮細(xì)胞將羊支氣管上皮細(xì)胞以5×105個細(xì)胞/孔鋪于六孔細(xì)胞培養(yǎng)板中。分為實驗組和對照組，每組各作3 個重復(fù)（A-LPS-1～A-LPS-3、D-LPS-1～D-LPS-3、NC-1～NC-3）。實驗組加入2 mL 含LPS 的DMEM 完全培養(yǎng)基使其終濃度為1 μg/mL，對照組加入等體積的DMEM 完全培養(yǎng)基，刺激細(xì)胞6 h 后，每孔加1 mL 裂解液裂解細(xì)胞后移至凍存管中，迅速放入液氮冷凍，郵寄至上海派森諾生物科技有限公司。并由其構(gòu)建各組細(xì)胞的cDNA 文庫及基于Illumina HiSeqTM6000 測序平臺，對構(gòu)建的文庫經(jīng)雙末端（Paired-end，PE）測序。

1.3 數(shù)據(jù)質(zhì)控為保證數(shù)據(jù)的可靠性對測序結(jié)果采用下述方式過濾，（1）采用Cutadapt 去除3'端帶接頭的序列；（2）去除低質(zhì)量reads[Q20（%）或Q30（%）低于90%]；（3）去除含N（即無法識別的堿基）比例大于5%的reads。同時統(tǒng)計有效測序量（Clean data）、高質(zhì)量序列reads（Clean reads No.）、Q30（堿基識別準(zhǔn)確率在99.9%以上的堿基總數(shù)）、N（%）（模糊堿基所占百分比）、Q20（%）（堿基識別準(zhǔn)確率在99%以上的堿基所占百分比）和Q30（%）（堿基識別準(zhǔn)確率在99.9%以上的堿基所占百分比）。

1.4 轉(zhuǎn)錄差異基因的篩選及聚類分析采用DESeq2（R 語言軟件包，用于測序數(shù)據(jù)的標(biāo)準(zhǔn)化、可視化和差異分析）篩選各組之間轉(zhuǎn)錄差異的基因（篩選標(biāo)準(zhǔn)log2（Fold Change）＞1，P＜0.05），利用ggplots2繪制轉(zhuǎn)錄差異基因的火山圖。利用Pheatmap 分析各組的組內(nèi)重復(fù)性，并繪制熱圖對各組轉(zhuǎn)錄差異基因并集和對樣品進(jìn)行雙向聚類分析。

1.5 轉(zhuǎn)錄差異基因的GO 功能注釋和KEGG 富集分析利用top GO 對轉(zhuǎn)錄差異基因進(jìn)行GO 聚類分析，分析時利用GO term 注釋的差異基因計算每個term 的基因列表和基因數(shù)目，篩選出與整個山羊參考基因組背景（Capra hircus.ARS1.dna.toplevel.fa）相比，差異基因顯著富集的GO term（顯著富集的標(biāo)準(zhǔn)為P＜0.05），從而確定差異基因行使的主要生物學(xué)功能。按照分子功能（Molecular function，MF）、生物過程（Biological process，BP）和細(xì)胞組分（Cellular components，CC）對轉(zhuǎn)錄差異基因的GO 功能注釋分析結(jié)果進(jìn)行GO 分類，將P值由小到大選出每組前10 的GO term。采用KEGG 富集分析轉(zhuǎn)錄差異基因（P＜0.05），并根據(jù)該結(jié)果，挑選P值最小即富集最顯著的前20 個信號通路。

1.6 轉(zhuǎn)錄差異基因的RT-qPCR 驗證從D-LPS 和A-LPS 兩組測序數(shù)據(jù)中各隨機選取共13 個轉(zhuǎn)錄差異基因進(jìn)行下述RT-qPCR 驗證。LPS 刺激細(xì)胞后6 h，收集各組細(xì)胞，采用TRIzol 法提取細(xì)胞樣品的RNA，反轉(zhuǎn)錄為cDNA 后作為模板，以GAPDH 作為內(nèi)參基因（引物序列見表1），經(jīng)RT-qPCR 驗證轉(zhuǎn)錄組測序結(jié)果的準(zhǔn)確性。根據(jù)2-ΔΔCt計算各基因的相對轉(zhuǎn)錄水平，采用GraphPad Prism8統(tǒng)計軟件作圖。

表1 RT-qPCR引物Table 1 RT-qPCR primers used in this study

2 結(jié) 果

2.1 各組細(xì)胞cDNA 文庫測序數(shù)據(jù)的整理與質(zhì)控結(jié)果對構(gòu)建的各組細(xì)胞的cDNA文庫（濃度是10 ng/μL～20 ng/μL），采用基于Illumina HiSeqTM6000 測序平臺的NGS 經(jīng)PE 測序。并對原始測序數(shù)據(jù)經(jīng)過濾后獲得有效測序量和高質(zhì)量序列堿基數(shù)。結(jié)果顯示平均有效測序量達(dá)182 億個，平均高質(zhì)量序列reads 數(shù)達(dá)1.22 億個且Q20 和Q30 值均大于90%，N 值均小于0.002%（表2），表明質(zhì)控后得到了高質(zhì)量的測序數(shù)據(jù)，證實了轉(zhuǎn)錄組測序數(shù)據(jù)的可靠性。

表2 測序數(shù)據(jù)統(tǒng)計分析Table 2 Statistical analysis of sequencing data

2.2 各組之間轉(zhuǎn)錄差異基因的篩選與聚類分析利用DESeq2 對各組轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行差異分析，設(shè)置篩選轉(zhuǎn)錄差異基因的標(biāo)準(zhǔn)為log2（Fold Change）≥1 和P＜0.05。結(jié)果顯示，D-LPS 中共有顯著轉(zhuǎn)錄差異基因86 個，其中轉(zhuǎn)錄顯著上調(diào)的基因62 個，轉(zhuǎn)錄顯著下調(diào)的基因24 個。A-LPS 中共有顯著轉(zhuǎn)錄差異基因56 個，其中轉(zhuǎn)錄顯著上調(diào)的基因27 個，轉(zhuǎn)錄顯著下調(diào)的基因29 個（圖1），兩組共同轉(zhuǎn)錄差異基因6 個，分別是LOC102172606、EFNA2、NCMAP、ADAMTS5、LRRIQ3、 ENSCHIG00000016421。 采 用Pheatmap 分析各組的組內(nèi)重復(fù)性，并對轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行轉(zhuǎn)錄差異基因的聚類分析，結(jié)果顯示，實驗組和對照組的組內(nèi)重復(fù)性均較好，且羊支氣管上皮細(xì)胞經(jīng)LPS 刺激前后基因的轉(zhuǎn)錄水平存在明顯的差異（圖2）。綜上表明，經(jīng)Pm LPS 刺激的羊支氣管上皮細(xì)胞中存在轉(zhuǎn)錄顯著差異的基因且各組組內(nèi)重復(fù)性良好，可以進(jìn)行后續(xù)的GO 功能注釋和KEGG 富集分析。

圖1 各組轉(zhuǎn)錄差異基因的火山圖分析Fig.1 Volcano map analysis of transcriptional differential genes between two groups

圖2 各組轉(zhuǎn)錄差異基因的熱圖分析Fig.2 Heat map analysis of transcriptional differential genes between two groups

2.3 轉(zhuǎn)錄差異基因的GO 功能注釋分析對D-LPS和A-LPS 兩組中轉(zhuǎn)錄差異基因進(jìn)行GO 功能注釋分析。結(jié)果顯示， D-LPS 組中共顯著富集到243 個GO term，包括146 個BP、44 個CC、53 個MF；A-LPS 組中共顯著富集到167 個GO term，包括109 個BP、19個CC、39 個MF。按照P值由小到大排列，取每組生物學(xué)功能前10的GO term進(jìn)行分析，結(jié)果顯示，在DLPS 中顯著富集的term 包括RAGE 受體結(jié)合（RAGE receptor binding，GO：0050786）、ATP 生物合成過程（ATP biosynthetic process，GO：0006754）、NADH脫氫酶活性（NADH dehydrogenase activity，GO：0003954）、嘌呤核糖核苷三磷酸生物合成過程（Purine ribonucleoside triphosphate biosynthetic process，GO：0009206）、氧化還原酶活性，作用于NAD（P）H、醌或類似化合物作為受體（oxidoreductase activity，acting on NAD（P）H，quinone or similar compound as acceptor GO：0016655）等（圖3A），主要與免疫和能量生成代謝相關(guān)，表明D-LPS 刺激細(xì)胞后，細(xì)胞的免疫應(yīng)答、細(xì)胞增殖分化等活動比較活躍；而在A-LPS 中富集的term 包括IL-11 受體結(jié)合（Interleukin-11 receptor binding，GO：0005142）、磷脂酶A2 活性（phospholipase A2 activity，GO：0004623）、GTP 酶活性（GTPase activity，GO：0003924）和尿嘧啶分解代謝過程（Uracil catabolic process，GO：0006212）等（圖3B），表明A-LPS 刺激羊支氣管上皮細(xì)胞后可能影響宿主細(xì)胞的代謝、免疫應(yīng)答等通路。

2.4 轉(zhuǎn)錄差異基因的KEGG 富集分析對轉(zhuǎn)錄差異基因進(jìn)行KEGG 富集分析，根據(jù)P值由小到大排列，每組共展示了前20 條的富集通路。結(jié)果顯示，D-LPS 組顯著富集到細(xì)胞氧化磷酸化、產(chǎn)熱、疾病以及ECM-受體相互作用途徑等通路（圖4A）；ALPS 組顯著富集到Ras 信號通路以及泛酸和輔酶A 生物合成等信號通路（圖4B）。兩組數(shù)據(jù)KEGG 分析均富集到了能量生成與細(xì)胞增殖分化相關(guān)通路，表明D-LPS 和A-LPS 均會影響宿主細(xì)胞的能量生成與細(xì)胞增殖。

圖4 各組轉(zhuǎn)錄差異基因的KEGG富集分析Fig.4 KEGG enrichment analysis of transcriptional differential genes between two groups

2.5 轉(zhuǎn)錄組測序結(jié)果的RT-qPCR 驗證為了驗證轉(zhuǎn)錄組測序數(shù)據(jù)的準(zhǔn)確性，本研究從D-LPS、ALPS 兩組中各隨機篩選出共13 個轉(zhuǎn)錄差異基因經(jīng)RT-qPCR 驗證并分別與轉(zhuǎn)錄組數(shù)據(jù)比較。結(jié)果顯示，11 個基因轉(zhuǎn)錄水平的上下調(diào)結(jié)果與轉(zhuǎn)錄組測序結(jié)果完全一致（圖5），表明利用轉(zhuǎn)錄組測序獲得的結(jié)果雖存在一定偏差但可靠性仍較高，可以較真實反映經(jīng)不同血清型Pm LPS 刺激后的羊支氣管上皮細(xì)胞轉(zhuǎn)錄組學(xué)的差異水平。

圖5 轉(zhuǎn)錄差異基因的RT-qPCR驗證Fig.5 RT-qPCR verification of transcriptional differential gene

3 討 論

Pm 是一種常見的動物性致病菌，根據(jù)其莢膜抗原分為5 種血清型（A、B、D、E 和F）[14]。在我國主要流行的是A、B、D 這3 種血清型。本實驗室前期從病羊中分離到的莢膜血清D 型Pm HN01 株可引起小鼠的肺泡損傷[15]，莢膜血清A 型Pm HN02 株經(jīng)PCR 鑒定發(fā)現(xiàn)攜帶fimA、toxA、exbB、exbD等12 個毒力基因[16]。后續(xù)本實驗室分別利用D 型和A 型Pm感染小鼠，結(jié)果顯示D 型菌株毒力強于A 型菌株（文章未發(fā)表）。LPS 是Pm 的主要毒力因子，Horadagoda 等將B 型Pm LPS 注射到水牛體內(nèi)可以引發(fā)牛的出血性敗血癥[17]。Periasamy 利用B 型Pm LPS 刺激能夠引起牛白細(xì)胞線粒體功能障礙并導(dǎo)致細(xì)胞死亡[18]，表明Pm LPS 具有內(nèi)毒素特性。由于目前關(guān)于Pm LPS 的相關(guān)研究較少，為確定Pm LPS 是否是造成D 型和A 型Pm 毒力不同的關(guān)鍵因子以及探究LPS對宿主的致病機制，本實驗分別利用D 型和A 型Pm LPS 刺激羊支氣管上皮細(xì)胞，6 h 后收獲各組細(xì)胞，裂解后進(jìn)行轉(zhuǎn)錄組測序，測序數(shù)據(jù)經(jīng)過濾后，最終從D-LPS 組篩選出上調(diào)的轉(zhuǎn)錄差異基因62 個，下調(diào)的轉(zhuǎn)錄差異基因24 個，A-LPS 組篩選出上調(diào)的轉(zhuǎn)錄差異基因27 個，下調(diào)的轉(zhuǎn)錄差異基因29 個，DLPS 組轉(zhuǎn)錄差異基因的數(shù)量多于A-LPS 組，推測可能與D 型菌株毒力強于A 型菌株有關(guān)。D-LPS 組轉(zhuǎn)錄差異基因的GO 功能注釋結(jié)果顯示，顯著富集且與免疫相關(guān)的生物學(xué)過程為RAGE 受體結(jié)合，該過程包含基因S100A8 和S100A12。S100A8 和S100A12 屬于S100家族的鈣結(jié)合蛋白，S100 家族與多種宿主免疫信號和細(xì)胞因子的產(chǎn)生和分泌有關(guān)，能引發(fā)全身炎癥反應(yīng)和繼發(fā)性器官衰竭[19]。革蘭氏陰性菌可引起患者刺激性休克，其體內(nèi)高水平的S100A8 可以激活RAGE 信號并導(dǎo)致炎癥損傷，這與本實驗結(jié)果一致，并且在革蘭氏陰性菌刺激過程中，S100A8 作為TLR4的配體，被強烈誘導(dǎo)表達(dá)[20]。不僅如此，S100A8通常與S100A9以異二聚體的形式存在，在體外，S100A8/A9 可以抑制中性粒細(xì)胞的氧化代謝，清除釋放的活性氧（ROS），減輕LPS 刺激小鼠肺和肝臟的氧化損傷[21]。并且S100A8 在氣道上皮細(xì)胞中可降低促炎細(xì)胞因子的分泌水平[22]。這或許是D-LPS 組轉(zhuǎn)錄差異基因中無氧化應(yīng)激相關(guān)基因和促炎細(xì)胞因子基因的原因。對A-LPS 組轉(zhuǎn)錄差異基因進(jìn)行GO 功能注釋分析顯示，顯著富集且與免疫相關(guān)的生物學(xué)過程為IL-11 與受體的結(jié)合，該過程涉及到細(xì)胞因子IL-11。IL-11是糖蛋白（GP）130細(xì)胞因子家族的成員，與其他促炎細(xì)胞因子不同，IL-11是一種抗炎細(xì)胞因子[23]，其抗炎作用的機制是抑制LPS 介導(dǎo)的NO、腫瘤壞死因子-α（TNF-α）、IL-1β、IL-12、IL-10、轉(zhuǎn)化生長因子-β（TGF-β）和IL-6的釋放，以及刺激可溶性的腫瘤壞死因子-αR1（TNF-αR1）直接拮抗巨噬細(xì)胞/單核細(xì)胞誘導(dǎo)的TNF-α 信號轉(zhuǎn)導(dǎo)[24]。在體內(nèi)和體外模型中，IL-11 也被證明通過激活內(nèi)皮細(xì)胞來減少炎癥對其的相關(guān)損傷，從而在免疫介導(dǎo)的損傷中起到保護(hù)宿主細(xì)胞的作用[25]。本實驗中是否因為IL-11 的顯著上調(diào)導(dǎo)致炎癥相關(guān)基因與通路均無顯著變化，后續(xù)還值得探究。Mao 的研究表明人類磷脂磷酸酶3（PLPP3）啟動子包含3 個功能性NF-κB 響應(yīng)元件，因此人類PLPP3 啟動子對NF-κB 通路的激活高度敏感[26]。多種細(xì)胞被刺激后的PLPP3 表達(dá)水平顯著上調(diào)[27-28]，這與本實驗中D-LPS 組的結(jié)果一致。

KEGG 富集分析結(jié)果顯示，在D-LPS 組中顯著富集到的通路有氧化磷酸化、產(chǎn)熱、ECM-受體相互作用、視黃醇代謝、鉑類耐藥等，涉及到能量生成與細(xì)胞增殖分化；在A-LPS 組中顯著富集到的通路有泛酸和輔酶A、Ras 信號通路等，同樣也涉及到能量生成與細(xì)胞增殖。本實驗中，D-LPS 組中與細(xì)胞凋亡相關(guān)的基因MAP2K6 顯著上調(diào)，與抑制細(xì)胞凋亡相關(guān)的基因BIRC5 顯著下調(diào)； A-LPS 組中與細(xì)胞凋亡相關(guān)的基因MXRA8 和BEX3 均顯著下調(diào)。Periasamy 的研究結(jié)果顯示，A、B 型Pm LPS 均能夠誘導(dǎo)牛白細(xì)胞（包括淋巴細(xì)胞、中性粒細(xì)胞、單核細(xì)胞和單核巨噬細(xì)胞）的凋亡，導(dǎo)致組織損傷[18]，與本實驗D-LPS 組的結(jié)果一致。據(jù)此推測D-LPS 會引起細(xì)胞凋亡，從而導(dǎo)致宿主損傷；但A-LPS 組中與細(xì)胞凋亡相關(guān)基因顯著下調(diào)是否和A 型菌株毒力較弱有關(guān)，還有待深入研究。GO 和KEGG 分析結(jié)果表明，兩組轉(zhuǎn)錄差異基因均與細(xì)胞免疫、細(xì)胞增殖分化、能量產(chǎn)生代謝等信號通路和功能有關(guān)。表明DLPS 和A-LPS 均能影響宿主細(xì)胞的免疫、增殖分化和能量的生成與代謝。

本研究選取了D-LPS、A-LPS 兩組中13 個轉(zhuǎn)錄差異基因進(jìn)行RT-qPCR 驗證，其中11 個基因的檢測結(jié)果與轉(zhuǎn)錄組測序結(jié)果一致，表明轉(zhuǎn)錄組測序結(jié)果基本可靠。并且RT-qPCR 結(jié)果驗證了與免疫相關(guān)的基因IL-11、S100A8 轉(zhuǎn)錄水平顯著上調(diào)，與氧化還原酶活性相關(guān)的基因LOXL4 轉(zhuǎn)錄水平顯著下調(diào)，與能量代謝相關(guān)的基因IGFBP1 轉(zhuǎn)錄水平顯著下調(diào)，進(jìn)一步表明Pm LPS 刺激羊支氣管上皮細(xì)胞會影響其免疫功能和能量的生成與代謝。

綜上所述，本實驗采用轉(zhuǎn)錄組測序技術(shù)篩選并分析了D-LPS 和A-LPS 刺激羊支氣管上皮細(xì)胞后的免疫、能量生成與代謝以及細(xì)胞凋亡的相關(guān)基因變化及其所富集的通路。為進(jìn)一步探究Pm LPS 的致病機理及Pm 疫苗的研發(fā)奠定了基礎(chǔ)。