丁恒武, 蔣 瀾, 章 勤, 王青青, 王 瑩, 董錦繡, 吳 璇, 易雙龍, 陳仁瑞, 闞顯照

(安徽師范大學(xué) 生命科學(xué)學(xué)院 生物信息研究所, 蕪湖 241000)

干擾素(interferon，IFN)是一類在特定誘導(dǎo)劑(細(xì)菌、病毒等)作用下產(chǎn)生的一類具有抗病毒、抗腫瘤、免疫調(diào)節(jié)等生物學(xué)功能的細(xì)胞因子[1-4]，由Isaacs等[5]于1957年利用雞胚絨毛尿囊膜研究流感病毒的干擾現(xiàn)象時發(fā)現(xiàn)。根據(jù)氨基酸的序列、干擾素的功能及相應(yīng)受體的結(jié)合，可將干擾素分為3種類型：Ⅰ型、Ⅱ型和Ⅲ型[6]。Ⅰ型干擾素包括IFN-α、IFN-β、IFN-ω、IFN-ε、IFN-τ、IFN-δ、IFN-κ及IFN-ξ，Ⅱ型和III型干擾素只有一種類型，分別為IFN-γ和IFN-λ[7]。其中，IFN-α(干擾素-α)由基因IFNA編碼，是脊椎動物細(xì)胞受到病毒侵染時，由白細(xì)胞產(chǎn)生的一種調(diào)節(jié)蛋白，具有廣譜的抗病毒作用[8]。各種致病性細(xì)菌或病毒對鳥類的生長繁殖影響很大，尤其是禽流感疫情帶來的影響，不僅對家禽養(yǎng)殖業(yè)造成巨大的沖擊，也對人類健康帶來了巨大的隱患。已有實驗表明，經(jīng)過IFN-α預(yù)處理過的細(xì)胞對抵抗H7N9等禽流感病毒有重要作用[9]。

IFNA為一種重要的免疫基因，在禽類中關(guān)于其免疫功能研究已有相關(guān)報道。但關(guān)于IFNA在禽類中的分子進(jìn)化研究的報道較少，因此開展禽類IFNA序列結(jié)構(gòu)特征、進(jìn)化選擇模式和系統(tǒng)發(fā)育分析的研究顯得非常重要。

本研究采用PCR的方法，對雁形目8種鳥類IFNA進(jìn)行擴(kuò)增和測序[10]，并聯(lián)合GenBank已釋放的雞雁小綱數(shù)據(jù)，對IFNA核苷酸序列和對應(yīng)的氨基酸序列進(jìn)行初步的生物信息學(xué)分析。本研究主要探討：1)IFNA的序列結(jié)構(gòu)特征；2)IFN-α的理化性質(zhì)；3)IFNA編碼區(qū)的進(jìn)化速率；4)基于IFNA核苷酸序列的系統(tǒng)發(fā)育關(guān)系。

1　材料與方法

1.1　材料

本研究實驗材料選取了雁形目鴨科鳥類8種，分別為赤嘴潛鴨(Nettarufina)、斑嘴鴨(Anaspoecilorhyncha)、赤麻鴨(Tadornaferruginea)、翹鼻麻鴨(Tadornatadorna)、大天鵝(Cygnuscygnus)、鴻雁(Ansercygnoides)、豆雁(Anserfabalis)、白額雁(Anseralbifrons)，樣本由安徽師范大學(xué)生命科學(xué)學(xué)院生物資源保護(hù)與開發(fā)省級重點實驗室保存(表 1)。

本研究還選取了GenBank已釋放的9種雞雁小綱鳥類的IFNA序列作為分析材料：綠頭鴨(Anasplatyrhynchos，AY879230)、灰雁(Anseranser，EU925650)、斑頭雁(Anserindicus，KF731865)、疣鼻棲鴨(Cairinamoschata，JF894229)、原雞(Gallusgallus，AM049251)、火雞(Meleagrisgallopavo，U28140)、藍(lán)孔雀(Pavocristatus，KJ001188)、中華鷓鴣(Francolinuspintadeanus，HM196761)和日本鵪鶉(Coturnixjaponica，AB154298)。

1.2　方法

1.2.1IFNA的序列擴(kuò)增及測序

