基于全基因組重測序及等位基因頻率解析柔嫩艾美耳球蟲對地克珠利耐藥分子機制的初步研究

畢菲菲 ， 郝振凱 ， 孫 霈 ， 于詠蘭 ， 劉賢勇 ， 索 勛

(中國農業大學動物醫學院 ， 北京 海淀100193)

艾美耳球蟲是一類寄生于禽類和哺乳細胞內的頂復門原蟲，是導致畜禽球蟲病的病原[1]。球蟲病對養殖業危害巨大，每年對養殖業造成巨額的經濟損失[2]。目前，防控該病主要依賴抗球蟲藥物和疫苗免疫[3]。我國每年用于防治雞球蟲病的藥物費用高達20億～30億人民幣[4]。目前由于抗球蟲藥的長期使用及不規范使用，使球蟲產生了耐藥性，但是耐藥性產生的具體機制尚不清楚。地克珠利是一種廣譜抗球蟲藥，屬于三嗪類，對雞、火雞、兔、羊和牛的不同艾美耳球蟲蟲種有很好的療效[5]。研究表明，每1 kg飼料中添加1 mg的地克珠利對感染雞的大部分艾美耳球蟲都有效[6]。對于地克珠利耐藥機制的研究，有學者比較了耐藥株與敏感株在同工酶上的差異，發現耐藥株有特異性酶條帶[7]； 而利用隨機引物擴增(Random amplified polymorphic DNA,RAPD)顯示了不同耐藥分離株特定基因的多樣性[8]。mRNA 差異顯示技術、SDS-PAGE等也被用于球蟲耐藥株與敏感株在基因轉錄和蛋白表達水平的差異研究[9-12]。Shirley 等于2000年采用遺傳學方法繪制的柔嫩艾美耳球蟲的遺傳連鎖圖，推測氟腺呤(Arprinocid)耐藥性相關的連鎖群可能位于第 1 號染色體上[13]。然而，至今還未找到艾美耳球蟲地克珠利耐藥性的靶基因，亦無法解析艾美耳球蟲耐藥性產生的分子機制。

本研究將結合艾美耳球蟲耐藥性表型特征、高通量測序和生物信息學分析，試圖鑒定球蟲針對地克珠利的耐藥性基因并解析產生耐藥性的分子機制，為建立球蟲耐藥性田間快速分子診斷方法提供技術支撐。

1 材料與方法

1.1 試驗蟲株 本研究從山東省、云南省、上海市和河南省的養雞場收集新鮮糞便樣品，在實驗室分離培養出柔嫩艾美耳球蟲后保存備用。

1.2 實驗動物 本實驗室用于球蟲傳代的AA肉雞，購自北京愛拔益加家禽育種有限公司。雞只飼養于無球蟲環境中，保障充足的水和飼料供給。提供的飼料和水均經過嚴格高溫處理(80 ℃處理2 h以上)，防止球蟲污染。本研究涉及的所有動物的飼養和試驗均嚴格按照中國農業大學實驗動物福利與動物實驗倫理委員會的相關規章制度進行。用于傳代或收集裂殖子的雞只均為AA肉雞。

1.3 主要試劑和器材 飽和食鹽水:將蒸餾水煮沸后加入食鹽不斷攪拌直至飽和，冷卻后倒入封閉容器中保存備用；2.5%重鉻酸鉀溶液:稱取2.5 g重鉻酸鉀粉末，將其加入到1 L蒸餾水中，完全溶解后倒入封閉容器中保存備用；PBS溶液:將商品化的PBS粉末加到1 L蒸餾水中，攪拌直至完全溶解，倒入封閉容器中保存備用；次氯酸鈉溶液:有效氯濃度>10%；過濾篩:80目、200目；錐形瓶:250 mL、500 mL；載玻片；以上試劑和器材均購自北京吉普騰生物技術有限公司。地克珠利：濃度為0.5%，購自上海獸醫研究所。巴斯德吸管，購自上海晶安生物科技有限公司。

1.4 主要儀器設備 大容量低速離心機(Anke LXJ-IIB型)和小容量離心機(Anke LXJ-40B型)，均購自上海安亭科學儀器廠；麥克馬斯特計數板，購自上海獸醫研究所；不銹鋼禽用隔離器(IPQ-3型)，購自蘇州市蘇杭實驗動物設備廠；電熱恒溫鼓風干燥箱(DGG-9070A型)，購自上海森信實驗儀器有限公司。

1.5 試驗方法

1.5.1 田間分離蟲株的采集及樣品中球蟲卵囊的觀察 選擇育雛雞舍采集糞樣，按照隨機取樣原則進行，若雞場規模較小，則在每個雞舍隨機取樣；若雞場規模較大，則隨機選擇雞舍進行采樣。采集的糞便樣品注意分類保存到自封袋或封閉容器中，加入2.5%的重鉻酸鉀溶液后運回實驗室進行后續處理。

