Gene 389 (2007) 87 – 95 www.elsevier.com/locate/gene

Genome reduction of the aphid endosymbiont Buchnera aphidicola in a recent evolutionary time scale Laura Gómez-Valero a , Francisco J. Silva a , Jean Christophe Simon b , Amparo Latorre a,⁎ a

Institut Cavanilles de Biodiversitat i Biologia Evolutiva and Departament de Genètica, Universitat de València, Apartat 22085, 46071 Valencia, Spain b INRA-Rennes, UMR BiO3P INRA/Agrocampus Rennes, BP 35327, 35653 Le Rheu Cedex, France Received 14 June 2006; received in revised form 27 September 2006; accepted 1 October 2006 Available online 14 October 2006 Received by F.G. Alvarez-Valin

Abstract Genome reduction, a typical feature of symbiotic bacteria, was analyzed in the last stages of evolution of Buchnera aphidicola, the primary aphid endosymbiont, in two neutrally evolving regions: the pseudogene cmk and an intergenic region. These two regions were examined in endosymbionts from several lineages of their aphid host Rhopalosiphum padi, and different species of the same genus, whose divergence times ranged from 0.62 to 19.51 million years. Estimates of nucleotide substitution rates were between 4.3 and 6.7 × 10− 9 substitution/site/year, with G or C nucleotides being substituted around four times more frequently than A or T. Two different types of indel events were detected, of which many were small (1–10 nt) but one was large (about 200 nucleotides).With respect to the large one and considering the proportion and size of the deletions and insertions, the reduction rate was 1.3 × 10− 8 lost nucleotides/site/year. We propose a stepwise scenario for the last stages of evolution in B. aphidicola: together with a very slow and gradual degradation, considerable indels would punctually emerge. The only restriction to large deletion fixation is that the lost fragment does not contain essential genes. © 2006 Elsevier B.V. All rights reserved. Keywords: DNA loss; Indels; Pseudogenes; Symbiosis; Molecular clock

1. Introduction The range of known bacterial genome sizes varies between 9.2 Mb in the soil-borne bacterium Myxococcus xanthus (Stêpkowski and Legocki, 2001) and 0.45 Mb in Buchnera aphidicola primary symbiont of the aphid Cinara cedri (Gil et al., 2002). It is well accepted that small-sized genomes derived from ancestors with larger genomes, as demonstrated for several bacterial endosymbionts of insects, such as Abbreviations: BAp, B. aphidicola from Acyrthosiphon pisum; BSg, B. aphidicola from Schizaphis graminum; BBp, B. aphidicola from Baizongia pistaciae; BRm, B. aphidicola from R. maidis; BRi, B. aphidicola from R. insertum; BRn, B. aphidicola from R. nymphaeae; BRc, B. aphidicola from R. cerasipholiae; IGR, intergenic region; A + T, molar fraction of adenine and thymine in DNA; G + C, molar fraction of guanine and cytosine in DNA; kb, kilobases; dNTP, deoxynucleoside triphosphate; Myr, million years. ⁎ Corresponding author. Tel.: +34 96 3543649; fax: +34 96 3543670. E-mail address: [email protected] (A. Latorre). 0378-1119/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.gene.2006.10.001

B. aphidicola (Shigenoubu et al., 2000; Tamas et al., 2002; van Ham et al., 2003), Blochmannia floridanus (Gil et al., 2003), B. pensilvannicus (Degnan et al., 2005), or Wigglesworthia glossinidia (Akman et al., 2002). Due to their strictly vertical transmission mode, all obligate intracellular bacteria frequently undergo bottlenecks leading to a very low effective population size, as compared to free-living bacteria. The small effective population size and the absence of recombination lead to the accumulation of slightly deleterious mutations in both essential and non-essential genes due to Muller's ratchet (Moran, 1996). The reduced strength or efficacy of selection gives rise to the loss of the expendable genes. Moreover, in endosymbionts the acquisition of genetic material via horizontal transfer from other bacterial species does not take place (Silva et al., 2003). The drop in DNA incorporation, along with the loss of genetic material, leads to a drastic reduction of the endosymbiont genome. A mutational deletion bias has been proposed as the main mechanism of genome shrinkage (Moran

88

L. Gómez-Valero et al. / Gene 389 (2007) 87–95

and Mira, 2001; Gregory, 2003), involving the continuous erosion of material that is not subject to selection (i.e. noncoding DNA; pseudogenes, etc.). This bias has been found both in eukaryotes, like Drosophila or mammals (Ophir and Graur, 1997; Petrov and Hartl, 1999), and in microorganisms, like parasites and symbionts (Andersson and Andersson, 1999, 2001; Moran and Mira, 2001; Wernegreen, 2002, 2005). A model system to study the process of genome reduction is the obligate association between aphids and their maternally transmitted intracellular symbiont B. aphidicola. The association is ancient, and the congruence between the phylogenetic trees of hosts and symbionts indicates a unique infection occurred over 100 million years (Myr) ago, followed by the coevolution of both partners (Moran et al., 1993). It has been estimated that all extant Buchnera diverged from free-living bacteria with genome sizes from 2.0 to 2.5 Mb, containing 1800 to 2500 genes (Moran and Mira, 2001; Silva et al., 2001). At the moment of the analysis, three complete genomes of B. aphidicola were available: B. aphidicola from Acyrthosiphon pisum (BAp) (Shigenobu et al., 2000), B. aphidicola from Schizaphis graminum (BSg) (Tamas et al., 2002), and B. aphidicola from Baizongia pistaciae (BBp) (van Ham et al., 2003). BAp and BSg, whose hosts belong to different tribes of the subfamily Aphidinae (Macrosiphini and Aphidini, respectively) are estimated to have diverged 50 to 70 Myr ago (Clark et al., 1999). On the other hand, the strain BBp, whose host belongs to the subfamily Pemphiginae, probably diverged from BAp and BSg 80 to 150 Myr ago (Moran et al., 1993). The genome size of these three B. aphidicola strains varies only from 616 to 641 Kb. However, the estimation of genome size of B. aphidicola from several aphid subfamilies showed differences of up to 200 Kb (Gil et al., 2002). It has been proposed that DNA loss must have occurred very quickly in the initial stages of the symbiotic process, due to the elimination of several contiguous genes in single large deletion events (Moran and Mira, 2001), as well as the simultaneous disintegration of many pseudogenes (Silva et al., 2001). However, the final steps must have been slow because most of the retained genes would be necessary, most of the repeat elements favouring recombination would have been lost, and an efficient recombination system would not exist (Mira et al., 2001; Frank et al., 2002; Rocha, 2003). Due to the uncertainties in estimated divergence times, the speed of the genome reduction process is difficult to assess. Thus, different analysis using the three B. aphidicola sequenced genomes gave diverse estimations (Mira et al., 2002; Tamas et al., 2002; GómezValero et al., 2004). Recently, Nilsson et al. (2005) have shown that large deletions can occur in a short evolutionary time scale in Salmonella enterica which is defective in mismatch repair, this being the case of most insect endosymbionts. The aim of this work was to characterize the degradation of neutral regions in the latter stages of B. aphidicola evolution, by analyzing nucleotide substitutions and indel rates at two relatively recent evolutionary time scales: within and between host species of the same genus. The system chosen was B. aphidicola from Rhopalosiphum padi, and related species of the same genus. This choice was made because recent molecular