根據(jù)GenBank已釋放的雁形目IFNA序列(GenBank登錄號：X84764)，使用在線引物設(shè)計工具Primer-BLAST(https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi?LINK_LOC=BlastHome)設(shè)計引物，上游引物為DNIFNA-F：5′-CCATGACCTGAAAGCGACGA-3′；下游引物為DNIFNA-R：5′-GGGCTCCGGTCAGTTCTTG-3′。

所選實驗材料總DNA的提取采用酚/氯仿抽提法[11]。

以8種雁形目鴨科鳥類的總DNA為模板，進(jìn)行PCR擴(kuò)增。PCR擴(kuò)增采用50 μL體系: ddH2O 25.94 μL，10×PCR buffer 5 μL，dNTPs(2.5 mmol/L) 6 μL，DMSO 1 μL，BSA 1 μL，Template 2 μL，DNIFNA-F(5 μmol/L) 4.4 μL，DNIFNA-R(5 μmol/L) 4.4 μL，Taqplus(5 μ/μL) 0.26 μL。PCR反應(yīng)體系為：95℃預(yù)變性5 min；94℃變性30 s，55℃退火30 s，72℃延伸1 min，共36個循環(huán)；最后72℃再延伸8 min。

將PCR產(chǎn)物進(jìn)行瓊脂糖凝膠電泳檢測，并將陽性擴(kuò)增的產(chǎn)物送至生工生物工程(上海)股份有限公司進(jìn)行測序。

表1 研究IFNA所用到的物種來源

1.2.2IFNA的核苷酸序列分析

