Complete sequence of the mitochondrial genome of Tetrahymena thermophila and comparative methods for identifying highly divergent genes

INTRODUCTION

The mitochondrial DNA of Tetrahymena pyriformis was among the first protist mitochondrial DNAs characterized (1,2). The Tetrahymena mitochondrial DNA is linear (∼47 kb in length). An examination of the mitochondrial genome of Tetrahymena reveals several unusual properties. Only seven different tRNAs are coded by the mitochondria. The remaining mitochondrial tRNAs are of nuclear origin and are imported into the mitochondria (3–7). The terminal ends are long inverted repeats that contain the large subunit ribosomal RNA (LSU rRNA) genes and the leucine tRNA genes (8,9). In the Tetrahymena mitochondrial DNA, both the LSU rRNA gene and the small subunit ribosomal RNA (SSU rRNA) genes are split. The nad1 gene is also split into two transcriptional units (7,10–12).

Most mitochondrial genomes contain genes that code for proteins involved in pathways of oxidative phosphorlyation and protein synthesis (13). In T.pyriformis there are 22 proteins with assigned function. Of these, 13 are involved in oxidative phosphorylation pathways, eight are ribosomal proteins and one is involved in cytochrome c biosynthesis. An additional 22 putative proteins are of unknown function (7).

Among ciliates, Tetrahymena thermophila and T.pyriformis are close evolutionary relatives, but within the more than 30 species comprising the T.pyriformis complex, T.thermophila and T.pyriformis are fairly distant relatives (14,15). Among all known organisms, Tetrahymena malaccensis is the most closely related species to T.thermophila. A detailed comparison of the mitochondrial genomes of T.thermophila and T.pyriformis allows the unequivocal identification of protein coding regions. The degree of divergence among these regions is indicative of the selective pressure acting on these various genes. In cases where the T.thermophila and T.pyriformis genome organization differs, we have examined the T.malaccensis mitochondrial sequences to gain further insight into the evolutionary histories of these mitochondrial genes.

Many ciliates, including T.thermophila and T.pyriformis, use a mitochondrial genetic code in which TAA is the sole termination codon, TGA codes for tryptophan and TAG is not used (13). The genes in the mitochondrial genome of T.pyriformis apparently use four initiation codons (ATA, ATT, TTG and GTG) in addition to the standard ATG (16). Polygenic transcriptional units are also indicated, but there is no evidence of RNA editing (16).

During the evolution of modern eukaryotes the mitochondrial genome has been dramatically modified in various different lineages. The vast majority of the genes present in the original endosymbiont have been lost from the mitochondrial genome. Most of these genes are simply gone from the cell, but a number of genes originally in the mitochondria have been transferred to the nucleus (17,18). In a few eukaryotic lineages, the mitochondria and the mitochondrial genome have been lost entirely (19,20). This process of elimination and transfer of mitochondrial genes has followed different evolutionary courses in different lineages.

The genes present in mitochondrial genomes fall into distinct categories: (i) rRNAs; (ii) tRNAs; (iii) ribosomal proteins; (iv) proteins involved in oxidative phosphorylation; (v) other proteins (not ribosomal proteins and not directly involved in oxidative phosphorylation); and (vi) open reading frames (ORFs) coding for putative proteins for which homologs of known function have not be identified (21). Among the putative mitochondrial protein coding genes, ∼50% are involved in oxidative phosphorylation, 30% are ribosomal proteins, 7% are other proteins and 13% are of unknown function. Thus, there are a limited set of genes that are expected among the mitochondrial ORFs, which makes their identification considerably easier.

The ciliates have among the most rapidly evolving mitochondrial genomes as evidenced by the large number of putative proteins for which function cannot be readily assigned. By comparison, in another protist, Reclinomonas americana, 92% of the mitochondrial genes can be assigned functions (22).