phylogenetic studies provided time divergence estimates between the two main mitochondrial DNA lineages (haplotypes I and II) of R. padi (Martínez-Torres et al., 1996; Simon et al., 1999; Delmotte et al., 2003). Mitochondrial and endosymbiontic genomes should have similar evolutionary histories because they are both maternally inherited, as shown in R. padi (Simon et al., 1996). Hence, time estimates obtained on mtDNA lineages should also apply to associated endosymbiont ones. Through the reconstruction of a phylogeny and the calibration of the divergence, in the present work we have been able to study genome degradation in B. aphidicola belonging to aphid host lineages of known ages. This study covers both from intraspecific and interspecific variation. 2. Material and methods 2.1. Aphid samples and DNA sequences of B. aphidicola B. aphidicola from thirty-seven R. padi strains from different geographical localities distributed worldwide were analyzed (Table 1), 16 of which were of mtDNA haplotype I (hI) and 21 of haplotype II (hII). Discrimination between these two haplotypes is detailed in Delmotte et al. (2003). Furthermore, B. aphidicola from four other aphid species of the genus Rhopalosiphum were studied, these being: Rhopalosiphum maidis (BRm), Rhopalosiphum insertum (BRi), Rhopalosiphum nymphaeae (BRn) and Rhopalosiphum cerasipholiae (BRc) (only for repA2 analysis, see below). These aphid species belong to the subfamily Aphidinae, tribe Aphidini. To calibrate the molecular clock in B. aphidicola we also included sequences of the endosymbiont of two aphid species deposited in databases: S. graminum (BSg) and A. pisum (BAp) belonging to the Aphidini and Macrosiphini tribes, respectively. Two neutrally evolving regions were studied: the pseudogene cmk and the intergenic region between the genes hupA and rpoC (IGR). The pseudogene cmk in the Aphidini lineage was assumed neutral, due to the low G + C content and small size of the cmk pseudogene in BSg (Tamas et al., 2002) and BRp (present work). IGR neutrality could be affected by the possible existence of regulatory elements, but as most of these elements have been lost in the extant B. aphidicola genomes, neutrality was assumed. We also worked with the plasmid gene repA2 to calibrate the molecular clock. This gene was chosen because it had previously been shown to behave as a good phylogenetic marker for the subfamily Aphidinae, with the predicted separation of the B. aphidicola strains from the tribes Macrosiphini and Aphidini (Silva et al., 1998). We obtained the repA2 sequences from BRp (one for each mtDNA haplotype), BRi and BRm. Sequences were obtained from the GenBank from BRc and BSg, which belong to the same tribe but to a different genus, and from BAp from the tribe Macrosiphini. 2.2. DNA isolation and sequencing Total DNA was extracted from aphids following the method by Latorre et al. (1986).

L. Gómez-Valero et al. / Gene 389 (2007) 87–95

89

Table 1 Name, code, geographical origin, mitochondrial DNA haplotype (mtDNA), and accession number of the species used Species

Code

Rhopalosiphum padi

(strains)

Outgroups Rhopalosiphum insertum Rhopalosiphum maidis Rhopalosiphum cerasipholiae Rhopalosiphum nymphaeae Schizaphis graminum Acyrthosiphon pisum

Geographical origin

mtDNA

BRp1 BRp2 BRp3 BRp4 BRp5 BRp6 BRp7 BRp8 BRp9 BRp10 BRp11 BRp12 BRp13 BRp14 BRp15 BRp16 BRp17 BRp18 BRp19 BRp20 BRp21 BRp22 BRp23 BRp24 BRp25 BRp26 BRp27 BRp28 BRp29 BRp30 BRp31 BRp32 BRp33 BRp34 BRp35 BRp36 BRp37

Rennes, France Rennes, France Rennes, France Rennes, France Perth, Australia Colmar, France Bristol, Great Britain Bristol, Great Britain Dundee, Scotland Kendenup, Australia Nairobi, Kenya Montpellier, France Rennes, France Opava, The CzechRepublic Alep, Syria Tunis, Tunisia St. Amand, France Rennes, France Rennes, France Guelph, Canada Windsor, Canada Colmar, France Colmar, France Colmar, France Göttingen, Germany CYMMIT, Mexico Rennes, France Rennes, France Rennes, France Rennes, France Rennes, France Rennes, France Rennes, France Rennes, France Rennes, France Belgrade, Serbia Kerguelen Islands (Indian Ocean)

BRi BRm BRc BRn BSg BAp

Rennes, France Rennes, France AlgonquinPark, Canada Rennes, France Database Database