選取已測序列和GenBank中下載的IFNA序列，使用MEGA7.0.14[12]分析IFNA的堿基組成，并根據(jù)堿基偏選公式計算: AT Skew=(A-T)/(A+T); GC Skew=(G-C)/(G+C)。通過EMBL開發(fā)的在線工具Cpgplot分析核苷酸序列的CpG island(http://www.ebi.ac.uk/Tools/seqstats/emboss_cpgplot/)。

1.2.3IFNA的編碼蛋白分析

IFNA編碼IFN-α，蛋白質(zhì)的理化性質(zhì)主要利用ExPASy在線軟件分析(http://web.expasy.org/protparam/)；通過Signal P4.1Server軟件預(yù)測IFN-α的信號肽(http://www.cbs.dtu.dk/services/SignalP/)；應(yīng)用TargetP 1.1工具對IFN-α的亞細(xì)胞定位進(jìn)行預(yù)測(http://www.cbs.dtu.dk/services/TargetP/)；通過ExPASy服務(wù)器的ProtScale程序中的K-D法，對IFN-α的疏水性/親水性進(jìn)行分析(http://web.expasy.org/protscale)；通過GOR IV在線分析平臺預(yù)測IFN-α的二級結(jié)構(gòu)(https://npsa-prabi.ibcp.fr)。

1.2.4IFNA的進(jìn)化速率分析

本研究分別選取線粒體編碼氧化磷酸化復(fù)合物I、IV、V 相關(guān)基因nad2、atp6、cytb與免疫基因IFNA的進(jìn)化速率分別從雞雁小綱(17個物種)和雁形目(12個物種)進(jìn)行比較。

本研究聯(lián)合測序結(jié)果和GenBank中已經(jīng)釋放的數(shù)據(jù)(表2)，采用PAML軟件[13]中的位點模型對所選基因進(jìn)行選擇壓力分析。此外本研究還采用了在線分析平臺Datamonkey(http://www.datamonkey.org)中的SLAC方法對所選基因進(jìn)行選擇壓力分析。

1.2.5系統(tǒng)發(fā)育樹的構(gòu)建

使用MAFFT version 7軟件[14]進(jìn)行多序列比對(http://mafft.cbrc.jp/alignment/server/)，使用DAMBE version 5.3.8軟件[15]對IFNA進(jìn)行了飽和替代性分析，并對將多序列比對所得的矩陣用NotePad 7.3.1對序列進(jìn)行拼接。

最大簡約法(MP)分析，使用PAUP*4.0b10軟件[16]。 采用TBR分支交換算法，自展檢驗1000次以確認(rèn)各分支的支持率得到穩(wěn)定的MP樹拓?fù)浣Y(jié)構(gòu)。最大似然法(ML)分析，使用RaxML GUI v.1.3.1軟件[17]。核苷酸替代模型采用GTRGAMMA，ML + rapid bootstrap，通過1000次自展檢驗評估得到穩(wěn)定的ML樹拓?fù)浣Y(jié)構(gòu)。貝葉斯法(BI)分析，使用ModelGenerator V. 0.851軟件分段對核苷酸序列計算最適堿基替代模型，以隨機(jī)樹為起始樹，4條馬爾可夫鏈(Markov chains)運行10 000 000代(每1000代抽樣1次)，最后舍棄前25%的樹，剩下的75%的樹被用于推測支持率大于50%的一致樹以及貝葉斯后驗概率(Bayesian posterior probability)[18]。

本研究選擇4種線粒體基因(nad2、cox1、atp6和cytb)、2種常染色體基因(IFNG和IL-2)和1種性染色體連鎖基因(IFNA，位于Z染色體)來構(gòu)建進(jìn)化樹。以企鵝目中的阿德利企鵝(Pygoscelisadeliae)作為外類群，內(nèi)類群分別來自雞雁小綱的雁形目和雞形目[19](表2)。

表2 本研究選擇的物種及基因的GenBank登錄號

* 表示該數(shù)據(jù)為本實驗數(shù)據(jù)

表3 IFNA編碼序列特征分析

2　結(jié)果與分析

2.1　核苷酸序列分析

本實驗IFNA測序結(jié)果已遞交公共基因數(shù)據(jù)庫NCBI(www.ncbi.nlm.nih.gov)，GenBank登錄號見表2。

如表3所示，12種雁形目鴨科鳥類CDS序列長度除了赤嘴潛鴨和大天鵝為579 bp外，其他均為576 bp，GC含量在67.01%～68.92%之間，起始和終止密碼子分別為ATG、TAA； 5種雞形目雉科鳥類CDS長度在579～594 bp之間，GC含量在61.34%～63.23%之間，起始和終止密碼子分別為ATG、TAG。雁形目和雞形目鳥類均存在明顯的正值A(chǔ)T skew與負(fù)值GC skew(圖 1)。

圖1 雞雁小綱IFNA核苷酸序列的堿基偏選

上圖為AT skew;下圖為GC skew

所選雞雁小綱17種鳥類中，除了藍(lán)孔雀和中華鷓鴣IFNA核苷酸序列中不存在CpG island，其他15種鳥類均存在不同長度的CpG island，其中火雞有2個CpG island(圖2)。

2.2　IFN-α蛋白質(zhì)性質(zhì)分析

2.2.1IFN-α氨基酸序列比對分析

對12種雁形目鴨科鳥類和5種雞形目雉科鳥類IFN-α氨基酸序列進(jìn)行比對分析。12種雁形目鴨科鳥類IFN-α氨基酸序列長度除了赤嘴潛鴨和大天鵝為191個氨基酸外，其他均為190個氨基酸；雞形目雉科鳥類中IFN-α氨基酸序列長度火雞為192個氨基酸，日本鵪鶉為197個氨基酸，其他均為193個氨基酸。雞雁小綱內(nèi)IFN-α氨基酸序列矩陣為197個氨基酸殘基(圖3)，其中66個氨基酸高度保守。

2.2.2IFN-α理化性質(zhì)分析

如表4顯示，在12種雁形目鴨科鳥類中，IFN-α蛋白分子質(zhì)量在21.6～21.8 ku之間，等電點在7.81～9.03之間，不穩(wěn)定指數(shù)在55.31～64.69之間，均為分泌通路；在5種雞形目雉科鳥類中，IFN-α蛋白分子質(zhì)量在21.9～22.3 ku之間，等電點在8.28～9.08之間，不穩(wěn)定指數(shù)在65.21～80.44，均為分泌通路；雞雁小綱17個物種的親水性平均系數(shù)均為負(fù)值，說明IFN-α蛋白是可溶性蛋白。

圖2 IFNA CpG island預(yù)測

A～C為鴻雁；D～F為原雞的IFNACpG island。圖中Percentage表示C%+G%，Obs/Exp表示觀測到的CG與CG%的比值，Putative Islands表示綜合預(yù)測的CpG island的位置

將所選IFN-α氨基酸序列用TargetP工具分析顯示，在12種雁形目鴨科鳥類中IFN-α除大天鵝、赤嘴潛鴨編碼的IFN-α信號肽剪切位點位于氨基酸序列29位與30位之間，其他物種蛋白信號肽剪切位點位于氨基酸序列28位與29位之間；在5種雞形目雉科鳥類中，原雞、藍(lán)孔雀、中華鷓鴣信號肽剪切位點位于氨基酸序列31位與32位之間，火雞信號肽剪切位點位于氨基酸序列30位與31位之間，日本鵪鶉信號肽剪切位點位于氨基酸序列25位與26位之間。圖4顯示的是鴻雁和原雞的信號肽預(yù)測。

圖3 干擾素-α氨基酸序列比對Fig 3 Amino acid sequence alignment of IFN-α

2.2.3IFN-α的二級結(jié)構(gòu)預(yù)測

IFN-α的二級結(jié)構(gòu)包括無規(guī)則卷曲、5個α-螺旋、伸展鏈。12種雁形目鳥類IFN-α無規(guī)則卷曲占52.15%～61.96%，α-螺旋占28.22%～38.04 %，伸展鏈占9.82%～12.27%；5種雞形目鳥類IFN-α無規(guī)則卷曲占53.89%～58.64%，α-螺旋占26.54%～37.65%，伸展鏈占6.17%～14.81%(表5)。

圖4 IFNA信號肽預(yù)測

上圖為鴻雁；下圖為原雞的信號肽預(yù)測。C值表示剪切位點分值，S值表示信號肽分值，Y值表示綜合剪切點分值

2.3 IFNA進(jìn)化速率分析

現(xiàn)在常用dN/dS描述基因的進(jìn)化速率，如表6所示。在雞雁小綱內(nèi)，17種鳥類基因atp6、cytb、nad2的dN/dS值均小于IFNA的dN/dS值，說明所選線粒體基因進(jìn)化速率小于免疫基因IFNA。在雁形目內(nèi)，12種鳥類所選的基因也出現(xiàn)類似的情況。使用PAML位點特異性模型M0(單一速率)和M3(離散的速率)的比較來評判位點間選擇壓力未檢測到正選擇位點；使用Datamonkey的SLAC算法，在雞雁小綱水平，免疫基因IFNA找到1個正選擇位點：Codon 16(P=0.0115)。

2.4　系統(tǒng)發(fā)育分析

經(jīng)過多序列比對，我們得到了以下序列矩陣：1)nad2編碼區(qū)長度為1041 bp；2)cytb編碼區(qū)長度為1143 bp；3)cox1編碼區(qū)長度為1551 bp；4)atp6編碼區(qū)長度為684 bp；5)IFNA編碼區(qū)長度為558 bp；6)IFNG編碼區(qū)長度為498 bp；7)IL-2編碼區(qū)長度為429 bp。