The most common means of assigning putative function to mitochondrial genes is to identify proteins of known function that have highly similar amino acid sequences (23). It is generally assumed that homologous proteins will have homologous functions. Our objective is to assign tentative functions to as many of the putative protein coding genes in the Tetrahymena mitochondrial genome as possible. We have examined a series of characteristics, in addition to sequence similarity, that may also indicate homology between proteins. These include physico-chemical properties such as theoretical pI, molecular weight and particularly the predicted transmembrane (TM) spanning regions.

MATERIALS AND METHODS

A B strain of T.thermophila (SB210, provided by E. Orias, University of California, Santa Barbara), T.malaccensis (MP75, obtained from ATCC) and T.pyriformis (GL, obtained from H. Buhse, University of Illinois, Chicago) were grown in PYG medium (24,25). Mitochondria were isolated from T.thermophila and total DNA was prepared according to Brunk and Hanawalt (2). The mitochondrial DNA was sheared to fragments of ∼3 kb in length using a shear point device (26). The ends of these fragments were filled in with T4 DNA polymerase and phosphorylated with polynucleotide kinase (27). Oligonucleotide adapters were ligated to the blunt ends of the mitochondrial DNA fragments and these fragments were directly cloned into pAMP 10 vectors (Life Technology, San Diego, CA). The sequences (∼500 bp) at each end of the cloned fragments were determined. These sequences (∼150 kb) were aligned using Phred, Phrap and Consed to form contigs (28–30). Clones that bridge the gaps between contigs were completely sequenced using specific primers to connect the contigs. The long inverted repeats in the terminal regions were independently PCR amplified using specific primers to achieve accurate consensus sequences for these regions. The final sequence has a >3-fold coverage of all regions and an error rate of 0.75 bases per 10 000 bases as estimated by the Consed program. The organization of the nad9 genes (tandemly repeated) and the atp9 gene in T.malaccensis were examined by cloning and sequencing these regions. The sequence of the T.thermophila mitochondrial genome was aligned with the sequence of the T.pyriformis mitochondrial genome, matching sequences of the ORFs at both the amino acid and nucleotide levels using Clustal X (31).

Table 1 shows the similarity metric (z score) used to measure the degree of sequence similarity between ORFs of mitochondrial genes in T.thermophila and T.pyriformis (32–34). The ORFs were translated into amino acid sequences, aligned and a normalized alignment score (NAS) was computed. The amino acid sequences were then randomized, aligned and additional NASs were computed. This process was repeated 1000 times and the mean and standard deviation of the resulting distribution was computed. The z score for two ORFs is the difference between their NAS value and the mean for the random distribution divided by the standard deviation of this distribution. This distribution approximates a Gumbel extreme value distribution (34). A z score >6 is expected by chance less than three times in 1000 trials, which is consistent with our random distributions. Thus, z scores >6 strongly indicate homology (34,35). A similar metric for sequence similarity was used in comparing nucleotide sequences.

The proportions of silent substitutions (Ks) and substitutions that change an amino acid (Ka) were calculated for each homolog in the T.thermophila and T.pyriformis mitochondrial genome and the Ks/Ka ratios were determined. The nine nucleotide changes possible for each codon were characterized as either silent or amino acid changing. These changes were averaged per nucleotide position. The actual nucleotide changes occurring between the sequences were characterized as well. When more that one nucleotide change occurred in a codon, all potential paths were considered on an unweighted basis and the changes were prorated per nucleotide position. The actual nucleotide changes per codon were normalized by the potential changes. These observed normalized nucleotide changes were corrected for multiple hits producing the Ks and Ka values (36,37). The Ks values and Ks/Ka ratios are shown in Table 1.

Each of the 45 ORFs in the T.thermophila mitochondrial genome was used as a query sequence with the NCBI BLAST program to produce a list of BLAST hits (38). ORFs with a list of BLAST hits containing numerous versions of the same protein (consensus protein) that had a low probability of occurring by chance were tentatively assigned the function of the consensus protein. Consensus proteins involved in oxidative phosphorylation and electron transport were compared to the homologs found in Bos taurus (cow), directly or via a chain of homology. Similarly, for the ribosomal proteins, T.thermophila was compared to Saccharomyces cerevisiae (yeast).