Degenerated primers were designed based on the genome sequence of flanking genes of the two regions of BAp and BSg. For cmk the primers were: rpsA1R (5′CTGAGAAAAGGTTCGGTATAGTC3′) and aroA1F (5′TTATAACGATCATCGCATGGC3′), and for the IGR the primers were: hupABu (5′ DTTAATTAATTGAGTTTTATTCAT3′) and rpoCBu (5′ ACWGGATATGCATATCAYAAARAAGG3′). Each PCR product was first sequenced with these primers in BRp. The first sequence obtained enabled us to design the following species-specific PCR primers: hupABuRp (5′TTAATTGAGTTTTATTCAT3′) and rpoBuRp (5′TATGCATATCAYAAARAACG3′) for the IGR, and cmkRp1F (5′GGCATGAGTGTTGGTATACGTAATC3′) and cmkRp1R (5′TCTCCTGGATTAATAGAAGAAGACG3′) for cmk. Two more primers were

I I I I I I I I I I I I I I I I II II II II II II II II II II II II II II II II II II II II II

Acc. No cmk

IGR

repA2

DQ105477 DQ105468 DQ105470 DQ105475 DQ105466 DQ105465 DQ105476 DQ105469 DQ105478 DQ105471 DQ105473 DQ105480 DQ105467 DQ105479 DQ105474 DQ105472 DQ105495 DQ105490 DQ105492 DQ105498 DQ105494 DQ105485 DQ105497 DQ105496 DQ105489 DQ105483 DQ105481 DQ105484 DQ105501 DQ105486 DQ105487 DQ105482 DQ105488 DQ105499 DQ105493 DQ105491 DQ105500

– – – – – DQ105510 DQ105511 DQ105513 – DQ105512 DQ105515 – – DQ105514 – – – – – DQ105509 – – DQ105508 – – DQ105505 DQ105506 DQ105507 – – – – – – – – –

– – – – DQ105461 – – – – – – – – – – – – – – – – – DQ105460 – – – – – DQ105459 – – – – – – – –

DQ105503 DQ105504 – DQ105502 – –

DQ105516 DQ105517 – DQ105518 – –

DQ105462 DQ105464 DQ105463 AJ006876 AJ006878

designed to obtain the sequence of cmk in BRn: cmkRnymF (5′ CCAAATTGTATTTCTAAAACTTTTCCATC) and cmkRnymR (5′TCAATAGTTGAGCAAAAGATTCATTC3′) and one more upstream to the gene serC, to obtain the sequence of cmk in BRm: serCRp (5′GGSCARTTTGCHGCTGTY CCWATR3′). In the case of the gene repA2 the sequences were obtained from a long PCR product amplified with the primers leuA.lo3 (5′ARACTWGCTTGWARWGCTTGTTCWCCATC3′) and ORF1up2 (5′GTWATGGTWATGTTTTCWGGWTA3′). PCR mixtures consisted of 1.5 U of Taq DNA polymerase (Promega), 200 μM of each dNTP, 300 nM of each primer and 10 ng of DNA template in a final volume of 50 μl. In the case of the PCR for cmk in BRn we carried out a PCR reaction, using

90

L. Gómez-Valero et al. / Gene 389 (2007) 87–95

Table 2 Values of indels, cost and ratio ts/tv (transitions/transversions) depending on different combination of GEP penalties, for a constant GOP penalty GOP

GEP

Indel cost according to the indel size

BRp–BRn

1 nt

2 nt

3 nt

10 nt

10*1 nt

ts

tv

ratio

cmk 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5

0.03 0.5 4 6 8 10 12 14

0.53 1 4.5 6.5 8.5 10.5 12.5 14.5

0.56 1.5 8.5 12.5 16.5 20.5 24.5 28.5

0.59 2 12.5 18.5 24.5 30.5 36.5 42.5

0.81 5.5 40.5 60.5 80.5 101 121 141

5.3 10 45 65 85 105 125 145

27 36 39 48 51 54 54 55

0 0 11 17 23 27 29 32

– – 3.5 2.8 2.2 2 1.9 1.7

IGR 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5

0.03 0.5 4 6 8 10 12 14

0.53 1 4.5 6.5 8.5 10.5 12.5 14.5

0.56 1.5 8.5 12.5 16.5 20.5 24.5 28.5

0.59 2 12.5 18.5 24.5 30.5 36.5 42.5

0.81 5.5 40.5 60.5 80.5 101 121 141

5.3 10 45 65 85 105 125 145

19 15 36 39 45 46 26 42

0 0 12 16 17 23 14 30

– – 3 2.4 2.6 2 1.8 1.4

*Ten indels of 1 nucleotide each one.

Ex Taq polymerase (Takara) and 500 μM of each dNTP. PCR amplification conditions were 94 °C for 2 min followed by 35 cycles of 94 °C for 30 s, 54 °C for 30 s, and 72 °C for 35 s, or 68 °C for 2 min, in the case of BRn. For those primers not yielding a PCR product when the above-cited temperature was used, the annealing temperature was gradually lowered until the desired product was obtained. Sequencing was conducted with an ABI3700 automated sequencer using the ABI PRISM Dye Terminator Cycle Sequencing Ready Reaction Kit (PerkinElmer) following the manufacturer's specifications. 2.3. Computer analysis of DNA sequences Multiple alignments were carried out with Clustal X (Thompson et al., 1997); however, as many indels were found in these neutral regions, they were difficult to align. Thus, we carried out a preliminary study to examine the most suitable values for the alignment parameters. Previous studies comparing sequences between BSg and BRp determined that at these divergence levels (less than 50 Myr), the transitional differences were more abundant than the transversional ones (Silva et al., 1998). We prepared a matrix for the alignments in which the transitions had an intermediate value (−1) between transversions (−9) and matches (10). Regarding the values of the GOP (the cost of inserting a new gap into a sequence) and the GEP (the cost of extending an existing gap), Morrison and Ellis (1997) proposed the strategy of logarithmically varying them, with a total of 64 possible combinations. This involves varying the GOP from 0.5 to 64 times the cost of a substitution (log2 = −1, 0, 1, 2, 3, 4, 5, 6) and the GEP from 0.031 to 8 times the cost of a substitution (log2 GEP = −5, −4, −3, −2, −1, 0, 1, 2, 3). Sanchis et al. (2001) showed that with only six of these combinations the strength of the phylogenetic signal could be detected. These six combinations