表6 基因進(jìn)化速率分析

上述7個序列矩陣的最適堿基替代模型：1)nad2最適模型為TrN+G；2)cytb最適模型為K81uf+G；3)cox1最適模型為GTR+I；4)atp6最適模型為HKY+G；5)IFNA最適模型為TrN；6)IFNG最適模型為HKY+I；7)IL-2最適模型為HKY+G。

本研究中系統(tǒng)發(fā)育樹圖的枝長和拓?fù)浣Y(jié)構(gòu)源于貝葉斯推斷法樹圖，并分別用MP、ML法數(shù)據(jù)標(biāo)注支持率，BI法數(shù)據(jù)標(biāo)注后驗概率。本研究先后構(gòu)建基于IFNA單基因的系統(tǒng)發(fā)育樹和基于上述7種基因聯(lián)合構(gòu)建的系統(tǒng)發(fā)育樹。

飽和替代性分析結(jié)果是Iss=6608

表4 干擾素-α的氨基酸理化性質(zhì)的比較

表5 IFN-α的二級結(jié)構(gòu)預(yù)測

圖5基于IFNA序列構(gòu)建的系統(tǒng)發(fā)育樹

Fig 5 Phylogenetic trees based on sequences ofIFNA

通過nad2、cytb、cox1、atp6、IFNA、IFNG和IL-2 7種基因聯(lián)合構(gòu)建的系統(tǒng)發(fā)育樹顯示(圖6)，雞形目3種鳥類聚成一大枝，ML和MP的自展檢驗值均為100，BI后驗概率為1，雞形目3種鳥構(gòu)成很強(qiáng)的單系；雁形目6種鳥類聚成一大枝，ML和MP的自展檢驗值均為100，BI后驗概率為0.99。雁形目鴨科不同屬的物種的系統(tǒng)發(fā)育關(guān)系較為穩(wěn)定，有著很高的支持率和后驗概率：Tadorna屬與Anas屬具為一支(ML和MP的自展檢驗值均為100，BI后驗概率為0.99)與Anser屬互為姐妹群(ML和MP的自展檢驗值均為100，BI后驗概率為0.99)，多基因聯(lián)合構(gòu)建的系統(tǒng)發(fā)育樹可以較好地反映所選物種進(jìn)化的歷史。