The T.thermophila ORFs were translated into amino acid sequences and aligned with the presumptive homolog using Clustal X (31). The aligned sequences were examined for the presence of sizable extensions at the N- or C-terminus and unaligned terminal additions >10 amino acids were trimmed. The z scores and trimming details for each homolog are shown in Table 2. The nad1 split gene was treated as a single protein. Protein properties, including the theoretical pI, ratio of negative to positive residues, aliphatic index and the grand average of hydropathicity (GRAVY) for all of the T.thermophila ORFs were calculated using ProtParam (39). TM of the T.thermophila ORFs were predicted using TMHMM version 2 (40). The protein properties and number of TM regions are shown in Table 2. The T.thermophila ORFs were also used as queries for PRODOM and pFAM protein domain databases (41).

RESULTS

An alignment of the T.thermophila and T.pyriformis mitochondrial genomes indicates that these genomes are virtually identical with respect to gene content and order (Fig. 1), except for the tandem duplication of nad9 in T.thermophila. Thus, we have labeled the T.thermophila ORFs with the same registered ymf designations used for the T.pyriformis genes (7).

The z scores calculated for each of the ORFs in T.thermophila and T.pyriformis (Table 1) all substantially exceed our criteria (z > 6) as an indicator of homology between the sequences. The RNA coding genes also show a high degree of sequence similarity as well.

Strong selective pressure on an ORF will favor silent substitutions (Ks) over nucleotide substitutions that lead to a change in an amino acid (Ka). The values for Ks (Table 1) indicate that, while there has been substantial nucleotide substitution in the T.thermophila and T.pyriformis mitochondrial genomes, the silent substitutions have not saturated the potential sites. The Ks/Ka ratio varies dramatically from gene to gene, providing an indication of the selective pressure on each protein.

As noted for the T.pyriformis mitochondrial genome, the T.thermophila mitochondrial genome appears to have a limited number of transcriptional units (7). To a first approximation there appear to be two transcriptional units, transcribed from a nearly central bi-directional promoter set (Fig. 1). The largest intergenic gap in the genome (424 bp) is at the position of the putative central set of bi-directional promoters. A short transcriptional unit, including the three genes transcribed in the 5′→3′ direction, interrupts the large 3′→5′ transcriptional unit.

The most apparent difference between the two genomes is the tandem duplication of the nad9 gene in the T.thermophila genome. Tetrahymena malaccensis, the closest known relative of T.thermophila, also has tandemly repeated nad9 genes. The tandemly repeated genes are highly similar within each species, but diverge considerably between the two species (Table 3a). Similarily, the long inverted repeats in the terminal regions are virtually identical within a species, but diverge significantly between species (Table 3b).

The mitochondrial genes in T.pyriformis apparently use four initiation codons (ATA, ATT, TTG and GTG) in addition to the standard ATG (16). The vast majority of the T.thermophila mitochondrial putative proteins appear to start at the same relative position as their T.pyriformis homolog and use the same initiation codon. Seven putative proteins appear to initiate at somewhat different positions in T.thermophila and T.pyriformis. The most plausible initiation codons for T.thermophila were assigned based on three criteria: (i) their proximity to the 5′ end of the ORF; (ii) their juxtaposition to the assigned T.pyriformis initiation codons; and (iii) the utilization of only initiation codons assigned in T.pyriformis. In five cases, proteins that initiate with unconventional codons in one Tetrahnymena species are initiated in the other species at an identical nucleotide position with ATG (Table 1). This strongly supports the correct identification of these unconventional initiation codons in these sequences.

By far the largest gene in T.thermophila or T.pyriformis is ymf77. It is highly improbable that a region of >4000 nt would be devoid of TAA stop codons or TAG unused codons in both species by chance. A comparison of the nucleotide substitutions between T.thermophila and T.pyriformis for ymf77 indicates that the selective pressure on this protein is probably not large (Ks/Ka = 1.61, Table 1).