should include three GOP values (representing strong, medium and weak gap penalties) versus two GEP values. With these six combinations we always obtained the same topology, but not the same alignment (data not shown). Moreover, we realized that when the GOP had a large or medium value, the alignment favoured the formation of one single large deletion instead of several small ones. However, it has been postulated that the indels in obligate intracellular bacteria are most frequently small in size (Andersson and Andersson, 1999; Tamas et al., 2002; van Ham et al., 2003). To avoid this problem, we decided to fix the GOP at a low value (0.05) and combine this value with the eight different values of GEP from the least (0.03) to the most restrictive (14) (Table 2). This meant that the gap penalty was mainly associated to the number of nucleotide sites in the indel rather than to the number of indels. For each combination, we estimated the ratio of transitional differences to transversional ones for BRp and BRn. As an indication of the expected values, we estimated this ratio in BRp and BSg (a more divergent species pair than the previous pair) in the third codon positions of the genes available in the databases. This value was 1.6. Thus, based on many evolutionary studies showing that this ratio decreases as the divergence time increases (Li, 1997), we expected higher values for the pair BRp–BRn. As a result of all the analysis carried out, the GEP values ranging from 6 to 12 were finally selected. 2.4. Phylogenetic reconstruction Phylogenetic reconstruction was done using the neighborjoining method implemented in MEGA3 (Kumar et al., 2004), and maximum-likelihood (ML) and maximum-parsimony (MP) methods implemented in PAUP⁎4.0b10 (Swofford, 2002). Modeltest, version 1.05 (Posada and Crandall, 1998) was used in combination with PAUP to find the best-fit model of DNA substitution for the entire data set. In the case of repA2 only the most conserved second codon position was used, because the BAp sequence which hosts belongs to a different clade (Macrosiphini) increases the divergence time. Heuristic searches were performed in all cases with random sequence addition (n = 500 and 1000 in ML and MP analyses, respectively) and TBR branch swapping. Bootstrap was conducted with all methods performing 1000 replicates in all cases. Once the topology was established, the nucleotide substitutions with a single parsimonious solution were located on the tree. Likewise, insertions can be distinguished from deletions by Table 3 Possible scenarios (from A to N) for the indel events: (X), presence of nucleotide/s (-), absence of nucleotide/s

BRp BRi BRm BRn

A

B

C

D

E

F

G

H

I

J

K

L

M

N

X X X D

X X X D

X X X D

X X X u

X X D

X X u

X X u

X X u

X X u

X X I

X u

X I

X I

X I

Each case is solved with the single most parsimonious solution (D, deletion; I, insertion; u, unsolved).

L. Gómez-Valero et al. / Gene 389 (2007) 87–95

91

Fig. 1. Phylogenetic trees of a) cmk and b) IGR sequences. Insertions (I), deletions (D), transitions (s) and transversions (v) are located above the corresponding branch. Bootstrap values are indicated next to the corresponding node.

superimposing the observed mutations onto the reconstructed phylogeny and assuming that the mutational events, which correspond to the smallest number of changes, represent the actual events (Andersson and Andersson, 1999). With this criterion we were able to resolve 8 of the 14 possible scenarios for the indel mutations (see Table 3). 2.5. Molecular clock test Molecular clock was tested by means of a likelihood ratio test using PAUP*4.0b10 (Swofford, 2002). Twice the difference between both Lns is χ2 distributed with n - 2 df, where n is the number of sequences. Tajima's relative test (Tajima, 1993) implemented in MEGA3 (Kumar et al., 2004) was also applied using selected representative sequences of the clusters obtained in the phylogenetic analyses. The divergence times were finally dated with a Bayesian statistical inference by using the Markov chain Monte Carlo integration (Drummond et al., 2002) and the Kingman coalescent model using the program Beast v1.0.3 (Drummond and Rambaut, 2003; http://evolve.zoo.ox.ac.uk/beast). 3. Results 3.1. Phylogenetic reconstruction and location of indels and nucleotide substitutions The cmk pseudogene was sequenced in the 37 R. padi strains, in BRi, BRm and BRn. In the case of R. padi all the

strains with mtDNA hI possessed the same sequence, whereas three variants were found in the strains with mtDNA hII. The size was smaller in BRp and BRi (around 600 bp) than in BRm and BRn (around 800 bp). Sequence alignments revealed a large indel towards the end of the sequence caused by a deletion that occurred in the ancestor of BRi and BRp. We excluded this gap from the initial analysis in order to estimate the effect of the small-size indels alone. After removing the large indel the size of cmk ranged from 548 to 552 bp. To reconstruct the phylogenetic relationships between Buchnera lineages from distinct Rhopalosiphum taxa, we used two alignments, the most and the least restrictive taking the GEP values 6 and 12 (Table 2). We only took representative sequences of R. padi that showed polymorphic sites in the previous analysis. The evolutionary model obtained with Modeltest was the same for the two alignments: K81 + uf

Table 4 Number and type of indel events classified by size Indels

Deletions

Insertions

Size range (nt)

Size range (nt)

Size range (nt)

1

2–5 6–10 I. nt 1

cmk 16 0 IGR 9 4 Total 25 4

0 3 3

16 37 53

2–5 6–10 I. nt 1 2–5 6–10 I. nt

10 0 6 2 16 2

0 2 2

10 25 35

6 0 3 2 9 2

0 1 1

The total number of involved nucleotides (I. nt) in each case is shown.

6 12 18

92

L. Gómez-Valero et al. / Gene 389 (2007) 87–95

Table 5 Percentage of each type of nucleotide substitution Nt change

cmk

IGR

A → G=T → C A → T=T → A A → C=T → G C → A=G → T C → G=G → C G → A=C → T GC → AT AT → GC

17.74 5.28 1.42 10.82 2.06 62.67 73.5 19.6

35.52 4.44 3.07 8.99 0 47.97 56.96 38.59

