staphylococcus aureus strain sa957 was obtained from chang gung children s hospital in taiwan . for genome sequencing , a mono - clone of sa957 was picked from lb agar and grown in liquid lb broth at 37 c overnight . the overnight culture was diluted 40-fold with fresh lb broth and incubated at 37 c to reach log phase . the genomic dna was then extracted by qiaamp dna mini kit ( qiagen , valencia , ca , usa ) and subjected to solid mate - pair library construction . the library was introduced to a solid 3 plus sequencer for 2  50 bp sequencing . cp003603.1 ) was used as the reference genome and the socs package ( ondov et al . 2008 ) was used to map the solid reads and to detect single - base substitutions . in detail , the solid reads were firstly filtered by the script filterreads.pl ( within socs package ) with the quality threshold set as q15 . according to the average error rate obtained from a pilot mapping , we decided to keep the first 42 bp of each read and ran the mapping again . the variant mode  short variants  was set for detecting valid sequence - space mismatches indicated by color - space mismatches . the results were then outputted as the coverage at each position and the valid sequence - space mismatches detected from each strand of the chromosome . only reads with unique hits against the genome were included in the downstream analysis so that spurious alignments were avoided . in addition , the positions with coverage smaller than 500 or larger than 15,000 were discarded because with such low or high coverage there was a higher likelihood of misalignment . if a given position had more than one type of substitution meeting the above criteria , only the major substitution was counted . to differentiate mutations from the sequencing errors , we defined firstly that a mutant allele must have at least five supporting reads on both the forward and reverse strands with the aim of reducing the likelihood that a mutant allele came exclusively from duplicate molecules amplified during library construction . then , the bernoulli test was performed for each variant site . in detail , the function binomdist in microsoft excel was used , where the number_trail was the total coverage of the site , the number_success was the number of reads supporting the mutant allele , and the prob_success was the presumed error rate . the chromosome of sa957 was 2,789,538 bp long . according to the cumulative gc skew and at skew , which were calculated using genskewapp.jar ( downloaded from http://genskew.csb.univie.ac.at , last accessed march 13 , 2013 ) , the oric and ter site were presumed to be at 0 and 1,414 kb , respectively . the substitutions in the second half of the chromosome were converted into their reverse complement so that the substitution pattern could be entirely described for the leading strand . the frequency of a certain substitution was defined to be the number of this type of substitutions divided by the total number substitutions . for example , f(c  a ) = number(c  a)/number(all ) . the trend for how each base changed as a percentage during short- and long - term evolution were then predicted based on the frequency of the 12 possible substitution types . in detail , 
 net frequency ( a ) = f(c  a ) + f(g  a ) + f(t  a )  f(a  c )  f(a  t ) f(a  g)net frequency ( t ) = f(a  t ) + f(c  t ) + f(g  t ) f(t  a )  f(t  c )  f(t  g)net frequency ( c ) = f(a  c ) + f(t  c ) + f(g c )  f(c  a )  f(c  t )  f(c  g)net frequency ( g ) = f(a  g ) + f(t  g ) + f(c  g ) f(g  c ) 
  net frequency ( a ) = f(c  a ) + f(g  a ) + f(t  a )  f(a  c )  f(a  t ) g ) net frequency ( t ) = f(a  t ) + f(c  t ) + f(g  t ) f(t  a )  f(t  c )  f(t  g ) net frequency ( c ) = f(a  c ) + f(t  c ) + f(g  c )  f(c f(c  g ) net frequency ( g ) = f(a  g ) + f(t  g ) + f(c  g ) f(g  c ) the socs package does not support gapped alignment and therefore we used shrimp2 ( david et al . 2011 ) to create a gapped alignment that would allow the detection of short indel events . strata  was adopted to make sure that only the highest scoring mapping for a given read was outputted . the score to extend a gap along the genome sequence was set at 2 and the score to extend a gap along the read sequence was set at 1 . indels that were located in the first / last 5 bp of the read were not considered to be valid ones and such hits were discarded from the alignment file using a home - made perl script . 2010 ) was used to left - align indels so that the same indel would not be represented as a different variant . the bambino program ( edmonson et al . 2011 ) was used to process the realigned file and detect the indels within it . the following options were set : the minimum number of reads passing all quality filters required at a variant site was set at two ; the required minimum number of unique read names supporting the mutant allele was set at two ; and other options were set at default . only indels that were supported by reads in both the forward and reverse direction were counted . then bernoulli test was performed on each indel site ( p < 0.001 ) . only the sites with significantly more reads supporting the indel mutant proceeded to the subsequent analysis . cp003603.1 , multi - locus - sequence - type st59 ) together with s. aureus strains fpr3757 ( nc_007793.1 , st8 ) , jkd6159 ( nc_017338.1 , st93 ) , lga251 ( nc_017349.1 , st425 ) , mrsa252 ( nc_002952.2 , st36 ) , mw2 ( nc_003923.1 , st1 ) , n315 ( nc_002745.2 , st5 ) , and rf122 ( nc_007622.1 , st151 ) were used to infer the genomic sequence of an ancestral s. aureus . the derived backbone regions , namely regions that are conserved in all these strains , were extracted from the output backbone file and were further aligned by mafft program ( katoh et al . the length of the concatenated conserved regions totaled 2,393 kb , which is approximately 85% of the whole genome of each strain . the concatenated alignment of the conserved regions was inputted into the baseml program from the paml package ( yang 2007 ) and the marginal reconstruction algorithm was used to infer the genome sequences of the ancestral s. aureus . a neighbor - joining tree was produced by mega5 ( tamura et al . 2011 ) to provide baseml with a tree structure file . the run mode was set at 0 , and the gtr nucleotide substitution model was adopted . clock was set to be 0 , indicating no clock and the rates were thus entirely free to vary from branch to branch . four categories for the discrete - gamma model were specified ( ncatg = 4 ) . the rateancestor was set to be 1 to force the program to infer the ancestral sequence . single - base substitutions were identified by comparing the ancestral sequences with that of sa957 using the mummer package ( kurtz et al . 2004 ) . in detail , the nucmer program of the package was used to generate a chromosomal alignment ; then the delta - filter program was used to filter the redundant alignment and the show - snps program was finally used to report the substitutions . a home - made perl script was written to make sure that only substitutions that were 5 bp apart were considered to be valid single - base substitutions . home - made perl scripts were written to extract the short indels from the alignments of the conserved regions . to determine the direction of the indel , namely whether it is an insertion or deletion , we firstly need to know the ancestral state of the indel . to do this , we defined that if a certain state was the same in more than five out of the eight clonal complexes , then that state was deemed the ancestral state . staphylococcus aureus strain sa957 was obtained from chang gung children s hospital in taiwan . for genome sequencing , a mono - clone of sa957 was picked from lb agar and grown in liquid lb broth at 37 c overnight . the overnight culture was diluted 40-fold with fresh lb broth and incubated at 37 c to reach log phase . the genomic dna was then extracted by qiaamp dna mini kit ( qiagen , valencia , ca , usa ) and subjected to solid mate - pair library construction . the library was introduced to a solid 3 plus sequencer for 2  50 bp sequencing . cp003603.1 ) was used as the reference genome and the socs package ( ondov et al . 2008 ) was used to map the solid reads and to detect single - base substitutions . in detail , the solid reads were firstly filtered by the script filterreads.pl ( within socs package ) with the quality threshold set as q15 . according to the average error rate obtained from a pilot mapping , we decided to keep the first 42 bp of each read and ran the mapping again . the variant mode  short variants  was set for detecting valid sequence - space mismatches indicated by color - space mismatches . the results were then outputted as the coverage at each position and the valid sequence - space mismatches detected from each strand of the chromosome . only reads with unique hits against the genome were included in the downstream analysis so that spurious alignments were avoided . in addition , the positions with coverage smaller than 500 or larger than 15,000 were discarded because with such low or high coverage there was a higher likelihood of misalignment . if a given position had more than one type of substitution meeting the above criteria , only the major substitution was counted . to differentiate mutations from the sequencing errors , we defined firstly that a mutant allele must have at least five supporting reads on both the forward and reverse strands with the aim of reducing the likelihood that a mutant allele came exclusively from duplicate molecules amplified during library construction . then , the bernoulli test was performed for each variant site . in detail , the function binomdist in microsoft excel was used , where the number_trail was the total coverage of the site , the number_success was the number of reads supporting the mutant allele , and the prob_success was the presumed error rate . the chromosome of sa957 was 2,789,538 bp long . according to the cumulative gc skew and at skew , which were calculated using genskewapp.jar ( downloaded from http://genskew.csb.univie.ac.at , last accessed march 13 , 2013 ) , the oric and ter site were presumed to be at 0 and 1,414 kb , respectively . the substitutions in the second half of the chromosome were converted into their reverse complement so that the substitution pattern could be entirely described for the leading strand . the frequency of a certain substitution was defined to be the number of this type of substitutions divided by the total number substitutions . for example , f(c  a ) = number(c  a)/number(all ) . the trend for how each base changed as a percentage during short- and long - term evolution were then predicted based on the frequency of the 12 possible substitution types . in detail , 
 net frequency ( a ) = f(c  a ) + f(g  a ) + f(t  a )  f(a  c )  f(a  t ) f(a  g)net frequency ( t ) = f(a  t ) + f(c  t ) + f(g  t ) f(t  a )  f(t  c )  f(t  g)net frequency ( c ) = f(a  c ) + f(t  c ) + f(g  c )  f(c  a )  f(c  t ) f(c  g)net frequency ( g ) = f(a  g ) + f(t  g ) + f(c  g ) f(g  c ) 
  net frequency ( a ) = f(c  a ) + f(g  a ) + f(t  a )  f(a  c )  f(a  t ) f(a  g ) net frequency ( t ) = f(a  t ) + f(c  t ) + f(g  t ) f(t  a )  f(t  c )  f(t  g ) net frequency ( c ) = f(a  c ) + f(t  c ) + f(g c )  f(c  a )  f(c  t )  f(c  g ) net frequency ( g ) = f(a  g ) + f(t  g ) + f(c  g ) the socs package does not support gapped alignment and therefore we used shrimp2 ( david et al . 2011 ) to create a gapped alignment that would allow the detection of short indel events . was adopted to make sure that only the highest scoring mapping for a given read was outputted . the score to extend a gap along the genome sequence was set at 2 and the score to extend a gap along the read sequence was set at 1 . indels that were located in the first / last 5 bp of the read were not considered to be valid ones and such hits were discarded from the alignment file using a home - made perl script . 2010 ) was used to left - align indels so that the same indel would not be represented as a different variant . the bambino program ( edmonson et al . 2011 ) was used to process the realigned file and detect the indels within it . the following options were set : the minimum number of reads passing all quality filters required at a variant site was set at two ; the required minimum number of unique read names supporting the mutant allele was set at two ; and other options were set at default . only indels that were supported by reads in both the forward and reverse direction were counted . then bernoulli test was performed on each indel site ( p < 0.001 ) . only the sites with cp003603.1 , multi - locus - sequence - type st59 ) together with s. aureus strains fpr3757 ( nc_007793.1 , st8 ) , jkd6159 ( nc_017338.1 , st93 ) , lga251 ( nc_017349.1 , st425 ) , mrsa252 ( nc_002952.2 , st36 ) , mw2 ( nc_003923.1 , st1 ) , n315 ( nc_002745.2 , st5 ) , and rf122 ( nc_007622.1 , st151 ) were used to infer the genomic sequence of an ancestral s. aureus . the derived backbone regions , namely regions that are conserved in all these strains , were extracted from the output backbone file and were further aligned by mafft program ( katoh et al . the length of the concatenated conserved regions totaled 2,393 kb , which is approximately 85% of the whole genome of each strain . the concatenated alignment of the conserved regions was inputted into the baseml program from the paml package ( yang 2007 ) and the marginal reconstruction algorithm was used to infer the genome sequences of the ancestral s. aureus . a neighbor - joining tree was produced by mega5 ( tamura et al . 2011 ) to provide baseml with a tree structure file . the run mode was set at 0 , and the gtr nucleotide substitution model was adopted . clock was set to be 0 , indicating no clock and the rates were thus entirely free to vary from branch to branch . four categories for the discrete - gamma model were specified ( ncatg = 4 ) . the rateancestor was set to be 1 to force the program to infer the ancestral sequence . single - base substitutions were identified by comparing the ancestral sequences with that of sa957 using the mummer package ( kurtz et al . in detail , the nucmer program of the package was used to generate a chromosomal alignment ; then the delta - filter program was used to filter the redundant alignment and the show - snps program was finally used to report the substitutions . a home - made perl script was written to make sure that only substitutions that were 5 bp apart were considered to be valid single - base substitutions . home - made perl scripts were written to extract the short indels from the alignments of the conserved regions . to determine the direction of the indel , namely whether it is an insertion or deletion , we firstly need to know the ancestral state of the indel . to do this , we defined that if a certain state was the same in more than five out of the eight clonal complexes , then that state was deemed the ancestral state .