將糞便樣品攪拌均勻后，用天平稱取2 g至100 mL燒杯中。用量筒量取30 mL飽和食鹽水，倒入盛有糞便樣品的燒杯中并攪拌均勻。靜置5～10 min后，用巴斯德吸管吸取上層液體并注入麥克馬斯特計數板中，然后放到顯微鏡下靜置3～5 min。在10倍物鏡下觀察和計數球蟲卵囊。

1.5.2 田間分離蟲株的分離和收集 對鏡檢下有球蟲的樣品，收集其中的卵囊。以混合卵囊接種2周齡左右AA肉雞。感染后7 d收集盲腸，將腸道剪碎后，以2.5%胰酶42 ℃消化30 min，紗布過濾腸道組織后，3 600 r/min離心6 min。將沉淀中的卵囊用飽和食鹽水漂浮法進行收集和孢子化。以多卵囊的方式從中分離出柔嫩艾美耳球蟲，即將卵囊稀釋成約1 μL的濃度，加1 μL稀釋好的卵囊液于載玻片上，在顯微鏡下觀察確認每滴液體中的卵囊是柔嫩艾美耳球蟲后，小心收集10個卵囊，并接種于2只1周齡左右的AA肉雞，收集后代卵囊，并對后代卵囊進行傳代擴增。

1.5.3 田間分離株的耐藥性測試

1.5.3.1 試驗分組及感染 試驗時所用雞只為13日齡AA肉雞，蟲株為分離獲得的田間蟲株，分組時每個蟲株均設立感染給藥組和感染不給藥組，每組5只雞，同時設立不感染不給藥組(10只雞)。取保存在重鉻酸鉀溶液中的孢子化卵囊，用PBS清洗除去重鉻酸鉀，各感染給藥組和感染不給藥組同一天接種不同分離蟲株的孢子化卵囊5×104個/只。試驗劑量采用地克珠利臨床推薦的劑量 1 mg/L，飲水服用，每1 L水中加入 200 μL 0.5%地克珠利溶液。感染給藥組從感染前1天飼喂含地克珠利藥物的飲用水直到試驗結束。分別稱量感染前、感染后第7天的體重。試驗過程中觀察雞只的狀態，如有無便血等，如雞只死亡，對其進行剖檢，觀察腸道病變情況，判斷雞只是否死于球蟲感染。感染后第7天將所有雞只剖殺，對腸道進行病變計分，然后將每只雞腸道分別收集到 50 mL 離心管中，加入重鉻酸鉀溶液，進行卵囊勻漿，計算每只雞腸道中球蟲卵囊的數量。

1.5.3.2 耐藥性評價指標 本研究選取的耐藥性評價指標為相對卵囊產量(Relative oocyst production,ROP)、最適抗球蟲活性百分率(Percent of optimum anticoccidial activity,POAA)、病變記分減少率(Reduction of lesion scores,RLS)和抗球蟲指數(Anticoccidia indexes,ACI)[14]。

相對卵囊產量(ROP)/%=(感染給藥組的平均卵囊產量/感染不給藥組的平均卵囊產量)×100%。ROP<15%，耐藥性為陰性；ROP≥15%，耐藥性為陽性。

最適抗球蟲活性百分率(POAA)/%=(感染給藥組GSR-感染不給藥組GSR)/(不感染不給藥組GSR-感染不給藥組GSR)×100%，其中增重和存活率(Growth and survival ratio,GSR)=籠末重/籠初重。POAA>50%，耐藥性為陰性；POAA≤50%，耐藥性為陽性。

病變記分減少率(RLS)/%=(感染不給藥組的平均病變記分-感染給藥組的平均病變記分)/感染不給藥組的平均病變記分×100%。RLS≥50%，耐藥性為陰性；RLS<50%，耐藥性為陽性。

抗球蟲指數(ACI)=(存活率+相對增重率)-(病變值+卵囊值)，其中存活率/%=(存活雞只數/開始雞只數)×100%；相對增重率/%=(感染不給藥組雞只的平均增重/不感染不給藥組雞只的平均增重)×100%；病變值=每組的平均盲腸病變記分×10。ROP為0～1%，卵囊值為0；ROP為26%～50%，卵囊值為10；ROP為51%～75%，卵囊值為20；ROP為76%～100%，卵囊值為40。ACI>160，耐藥性為陰性;ACI≤160，耐藥性為陽性。

對雞球蟲耐藥性進行評價時，結合上面提到的4項指標進行綜合判定：4項指標均為耐藥性陰性時，判定為對藥物敏感；1項指標為陽性時，判定為輕度耐藥；2項指標為陽性時，判定為中度耐藥；3項或4項指標為陽性時，判定為完全耐藥。