(unequal frequencies). The phylogenetic reconstructions gave similar topologies, whatever the method used: NJ, MP and ML (Fig. 1a). The closest related species to BRp was BRi, followed by BRm and finally by BRn. The IGR was amplified in B. aphidicola from the same species but in a restricted number of R. padi strains considering the low variability, previously found in cmk (Table 1). The IGR size ranged from 463 to 499 bp. As in the case of cmk, we used the most and the least restrictive alignments. In this case, the evolutionary model obtained with Modeltest was different for each alignment, TVM + I for the alignment with a GEP of 6, and K81 + uf for the alignment with a GEP of 12. As for cmk we only took one representative of the samples showing some sequence differences (BRp 11, 26 and 29 in the phylogenetic tree). The topology of the tree was the same with the three methods used (Fig. 1b). To differentiate insertion from deletion events, a descriptive table was drawn up for each region with the common indels to the four alignments indicated in Material and methods (data not shown and available upon request). As both the position and the size of each indel varied slightly between the different alignments, we took the common number of nucleotides for all the alignments as the final size (Table 4). It is remarkable that one third of the indels (34%) was part of a mononucleotide tail. Deletions prevailed over insertions in the two regions. In the case of cmk the size of the deletion event was always 1 nt, whereas deletions greater than 1 nt were found in the IGR, but their sizes never exceeded 10 nt. The phylogenetic reconstruction described above allowed us to locate the indels and nucleotide substitution events on the tree with a single parsimonious solution (Fig. 1). We only considered those events shared by the four alignments, be they

either substitutions or indels. We did not locate indel events in the BRn branch because this strain was used as an outgroup. Regarding nucleotide substitutions, we studied the possible existence of some bias in the changes. First, we corrected the detected substitutions according to the initial composition of the sequences, estimating the number of directional nucleotide changes in each region and dividing each value by the frequency of the original nucleotide. Next, we obtained the relative frequency of each class of nucleotide substitution ( fij) (from Tamura and Nei, 1993). Finally, we grouped the changes that could not be differentiated because the mutated and complementary strands cannot be distinguished ( fAG = fTC; fAT = fTA; fAC = fTG; fCA = fGT; fGC = fCG; fGA = fCT). As shown in Table 5, transitional substitutions overcome transversions in both regions. In addition, G or C nucleotide substitutions to A or T were more frequently detected than the opposite ones in both regions. However, the bias was higher in cmk (73.5% changes from G or C to A or T versus 19.6% from A or T to G or C) than in IGR (59.96% versus 38.59%, respectively). 3.2. Calibrating the molecular clock A repA2 phylogeny was reconstructed using distance, likelihood and parsimony methods, obtaining the same topology (Fig. 2). The molecular clock hypothesis could not be rejected ( p b 0.05), and the relative rate tests carried out for all combinations between all species were nonsignificant, with only one exception corresponding to the triad BRp, BSg and BAp (p = 0.027). The clock was calibrated using the divergence time between the tribes Macrosiphini and Aphidini estimated according to fossil record as 50 to 70 Myr (Clark et al., 1999). In order to date the time at which splits occurred, a Bayesian statistical inference was carried out using the BEAST program and assuming constant rates. The advantage of this program is that the split events can be dated probabilistically within a range, which is a more realistic approximation to the actual date. Results showed that evolutionary events under investigation occurred from around 600,000 years to less than 20 Myr (Table 6). 3.3. Rates of indel and nucleotide substitution We estimated the rates of indel and nucleotide substitution using the mean time of each split event. When the number of

Fig. 2. Topology of repA2 genes inferred from the neighbor-joining method. Bootstrap values are indicated next to the corresponding node. The nodes are named from A to F. See Table 6 for the dating of the split events.

L. Gómez-Valero et al. / Gene 389 (2007) 87–95 Table 6 Dates of split events corresponding to the nodes from Fig. 2

93

4. Discussion

NODE

Mean time (Myr)

hpdLower (Myr)

hpdUpper (Myr)

ESS

A B C D E F

59.42 19.51 14.88 11.21 4.77 0.62

50.01 12.41 8.57 6.39 2.27 0.07

68.08 27.64 21.54 16.64 7.41 1.39

8085.454 5456.37 5427.05 6057.6 6411.72 8623.91

For each node the mean time and the lower and the upper time limits are given. The ESS is the effective sample sizes.

gaps was approximately proportional to the length (evolutionary time) of each branch of the tree, the evolutionary rate of insertion and deletion was assumed to be more or less constant during the last stages of B. aphidicola evolution in neutral regions. This result was obtained in all cases except for the branch leading to BRm in the cmk phylogeny, where the number of gaps was lower than expected according to the divergence time. One possibility for this result is that, at least during some period of time, cmk was evolving as a gene in this lineage. This hypothesis was supported by a higher GC content of the BRm cmk region than the other B. aphidicola strains (21% versus 16–17%) and by a higher degree of similarity with functional cmk genes. These results would indicate a convergent gene loss of cmk in B. aphidicola lineages, similar to that demonstrated in other genes losses (Gómez-Valero et al., 2004). Thus, a cmk gene probably existed in B. aphidicola at the moment of the divergence of R. maidis, and was inactivated in the early stage of lineage leading to R. padi/R. insertum but after the split of the R. maidis lineage. Therefore, we decided to remove this sequence on calculating the indel rate in cmk. Indel rate was calculated using the Saitou method (Saitou and Ueda, 1994). In this method, the total number of gaps is divided by the total time (Myr) according to the addition of the branch lengths, and by the average length of the compared nucleotide sequences. The number of events was 1.0 × 10− 9 and 1.2 × 10− 9 per site and year in IGR and cmk, respectively. When the total number of inserted or deleted nucleotides was considered, a bias was detected towards the loss of DNA at a rate of 7.7 × 10− 10 and 2.7 × 10− 10 lost nucleotides/site/year in IGR and cmk, respectively. This estimation did not take into account the large indel (about 200 nt), which took place in cmk between the split of BRm and the ancestor of BRp and BRi (4.77– 14.88 Myr). Were it to be included in the analysis, the deletion rate in cmk would be much higher (1.3 × 10− 8 lost nucleotides/ site/year). To calculate the substitution rate, we estimated the number of substitutions per site in each branch using MEGA3 (Kumar et al., 2004). Kimura 2P was the chosen model with the pairwise deletion option, and the distance tree reconstruction was done using the neighbor-joining method. Once these values were obtained, they were divided by the total time considered, similar to what was done to calculate the indel rates. Finally, we obtained an average value from the four alignments for each region. The overall results were 4.3 × 10− 9 and 6.7 × 10− 9 substitutions/site/year for IGR and cmk, respectively.