Although most of the intergenic regions in T.thermophila and T.pyriformis are of similar length, the intergenic region following the T.pyriformis atp9 gene is exceptionally large (95 bp) compared to the intergenic in T.thermophila (14 bp). The atp9 gene and its flanking regions from T.malaccensis were sequenced and compared with T.thermophila and T.pyriformis. This comparison indicates that the T.malaccensis and T.thermophila atp9 genes and intergenic regions are virtually identical.

Our objective is to identify the probable function of each of the Tetrahymena genes based on their sequence similarity to proteins of known function. BLAST analysis provided initial identification of the Tetrahymena genes. Thirteen of the T.thermophila mitochondrial ORFs can be immediately assigned function by their similarity to standard proteins of B.taurus (Atp9, Cob, Cox1, Cox2, Nad1, Nad2, Nad3, Nad4, Nad5, Nad7 and Nad10) or S.cerevisiae (Rpl2 and Rps12) (Table 2). Two additional T.thermophila mitochondrial ORFs (Nad9_1 and Nad9_2) are linked via a chain of homology to B.taurus proteins and five further ORFs (Rpl14, Rpl16, Rps13, Rps14 and Rps19) are linked to S.cerevisiae proteins (Table 2). In this manner, the functions for 21 of the 45 T.thermophila ORFs can be reliably assigned by chains of sequence similarity comparison.

The remainder of the T.thermophila ORFs have BLAST hit lists that include a wide variety of proteins having different functions, with no clear consensus protein indicated. At this point the physico-chemical parameters for these protein sequences were examined for clues to their identity. The theoretical pI, ratio of negative to positive residues, aliphatic index and the GRAVY and TM were calculated for the ORFs with assigned function (Table 2).

TM are not present in the ribosomal proteins of Tetrahymena or in virtually any ribosomal proteins found in GenBank. Thus, the presence of TM strongly indicates that the ORF is not a ribosomal protein. Alternatively, most of the proteins involved in oxidative phosphorylation and electron transport have TM and many have multiple TM. However, some proteins of the oxidative phosphorylation NADH Fo complex (Nad5, Nad9 and Nad10) do not have TM. Generally, the theoretical pI for ribosomal proteins is high and rarely below 10, while the theoretical pI for proteins involved in oxidative phosphorylation and electron transport is low, usually below 8 and never above 10. The ratio of negative to positive residues for ribosomal proteins is low, rarely above 0.5, while this ratio is usually above 1.0 for proteins involved in oxidative phosphorylation and electron transport. Among the ribosomal proteins the aliphatic index is never above 100, while the aliphatic index for proteins involved in oxidative phosphorylation and electron transport is seldom below 100. The GRAVY for ribosomal proteins is always negative, averaging –0.51, while the GRAVY for proteins involved in oxidative phosphorylation and electron transport is seldom negative, averaging 0.49. Although these parameters are not totally independent of one another, taken together they present a valuable signature for ribosomal proteins and proteins involved in oxidative phosphorylation and electron transport.

The physico-chemical parameters for Ymf60 indicate that it is potentially a ribosomal protein (Table 2). It has no TM and all of the other parameters match those of ribosomes. Among the weak BLAST hits for Ymf60 were two hits for the large subunit ribosomal protein 6. We compared the sequence of Ymf60 to virtually all of the Rpl6 sequences in GenBank. The Rpl6 sequence from Thermatoga maritima was found to have sufficient similarity to Ymf60 to indicate homology (Fig. 2a), and a chain of homology linked Ymf60 to Rpl6 of S.cerevisiae (Table 2). Thus, Rpl6 function can be assigned to Ymf60 with confidence.

On the basis of the additional physico-chemical parameters, Ym62 and Ymf58 appear to be proteins involved in oxidative phosphorylation and electron transport. They have five and three TM respectively and they each have low pIs, ratios of negative to positive residues above 1.0, aliphatic indices well above 100 and positive GRAVYs. The TM profiles of Ymf62 and Ymf58 were compared with the TM profiles of representatives of all of the proteins involved in oxidative phosphorylation and electron transport whose function had not already been identified within the T.thermophila mitochondrial genome.