1.5.4 田間分離株第2代裂殖子基因組DNA的提取 取新鮮柔嫩艾美耳球蟲分離株分別接種1月齡左右的AA肉雞，每只肉雞接種2×105個卵囊，每株接種3只雞。接種后第120 小時，剖殺雞只，取盲腸，輕輕移除盲腸內容物并清洗，將盲腸剪碎后，向其中加入42 ℃溫熱的消化液，置于恒溫攪拌器上，42 ℃消化30 min。取消化后懸液，用4層紗布過濾。濾液轉入50 mL離心管中，3 600 r/min 離心6 min。離心后，小心棄除一半的上層消化液，用1 mL加樣槍，小心吸除漂浮于白色沉淀上的絮狀組織和消化液，將下層的白色沉淀轉移至1.5 mL EP管中，加入1 mL PBS，混勻，3 000 r/min 離心6 min。此時沉淀分為3層，最下層為白色，即柔嫩艾美耳球蟲第2代裂殖子，用PBS重懸，將其轉入新的1.5 mL EP管中，并用血細胞計數板計數。每個樣品取約5×108個第2代裂殖子，加入1 mL PBS，混勻，3 600 r/min離心5 min，棄上清。基因組DNA提取使用Blood & Cell Culture DNA Midi Kit，嚴格按照試劑盒說明書進行操作。DENOVIX 超微量分光光度計測量每個樣品的基因組DNA的濃度和純度。

1.5.5 高通量全基因組重測序及生物信息學分析 田間分離柔嫩艾美耳球蟲蟲株基因組DNA經分光光度計測量合格后，送至安諾優達基因科技(北京)有限公司進行 Illumina小片段文庫構建及全基因組重測序。單個樣品測序數據量大于5 Gb(大于100×平均基因組覆蓋度)。

1.5.6 柔嫩艾美耳球蟲耐藥性相關基因的分析 將田間分離株全基因組重測序得到的數據進行過濾，再以Fastqc軟件進行質檢，符合Q30 > 80%后即認為Clean data，可用于后續分析。首先構建柔嫩艾美耳球蟲基因組的Index，再用BWA軟件分別將Clean reads比對到新版參考基因組(豪頓株，ETH)上，生成SAM文件，然后用SAMtools將SAM文件進行排序并轉換成BAM文件，隨后借助GATK v4.011軟件中的MarkDuplicates工具將由建庫過程中聚合酶鏈式反應(Polymerase chain reaction,PCR)產生的復制情況進行標記。因本次分析需要將耐藥株和敏感株視為2個由不同數目個體組成的微生物群體，而分析2個群體之間的差異。在最后Call SNP時需用不同的參數進行。Freebayes v1.0[15]具體參數：-F 0.01 -C 3 --pooled-continuous --min-mapping-quality 30 --min-base-quality 20 --min-coverage 10。將生成的VCF文件以VCFtools[16]提取SNP，并用Python腳本進行質量過濾，過濾標準：QUAL<60.0,DP<10。

ΔSNP-index計算公式：單點SNP-index=突變基因型的覆蓋度/總覆蓋度，ΔSNP-index=SNP-index(耐藥)-SNP-index(敏感)。本試驗利用Python編寫腳本程序，計算耐藥株與敏感株基因組之間每50 kb窗口的ΔSNP-index，并以20 kb作為滑窗步長。最后以R語言qqman語言包進行可視化作圖。

2 結果

2.1 田間分離柔嫩艾美耳球蟲蟲株的分離 本試驗共獲得4株田間柔嫩艾美耳球蟲蟲株，分別來自山東、云南、上海和河南，依次命名為Et-SD、Et-YN、Et-SH和Et-HN。

2.2 田間分離的柔嫩艾美耳球蟲蟲株的耐藥性測試 本研究對4株田間分離蟲株進行了耐藥性相關試驗，結果如表1所示，通過4項指標的評價，本研究分離到的4株田間蟲株中，河南蟲株4項指標耐藥表型為陽性，表明河南地方株對地克珠利是完全耐藥的，而山東株、云南株、上海株的4項指標耐藥表型為陰性，說明這3株地方株對地克珠利具有敏感性。

表1 不同分離蟲株對地克珠利的耐藥性Table 1 Resistance of different Eimeria tenella strains when treating with diclazuril

2.3 田間分離蟲株第2代裂殖子基因組DNA的提取 本試驗獲得4株分離蟲株第2代裂殖子的高質量DNA(表2)，并送樣至公司進行全基因組重測序，另一部分存儲于-80 ℃備用。