The rate of genome shrinkage is difficult to estimate because divergence rates are not available for many systems. Because symbionts have co-evolved with their aphid hosts, the divergence dates of the aphids are parallel to those of their endosymbionts, facilitating the use of calibrated molecular clocks. Analysing genome reduction in B. aphidicola led to several conclusions. First, based on a divergence date of 50 to 100 Myr for BAp and BSg, and the observed loss of 14 genes in the two lineages, one gene loss for every 5 to 10 Myr was estimated (Tamas et al., 2002). Second, the comparison of genome sizes of the three B. aphidicola sequenced genomes (BAp, BSg and BBp), led to the proposal that the rate at which sequences in the B. aphidicola genome had more recently eroded was as low as 1 nt per 10,000 years (Mira et al., 2002). Finally, Gómez-Valero et al. (2004) proposed that DNA loss was taking place at a higher rate (between 1 × 10− 8 and 5 × 10− 8 nucleotides/site/year) on analyzing the fate of the DNA from those genes that were lost throughout the evolution of the genome in the three B. aphidicola strains (so becoming part of the non-functional DNA). Genome reduction is an ongoing process in B. aphidicola as shown by the extremely small genome size of B. aphidicola from C. cedri (416 Kb), which is the smallest microbial genome sequenced so far (Pérez-Brocal et al., 2006). The comparison of BCc with the previously sequenced genomes revealed that more than 120 genes have been completely lost, and only three pseudogenes are still present. Thus, in order to gain in-depth knowledge of the genome reduction process we need to study recent evolutionary times where we can detect and measure mutational events. Insertion and deletion rates are typically calculated according to the nucleotide substitutions because the times of divergence between species are unknown (Ophir and Graur, 1997; Petrov and Hartl, 1999; Bensasson et al., 2001). In the present work, and due to the co-evolution of B. aphidicola with its aphid host, we have used the divergence time of the tribes Aphidini and Macrosiphini (Clark et al., 1999), as a way to calibrate the molecular clock. Thus, we have been able to estimate the age of the lineages of the genus Rhopalosiphum placed on the tree. It is worth mentioning, that the estimates for time divergences between the two BRp lineages (node F in Fig. 2) were similar to those obtained in previous studies for R. padi mtDNA hI and hII, assuming the rate estimated for Drosophila mtDNA (Martínez-Torres et al., 1996, 0.46–1.40 Myr). In both regions, transitional substitutions prevailed over transversional ones. However, the rate of nucleotide substitution was slightly different with a higher rate in the case of cmk than in IGR (6.65 × 10− 9 and 4.3 × 10− 9 substitution/site/year, respectively), and likewise a stronger bias of GC to AT versus AT to GC substitutions in the pseudogene (73.5 versus 19.6) than in the IGR (56.96 versus 38.59). The weak bias in IGR could be explained by the presence of two GC-rich functional invariable motifs. We searched for G or C nucleotide sites with a putative function in the IGR and found 3 C nucleotides corresponding to the complementary Shine–Dalgarno sequence

94

L. Gómez-Valero et al. / Gene 389 (2007) 87–95

of the hupA gene and 10 G or C nucleotides located in a putative transcription terminator hairpin of the rpoC gene. These invariant nucleotides represent 20% of the total G + C for this intergenic region, and are responsible for the bias of the IGR. In the case of cmk, the higher relative substitution frequency of GC to AT versus AT to GC pairs is very similar to the data obtained for Rickettsia pseudogenes (70.4 versus 22.9) (Andersson and Andersson, 1999). The most probable explanation for this bias towards A + T in the B. aphidicola genome is the elimination of genes encoding DNA repair enzymes (Wernegreen, 2005), but see Rocha and Danchin (2002) for a more general discussion of the A + T bias in parasitic and symbiont sequences. In particular, the incorporation of uracil into DNA, due either to replication error or to C → U deamination, leads to mutational pressure towards A + T if it is not prevented or corrected (Glass et al., 2000). Regarding the indels, our results showed the predominance of deletions over insertions in accordance with the previously shown existence of a mutational deletion bias (Gregory, 2003). We also showed a high prevalence of single nucleotide gaps, which seems to be a common phenomenon both in eukaryotic and prokaryotic nucleotide evolutions (Mira et al., 2001; Gregory, 2003). The average size of the deletion and insertion events is an important factor affecting the reduction rate. This has been demonstrated in several organisms displaying different rates as a result of the difference in the size of the events rather than the difference in their number. For instance, this is the case of the different rates of DNA loss between grasshoppers, Podisma and Italopodisma (Bensasson, 2001). The big indel found in cmk indicates that deletions of considerable size may still occur during the latter period of B. aphidicola evolution. The rate obtained in the present work was 2.7 and 7.7 × 10− 10 lost nucleotides/site/year in cmk and IGR, respectively. Such a rate is not fast enough to explain the DNA lost by B. aphidicola during its adaptation to the intracellular life since the last common ancestor. In fact, this rate is lower than that of 2.9 × 10− 8 lost nucleotides/site/year estimated previously (Gómez-Valero et al., 2004). The difference in the rate of nucleotide loss is presumably a result of the periods covered by the two studies, which was 100 MY in the previous one, and less than 20 MY in the present work. This would indicate that the loss rate was slower in recent periods of evolution, which would be in agreement with the model of gradual genome disintegration postulated previously for the B. aphidicola evolution (Silva et al., 2001). However, if we had included the large indel detected in cmk, the deletion rate estimation would have been two orders of magnitude higher (1.3 × 10− 8 lost nucleotides/site/ year). Thus, the gradual process of genome reduction can be sporadically altered by more drastic losses. Either slippage replication or unequal crossing over recombination provoked by flanking repeats could be responsible for the small indels detected in both regions. However, the replication slippage does not produce deletions on the order of several hundred base pairs (Gregory, 2003). Thus, the large indel detected in cmk probably originated through unequal crossing over that also involves short flanking repeats (Bzymek and Lovett, 2001). In fact, the necessary components for large