ORF Ymf62 has a TM profile very similar to that of Nad6 from B.taurus (Fig. 3A). Each sequence has four TM helices followed by a non-TM region of ∼40 amino acids followed by a fifth TM region. A survey of virtually all of the Nad6 proteins in GenBank indicates that the Nad6 protein of Porphyra purpurea has the greatest sequence similarity to the T.thermophila Ymf62 protein (Fig. 2b). A chain of homology was established linking Ymf62 to B.taurus Nad6 (Table 2). The combination of TM profile and sequence similarity allows us to assign Nad6 function to the T.thermophila Ymf62 protein.

The TM profile of Ymf58 is very similar to the TM profile of B.taurus Nad4L, both for the TM regions and the inter-TM spacings (Fig. 3B). A survey of virtually all of the Nad4L proteins in GenBank indicates that Nad4L of Nephroselmis olivacea has the greatest sequence similarity to the T.thermophila Ymf58 protein (Fig. 2c). Although the sequence similarity is only moderate (5.47), coupled with the similarity of TM profiles, homology is highly probable. A chain of homology extending from T.thermophila to the B.taurus Nad4L sequence via N.olivacea and Penaeus monodon can be readily established (Table 2). The function of Nad4L can be assigned to the T.thermophila Ymf58.

The ORF identified as yejR by Burger et al. (7) does not have an apparent homolog in either S.cerevisiae or B.taurus. However, this ORF is clearly homologous to a Triticum aestivum (wheat) gene involved in the biogenesis of c-type cytochromes via a homology chain through Pseudomonas putida (Table 2). The T.aestivum protein is a member of a family that consists of various proteins involved in cytochrome c assembly from mitochondria and bacteria (41). Thus, the T.thermophila YejR protein is assigned a similar function. This protein family differs from the heme lyase found in yeast (42).

Virtually all of the identified ORFs have reasonably strong hits with the appropriate domains in the PRODOM and pFAM databases. This strengthens our confidence that the appropriate function has been assigned to these ORFs. Most of the ORFs without assigned functions do not have any hits with domains in the PRODOM and pFAM databases. The few hits that do occur for these ORFs have not led to the identification of proteins with reasonably high sequence similarity.

One of the T.pyriformis ORFs, which was identified as Rps3 by Burger et al. (7), does not have significant sequence similarity with any of the Rps3 genes found in GenBank. No protein has been found that has significant sequence similarity to this ORF. Most of the physico-chemical parameters suggest that this protein is potentially a ribosomal protein. However, it has a predicted TM, which is not characteristic of ribosomal proteins (Table 4).

The T.thermophila mitochondrial genome has 19 ORFs, in addition to Rps3, that do not have sufficient sequence similarity to proteins of known function to permit the confident assignment of function. On the basis of their physico-chemical parameters these ORFs can be grouped into two broad classes: putative ribosomal proteins and putative non-ribosomal proteins. Table 4 lists the parameters of these ORFs along with the mean values of these parameters from the ORFs for which function has been assigned.

DISCUSSION

The pattern of nucleotide substitution between homologous genes gives a good indication of the selective pressure on a coding region. In general, mitochondrial proteins for which function can be assigned have higher Ks/Ka ratios than proteins of unknown function, whereas proteins involved in oxidative phosphorylation have among the highest Ks/Ka ratios. Genes atp9 and nad5 are exceptions. Their Ks/Ka ratios indicate less selective pressure than expected for proteins with assigned function. The divergence at the C-terminus of atp9 may account for its low Ks/Ka ratio, but the reason for the low Ks/Ka ratio of nad5 is unclear. The proteins for which function cannot be readily assigned have generally lower Ks/Ka ratios, with some exceptions like Ymf56, Ymf65 and Ymf68, which have Ks/Ka ratios >10. The largest ORF, ymf77, has a very low Ks/Ka ratio (1.61) suggesting reduced selective pressure, which may contribute to the divergence of this protein and the difficulty in matching it with a protein of known function.

