geminiviruses are emerging as one of the most economically important plant viruses in agricultural production ( abou - jawdah et al . , 2006 ; briddon et al . , 2008 ; zhou et al . ,
begomovirus is the largest genus of the family of geminiviridae and is phylogenetically and geographically divided into two groups ; the old world viruses and the new world viruses .
the new world begomovirus consists of two viral genomes , dna - a and dna - b , while most of the old world begomovirus just has one partite dna - a ( briddon et al . , 2008 ) .
about a decade ago , a satellite molecule called dna- was found to associate with some of the old world geminivirus ( saunders et al .
dna- has a genome approximately 1.31.5  kb long , and depends on the helper virus dna - a for its replication , movement , and transmission ( saunders et al . , 2000 ;
it is grouped into sub - viral agents by the international committee on taxonomy of viruses ( ictv ) .
the most typical plant symptoms caused by geminivirus are due to an association of dna- with dna - a , whereas dna - a alone does not lead to severe damage to crops ( cui et al . , 2004 ; briddon et al . ,
c1 gene encoded by dna- were found to suppress host defense systems ( cui et al . , 2005 ) and modulate host development ( yang et al . , 2008 ) , and was believed to be one of the determining factors for geminivirus - induced disease symptom development ( briddon et al . , 2008 ) .
dna- has not been found in the new world ( north american and south american continents ) and is believed to be associated with old world begomoviruses after the geographical divergence of  old  and  new  continents ( mansoor et al . , 2003b ) .
although dna- has relatively a large range of its selection on different species of the helper virus dna - a ( mansoor et al . , 2003a ) , it is proposed to co - evolve with the dna - a component ( briddon et al . , 2008 ) .
recombination plays an important role in geminivirus ( lefeuvre et al . , 2009 ) and dna- evolution ( amin et al . , 2006 ; lefeuvre et al . , 2007 ) .
a fragment of dna- genome infecting tomato was reported to migrate to cotton via recombination with other adaptive dna- molecules ( amin et al . ,
2006 ) , indicating the role of a recombination event in evolution of dna- molecules . because of the important role of recombination in dna- evolution , analysis on recombination events of dna- becomes specially important for understanding this viral evolution and disease epidemic as well as development of potential control strategies . in this paper
, we apply a statistical phylogenetic analysis using a bayesian stochastic method to infer changes in phylogeny along multiple sequence alignments while accounting for rate heterogeneity developed by webb et al .
( 2009 ) to estimate potential recombination spots of dna-. it is the first time that this statistic method has been used on dna- recombination study and give a clear recombination history of dna- recombination . in order to confirm our results , we also apply a statistical phylogenetic method developed by martin et al .
one strain of tomato leaf curl maharashtra betasatellite ( tolcmab ) has a recombination pattern and is possibly recombinant molecule between two strains from two distinct species , papaya leaf curl betasatellite ( palcub ) and tomato leaf curl betasatellite ( tolcb ) , palcub-[in : chi:05 ] ( major parent ) and tolcb-[in : cp:04 ] ( minor parent ) .
a proposed taxonomy of dna- using 78% nucleotide sequence identity as demarcation threshold was accepted and widely used for distinguishing species from strains of dna- ( briddon et al . , 2008 ) .
a recent report showed that different species of dna- associated with tolcd in india are geographically isolated and distributed ( sivalingam et al . , 2010 ) .
the dna- molecules in southern and central india are more closely related to each other than those in northern india . to observe potential recombination events among these geographically related dna- species
, we chose four strains from four distinct species of dna- associated with tolcd in india . among the four strains , tolcbdb-[in;luk;05 ] ( taxon-0 ) and
tolcb-[pk;ryk;97 ] ( taxon-1 ) are from northern india , while palcub-[in;chi;05 ] ( taxon-2 ) and tolcmab-[in;pun;04 ] ( taxon-3 ) are from southern india . in the same report as well as another report ( mazhar et al . , 2009 ) ,
species of tolcbdb and tolcb are closely related in phylogenetic tree , while palcub and tolcmab are sisters ( neighbors ) .
another tolcd associated dna- from indonesia ( taxon-4 ) was chosen as an out group .
other five species of non - tolcd related dna- from eastern asia and southeastern asia ( taxa-5 , 6 , 7 , 8 , and 9 ) were also chosen for the out group .
a proposed taxonomy of dna- using 78% nucleotide sequence identity as demarcation threshold was accepted and widely used for distinguishing species from strains of dna- ( briddon et al . , 2008 ) .
a recent report showed that different species of dna- associated with tolcd in india are geographically isolated and distributed ( sivalingam et al . , 2010 ) .
the dna- molecules in southern and central india are more closely related to each other than those in northern india . to observe potential recombination events among these geographically related dna- species , we chose four strains from four distinct species of dna- associated with tolcd in india . among the four strains , tolcbdb-[in;luk;05 ] ( taxon-0 ) and
tolcb-[pk;ryk;97 ] ( taxon-1 ) are from northern india , while palcub-[in;chi;05 ] ( taxon-2 ) and tolcmab-[in;pun;04 ] ( taxon-3 ) are from southern india . in the same report as well as another report ( mazhar et al . , 2009 ) ,
species of tolcbdb and tolcb are closely related in phylogenetic tree , while palcub and tolcmab are sisters ( neighbors ) .
another tolcd associated dna- from indonesia ( taxon-4 ) was chosen as an out group .
other five species of non - tolcd related dna- from eastern asia and southeastern asia ( taxa-5 , 6 , 7 , 8 , and 9 ) were also chosen for the out group .
first , a data set of 10 dna- genome sequences in fasta format was aligned using clustalw - multialign software with the following parameters : ( gap opening penalty 10.0 , gap extension penalty 0.2 , gap separation penalty range 8 , dna weight matrix : iub ) ( thompson et al . , 1994 ) .
to analyze recombination for dna- from geminiviruses , we used the software package from webb et al .
they applied a hidden markov model ( hmm ) to infer changes in phylogeny along multiple sequence alignments while accounting for rate heterogeneity . under the hmm ,
the hidden states are all possible unrooted tree topologies with the number of leaves n fixed along each site .
the observed state space is { a , c , g , t , }. under the evolutionary model , the evolution of homologous dna / rna sequences ( or protein - coding sequences where the state space is of size 61 ) can be described by continuous time markov chains on a phylogenetic tree . a continuous time markov chain is characterized by a substitution rate matrix , and the phylogenetic tree summarizes the relationships between the species in terms of edge lengths ( times since divergence ) and common ancestors .
the dna sequences are only observed in the leaves , and information on the phylogenetic tree , substitution events ( time and type ) and edge lengths is missing .
the transition matrix p(t ) for a continuous time markov process can be written as exp(qt ) , where q is a parametrized substitution rate matrix which determines the markov process . in this method
the evolutionary model was set as hasegawa  kishino  yano ( hky ) model ( hasegawa et al . , 1985 ) .
the rate matrix q under hky model is written as the following : let   =  { a , c , g , t } and let a , a    , a a  =  1 , denote the stationary distribution of the markov chain .
hky model has substitution rate matrix : where the diagonal elements are such that each row sums to 0 and the two unknown parameters are  and . the software from webb et al .
( 2009 ) estimates the posterior distribution using monte carlo markov chain ( mcmc ) method under the hmm and then it outputs each tree topology with its posterior probability along each site ( see webb et al . , 2009 for details ) .
we have used hky model for phylogenetic analysis on our data sets in this paper , since the hmm software in webb et al .
( 2009 ) uses hky model . also note that we have used the generalized time reversible ( gtr )  + gamma  + invariant model , which is within the 95% confidence interval computed via akaike 's information criteria ( aic ) in the software jmodeltest ( guindon and gascuel , 2003 ; posada , 2008 ) , to reconstructing a ml tree and the ml tree under the gtr  + gamma  + invariant model has the same tree topology as the ml tree under hky model in figure 5 as well as the consensus tree under hky model in figure 4 .
the generated alignment file in phylip format was put in to the hmm software ( webb et al . ,
java -jar st-hmm.jar  with the following parameter ( iterations : 50000 , burn - in : 25000 , rates : 0.001 , 0.003 , 0.01 , 0.03 , 0.1 , 0.3 , 1.0 , 3.0 , 10.0 , 100.0 , lambda : 5 , kappa : 2.0 , tuningpar 0.4 ) .
java -jar sthmmposterior.jar  was used to summarize the posterior distribution , and trees with posterior probability above 0.05 were selected using the command 
the region 11000 nucleotide ( nt ) was found to have a clear pattern of recombination , while the region 10001505  nt seems to have a massive pattern of tree probability . in order to apply phylogenetic analysis to the sequences of 11505 and 10001505  nt of the 10 viral sequences after aligning with the clustalw - multialign software into nexus format , we estimated the posterior distribution under the generalized time reversible ( gtr )  +  model and hky model , and we estimated the maximum likelihood estimators .
first we applied a software mrbayes ( bronquist and huelsenbeck , 2003 ) to analyze the split of different taxa on the most consensus tree under the gtr  +  and hky models .
647300 generations were sampled for 11505  nt alignment , while 3600000 generations were sampled for 10001505  nt alignment .
we followed the recommendation of mrbayes which suggests running the chains until the standard deviation of the chains split frequencies is less than 0.01 .
in addition , to verify our results we applied the software rdp3 ( martin et al . , 2005b ) to the same data sets . sequence alignment in phylip format was used as input for rdp3 .
parameters were set to default used by rdp3 . in the software rdp3 they have implemented several different methods to find recombination sites , rdp ( martin and rybicki , 2000 ) , geneconv ( padidam et al . , 1999 ) , bootscan ( martin et al . , 2005a ) , maxchi ( smith , 1992 ) , chimaera ( posada and crandall , 2001 ) , siscan ( gibbs et al . , 2000 ) , and 3seq ( boni et al . , 2007 ) .
the software rdp takes basically three steps : first they discard non - informative sites from the input data sets and then for every triplet of taxa { a , b , c  } , from the data set , choose the sister a and b. second , they use a window of user - defined width moved among the aligned sub - sequences one nucleotide at a time and take an average percentage identifying each of the three possible sequence pairs among { a , b , c } at the each position .
third , the probability that the nucleotide arrangement in the identified region that results in a , b appearing more closely related to c may have occurred by chance is computed using a binomial distribution . the software geneconv is based on an earlier statistical approach for detecting gene conversion ( sawyer , 1989 ) .
they use the term fragment for an aligned or homologous pair of segments in the input alignment . in the process ,
the highest - scoring fragments in the given alignment are listed and assigned p - values based on the assumption of a random distribution of polymorphic sites .
they assign scores as follows : first , all sites that are monomorphic in the alignment are discarded so that only polymorphic sites are considered . secondly , for a given pair of sequences , matching bases are scored as + 1 and mismatches as m , where m depends on the pair of sequences .
fragments are assigned p - values similar to the blast procedure ( altschul et al . , 1990 ; karlin and altschul , 1993 ) .
this p - value is an approximation of the proportion of permutations of the polymorphic sites for which that pair of sequences has some fragment with the observed score or larger ( sawyer , 1989 ) .
the software bootscan takes two phases :  scanning phase  and  detection phase .
 in  scanning phase  first they discard non - informative sites from the input data sets and in each window of user - defined width move among the given aligned sequences .
it makes bootstrap samples and compute rooted upgmas by definition rooted or mid - pointed neighbor - joining ( nj ) trees . in
 detection phase  every combination of triplets is individually examined for bootstrap evidence that one of the sequences may be alternatively more closely related to each of the other two sequences at different positions along its length .
the probability that the pattern of sites within a potential recombinant region could have occurred by a chance distribution of mutations is approximated using a bonferroni corrected version of the binomial distribution .
the software maxchi considers only polymorphic sites : for a given position of the moving window on the input sequence alignment and for a given pair of sequences , a chi - square statistic is computed to compare two proportions : the proportion of sites at which the sequences agree in the left half - window and the proportion of sites at which the sequences agree in the right half - window .
discordance between these two proportions may reflect a recombination event in the history of the two sequences .
the maximum chi - square over all sequence pairs is recorded as a summary of the evidence for recombination at the window center .
significance of observed chi - square statistics is assessed by a monte carlo permutation test .
the software chimaera is also a modification of maynard smith 's maximum  method ( wiuf et al . , 2001 ) with only variable sites .
the p - value equals the number of times the original statistic is smaller than the statistic from permuted alignments divided by the number of permutations .
for all calculations , a sliding window was used , with the width of the window set to the number of polymorphic sites divided by 1.5 .
the software siscan uses a similar idea as algorithms implemented in maxchi and chimaera , but instead of using contingency tables they use gaussian distribution and use z - score to compute the p - value .
the software 3seq is similar to rdp discards non - informative sites from the input data sets and then for every triplet of taxa { a , b , c } , from the data set , it chooses the sister a and b : two parent sequences that may have recombined , with one or two breakpoints , to form the third sequence ( the child sequence ) .
excess similarity of the child sequence to a candidate recombinant of the parents is a sign of recombination ; they take the maximum value of this excess similarity as the test statistic .
then they rapidly calculate the distribution of the excess similarity and using this method they estimate the p - value .
the most consensus trees found with the 11505 and 10001505  nt alignment were the same as the most dominant tree found with the hmm software ( the pink tree in figure 1 ) .
this is the tree topology written in pink ( series 14 ) in figure 2 .
this is the most likely tree topology from position 1 to 140 and position 300 to 1000 .
( 2005b ) indicate a potential recombination event among taxa 0 , 1 , 2 , and 3 in the red rectangle .
also the ml tree estimated by the software phyml has the same tree topology under hky model as well as the consensus tree estimated by the software mrbayes under hky and gtr  + . then we estimated the maximum likelihood ( ml ) tree from the whole alignment ( including position 1 through position 1505 ) .
next we infer phylogenetic tree using maximum likelihoods method , using phyml v3.0 software ( guindon and gascuel , 2003 ) , with all settings default , namely the evolutionary model is hky model , the tree topology search operation method is nearest neighbor interchange ( nni ) , and the starting tree was computed using bionj ( gascuel , 1997 ) . to analyze the splits of different taxa on the ml tree we applied bootstrapping on the columns of each alignment with the bootstrap sample size 1000 .
the ml tree found with the 11505  nt alignment was the same as the most dominant tree found with the hmm software ( the pink tree in figure 1 ) . from position 1
to position 141 and from position 312 to position 1000 , the tree topology in figure 1 has almost probability 1.0 ( see figure 2 ) .
note that the estimated ml tree and the estimated consensus tree reconstructed with the whole sequences from an estimated posterior distribution have the same tree topology .
however , from position 141 to position 311 in the alignment , the tree topology in figure 3 has almost probability 1.0 ( see figure 2 ) . the robinson
 foulds ( rf ) distance ( robinson and foulds , 1981 ) between the tree topology in figure 3 and tree topology in figure 1 is 6 .
note that the largest possible rf distance for trees with n taxa is 2n   6 which is 14 in our case ( the normalized rf distance between these tree topologies is 0.43 ) .
thus we do not think this happened because of the low support of a split but this seems to indicate strongly that around position 142 and position 311 there are possible recombination sites .
the figure shows an estimated probability of each tree topology along each site computed using the software from webb et al .
the label of  series i  for i  =  1,  ,17 in the figure represents each different tree topology .
the y - axis represents the probability for each tree topology and the x - axis represents position number .
the tree written in pink is in figure 1 and the tree written in the dark blue dominating from position 140 to 300 is in figure 3 .
( 2005b ) indicate a potential recombination event among taxa 0 , 1 , 2 , and 3 in the red rectangle . in order to compute the support for each split
we have also computed the consensus tree using the software mrbayes ( figure 4 ) and the ml tree using phyml ( figure 5 ) . for the consensus tree we used the posterior distribution and for the ml tree we use the bootstrap with the sample size 1000 to compute the support for each split .
they have the same tree topology as the tree in figure 1 and the support for each split in the ml tree and the consensus tree has very high probability .
especially , the probability of each split on the consensus tree estimated with the whole sequences under hky is 1.0 ( 100% ) .
( even though one of the splits on the ml tree reconstructed with the whole sequences under hky has about 90% of its support all other splits have strong support ; figure 5 . ) . the consensus tree estimated by the software mrbayes under hky from the whole alignment ( including position 1 through position 1505 ) .
the consensus tree estimated under the gtr +  also has the same tree topology but it has smaller probabilities of some splits . note that the tree topology of the consensus tree is the same as the tree topology of the ml tree in figure 5 and the tree topology in figure 1 . the ml tree estimated by the software phyml under hky model from the whole alignment ( including position 1 through position 1505 ) .
the number in each split represents the probability of the split estimated by bootstrapping with the bootstrap sample size 1000 .
note that the tree topology of the ml tree is the same as the tree topology of the consensus tree in figure 4 and the tree topology in figure 1 .
the mutation rates along each site are also estimated by the software from webb et al .
( 2009 ) and it seems that the mutation rates are between 0.1 and 0.3 ( figure 6 ) .
the figure shows an estimated probability of each mutation rate along each site computed using the software from webb et al .
the y - axis represents the probability for each mutation rate and the x - axis represents position number .
rdp3 estimated a similar recombination event , where a small genome fragment of tolcmab-[in;pun;04 ] ( taxon-3 ) ( position 142311 in alignment ) is migrated from tolcb-[pk;ryk;97 ] ( taxon-1 ) , as circled by red rectangle in figure 7 .
rdp3 uses multiple methods for recombination estimation , and the average p - value from different methods are listed below ( table 2 ) .
red rectangles indicate the same event inferred by a hmm method from webb et al .
2005b ) inferring the recombination event between tolcmab-[in;pun;04 ] and tolcb-[pk;ryk;97 ] at position of 142311 .
we used rdp ( martin and rybicki , 2000 ) , geneconv ( padidam et al . , 1999 ) , bootscan ( martin et al . , 2005a ) , maxchi ( smith , 1992 ) , chimaera ( posada and crandall , 2001 ) , and 3seq ( boni et al . , 2007 ) .
we set parameters for each software as follows ; 
 rdp 
 reference sequence no : window size : 30 ; detect recombination between sequence identity : 0100% ; 
 geneconv : 
 sequence option : treat blocs as one polymorphism ; g - scale : 1 ; maximum number of global fragments listed per sequence pair : 2000 ; maximum number of pairwise fragments listed per sequence pair : 0 ; minimum aligned fragment length : 1 ; minimum polymorphisms in fragments : 2 ; minimum pairwise fragments score : 2 ; maximum number overlapping fragments : 1 ; 
 bootscan : 
 window size : 200 ; step size : 20 ; use distances ; number of bootstraps replicates : 100 ; random number seed : 3 ; cutoff percentage : 70 ; transversion rate ratio : 0.5 ; coefficient of variation : 1 ; 
 maxchi : 
 window size : 70 ; gaps : no ; 
 chimaera : 
 window size : 60 ; and 
 3seq : 
 sequences are circular ; highest acceptable p - value : 0.05 ; bonferroni correction ; number of permutations : 0 ; use seqgen parametric simulations .
we first reported a potential recombination event between taxa 1 , 2 , and 3 , indicating that the strain tolcmab-(in;pun;04 ] ( taxon-3 ) from tolcmab is a recombinant of two strains from two different species , tolcb-[pk;ryk;97 ] ( taxon-1 ) and palcub-[in;chi;05 ] ( taxon-2 ) .
as one study reported , tolcmab-[in;pun;04 ] ( taxon-3 ) and palcub-[in;chi;05 ] ( taxon-2 ) are closely related in their phylogeny compared to other species ( sivalingam et al . , 2010 ) .
our study showed that tolcmab-[in;pun;04 ] ( taxon-3 ) shares sequence identity mainly with palcub-[in;chi;05 ] ( taxon-2 ) , while a small portion of its genome ( position 141312  nt in the alignment ) is potentially migrated from another species , tolcb-[pk;ryk;97 ] ( taxon-1 ) .
our results indicate a recombination event happened between a northern india dna- strain tolcb-[pk;ryk;97 ] ( taxon-1 ) and a southern india dna- strain palcub-[in;chi;05 ] ( taxon-2 ) , resulting a new strain tolcmab-[in;pun;04 ] ( taxon-3 ) which was found in southern india .
different geographic locations provide different physiology of host , weather conditions , helper viruses , and so on .
the phylogenetic relationship among tolcb-[pk;ryk;97 ] ( taxon-1 ) , palcub-[in;chi;05 ] ( taxon-2 ) , and tolcmab-[in;pun;04 ] ( taxon-3 ) coincides with their distinct geographic relationship , suggesting that different genetic information on the viral genome from northern india or southern india may already adapt to their geographic distribution ( figure 8) .
however , although the recombination event lead to the possible emergence of a new strain in a different epidemic location in india , it still has a stronger relationship within its parents geographically and phlegmatically than other strains which are epidemic in other asian countries .
the geographic distribution of four betasatellites , tolcbdb-[in;luk;05 ] , tolcb-[pk;ryk;97 ] , palcub-[in;chi;05 ] , and tolcmab-[in;pun;04 ] , associated with tolcd in the india sub - continent .
c1 protein , product of the c1 gene , can alter leaf development and suppress plant defense systems during infection ( cui et al .
the recombination happened in approximate 100220  nt of the genome ( 141312 in alignment ) , which partially covers the c - terminal of c1 gene on the betasatellite .
tolcb-[pk;ryk;97 ] ( taxon-1 ) c1 has 118 amino acids , while c1 of palcub-[in;chi;05 ] ( taxon-2 ) has 122 amino acids .
the recombination event leads to a c1 protein of tolcmab-[in;pun;04 ] ( taxon-3 ) with 118 amino acids , missing the 6 amino acids from major partent palcub-[in;chi;05 ] ( taxon-2 ) on the c - terminal of c1 , instead having 2 amino acids from c - terminal of c1 on minor parent tolcb-[pk;ryk;97 ] ( taxon-1 ) .
although functions of different domains of c1 were unknown , the recombination on c - terminal of c1 might modulate its function involving in virus 
dna- was known to be capable to adapt to a new helper virus from distinct geographic location by modifying its genome ( nawaz - ul rehman et al . , 2009 ) .
the genetic modification on this southern indian dna- strain tolcmab-[in;pun;04 ] ( taxon-3 ) via a recombination event may contribute to the fitness of this dna- strain on its host .
the advantage of our study is that estimating of phylogenetic tree through out the alignment by hmm method provide a clear history of dna- recombination .
it is the first time that researches on dna- recombination use such statistic method and give this clear recombination history .
our study also provides a way to understand dna virus evolution through recombination events . from our results , it is likely that the specie of tolcmab is a result of recombination from two different species , namely tolcb and palcub .
such recombination event contributed to the occurrence of new dna- species as well as the evolution of dna-. by providing the recombination history together with geographic information , we could link the phylogeny information to the geographic information of dna- strains , thus help us understand evolution and epidemic of the virus .
the authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest .
 although exchange of genetic information by recombination plays an important role in the evolution of viruses , it is not clear how it generates diversity . 
 understanding recombination events helps with the study of the evolution of new virus strains or new viruses . 
 geminiviruses are plant viruses which have ambisense single - stranded circular dna genomes and are one of the most economically important plant viruses in agricultural production . 
 small circular single - stranded dna satellites , termed dna- , have recently been found to be associated with some geminivirus infections . in this paper 
 we analyze several dna- sequences of geminiviruses for recombination events using phylogenetic and statistical analysis and we find that one strain from tolcmab has a recombination pattern and is a recombinant molecule between two strains from two species , palcub-[in : chi:05 ] ( major parent ) and tolcb-[in : cp:04 ] ( minor parent ) . 
 we propose that this recombination event contributed to the evolution of the strain of tolcmab in south india . the hidden markov chain ( hmm ) 
 method developed by webb et al . 
 ( 2009 ) estimating phylogenetic tree through out the whole alignment provide us a recombination history of these dna- strains . 
 it is the first time that this statistic method has been used on dna- recombination study and give a clear recombination history of dna- recombination . 
 
Introduction
Data set
Materials and Methods
Results
Conclusion
Discussion
Conflict of Interest Statement