indels to arise (greater than 50 nucleotides) are present in the extant genomes of B. aphidicola: short stretches of identity–up to ten bases–closely spaced. Moreover, the completed sequences of B. aphidicola genomes showed that the lost genes in BSg or BAp after the divergence of both lineages had significantly more repeats of this type than the conserved ones (Rocha, 2003). Furthermore, these repeats will arise more frequently in genomes with a high A + T content, as occurs in reduced genomes like B. aphidicola. Finally, as B. aphidicola lacks a functional recA gene, deletion formation must occur via RecA-independent mechanisms. Nilsson et al. (2005) demonstrated that the absence of recA does not reduce the deletion rate. Therefore, the only restriction to large deletion fixation is that the lost fragment does not contain essential genes. The essential information is completely dependent on the environment, which in the case of B. aphidicola depends mainly on the host diet. Thus, changes in the lifestyle of the host will change the requirements of B. aphidicola. In summary, we propose a stepwise scenario for the last stages of evolution of Buchnera, and probably other intracellular bacteria. A scenario in which, together with a gradual and very slow degradation, considerable deletions would also emerge punctually producing faster changes in genome size. The reductive process of endosymbiotic genomes might, in theory, lead to the loss of all genes except those that are essential to maintaining host-bacterial interaction, otherwise it might eventually lead these bacteria towards extinction. Acknowledgments This work was supported by grant BMC2003–00305 from the Ministerio de Ciencia y Tecnología (Spain) and grant Grupos03/204 from Generalitat Valenciana (Spain). L G-V was funded by a pre-doctoral fellowship from Generalitat Valenciana (Spain). We would like to thank B. Wrobel for support with the program BEAST v1.0.3; Dr. J.M. Michelena and Dr. N. Pérez-Hidalgo for providing some species used in this work, Dr. D. Posada for his help with the analysis of the molecular clock and Dr. Moya for his valuable comments and suggestions. We also acknowledge the SCSIE from the Universitat de València for sequencing support. References Akman, L., et al., 2002. Genome sequence of the endocellular obligate symbiont of tsetse flies, Wigglesworthia glossinidia. Nat. Genet. 32, 402–407. Andersson, J.O., Andersson, S.G.E., 1999. Genome degradation is an ongoing process in Rickettsia. Mol. Biol. Evol. 16, 1178–1191. Andersson, J.O., Andersson, S.G.E., 2001. Pseudogenes, junk DNA, and the dynamics of Rickettsia genomes. Mol. Biol. Evol. 18, 829–839. Bensasson, D., Petrov, D.A., Zhang, D.X., Hartl, D.L., Hewitt, G.M., 2001. Genomic gigantism: DNA loss is slow in mountain grasshoppers. Mol. Biol. Evol. 18, 246–253. Bzymek, M., Lovett, S.T., 2001. Evidence for two mechanisms of palindromestimulated deletion in Escherichia coli: single-strand annealing and replication slipped mispairing. Genetics 158, 527–540. Clark, M.A., Moran, N.A., Baumann, P., 1999. Sequence evolution in bacterial endosymbionts having extreme base compositions. Mol. Biol. Evol. 16, 1586–1598.

L. Gómez-Valero et al. / Gene 389 (2007) 87–95 Degnan, P.H., Lazarus, A.B., Wernegreen, J.J., 2005. Genome sequence of Blochmannia pennsylvanicus indicates parallel evolutionary trends among bacterial mutualists of insects. Genome Res. 15, 1023–1033. Delmotte, F., et al., 2003. Phylogenetic evidence for hybrid origins of asexual lineage in an aphid species. Evolution 57, 1291–1303. Drummond, A.J., Rambaut, A., 2003. BEAST v.1.0. Available from http:// evolve.zoo.ox.ac.uk/beast/. Drummond, A.J., Nicholls, G.K., Rodrigo, A.G., Solomon, W., 2002. Estimating mutation parameters, population history and genealogy simultaneously from temporally spaced sequence data. Genetics 161, 1307–1320. Frank, A.C., Amiri, H., Andersson, S.G.E., 2002. Genome deterioration: loss of repeated sequences and accumulation of junk DNA. Genetica 115, 1–12. Gil, R., Sabater-Munoz, B., Latorre, A., Silva, F.J., Moya, A., 2002. Extreme genome reduction in Buchnera spp.: toward the minimal genome needed for symbiotic life. Proc. Natl. Acad. Sci. U. S. A. 99, 4454–4458. Gil, R., et al., 2003. The genome sequence of Blochmannia floridanus: comparative analysis of reduced genomes. Proc. Natl. Acad. Sci. U. S. A. 100, 9388–9393. Glass, J.I., Lefkowitz, E.J., Glass, J.S., Heiner, C.R., Chen, E.Y., Cassell, G.H., 2000. The complete sequence of the mucosal pathogen Ureaplasma urealyticum. Nature 407, 757–762. Gómez-Valero, L., Latorre, A., Silva, F.J., 2004. The evolutionary fate of nonfunctional DNA in the bacterial endosymbiont Buchnera aphidicola. Mol. Biol. Evol. 21, 2172–2181. Gregory, T.R., 2003. Is small indel bias a determinant of genome size? Trends Genet. 19, 485–488. Kumar, S., Tamura, K., Nei, M., 2004. MEGA3: integrated software for molecular evolutionary genetics analysis and sequence alignment. Brief. Bioinform. 5, 150–163. Latorre, A., Moya, A., Ayala, F.J., 1986. Evolution of mitochondrial DNA in Drosophila subobscura. Proc. Natl. Acad. Sci. U. S. A. 83, 8649–8653. Li, W.-H., 1997. Molecular Evolution. Sinauer Associates, Massachusetts. Martinez-Torres, D., Simon, J.C., Fereres, A., Moya, A., 1996. Genetic variation in natural populations of the aphid Rhopalosiphum padi as revealed by maternally inherited markers. Mol. Ecol. 5, 659–669. Mira, A., Ochman, H., Moran, N.A., 2001. Deletional bias and the evolution of bacterial genomes. Trends Genet. 17, 589–596. Mira, A., Klasson, L., Andersson, S.G.E., 2002. Microbial genome evolution: sources of variability. Curr. Opin. Microbiol. 5, 506–512. Moran, N.A., 1996. Accelerated evolution and Muller's rachet in endosymbiotic bacteria. Proc. Natl. Acad. Sci. U. S. A. 93, 2873–2878. Moran, N.A., Mira, A., 2001. The process of genome shrinkage in the obligate symbiont Buchnera aphidicola. Genom. Biol. 2 RESEARCH0054. Moran, N.A., Munson, M.A., Baumann, P., Ishikawa, H., 1993. A molecular clock in endosymbiotic bacteria is calibrated using the insect hosts. Proc. R. Soc. Lond., B Biol. Sci. 253, 167–171. Morrison, D.A., Ellis, J.T., 1997. Effects of nucleotide sequence alignment on phylogeny estimation: a case study of 18S rDNAs of Apicomplexa. Mol. Biol. Evol. 14, 428–441. Nilsson, A.I., Koskiniemi, S., Eriksson, S., Kugelberd, E., Hinton, J.C.D., Andersson, D.I., 2005. Bacterial genome size reduction by experimental evolution. Proc. Natl. Acad. Sci. U. S. A. 102, 12112–12116. Ophir, R., Graur, D., 1997. Patterns and rates of indel evolution in processed pseudogenes from humans and murids. Gene 205, 191–202. Pérez-Brocal, V., et al., 2006. A small microbial genome: the end of a long symbiotic relationship? Science 314, 312–313.