表2 不同分離蟲株第2代裂殖子DNA的提取結果Table 2 The quality of total DNA extracted from the 2nd generation merozoites of different isolates

2.4 柔嫩艾美耳球蟲全基因組重測序及ΔSNP-index分析 本研究從田間分離到的3株敏感株和1株耐藥株的第2代裂殖子DNA送樣至公司進行全基因組重測序后，測序結果與柔嫩艾美耳球蟲參考基因組(豪頓株,ETH)進行比對，結果如表3所示，田間分離株與ETH株全基因組比對結果均大于90%，每個樣品與參考蟲株比較，提取出的SNP個數在40 000個左右。對于耐藥和敏感種群之間的高質量SNP(n=111 694)，通過計算核苷酸多樣性 Theta Pi (θπ) 沿基因組具有20 kb 窗口來估計遺傳多樣性。在第4、10、12、13、14號染色體上可以觀察到高水平的核苷酸多態性，相比之下，第1、2、3、7、8、9號染色體的多態性比較低(圖1)。將3株敏感株和1株耐藥株全基因組重測序數據進行群體遺傳學分析，以50 kb的滑窗統計其ΔSNP-index(耐藥/敏感)并進行分析(圖2、圖3)，結果顯示，經ΔSNP-index分析有1個窗口的ΔSNP-index值等于1，具體查看內部情況后，可見這些區域測序覆蓋度較低，區域內僅含單個SNP位點造成的假陽性位點。取ΔSNP-index值大于0.9的基因組區間，且核苷酸多態性降低的第1號和第6號染色體區域，在第1號染色體上的83 000～83 500 bp、第6號染色體的280 000～285 000 bp顯示出選擇信號，將這2個候選區域根據基因組注釋信息對這些區域內所包含的基因進行分析，可見這些區域中主要為串聯的表面抗原(Surface antigen,SAG)基因，耐藥株在其基因第36位點的堿基由C突變為A，氨基酸在第12位點由天冬氨酸變為谷氨酸GAC-GAA(D12E)(圖4)。

表3 敏感株與耐藥株全基因組比對的結果Table 3 Results of genome comparison between sensitive and resistant strains

圖1 不同染色體上多態性(θπ)的分布Fig.1 Distribution of the genetic diversity (θπ) on each chromosome

圖2 不同染色體上ΔSNP-index的分布Fig.2 Distribution of the ΔSNP-index on each chromosome

3 討論

據國內外文獻報道，目前雞球蟲的地克珠利耐藥性已經非常普遍[17]，其對耐藥性產生的分子機制研究是非常必要的，同時也可以為新藥的研發提供靶位點，并為田間快速診斷地克珠利耐藥性提供分子標記。李莎通過實時熒光定量PCR (Quantitative real-time PCR,qPCR)發現添加地克珠利藥物的蟲株的乳酸脫氫酶的 mRNA 水平顯著下降[18]；Wang等采用熒光定量 PCR 和蛋白質印跡等技術研究地克珠利對柔嫩艾美耳球蟲第2代裂殖子甘油醛-3-磷酸脫氫酶 mRNA 和蛋白表達量的影響，發現其 mRNA 表達水平顯著提高[19]；Shen等發現，地克珠利治療組熱休克蛋白90 mRNA的表達水平比對照組下降了29.7%[20]，推測地克珠利可能通過作用于乳酸脫氫酶、甘油醛-3-磷酸脫氫酶和熱休克蛋白90發揮抗球蟲作用；Zhou等發現用地克珠利處理后，烯醇化酶的mRNA和蛋白表達水平分別降低了36.30%和40.36%[21]。遺憾的是，這些研究限于方法的不足，未能直接對基因突變是否為造成球蟲的地克珠利耐藥性產生的原因進行驗證。

圖3 第1號和第6號染色體上ΔSNP-index的分布Fig.3 Distribution of the ΔSNP-index on 1 and 6 chromosome

圖4 敏感株和耐藥株在SAG基因上的突變位置信息Fig.4 Information on the mutation positions of sensitive and resistant strains in the SAG gene

隨著測序技術的不斷發展，全基因測序技術為尋找耐藥蟲株與敏感蟲株間差異表達的基因提供了極大的方便。通過全基因組測序可以找到大量差異的基因。本研究共對3株敏感株和1株耐藥株進行了全基因組重測序，通過比較等位基因頻率，首次發現這些區域中主要為串聯的SAG基因。Xie等也通過轉錄組比較發現，SAG基因在地克珠利耐藥蟲株中顯著高表達[22]，與本研究結果相一致。因此，我們猜測通過本研究計算找到的敏感株和耐藥株在SAG基因上存在的非同義位點突變很可能與柔嫩艾美耳球蟲的地克珠利耐藥性相關，后續還需通過基因編輯的方法對此基因進行功能驗證。