3　討論

雞雁小綱IFNA核苷酸序列較高的負(fù)值GC skew。其原因是IFNA編碼序列中高比例的C堿基，其中雁形目鴨科鳥類C>49.0%，雞形目雉科鳥類C>41.8%。通過在線工具Cpgplot證實了我們的發(fā)現(xiàn)。CpG island說明這段序列是含有豐富的5-甲基胞嘧啶，即可能存在DNA甲基化。

圖 6 基于多基因序列構(gòu)建的系統(tǒng)發(fā)育樹

干擾素-α蛋白分子質(zhì)量在21.6～22.3 ku之間，等電點在7.81～9.08之間，不穩(wěn)定指數(shù)在55.31～80.44之間，均為親水性分泌蛋白，存在信號肽及其切割位點，根據(jù)Guruprasad方法[20]表明當(dāng)不穩(wěn)定系數(shù)分值大于40時，預(yù)測蛋白質(zhì)不穩(wěn)定，故干擾素-α蛋白為不穩(wěn)定蛋白。本研究發(fā)現(xiàn)，所選物種中雁形目鴨科鳥類干擾素-α蛋白的分子質(zhì)量和不穩(wěn)定系數(shù)普遍小于雞形目稚科干擾素-α蛋白的分子質(zhì)量和不穩(wěn)定系數(shù)。

非同義替換(dN)指基因編碼的氨基酸從一種替換成另一種的現(xiàn)象；同義替換(dS)指基因編碼的氨基酸不改變的替換。如果dN/dS=1，則基因受中性進(jìn)化，dN/dS<1時，則基因受純化選擇，dN/dS>1，則認(rèn)為基因受正選擇，本研究中所選的基因均為dN/dS<1，故它們都受純化選擇作用。通常情況下，線粒體基因進(jìn)化速率要快于核基因[21]。在本研究中分析dN/dS時，發(fā)現(xiàn)所選物種無論是在雞雁小綱內(nèi)還是在雁形目內(nèi)，它們細(xì)胞核上免疫相關(guān)基因IFNA的進(jìn)化速率都快于線粒體基因nad2、cytb及atp6。并且在雞雁小綱中IFNA發(fā)現(xiàn)了1個正選擇位點，說明在雁形目和雞形目兩大類群的分歧過程中，IFNA可能經(jīng)歷了功能適應(yīng)性進(jìn)化。

由于線粒體基因嚴(yán)格的母系遺傳、進(jìn)化速率穩(wěn)定等優(yōu)點[21-22]，使得線粒體基因成為研究進(jìn)化的好材料，將線粒體基因加入到Z染色體連鎖基因或常染色體基因構(gòu)建聯(lián)合系統(tǒng)發(fā)育樹，結(jié)果將會更加可靠[19]。通過構(gòu)建系統(tǒng)發(fā)育樹，本研究發(fā)現(xiàn)所選物種中，IFNA單基因不能構(gòu)建穩(wěn)定的系統(tǒng)發(fā)育樹，當(dāng)Z染色體連鎖基因IFNA與常染色體基因IFNG、IL-2以及線粒體基因nad2、cox1、atp6和cytb聯(lián)合構(gòu)建系統(tǒng)發(fā)育樹時可以得到穩(wěn)定的、支持率高和后驗概率高的系統(tǒng)發(fā)育樹，這說明IFNA基因攜帶系統(tǒng)發(fā)育信息。

本次研究通過PCR測序并聯(lián)合已釋放的數(shù)據(jù)，通過生物信息學(xué)分析，探討了雞雁小綱IFNA的序列特征、進(jìn)化速率、系統(tǒng)發(fā)育意義以及其編碼蛋白IFN-α的性質(zhì)。為進(jìn)一步探討鳥類IFNA的分子進(jìn)化機(jī)制和生物學(xué)功能奠定了基礎(chǔ)。

[1]LITVIN O, SCHWARTZ S, WAN Z, et al. Interferon α/β enhances the cytotoxic response of MEK inhibition in melanoma[J]. Molecular Cell, 2015, 57(5): 784-796.