95

Petrov, D.A., Hartl, D.L., 1999. Patterns of nucleotide substitution in Drosophila and mammalian genomes. Proc. Natl. Acad. Sci. U. S. A. 96, 1475–1479. Posada, D., Crandall, K.A., 1998. MODELTEST: testing the model of DNA substitution. Bioinformatics 14, 817–818. Rocha, E.P.C., 2003. An appraisal of the potential for illegitimate recombination in bacterial genomes and its consequences: from duplications to genome reduction. Genome Res. 13, 1123–1132. Rocha, E., Danchin, A., 2002. Base composition bias might result from competition for metabolic resources. Trends Genet. 18, 291–294. Saitou, N., Ueda, S., 1994. Evolutionary rates of insertion and deletion in noncoding nucleotide-sequences of primates. Mol. Biol. Evol. 11, 504–512. Sanchis, A., Michelena, J.M., Latorre, A., Quicke, D.L.J., Gardenfors, U., Belshaw, R., 2001. The phylogenetic analysis of variable-length sequence data: elongation factor-1 alpha introns in European populations of the parasitoid wasp genus Pauesia (Hymenoptera: Braconidae: Aphidiinae). Mol. Biol. Evol. 18, 1117–1131. Shigenobu, S., Watanabe, H., Hattori, M., Sakaki, Y., Ishikawa, H., 2000. Genome sequence of the endocellular bacterial symbiont of aphids Buchnera sp. APS. Nature 407, 81–86. Silva, F.J., van Ham, R.C., Sabater, B., Latorre, A., 1998. Structure and evolution of the leucine plasmids carried by the endosymbiont (Buchnera aphidicola) from aphids of the family Aphididae. FEMS Microbiol. Lett. 168, 43–49. Silva, F.J., Latorre, A., Moya, A., 2001. Genome size reduction through multiple events of gene disintegration in Buchnera APS. Trends Genet. 17, 615–618. Silva, F.J., Latorre, A., Moya, A., 2003. Why are the genomes of endosymbiotic bacteria so stable? Trends Genet. 19, 176–180. Simon, J.C., Martinez-Torres, D., Latorre, A., Moya, A., Hebert, P.D.N., 1996. Molecular characterization of cyclic and obligate parthenogens in the aphid Rhopalosiphum padi (L). Proc. R. Soc. Lond., B Biol. Sci. 263, 481–486. Simon, J.C., Leterme, N., Latorre, A., 1999. Molecular markers linked to breeding system differences in segregating and natural populations of the cereal aphid Rhopalosiphum padi L. Mol. Ecol. 8, 965–973. Stêpkowski, T., Legocki, A.B., 2001. Reduction of bacterial genome size and expansion resulting from obligate intracellular lifestyle and adaptation to soil habitat. Acta Biochim. Pol. 48, 367–381. Swofford, D.L., 2002. PAUP* Phylogenetic Analysis Using Parsimony (And Other Methods). Version 4. MA Sinauer Associates, Sunderland. Tajima, F., 1993. Simple methods for testing the molecular evolutionary clock hypothesis. Genetics 135, 599–607. Tamas, I., et al., 2002. 50 million years of genomic stasis in endosymbiotic bacteria. Science 296, 2376–2379. Tamura, K., Nei, M., 1993. Estimation of the number of nucleotide substitutions in the control region of mitochondrial-DNA in humans and chimpanzees. Mol. Biol. Evol. 10, 512–526. Thompson, J.D., Gibson, T.J., Plewniak, F., Jeanmougin, F., Higgins, D.G., 1997. The CLUSTAL_X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Res. 25, 4876–4882. van Ham, R.C., et al., 2003. Reductive genome evolution in Buchnera aphidicola. Proc. Natl. Acad. Sci. U. S. A. 100, 581–586. Wernegreen, J.J., 2002. Genome evolution in bacterial endosymbionts of insects. Nat. Rev. Genet. 3, 850–861. Wernegreen, J.J., 2005. For better or worse: genomic consequences of intracellular mutualism and parasitism. Curr. Opin. Genet. Dev. 15, 572–583.