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 .
 mutation and selection are both thought to impact significantly the nucleotide composition of bacterial genomes . 
 earlier studies have compared closely related strains to obtain mutation patterns based on the hypothesis that these bacterial strains had diverged so recently that selection will not have had enough time to play its role . in this study 
 , we used a solid autosequencer that was based on a dual - base encoding scheme to sequence the genome of staphylococcus aureus with a mapping coverage of over 5,000. by directly counting the variation obtained from these ultradeep sequencing reads , we found that a  g was the predominant single - base substitution and 1 bp deletions were the major small indel . 
 these patterns are completely different from those obtained by comparison of closely related s. aureus strains , where c  t accounted for a larger proportion of mutations and deletions were shown to occur at an almost equal frequency to insertion . 
 these findings suggest that the genomic differences between closely related bacterial strains have already undergone selection and are therefore not representative of spontaneous mutation . 
Materials and Methods
Bacterial Strains and Genome Sequencing
Detection of Single-Base Substitution within SOLiD Data
Detection of Indels within SOLiD Data
Detection of Variation between Ancestral S. 