[2]MORO K, KABATA H, TANABE M, et al. Interferon and IL-27 antagonize the function of group 2 innate lymphoid cells and type 2 innate immune responses[J]. Nature Immunology, 2015, 17(1): 76-86.

[3]NICE T J, BALDRIDGE M T, MCCUNE B T, et al. Interferon-λ cures persistent murine norovirus infection in the absence of adaptive immunity[J]. Science, 2015, 347(6219): 269-273.

[4]YAMADA T, HORIMOTO H, KAMEYAMA T, et al. Constitutive aryl hydrocarbon receptor signaling constrains type I interferon-mediated antiviral innate defense[J]. Nature Immunology, 2016, 17(6): 687-694.

[5]ISAACS A, LINDENMANN J. Virus interference. I. The interferon[J]. Proceedings of the Royal Society of London, 1957, 147(927): 258-267.

[6]SHEPPARD P, KINDSVOGEL W, XU W, et al. IL-28, IL-29 and their class II cytokine receptor IL-28R[J]. Nature Immunology, 2003, 4(1): 63-68.

[7]XU L, YANG L, LIU W. Distinct evolution process among type I interferon in mammals[J]. Protein Cell, 2013, 4(5): 383-392.

[8]HOLKO M, WILLIAMS B R. Functional annotation of IFN-alpha-stimulated gene expression profiles from sensitive and resistant renal cell carcinoma cell lines[J]. Journal of Interferon & Cytokine Research, 2006, 26(8): 534-547.

[9]ARILAHTI V, MKELS M, TYNELL J, et al. Novel avian influenza A(H7N9) virus induces impaired interferon responses in human dendritic cells[J]. PLos One, 2014, 9(5): e96350.

[10]LOWENTHAL J W, STAEHELI P, SCHULTZ U, et al. Nomenclature of avian interferon proteins[J]. Journal of Interferon & Cytokine Research the Official Journal of the International Society for Interferon & Cytokine Research, 2001, 21(7): 547-549.

[11]SAMBROOK J, RUSSELL D W. Molecular cloning: a laboratory manual. 3rd edn[M]. New York:Cold Spring Harbor Laboratory Press, 2001.

[12]KUMAR S, STECHER G, TAMURA K. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets[J]. Mol Biol Evol, 2016, 33(7): 1870-1874.

[13]YANG Z. PAML: a program package for phylogenetic analysis by maximum likelihood[J]. Bioinformatics, 1997, 13(5): 555-556.

[14]KATOH K, STANDLEY D M. MAFFT multiple sequence alignment software version 7: improvements in performance and usability[J]. Molecular Biology & Evolution, 2013, 30(4): 772-780.

[15]XIA X. DAMBE5: a comprehensive software package for data analysis in molecular biology and evolution[J]. Molecular Biology & Evolution, 2013, 30(7): 1720-1728.

[16]SWOFFORD D L. Phylogenetic analysis using parsimony[J]. Sunderland, MA: Sinauer Associates Inc, software, 2002.

[17]SILVESTRO D, MICHALAK I. raxmlGUI: a graphical front-end for RAxML[J]. Organisms Diversity & Evolution, 2012, 12(4): 335-337.

[18]JIANG L, CHEN J, WANG P, et al. Correction: the mitochondrial genomes ofAquilafasciataandButeolagopus(Aves, Accipitriformes): sequence, structure and phylogenetic analyses[J]. Plos One, 2015, 10(10): e0141037.

[19]CORL A, ELLEGREN H. Sampling strategies for species trees: the effects on phylogenetic inference of the number of genes, number of individuals, and whether loci are mitochondrial, sex-linked, or autosomal[J]. Molecular Phylogenetics & Evolution, 2013, 67(2): 358-366.

[20]GURUPRASAD K, REDDY B V, PANDIT M W. Correlation between stability of a protein and its dipeptide composition: a novel approach for predicting in vivo stability of a protein from its primary sequence[J]. Protein Engineering, 1990, 4(2): 155-161.

[21]BOORE J L. Animal mitochondrial genomes[J]. Nucleic Acids Research, 1999, 27(8): 1767-1780.

[22]SOLIGNAC M, MONNEROT M, MOUNOLOU J C. Mitochondrial DNA heteroplasmy inDrosophilamauritiana[J]. Proceedings of the National Academy of Sciences of the United States of America, 1983, 80(22): 6942-6946.