The mitochondrial genome is very compact with very short intergenic regions (3.8% of the genome), only three of which are longer than 63 bp. These intergenic regions appear to be too small to accommodate promoters. This is consistent with the transcript mapping of the T.pyriformis mitochondrial genes which suggested multi-gene transcripts (16). A working hypothesis is that the entire genome is transcribed from a central bi-directional promoter, which would account for the transcription of all of the genes except the three gene cluster in the 5′ portion of the genome. These three genes could be transcribed from a unidirectional promoter in the 5′→3′ direction (Fig. 1).

The three long intergenic regions in T.thermophila have sequence z scores indicating that these regions are under strong selective pressures (data not shown). The largest intergenic region (493 bp), between ymf77 and cob, corresponds to the site of the putative central bi-directional promoter set. This position also corresponds to the origin of DNA replication identified in electron micrographs of partially replicated mitochondrial molecules in T.pyriformis (43). The sequence of the second largest intergenic region (261 bp) is highly conserved relative to the T.pyriformis homolog. It is located upstream from and immediately adjacent to the three gene cluster transcribed in the 5′→3′ direction. This may be the site of a unidirectional promoter responsible for transcribing this three gene cluster.

The most notable difference between the T.thermophila and T.pyriformis mitochondrial genomes is the tandem duplication of the nad9 gene in T.thermophila. A similar duplication occurs in T.malaccensis and the divergence between the genes from different species (orthologous) is much greater than the divergence among genes within the species (paralogous) indicating concerted evolution (Table 3a). Concerted evolution also appears to be operating on the terminally repeated regions (Table 3b). The regions under concerted evolution include the two split portions of the LSU rRNA gene, the leucine tRNA gene interrupting the LSU rRNA sequence and the two short (3 and 9 bp) intergenic regions flanking the tRNA gene. The telomeres at both ends of the mitochondrial DNA molecule (53 bp repeats) are also identical although they start at different positions in the telomere repeat sequence. Morin and Cech have proposed that unequal crossing over and terminal hybridization may be part of a process maintaining the telomeres (11,44). Whether this is the mechanism by which telomeres are regenerated or not, it is not a plausible mechanism for the concerted evolution of the inverted terminal repeats as it would be expected to homogenize the terminal tRNA genes. Concerted evolution of paralogs is relatively common (45–47).

The sequences for the C-terminus of the Atp9 protein and the intergenic region following the atp9 gene from T.thermophila and T.pyriformis have diverged significantly. The T.malaccensis and T.thermophila atp9 gene sequences are very similar at the C-termini and they are more similar to the Atp9 proteins from other organisms than is the T.pyriformis sequence. The most probable evolutionary history is that a DNA sequence was inserted into the C-terminus of the atp9 gene in the lineage leading to T.pyriformis, creating a slightly different C-terminal sequence and substantially lengthening the following intergenic region.

One of the putative mitochondrial proteins in Tetrahymena, Ymf77, presents a particular challenge. This gene has 1321 amino acids (almost 9% of the Tetrahymena mitochondrial genome) and apparently has 15 TM. It is conserved in both T.thermophila and T.pyriformis, although not found in Paramecium aurelia (48). In spite of its size and apparent function, BLASTP does not yield any reasonable protein candidates. Additional mitochondrial sequences from ciliates intermediate between Tetrahymena and Paramecium may offer additional versions of Ymf77 that give clues to its function.

The proteins found in mitochodrial genomes are relatively limited in potential function, which simplifies the identification of mitochondrial ORFs (13). In most mitochondrial genomes all of the genes are identified, even R.americana with 97 mitochondrial genes has virtually all of its genes functionally identified (22). In spite of this, almost half of the putative proteins in Tetrahymena remain unidentified. Clearly the mitochondrial genomes of the ciliates have diverged dramatically from other mitochondrial genomes, however, the evolutionary mechanisms responsible for this divergence are not immediately obvious.

Although sequence similarity remains the primary means of establishing homology from which function can be assigned, physico-chemical parameters provide a powerful additional indicator of potential homology. These parameters can be combined with nominal sequence similarity to establish homology. They are also valuable in clustering unidentified proteins into specific groups.

DDBJ/EMBL/GenBank accession no. AF396436