-
作为一类知名的昆虫致病菌,球孢白僵菌Beauveria bassiana由于其广泛的宿主范围(超过700种昆虫[1])及感染不同昆虫的毒力多样性,被世界各地广泛用于害虫生物防治[2],在农林业生产中展现出了巨大的价值.其中,作为全球最大的害虫防治项目,我国利用球孢白僵菌防治亚洲玉米螟Ostrinia furnacalis和马尾松毛虫Dendrolimus punctatus[3-4],取得了显著成效.此外,球孢白僵菌能够产生大量的次生代谢物,这些代谢物在工业上也具有巨大的应用价值.
球孢白僵菌ARSEF 2860的全基因组测序计划由中国科学院上海生命科学研究院和浙江大学的科学家在2012年完成[5].通过第二代测序技术,获得了76.6倍覆盖度的基因组测序数据,并组装成33.7 Mb基因组大小的242条scaffolds数据,以及10 364个蛋白质编码基因被注释.通过interproscan分析,7 283个基因被归类为3 002个基因家族,这表明球孢白僵菌存在着普遍的基因家族膨胀现象[5].许多物种特异的毒力基因,诸如编码一类小型富含半胱氨酸分泌蛋白的SSCPs基因,在球孢白僵菌基因组中被鉴定[5].上述2种特征被认为与其宿主范围和致病策略有关.此外,球孢白僵菌具有特殊的表达转录因子(TFs)用于调控和激活特定基因,以适应不同的昆虫宿主[5].
生物体内,在遗传物质(DNA或RNA)编码蛋白质(protein)的翻译过程中,遗传密码(genetic code)是一套起着关键作用的规则.标准遗传密码包含61个编码、20种氨基酸的密码子和3个终止密码子.由于除蛋氨酸Methionine和色氨酸Tryptophan外,其余氨基酸都可以由2个其以上的密码子进行编码,使得密码子具有简并性.这些简并的同义密码子synonymous codons虽然编码同一个氨基酸,但在不同物种中使用的频率却不尽相同[6],由此产生的现象称之为密码子使用偏好性(Codon Usage Bias).密码子偏好性广泛存在于生物体中,并受核苷酸组成[7]、翻译过程[8]、tRNA丰度[9]、基因功能[10]、mRNA结构[11]、基因长度[12]、蛋白质疏水结构[13]、环境温度[14]等多种因素的影响,尤其是基因突变与自然选择压力之间的平衡决定了密码子的使用偏好[15].因此,针对遗传密码子使用偏好的分析将能更好地理解不同生物体在分子进化、环境适应、基因组特征等方面的异同,是基因组学中的一个重要课题.
考虑到球孢白僵菌巨大的经济价值及其宿主范围广的特性,本文将通过生物信息学方法分析球孢白僵菌基因组数据,确定其密码子使用模式,明确其密码子偏好性,并描述影响其密码子使用偏好的主要因素,从而了解其基因组特征,为球孢白僵菌的真菌-昆虫互作研究奠定科学基础.
全文HTML
-
本文采用的球孢白僵菌ARSEF 2860基因组数据(序列号NZ_ADAH00000000)及10 364个蛋白质编码基因序列,均下载自GenBank数据库(http://www.ncbi.nlm.nih.gov/genbank).
-
运用自编PERL脚本(可从网址https://github.com/hxiang1019/calc_GC_content.git下载),计算上述10 364个球孢白僵菌编码基因及其密码子第一、第二、第三位碱基上的GC含量,分别记为GCcds,P1,P2,P3.将各基因P1和P2的平均值P12作为纵坐标,以P3作为横坐标,进行中性绘图分析(Neutrality analysis),以确定密码子3位置间核苷酸组成相关性,进而检测突变选择平衡是否显著影响密码子使用偏好[16].中性图中每个点代表一个编码基因,当该基因的P12与P3显著相关时,表明其变异更多是受到自身遗传漂变的非定向突变压力作用,而很少受到或不受到外界选择下定向突变压力的作用,因而成中性,该点应位于对角线上;而如果相关性不显著,该点应位于水平线上,表明其变异受到很大程度的定向突变压力作用.中性绘图分析回归线越呈水平状态,其斜率越接近于0,表明P12与P3越不相关,密码子使用偏好就更多地受定向突变压力影响;反之,斜率越接近于1,表明P12与P3相关,密码子偏好就更多地受遗传漂变作用[17].
此外,本研究也计算了10 364个基因在密码子第三位上的核苷酸组成,记为A3,U3,C3,G3,进而以AU偏量[A3/(A3+U3)]和GC偏量[G3/(G3+C3)]作为横纵坐标,进行PR2分析(Parity rule 2 analysis),以检测密码子在嘌呤和嘧啶使用上的偏好[18].在PR2分析中,每个点代表 1个编码基因,若该基因第三位密码子使用A多于U,则位于图中上半部分;若使用G多于C,则位于图中右半部分.而若某一基因的密码子偏好性仅受到核苷酸组成这1个因素的影响,那该基因的密码子第三位碱基上G3和C3及A3和U3的使用频率应该一致,该基因应处于中间位置[18].
-
运用CodonW软件(John Peden,http://www.molbiol.ox.ac.uk/cu,版本号1.4.2),计算10 364个球孢白僵菌编码基因的各项密码子指数,如密码子适应指数(Codon Adaptation Index,CAI)、密码子偏好指数(Codon Bias Index,CBI)、最优密码子使用频率(Frequency of optimal codon,Fop)、有效密码子数(Effective Number of Codon,ENc)、同义密码子第三位GC含量(GC3s)、氨基酸疏水性指数(Gravy)、氨基酸芳香性指数(Aromo)等.其中,密码子适应指数CAI一般用于反映基因的表达情况,其数值介于0到1间,该值越大表明基因的偏好性越强,表达水平也越高[19].密码子偏好指数CBI也能衡量基因的偏好性和表达程度,并能反映出基因中最优密码子的组成情况,其值介于0到1间时,数值越大基因偏好性越强;而当CBI值小于0时,则表示最优密码子使用频率比平均使用频率还要少,表现出弱偏好性[20].有效密码子数ENc是分析密码子使用偏好最具参考意义的指数,其值通常与CAI和CBI相关,ENc值越接近20,说明基因偏性越强,表达程度越高;越接近61则刚好相反[21].
此外,我们还以ENc值为纵坐标、GC3s值为横坐标,绘制ENc图(ENc-plot),用于探测各基因的密码子使用是否只受其核苷酸组成压力的影响,还是有其他因素在起作用[21].在ENc图中,曲线ENc=2+GC3s+29/[GC3s2+(1-GC3s)2]表示密码子偏好性仅由核苷酸组成所决定的标准情况.若某基因沿着或靠近标准曲线分布,表明该基因的密码子偏好性仅仅受GC含量的影响;若远离标准曲线,则表明该基因还受到其他诸如选择压力等因素的影响[21].
-
我们使用CodonW软件,计算了球孢白僵菌基因组的相对同义密码子使用值(relative synonymous codon usage,RSCU).作为密码子偏好性研究的一个重要指数,RSCU值是密码子实际观测值与理论观测值的比值.当RSCU为1时,表示实际观测值与理论值相同,即各同义密码子使用频率相同,密码子使用无偏好.当RSCU大于1时,表示该密码子使用频率高于理论值;反之则使用频率低于理论值.
在RSCU的基础上,我们将球孢白僵菌编码基因中CAI值居于最高、最低的5%样本取出生成2个数据集,分别作为高、低表达组,再利用卡方检验(chi-squared test)比较2组对应的RSCU值,若某密码子在高表达组的使用频率显著地(p-value<0.05)高于低表达组,则确定其为最优密码子(optimal codon)[22],并将RSCU值大于1的密码子认定为偏好密码子(preferred codon),而RSCU值小于0.10的密码子则认定为使用频率低的稀有密码子(rare codon).
基于球孢白僵菌各编码基因的RSCU值,利用CodonW软件进行对应分析(Correspondence analysis,COA)[23],并比较其基因组内部的变异情况,进而分析影响密码子偏好的因素.对应分析是一种多元统计方法,用于调查基因组内各基因间RSCU值的差异[23].通过多维空间,对应分析不仅能展示各基因的分布,还能反映影响密码子偏好性各种因素间的主次关系,从而探测影响密码子偏好的因素是否多元化[13].
文中相关性分析、方差分析、显著性检验等数据统计分析及绘图由SPSS 18.0和Microsoft Excel 2010完成.
1.1. 基因组及编码基因序列数据
1.2. 核苷酸组成及中性绘图和PR2分析
1.3. 密码子指数计算及ENc绘图分析
1.4. RSCU计算及对应分析
-
球孢白僵菌基因组的多项密码子特征,包括基因组GC含量(GCgenome)、编码基因GC含量(GCcds)、密码子第一、二、三位GC含量(P1,P2,P3)、以及密码子适应指数(CAI)、密码子偏好指数(CBI)、最优密码子使用频率(Fop)、有效密码子数(ENc)等被计算(表 1).其中,球孢白僵菌CAI平均值(0.079)和CBI平均值(0.018)接近于0,ENc平均值(48.82)较大,表明球孢白僵菌密码子使用偏好较弱,基因表达水平偏低.
核苷酸组成是影响密码子使用偏好的1个重要因素,球孢白僵菌的编码基因GC含量GCcds(57.5%)比其基因组GC含量GCgenome(51.5%)更高(表 1),表明编码区比非编码的序列更为稳定.相关性分析显示球孢白僵菌的GCcds不仅与其P1(r=0.623,p-value<0.01)和P2(r=0.480,p-value<0.01)显著相关,还与P3(r=0.869,p-value<0.01)显著相关,表明球孢白僵菌编码基因GC含量受到3个密码子位置的影响.有效密码子数ENc和密码子3个位置间存在相关性,P1(r=0.343,p-value<0.01),P2(r=0.130,0.01<p-value<0.05),P3(r=0.812,p-value<0.01),表明核苷酸组成可能是影响球孢白僵菌密码子使用偏好的因素.
-
为了检验球孢白僵菌密码子使用偏好是否仅受核苷酸组成的影响,课题组绘制了ENc图(图 1).在图 1中,除少部分基因沿着标准曲线分布外,大部分球孢白僵菌基因都位于标准曲线下方,表明核苷酸组成不是影响球孢白僵菌密码子偏好性的唯一因素,还存在其他因素的作用.某些基因ENc值等于61,位于图 1中最顶端,表明它们编码氨基酸的所有61个密码子都在使用,这些基因无密码子偏好性.此外,由于绝大多数球孢白僵菌基因的同义密码子第三位GC含量(GC3s)超过了50%,在ENc图中表现出了基因组右倾的情况(图 1).
-
为了调查球孢白僵菌基因组的密码子使用情况,课题组计算了球孢白僵菌的相对同义密码子使用值RSCU(表 2).
调查RSCU大于1的使用频率相对高的偏好密码子(表 2,下划线表示),发现球孢白僵菌基因组强烈偏好核苷酸G/C结尾密码子(表 3,25个偏好密码子中有23个都是G/C结尾的),尤其是C结尾的(表 3,25个偏好密码子中有16个是C结尾的).卡方鉴定出的最优密码子(表 2,符号*表示)也表现出同样的情况,22个最优密码子中有21个都是G/C结尾的、15个是C结尾的(表 3),其中CUC,GUC,UCC,CCC,ACC,CGC,GGC,GCC这8个RSCU大于2的最优密码子全是C结尾的.此外,3个稀有密码子UUA,AUA,GUA也被鉴定(表 2,符号-表示),它们以A结尾.
球孢白僵菌基因组G/C结尾密码子的偏好性表现与其基因组GC含量高(51.5%)有关,尤其与第三位密码子GC含量P3非常高(66.7%)(表 1)密切相关,也表明了核苷酸组成确实是影响球孢白僵菌密码子使用偏好的一个重要因素.
-
既然球孢白僵菌的最优密码子和偏好密码子使用C结尾远多于G结尾(表 3),为此,利用PR2图探讨密码子在嘌呤(A、G)和嘧啶(U、C)使用上的偏好性(图 2).如图 2所示,大多数基因位于图像的左下方.计算球孢白僵菌的GC偏量[G3/(G3+C3)]和AU偏量[A3/(A3+U3)]分别为0.436和0.439,表明密码子第三位偏好U高于A,C高于G,即嘧啶高于嘌呤.
球孢白僵菌密码子中嘌呤和嘧啶的使用频率具有不对称性,表明除核苷酸组成外还有其他因素影响其密码子使用偏好,这与ENc图(图 1)的结果相呼应.此外,个别基因位于PR2图的横坐标轴上(图 2),它们被认为不使用A结尾密码子.
-
既然球孢白僵菌密码子偏好性不只受到核苷酸组成的作用,因此通过中性绘图分析(图 3),检测突变选择压力是否为影响密码子使用偏好的一个重要因素.如中性分析图所示(图 3),球孢白僵菌基因组的斜率为0.117,趋近于0,其P12与P3相关系数为0.105,双尾检验未达到显著水平(p-value>0.05),表明P12与P3不相关,认为球孢白僵菌密码子偏好性受到外界自然选择下的定向突变压力作用多于自身遗传漂变下的非定向压力作用.
-
除上述核苷酸组成和突变选择压力外,课题组利用对应分析检测球孢白僵菌密码子偏好性是否还受到其他因素的作用,以及这些因素的主次关系.除对应分析图所展示的第一、二轴外(图 4),球孢白僵菌的前四轴分别占据了总变异度的25.13%,11.46%,4.10%,3.03%,表明球孢白僵菌密码子偏好性以第一轴(约占1/4)和第二轴(约占1/10)为主.因此,将前4轴与各项密码子指数进行相关性分析(表 4),发现第一轴与代表核苷酸组成的GCcds(r=-0.666,p-values<0.01)和GC3s(r=-0.915,p-values<0.01),以及代表基因功能的氨基酸疏水性指数Gravy(r=-0.245,0.01<p-values<0.05)和芳香性指数Aromo(r=-0.203,0.01<p-values<0.05)均显著相关,表明核苷酸组成和基因功能是影响球孢白僵菌密码子偏好性的最主要因素.同时,第二轴与CAI,CBI,Fop,ENc等密码子指数显著相关(p-values<0.01),表明也存在着其他因素对球孢白僵菌密码子使用偏好的影响.这与上述分析结果一致.
2.1. 密码子指数分析
2.2. ENc绘图分析
2.3. RSCU分析
2.4. PR2绘图分析
2.5. 中性绘图分析
2.6. 对应分析
-
密码子使用偏好是各物种在长期进化过程中形成的特征,是基因组及进化模式的特有反映[24],其相关研究已在诸多生物体中被报道.利用已公布的球孢白僵菌基因组及10 364条编码基因序列,课题组开展了球孢白僵菌密码子偏好性研究.通过密码子指数CAI,CBI,ENc的计算,发现球孢白僵菌密码子使用偏性较弱,基因表达水平也较低.尽管如此,根据RSCU值鉴定最优密码子,发现球孢白僵菌基因组使用G/C结尾密码子的频率远远大于A/U结尾的,几乎没有A/U结尾最优密码子的存在(表 2和表 3).球孢白僵菌这种G/C偏好模式与其基因组的高G/C含量一致,尤其与其第三位密码子G/C含量(P3)高达66.66%有关(表 1).大部分物种都表现出这种模式,即基因组G/C含量高的物种倾向于使用G/C丰富的密码子,而A/U高的偏好使用A/U丰富的密码子[25].因此,影响球孢白僵菌密码子偏好性的一个重要因素就是其核苷酸组成,而球孢白僵菌对应分析(图 4)及相关性分析(表 4)结果也支持了该结论.
对密码子第三位的碱基进行PR2分析,发现球孢白僵菌的偏好性还表现在倾向于使用嘧啶结尾的密码子,即U使用频率高于A,C使用频率高于G(图 2).因此,尽管核苷酸组成被认定为影响球孢白僵菌密码子的一个重要因素,但其密码子在嘌呤和嘧啶上的不对称性表明了还存在着其他因素的作用[18].而ENc图(图 1)和对应分析(表 4)的结果也支持了该结论.
由于Necsulea等[26]认为突变选择压力是引起密码子不对称性的最大因素,为此课题组对球孢白僵菌基因组进行了中性绘图分析,结果表明球孢白僵菌密码子受到外界自然选择下的定向突变压力作用(图 3),这与一类营寄生生活的真菌是一样的,它们的密码子也受到强烈的定向突变压力作用[27].因此,认为球孢白僵菌密码子受到的定向突变压力可能与其寄生生活有关,压力来自于其宿主.
为了维持基因稳定,大多数的序列突变都是同义突变.然而,在外界选择的定向突变压力下,P12与P3不相关、非同义突变增多、基因变化加快、毒力基因多样化,致使球孢白僵菌可以更快地调节自己,以适应其宿主范围广泛的寄生生活方式. 表 4中第一轴与代表基因功能的氨基酸疏水性指数Gravy和芳香性指数Aromo的相关性,揭示了影响球孢白僵菌密码子使用偏好的另一个主要因素是基因功能,这从一个侧面印证了定向突变压力下的基因功能变化确实是其基因组的一个重要特征.此外,球孢白僵菌基因组中存在着大量的转座元件[5],可能也是在定向突变压力下加快自身基因组进化和遗传变异的一种方式,以适应其在700多种昆虫宿主中的寄生.当然,更可靠的结论需要通过更多的关于球孢白僵菌及近源物种基因组的诸如Ka/Ks计算、正选择检验、系统进化树比较等大量分子进化研究来完成.