A Different Approach to Anti-Aging

In this fast-paced talk by Dr. Aubrey de Grey, we hear about a breathtakingly radical approach to forestalling aging, which (in a nutshell) involves moving all remaining mitochondrial genes into the nuclear DNA, so that mutations in mitochondrial DNA, per se, are rendered irrelevant. This technique of "obviation" of mitochondrial-DNA breakdown isn't a new idea (it's been around for at least 30 years). What's new is that, technologically, we're in a position to make it happen.

Most mitochondrial genes are, of course, already in the nucleus. The majority of scientists accept that mitochondria got their start as bacterial endosymbionts; a long-ago ancestor of today's alphaproteobacteria took up residency in an anaerobe. The anaerobe provided the invading bacterium with a nutrient-rich environment in which to live, while the bacterium provided oxygen-detoxification services (and a lot of adenosine triphosphate) to the host. Most likely, the invading bacterium had around 1,500 genes. Over time, ~500 redundant genes were lost and the remaining 1,000 or so migrated to the host cell's nuclear DNA (a much safer environment for DNA than the mitochondrion), leading to the present-day situation where (human) mitochondria have an extremely small circular chromosome encoding just 13 proteins. But we know mitochondria actually contain around 1,000 different proteins, most of which are encoded in nuclear genes.

The majority of mitochondrial genes (in the nucleus) encode proteins that are made in cyotplasm and imported into the mitochondrion. For import, proteins must be in an unfolded state. Some proteins (the most hydrophobic ones) are actually made on the surface of the mitochondrion and slurped into the interior of the mitochondrion as they're being made. Folding of the proteins then takes place inside the organelle.

Can the remaining 13 mitochondrial protein genes be moved to the nucleus? If we succeed in doing that, will cells live longer? What technical obstacles remain? What progress has been made? These and other questions are addressed in Aubrey de Grey's talk, which is well worth a listen.

Thermometer Genes

Heat shock proteins are an interesting class of proteins that provide "damage control" for enzymes when temperatures rise to the point where proteins start to unfold and refold improperly. Protein 3-dimensional structure is critical to proper enzyme function, and it doesn't take much thermal jostling to mess up a protein's structure. Therefore it's not surprising cells have their own miniature repair factories for refolding heat-misfolded proteins.

Collectively, heat shock proteins are part of a group of proteins known as chaperones, some of which are bonafide refoldases and others of which aid proteins in other ways. (For example, the ClpB protein rescues proteins from an aggregated state.)

GroEL is a so-called type I chaperonin involved in protein folding, assembly, and transport. Like many heat shock proteins, it's over-expressed at high temperatures and plays a critical role in growth and survival at non-permissive temperatures. Because of its importance in many cellular processes, GroEL is ubiquitous in bacteria, with most species having a single GroEL gene, but with about 30% of  genomes having two or more GroEL copies.

GroEL mRNA of Mycobacterium intracellulare can fold into the advanced low-energy secondary structure shown here.

The question of how cells up-regulate heat shock proteins during times of thermal stress is still largely open, although we know in some cases specific transcription regulator proteins are involved (but then the question becomes: how do the regulators know to up-regulate in times of heat stress?). The answer might not be that difficult. The messenger RNAs encoding GroEL and other heat shock proteins contain a great deal of secondary structure (that is to say, the RNA folds back on itself to form thermally sensitive structures). Recently, Wan et al. surveyed RNA thermal sensitivity in yeast and found thousands of so-called "mRNA thermometers": RNA molecules that unfold in response to heat. About three quarters of yeast RNA is thermo-stable at 37 degrees C, while around 55% of RNAs are unfolded at 55 degrees. In the folded state, RNA probably requires the help of helicases or other "helpers" to unfold, but at a high enough temperature, the molecules unfold by themselves and become eligible for translation by ribosomes.

To investigate the possible role of mRNA secondary structure in GroEL regulation, I wrote scripts that check a gene for all occurrences of length-9 (so-called "9-mer") nucleotide sequences that have a corresponding reverse-complement sequence in the same gene. When I checked the GroEL gene of Mycobacterium tuberculosis (Erdman strain), I found 14 pairs of complementarity 9-mers, representing regions of the gene that could, in theory, cause secondary structure to form in mRNA. A check of the sister organism M. intracellulare (whose GroEL gene is 80% identical to the M. tuberculosis version) showed 22 such complementary pairs.

Interestingly, the mutational differences between GroEL in M. tuberculosis and M. intracellulare do not appear to be randomly distributed along the gene. In M. tuberculosis, mutations occur at a rate of  0.12698 substitutions per site inside 9-mer regions (putative stems) versus a rate of 0.19722 for non-9-mer regions, indicating that (perhaps) selection pressure is different for self-complementing regions than for other regions. I found much the same thing in M. intracellulare, where the mutation rate was 0.16414 inside 9-mers and 0.19347 elsewhere.

Tending to confirm that selection pressure is different for the "secondary structure" regions versus other regions is the (surprising) finding that in M. tuberculosis complementary regions, the ratio of non-synonymous to synonymous mutations (Kn/Ks) is 0.526, versus 0.950 for other regions. In M. intracellulare, likewise, Kn/Ks is less in complementing regions (0.635) than in non-complementing regions (0.975).

To check whether these observations apply only to Mycobacterium or might be more widely applicable, I took a look at GroEL genes in Clostridium acetobutylicum strain ATCC 824 and Clostridium lentocellum strain DSM 5427. The Clostridia are phylogenetically quite distant from Mycobacteria (as confirmed by the fact that their GroEL genes share only 50% nucleotide sequence identity). A total of 462 mutations separated the two Clostridial genes. But again, the mutations segregated non-randomly according to whether they occurred in putative regions of complementarity (secondary structure) as opposed to non-complementing regions. In C. acetobutylicum the mutation rate in complementing regions was 0.23015 substitutions per site (29/126 bases) versus 0.28618 (433/1513 bases) for non-complementing regions, while in C. lentocellum the rates were
0.19047 (24/126) substitutions per site vs. 0.28949 (438/1513). For C. acetobutylicum the Kn/Ks ratios were 0.277 in 9-mers and 1.466 otherwise. For C. lentocellum, Kn/Ks was 0.538 in 9-mers and 1.415 outside 9-mers, tending to confirm that selection pressures are different in stems than in loops.

Bottom line, the data are consistent with a scenario in which secondary structure of GroEL mRNA (and/or ssDNA) plays a role in heat-activation of the gene, such that when temperatures exceed the melting point of secondary structures, the gene is eligible for transcription and/or translation. The gene is, in effect, its own thermometer.

How Do Bacteria Survive Radiation Damage?

Secondary structure of the Trad_1400 gene (encoding a MutT hydrolase) in Truepera radiovictrix.

In the 1950s, a tin of meat was exposed to a dose of radiation that was thought to be capable of killing all known forms of life, but the meat subsequently spoiled, and Deinococcus radiodurans (dubbed Conan the Bacterium by some) was isolated from it. Various members of the Deinococcus-Thermus group have shown themselves to be incredibly hardy, able to survive extremes of temperature and doses of radiation that, frankly, they shouldn't be able to survive. (If any group of bacteria were able to survive space travel, this group surely could.)

Members of the Deinococcus group probably learned their DNA repair tricks very early in the history of terrestrial life, before there was sufficient oxygen in the atmosphere to support an ozone layer. In those times (before about one billion years ago), ultaviolet light from the sun would have been strong enough to sterilize almost any exposed surface. UV radiation, when it's strong enough, causes single- and double-strand breaks in DNA (just as ionizing radiation from radioisotopes or cosmic rays will). How Deinococcus manages to survive such radiation is still something of a mystery, although several repair modalities have been elucidated. We know that these bacteria have high copy numbers of their genetic material, and this no doubt facilitates repair. Still, double-stranded DNA breaks, in most organisms, are quickly fatal if they accumulate.

It turns out, DNA from Deinococcus-group bacteria is unusually rich in internal (intra-strand) complementarity, which means single strands of DNA are capable (in theory) of folding back on themselves to form elaborate secondary structures of high thermal stability. One such structure, for the Trad_1400 gene of Truepera radiovictrix (a radiation-tolerant member of the Deinococcus group), is shown above. This particular structure has a 37°C Gibbs free energy of minus-71.47 kcal/mol (meaning the structure is more likely to form than randomly coiled ssDNA) and a Tm (melting temperature) of 63.1°C, meaning it should be thermo-stable to around 145°F. Almost the entire gene folds back on itself; the only portion that doesn't self-anneal is the flat line on the bottom containing 29 bases.

If the (separated) strands of Truepera DNA can assume stable self-annealed structures of this type, it would go a long way toward explaining how the organism could survive double-stranded breaks. Fire a random bullet at the DNA and you're bound to hit secondary structure, not canonical (B-form) duplex DNA. A double-strand break in a stem structure might liberate a stem/loop from one strand, but the other strand could unfold to form a template for immediate repair of the damaged strand. Something like this is probably going on in radiation-resistant Deinococcus members, which have evolved to allow more than the usual secondary structure in their DNA.

March of the Microbes

It can be hard to convey to non-biologists, who may think of germs as pests (an unseen, unclean residue to be neutralized with disinfectants), the grandeur and miraculous intricacy of microbial life. After more than three billion years of evolution, microbes are still, numerically speaking, the principal inhabitants of this planet; most of life on earth is invisible to the naked eye. But sightings of microbial life are all around us if we know where to look. In March of the Microbes: Sighting the Unseen (2010, Harvard University Press), Professor John L. Ingraham shows those of us with eyes to see just how and where to look.

It's impossible for me to come at this book with anything approaching objectivity, because I know the author well: I was Dr. Ingraham's student at U.C. Davis in the 1970s and did my graduate studies under his kind, insightful direction. And it was thanks largely to him that I developed what has turned out to be a deep, abiding, lifelong appreciation for the stunning ingeniousness (for there's no other way to put it) of microbial life.

The Ingraham lab, in those days, specialized in the study of pyrimidine metabolic pathways in Salmonella, but Dr. Ingraham's interests were broad and his understanding of microbial ecology impressively deep. Before undertaking some of the first studies of cold-sensitive mutants of bacteria, Ingraham had made important contributions to the understanding of fusel oil production by enologically important yeasts. But Ingraham was and is also an avid naturalist, as likely to identify a particular oak species as to diagnose the fungal disease causing its slow death. Serve him rainbow trout and he'll likely comment on the multiple layers of guanine crystals in the fish's shimmering scales before reminding you (as he does on page 24 of March of the Microbes) that the lack of a "fishy" smell in rainbow trout is due to the absence of trimethyl amine oxide, an osmolyte that would otherwise be converted by bacteria to oh-so-smelly trimethyl amine (TMA).

Go with Ingraham on a walk in the woods and he's sure to point out (as on p. 76 of MoM) that the age of a lichen on a boulder can be estimated from its size, based on its half-millimeter-a-year growth rate; and don't be surprised if he can name the particular cyanobacterium that gives the lichen its greenish hue. You'll forgive him, of course, when he uproots a small plant to show you the Rhizobium-fostered nodules on its roots, nodules that (when sliced open with a pocket knife) are red inside, crimson with the plant hemoglobin (p. 118) that only a symbiotically infected legume can produce. If you're worried about the safety of drinking water from that nearby mountain stream, he'll explain that you'd have to drink 250 gallons of such water before you're likely to imbibe a cell of Giardia (p. 298).

Reading March of the Microbes, I felt like I was on an extended nature walk with a world-class naturalist, except that this "walk" knows no physical bounds, for it extends from Sierra Nevada ski slopes, where slope operators mix dry powder made of Pseudomonas syringae cells (sold as Snomax) into the snow-machine slurry to facilitate ice-crystal nucleation (p. 263), all the way to miles-deep thermal chimneys on the ocean floor, where Pyrolobus furnarii grows in total darkness at 3000-psi pressure and 113-degree-Celsius temperatures (p.186); and from there to the inside of the 22-gallon rumen of a cow (p. 80), then to the 200-microliter hindgut of the termite, where, again (as in the cow), microbes do the vital work of digesting cellulose, something no higher life form can do.

The interconnectedness of life on earth is a recurring theme in March of the Microbes, as is the sheer adaptability of microbial life. We learn that only microbes can harness molecular nitrogen directly from the atmosphere; only microbes can return it. Only microbes can respire hydrogen sulfide and nitrate to sulfuric acid and nitrite, carving out caverns (such as the famous Carlsbad Caverns) in total darkness. Only microbes can utilize cellulose as a carbon source. Only microbes can survive temperatures well below freezing or significantly above boiling, or live through prolonged exposure to concentrated radiation in a nuclear reactor. Microbes are responsible for the world's great sulfur deposits; the calcium carbonate in the white cliffs of Dover; the 100-meter-thick diatomaceous earth pits in Clark County, Nevada; the massive saltpeter deposits (formed by decomposition of guano) in the coastal ranges of Chile's Atacama Desert. Only Archaea (prokaryotic extremophiles that were once thought to be bacteria) can make methane. And only a bacterium (Alcaligenes eutrophus, for example; there are others) can store 80 percent of its body weight as plastic (polyhydroxybutyrate).

Most bacteria are harmless to humans and other life forms, but March of the Microbes doesn't shy away from pointing out the exceptions. We learn about the (fungal) origins and history of potato blight, the history of syphilis and smallpox (one being a gift of the New World to the Old; the other, vice versa); the curious mode of action of botulism poisoning, and tetanus; the bacterial origins of peptic ulcers; the recent rise of the oddly pathogenic E. coli O157:H7; and the many intricate and interesting ways in which antibiotic resistant bacteria develop drug resistance. An entire chapter is devoted to viruses, which are rather harshly (but deservedly, some would say) villified, despite their contribution to evolution.

If one organism is glaringly absent from March of the Microbes, it would be (drum roll, please) Amoeba proteus, the prototype of the "lowly amoeba," which of course is not so lowly when you consider that the amoeba genome contains 100  times the DNA of a human genome. The freshwater amoeba is, conceptually at least, the progenitor of the modern-day phagocyte (white blood cell); many consider it the "nursery" organism on which Mycobacteria and other pathogens practiced and perfected their anti-phagocytic skills before adapting to warm-blooded animals.

Aside from the lack of an amoeba sighting, March of the Microbes lacks for very little, combining (as it does) refreshingly edifying bits of natural history and ecology with the occasional dash of epidemiology and immunology, plus just the right amount of common-sense chemistry, to provide a deeply satisfying view into the astonishingly diverse, always-full-of -surprises microbial world. It's the kind of book you wish you could never reach the end of, a 306-page Aha moment; a delight from start to finish.

A Manganese Catalase Fusion Protein

In science, it often happens that finding the answer to a particular mystery only leads to further questions. That's certainly the case with the non-heme/manganese-based catalases in bacteria (which I talked about in a previous post). Regardless of origin, catalases faclitate the breakdown of hydrogen peroxide to water and oxygen. The nearly universal heme-based catalase (found in almost every living thing) comes in large-subunit and small-subunit varieties but is always fairly big (726 amino acids, in E. coli). Manganese catalases, by contrast, are always fairly small: around 276 amino acids. But there's one group of manganese catalases that comes in at 416 to 431 amino acids, more than 50% larger than the "typical" manganese catalase. I wanted to know why. Why are these catalases bigger? 

Finding the answer wasn't hard (I got lucky). But the answer only leads to more questions.

It turns out very few organisms manufacture the "big" manganese catalase. The organisms in question belong to just 3 genera: Rhizobium, Bradyrhizobium, and Rhodopseudomonas. (The closely related Mesorhizobium has a Mn catalase, but it's the "normal" size Mn catalase. Meanwhile, its cousin, Sinorhizobium, has no Mn catalase.)

The first 280 or so amino acids of the Mn catalases made by Rhizobium, Bradyrhizobium, and Rhodopseudomonas align quite well with the normal-size Mn catalases made by other organisms, the only difference being the 140-amino-acid trailer on the end of the "long versions." I used the alignment editor in Mega6 (tantamount to Notepad) to Cut the trailer portion out of one of the sequences and Paste it into the BLAST search field at UniProt.org. A search against the protein database revealed something quite interesting: The trailer portion of the Bradyrhizobium Mn catalase is a 45%-identity match for the enteric-bacteria yciF gene (which is quite a good match considering the phylogenetic distance between E. coli and Bradyrhizobium).

The E. coli yciF gene (top) is a 58% match for the Bradyrhizobium yciF gene (middle), which in turn is a 64% match for the trailer portion of katN (manganese catalase gene) of Bradyrhizobium, bottom.

Further investigation left little doubt that the "big" Mn catalases are indeed fusion proteins: yciF fused to the C-terminal end of the Mn catalase.

But: What in the world is yciF? It turns out to be a very widely distributed, highly conserved protein of uncertain function. The protein has been purified and its crystal structure determined at 2.0 Å resolution, but its function is still uncertain. What we do know is that it is stress-inducible (along with other yci-series genes) and seems to bind a metal ligand, probably iron; and it shares structural features with rubrerythrin, a non-heme iron protein implicated in oxidative stress protection in anaerobic bacteria and archaea. In E. coli strains that have Mn catalase, the yciF gene occurs two genes upstream (on the 5' side) of the catalase (katN) gene.

Because yciF is more widely distributed (and more highly conserved) than manganese catalase, and because most Mn catalase producers (including those with the fusion enzyme) have an additional, separate copy of yciF, it seems likely (to me, anyway) that the fusion protein was created by chance in the common ancestor of the Rhizobiales when the original katN gene was laid down by a phage or other mobile genetic element. (In E. coli, katN often occurs near phage genes.) Sinorhizobium lost the combo gene entirely, while Mesorhizobium (which makes a small Mn catalase) either obtained katN on its own or lost the 3' trailer from its fusion protein over time.

Since Bradyrhizobium (and the others) already have a separate yciF gene, it's a mystery why the trailer portion of the fusion gene continues to exist. It might very well provide a favorable enhancement of katN function somehow (maybe exploiting iron in an auxiliary catalytic center). If the trailer's doing nothing useful, it should have disappeared over time. (Maybe it did disappear from other Rhizobiales members, and just hasn't disappeared yet in the three genera that still make the "big" enzyme.) I have a feeling the trailer piece does do something useful. Unfortunately, no one has characterized the Braydrhizobium Mn catalase yet, experimentally. We'll probably have to wait until that happens to find out what the "big" enzyme is capable of.

The Other Catalase

Microbiology students are taught from Day One that strict anaerobes (organisms that are killed by exposure to oxygen) lack the enzyme catalase, which breaks down hydrogen peroxide to water and O₂. You're already familiar with this enzyme if you've ever poured peroxide on a wound and seen it get foamy (blood is rich in catalase) or if you've used a peroxide bleach on your teeth, which causes saliva to become thick and gummy with microscopic bubbles. Catalase is a nearly universal enzyme, but strict anaerobes lack it (seemingly), because they are so seldom exposed to oxygen or its metabolites.

An unexpected result of genome data mining is the finding that non-heme (and thus non-iron-containing) catalase exists in a wide variety of bacteria once thought to contain no catalases. Manganese-containing catalase was first described in Lactobacillus, an organism that produces no cytochromes (and no porphyrins). The same type of enzyme was later experimentally verified in a thermophilic archeon, Pyrobaculum. It now turns out that many strict anaerobes previously believed to be catalase-free (such as most Clostridium species) contain this enzyme.

Phylogenetic distribution of manganese catalase (click to enlarge). At the top are the spore-forming Bacillus, with non-spore-forming Firmicutes (Aerococcus et al.) in a sub-branch below, and the cyanobacteria (Nostoc and Microcoleus) as an out-group to the Bacillus group. One cyanobacterial species, Cyanothece sp. PCC 7424 (a rice-field isolate), occurs as an out-node amongst the Gammaproteobacteria, indicating possible horizontal gene transfer. Archaeal organisms (Halalkalicoccus etc.) with this enzyme tend to be salt-lovers, although Pyrobaculum (not shown), which thrives in 2% salt, also has it. The phylo-tree was constructed from protein sequence alignments using Mega6 freeware. Node assignments were tested with 500 bootstraps.

Why would an anaerobe need catalase? Short answer: because small oxygenated molecules (like hydrogen peroxide and superoxide anion) are damaging to DNA, typically causing guanine to become oxidized to 8-oxo-guanine (which mispairs with adenine). Also, molecular oxygen irreversibly poisons the nitrogenase enzyme that many Clostridium members (and others) rely on to metabolize atmospheric nitrogen.

You would think that if molecular oxygen is a product of catalase action (and damages nitrogenase), Clostridium would not want to have catalase around in its cytoplasm. But Clostridium doesn't localize its catalase to the cytoplasm. It's a spore-surface enzyme. And this is key to understanding its distribution in nature.

Manganese catalase is rather sparsely distributed. It occurs just in certain taxonomic groups and certain species (see illustration, above). The enzyme seems to have been invented by the spore-forming Firmicutes (Bacillus and Clostridium) as a spore-coat enzyme. Many phylogenetically younger Firmicutes that have lost the ability to form spores also have the enzyme, as do some (not all) members of the Rhizobiales (not shown above), probably as an adaptation to low-iron niches. (The canonical heme-containing form of catalase requires iron.) There is evidence, in fact, that Lactobacillus leads a completely iron-free existence.

Where else do we find manganese catalase? Some (but not all) cyanobacteria have it. This is noteworthy in that the cyanobacteria form a specialized, environmentally hardened sessile cell called an akinete. (It's also noteworthy that both cyanobacteria and certain Clostridium members have nitrogenase.) The cyanobacteria that have the manganese catalase are primarily land-dwellers, however, not marine bacteria. This includes Anabaena (a moss symbiont and rice-paddy dweller), Microcoleus (which occurs in arid soils), and certain Nostoc members (which occur on rocks, in lichens), plus Cyanothece (from rice paddies).

Curiously, very few marine organisms have manganese catalase, the chief exceptions being Pirellula and Planctomyces. (And again, interestingly, Pirellula has a specialized sessile form as well as a motile form.) Many halophilic (salt-loving) archeons have the enzyme; surely those qualify as marine organisms? Not really. Halococcus, Natrinema, etc. are not planktonic; open-ocean waters have too little salt. (These organisms require upwards of 30% salinity, much stronger than the 3.5% salinity of sea water.) The salt-loving archeons live in drying-up seas (like the Dead Sea or Great Salt Lake) and at the edges of beaches, where salt concentrations skyrocket. These are, in effect, "terrestrial marine organisms," if that makes sense.

You can also find manganese catalase in some members of the Gammaproteobacteria (namely certain E. coli strains and some Pseudomonads), though curiously not in Shigella, Yersinia, or the Vibrio family (which is largely marine). The spotty nature of the enzyme's distribution among the enterics and Pseudomonads (which have no specialized sessile form) speaks to a possible horizontal gene transfer scenario.

There is little reason to believe manganese catalase is primordial. It certainly didn't come from the sea. In the first place, the enzyme is absent from Pelagibacter, Vibrio, and other important marine organisms. Secondly, sea water contains surprisingly little manganese (less than a part per billion). Marine organisms that do have catalase tend to have the heme-containing version of the enzyme (which makes sense, in that the oldest photosynthetic organisms long go mastered the art of porphyrin synthesis). Manganese catalase is a terrestrial adaptation, primarily of spore-, heterocyst-, and akinete-formers, with others obtaining the enzyme by lateral gene transfer. It's interesting that the Gammaproteobacteria that have manganese catalase (some enterics and some Pseudomonads) are opportunistic pathogens. They probably find the enzyme useful in combating the respiratory burst of phagocytes.

Odd Structures in a Prophage Genome

Contrary to what it looks like, this is not an arrival-gate diagram for an eastern European airport. It's a diagram of secondary structure of a portion of the mRNA for the yoqJ gene in B. subtilis phage SPBc2. Click to enlarge.

The common soil bacterium Bacillus subtilis has long been known to harbor an inducible prophage (virus) called SPBc2. Normally, SPBc2 is dormant in the Bacillus genome, but heat or treatment with mitomycin C will cause the phage to express itself and start a lytic infection cycle. Like many bacterial viruses, phage SPBc2 has a sizable double-stranded DNA genome, consisting (in this case) of 134,416 base pairs encoding 185 genes. It's one of the strangest collection of genes you'll ever want to meet.

The SPBc2 genome is strange, first of all, in its base composition. The following stats come from analysis of the coding regions of the genome:

Base	Abundance
A	0.3670
T	0.2825
G	0.1983
C	0.1506

Notice how the purines, A and G (adenine and guanine) are about 30% more abundant than the pyrimidines, C and T (cyotsine and thymine). This is a stark violation of Chargaff's Second Parity Rule (which says A=T and G=C not only for double-stranded DNA but for a single strand, as well).

When purine abundance is tallied on a per-gene basis, you get the following histogram:

Purine abundance for N=185 protein-coding genes in B. subtilis phage SPBc2. Only 8 genes contain less than 50% purine bases.

All but 8 (out of 185) genes have a message-strand purine content above 50%.

And if you look at the gene sequence data, the genes even look funny to the naked eye, containing, as they do, long runs of bases, with funny repeats, like TTGAAAGGAAAAAAAGACGGCCTAAATAAA. One gets the feeling, when  examining the sequence data, that one is looking at tRNA or rRNA data rather than protein-coding sequences. And yet, they really are protein-coding sequences.

But it turns out there's a lot of secondary-structure info in the sequences. (See the picture at the top of this post.) Many of the genes contain self-complementing intragenic regions. I decided to verify this with some custom scripts. First, I had scripts look inside each gene for length-10 nucleotide sequences with a matching reverse-complement sequence further downstream (in the same gene). I expected to find one or two (or ten) hits, based on the fact that any given length-10 sequence of the bases A, G, C, and T can turn up at random in about one in every million base pairs. (Each base is worth about two bits, so a length-10 DNA sequence can encode as much info as a length-20 binary sequence; 2-to-the-20th is a little over a million.) What I found was an astounding 419 complementary length-10 sequences within 89 genes.

When I upped the sequence length to 12, I still found 64 intragenic complement sequences in 38 genes. A length-12 match should happen randomly about once every 16 million base pairs. The SPBc2 genome (as I said) is 134,416 base pairs long.

The high degree of internal complementarity means that when the phage's DNA strands are separated during transcription or replication, they probably assume very particular 3-dimensional structures, and also, the messenger RNA made from the phage's genes are probably rich in secondary structure. Exactly why those structures are needed is anyone's guess. They may protect the virus from restriction nucleases (although frankly this chore is already taken care of by the prophage's own methylases). Or they may attract ribosomes in some special way (many of the genes form tRNA-lookalike structures). They may help with virion packaging. Or they may do nothing special at all. I kind of doubt the latter possibility, though.

A Gene Heatmap

Lately I've been using the great tools at genomevolution.org plus custom Canvas API scripts to render colorful heatmaps of aligned genes from phylogeneticaly diverse microorganisms. The following graphic is one such.

Glyceraldehyde-3-phosphate dehydrogenase genes of N=134 bacterial species, arranged in order of gene G+C content (high GC at the top). Hot colors are G and C. Cool colors are A and T.

What are we looking at? This is actually a composite rendering of the glyceraldehyde-3-phosphate dehydrogenase genes (DNA sequence info) from 134 bacterial species. Each gene is painted left-to-right (5' to 3') in a strip 4 pixels tall, with hot colors assigned to DNA bases G and C (guanine and cytosine), and cool colors assigned to bases A and T (adenine and thymine). Wherever there's a G or C, it gets painted red or red-orange. Wherever there's A or T, blue or blue-green. Same gene, 134 versions, varying significantly in G+C content. (The gene GC content ranges from a maximum of 69.2% at the top to 29.4% at the bottom.)

Why glyceraldehyde-3-phosphate dehydrogenase (GAPDH)? No real reason, except that it's a fairly universal (indeed, quite ancient) metabolic enzyme, reasonably compact (making possible a rendering that's not super-wide, as it would be for a larger gene), well-delineated genetically (not a fusion protein or an enzyme with multiple isoforms), and probably representative of a good many core metabolic enzymes. This is the enzyme that catalyzes the sixth step of glycolysis (sugar-breakdown). You may recall from Biochem 101 that the breakdown of glucose proceeds by splitting the twice phosphorylated molecule into two 3-carbon pieces. The triose phosphates in turn get phosphorylated by GAPDH before they transfer a phosphate to ADP to yield ATP, the 5-hour energy drink of all cells everywhere.

Alignment of genes was done via ClustalW in MEGA6 freeware. Rendering of the alignment FASTA file took about two seconds, in the browser, using 133 lines of custom JavaScript.

Not All Organisms Use Amino Acids the Same Ways

Organisms vary greatly in the GC (guanine plus cytosine) content of their DNA, and yet all organisms can still make ribosomal proteins, DNA and RNA polymerases, and the various other essential proteins of life, no matter what their DNA vocabulary limitations might be. A high-GC organism like Streptomyces can make a given enzyme (DNA polymerase, say) using mostly G and C bases in its DNA, but a low-GC organism like Clostridium botulinum can also make the same kind of enzyme, even though it uses mostly A and T in its DNA. How is this possible?

It's possible in part because of the many synonyms for amino acids available in the genetic code. But it's a mistake to think the same amino acids are used in equal numbers by high-GC organisms and low-GC organisms. Organisms at opposite ends of the GC spectrum use different amino acids.

I was curious to see which amino acids correlate most strongly with genomic GC, so I gathered codon usage statistics for 109 organisms of widely varying genomic GC content and used JavaScript to calculate Pearson correlation coefficients for all 20 amino acids with respect to  GC content. The results are shown in the following table.

TABLE 1. Correlation (r) between amino acid usage and genome GC content (N=109 organisms). 

Code	Amino Acid	r
A	Alanine (Ala)	0.9634
R	Arginine (Arg)	0.9495
G	Glycine (Gly)	0.9472
P	Proline (Pro)	0.9436
V	Valine (Val)	0.7725
W	Tryptophan (Trp)	0.7497
H	Histidine (His)	0.4660
L	Leucine (Leu)	0.3364
D	Aspartic Acid (Asp)	0.3347
T	Threonine (Thr)	0.3099
C	Cysteine (Cys)	-0.1280
Q	Glutamine (Gln)	-0.1668
M	Methionine (Met)	-0.2863
E	Glutamic Acid (Glu)	-0.4621
S	Serine (Ser)	-0.6831
F	Phenylalanine (Phe)	-0.8550
Y	Tyrosine (Tyr)	-0.8983
K	Lysine (Lys)	-0.9389
N	Asparagine (Asn)	-0.9391
I	Isoleucine (Ile)	-0.9558

Ten amino acids correlate positively with GC and ten correlate negatively. Alanine and arginine have the strongest positive correlation with GC, while isoleucine and asparagine have the strongest negative correlation with genomic GC content. (But note that these data apply only to the 109 organisms studied. For the complete list of 109 organisms, see this post.)

If you were to extract all the amino acids out of Clostridium botulinum (28% GC), you would get far more lysine than alanine. Conversely, if you were to hydrolyze all the proteins in Streptomyces griseus (GC 72%), you would find far more alanine than lysine.

Interestingly, serine has six synonymous codons (AGT, AGC, CTA, CTG, CTC, CTT) and can just as easily be specified with G and C as with A and T; so overall, you'd expect little correlation with genomic GC. And yet serine use correlates strongly with low GC. In a sense, this is not surprising. Certain low-GC organisms (like Streptococcus) are known to produce serine-rich cell-coat proteins, some of which are important determinants of pathogenicity. But it may simply be that the high utilization of serine in low-GC organisms is related to one-carbon chemistry. Serine, after all, is the source of the methyl group that, by way of methylenetetrahydrofolate, converts dUMP to TMP (thymidine monophosphate, a DNA precursor). Any organism whose DNA is unusually rich in thymine (low in GC) will almost certainly be processing large quantities of serine. Serine is also a carbon source in the biosynthetic pathways for cysteine and methionine, both of which, like serine itself, are negatively correlated with genomic GC content.

Comparative Genomics Using the Canvas API

An interesting and somewhat mysterious aspect of biodiversity is that the relative proportions of the bases in DNA can take on wildly different values in different organisms even though they're making many of the same proteins. I'm referring to the fact that the G+C (guanine plus cytosine) content of DNA can vary from more than 70% (e.g., Streptomyces species) to less than 30% for certain bacterial endosymbionts (and even for some free-living bacteria, such as Clostridium botulinum). Remarkably, the DNA of the tiny bacterium Buchnera aphidicola (which is distantly related to E. coli but entered into a symbiotic partnership with the aphid around 200 million years ago) has a GC content of only 26%, making its DNA look almost like a two-letter code (A and T, with the occasional G or C).

Many genome-reduced endosymbionts have lost most of their DNA repair enzymes (this is true of mitochondria, incidentally, which are thought to have arisen from a symbiosis with an ancient ancestor of today's Alphabproteobacteria), and this loss of repair capability could well explain much of the GC-to-AT shift seen in endosymbiont genomes. (Left on its own, DNA tends to accumulate 8-oxo-guanine residues, which incorrectly pair with thymine and result in GC-to-TA transversions at DNA replication time.) Whatever the cause(s), GC reduction has left many organisms with severely A- and T-enriched DNA. But the organisms in question are still able to encode perfectly functional proteins, even with a limited nucleic-acid vocabulary.

I decided it might be interesting to try to visualize "GC diversity" by creating a heat map of G and C (in hot colors) and A and T (in cool colors) for certain genes that occur across all bacteria. For example, the following graphic is a heat map of GC and AT usage in the gene for thymidine kinase as it occurs in 61 different organisms with genomic G+C ranging from just over 70% to just under 30%.

Thymidine kinase gene for 61 bacteria of widely varying genomic GC percentages. Guanine and cytosine are represented by hot colors, adenine and thymine by cool colors. Alignments were done with MEGA6 software via ClustalW using a relatively permissive gap-opening penalty of 9 and a gap-extension penalty of 5. The heat map was created with JavaScript using the Canvas API.

To create this map, I obtained DNA sequences for thymidine kinase from 61 organisms (using the excellent online tools at UniProt.org and genomevolution.org), then aligned them in MEGA6 and drew colors (red or red-orange for G or C, blue or blue-green for A or T) corresponding to the bases, using the Canvas API. What you're looking at are 61 rows of data (one row per gene; which is also to say, per organism). Gaps created during alignment are shown in grey.

Several things are apparent from this graphic, aside from the obvious fact that GC usage (indicated by red and orange) tends to be high in the genes for organisms like Brachybacterium faecium (top line, GC 72%) and low for bottom-of-the-chart organisms like Clostridium perfringens B strain ATCC 3626 (28.7% GC) and Ureaplasma urealyticum (25.9%). First, high-GC genes tend to be somewhat longer. (Indeed, in the upper right you can see that the longest genes opened up a sizable alignment gap that extends down the whole graphic.) Also, the genes differ substantially in their leader and trailer sequences, although I think what we're really seeing here is (at least in part) inaccurate annotation of start and stop codons. What's interesting to me is the way certain GC regions trail all the way down to the bottom of the graph while others fade to blue. I think it could be argued that the nucleotides represented by the red blips in the final few lines of the graph, at the bottom, are positions in the gene that are under strong functional constraints. It would be interesting to test those positions for evidence of selection pressure. It could be argued that all areas of low selection pressure have turned blue by the time you get to the bottom of the graph.

I think it's also interesting that several red-orange clusters just to the right of the midpoint, about three-quarters of the way down the graph, all but disappear just as a new red-orange zone begins to appear on the left around 80% of the way down. It's as if certain GC-to-AT mutations in one protein domain led to AT-to-GC mutations in another domain upstream. (The 5' side of the graph is on the left, 3' on the right.)

Getting this many genes (from such divergent sources) to align is not easy. You can do it, though, by setting the gap-open penalty in ClustalW to a low value and aligning genes in small batches, row by row, if need be.

As a technical aside: I first tried creating this graph using SVG (Scalable Vector Graphics), but the burden of creating a separate DOM node for every pixel (which is what it amounts to) was way too much for the browser to handle (Firefox choked, as did Chrome), so I quickly switched to the Canvas API, which puts no heavy DOM burdens on the browser and can convert a FASTA-formatted alignment file to a nice picture in about two seconds.

For what it's worth, here are the names of the organisms whose genes appeared in the above graphic, arranged in order of GC content of the kinase gene only (not whole-genome GC). High-GC organisms are listed first:

Blastococcus saxobsidens strain DD2
Geodermatophilus obscurus strain DSM 43160
Streptomyces cf. griseus strain XylebKG-1
Streptosporangium roseum strain DSM 43021
Kribbella flavida strain DSM 17836
Brachybacterium faecium strain DSM 4810
Deinococcus radiodurans strain R1
Gordonia bronchialis strain DSM 43247
Mesorhizobium australicum strain WSM2073
Turneriella parva strain DSM 21527
Mesorhizobium ciceri biovar biserrulae strain WSM1271
Propionibacterium acnes TypeIA2 strain P.acn33
Novosphingobium aromaticivorans strain DSM 12444
Geobacillus kaustophilus strain HTA426
Geobacillus thermoleovorans strain CCB_US3_UF5
Halogeometricum borinquense DSM 11551
Alistipes finegoldii strain DSM 17242
Agrobacterium radiobacter strain K84
Rhizobium tropici strain CIAT 899
Aeromonas hydrophila strain ML09-119
Rhodopirellula baltica SH strain 1
Porphyromonas gingivalis strain ATCC 33277
Dyadobacter fermentans strain DSM 18053
Enterobacter cloacae strain SCF1
Paenibacillus polymyxa strain M1
Parabacteroides distasonis strain ATCC 8503
Bacillus amyloliquefaciens strain Y2
Vibrio cholerae strain BX 330286
Erwinia amylovora strain ATCC 49946
Bacillus subtilis BEST7613 strain PCC 6803
Gramella forsetii strain KT0803
Prevotella copri strain DSM 18205
Anaerolinea thermophila strain UNI-1
Klebsiella oxytoca strain 10-5243
Bacteroides ovatus strain 3_8_47FAA
Aggregatibacter phage S1249
Escherichia coli B strain REL606
Shigella boydii strain Sb227
Psychroflexus torquis strain ATCC 700755
Bacteroides dorei strain 5_1_36/D4
Leuconostoc gasicomitatum LMG 18811 strain type LMG 18811
Mycoplasma gallisepticum strain F
Myroides odoratimimus strain CIP 101113
Oceanobacillus kimchii strain X50
Chryseobacterium gleum strain ATCC 35910
Carboxydothermus hydrogenoformans strain Z-2901
Yersinia pestis D106004
Bacillus thuringiensis serovar andalousiensis strain BGSC 4AW1
Bacillus anthracis strain CDC 684
Coprobacillus sp. strain 8_2_54BFAA
Proteus mirabilis strain HI4320
Staphylococcus aureus strain 04-02981
Aerococcus urinae strain ACS-120-V-Col10a
Lactobacillus acidophilus strain 30SC
Lactobacillus reuteri strain MM4-1A
Lactococcus lactis subsp. cremoris strain A76
Haemophilus ducreyi strain 35000HP
Ureaplasma urealyticum serovar 5 strain ATCC 27817
Clostridium botulinum A strain Hall
Streptococcus agalactiae strain 2603V/R
Clostridium perfringens B strain ATCC 3626 

An Evolutionary Debate that Misses the Point

One of many hotly debated topics in evolutionary biology is how codon usage bias (preference of an organism for certain codons, when other, synonymous variants are available) relates to transfer-RNA abundance. It's clear the two are related; no one disagrees on that. The question is whether codon usage bias is an outcome of tRNA abundance ratios, or the reverse.

What I think most people are missing in this argument is that the whole discussion might very well be mooted by a huge factor in tRNA evolution that no one seems to be taking into account. I'm talking about the fact that tRNAs are insertion targets for various kinds of mobile genetic elements, from phages to plasmid-borne genomic islands to transposons. As Jörg Hacker and Elisabeth Carniel point out in "Ecological fitness, genomic islands and bacterial pathogenicity" (EMBO Reports, 2001):

Genomic islands are part of the flexible bacterial gene pool and are somewhere between 10 and 100 kilobases (kb) in length. They frequently harbor phage- and/or plasmid-derived sequences, including transfer genes or integrases and IS elements. These particular blocks of DNA are most often inserted into tRNA genes and may be unstable.

(Emphasis added.) Transfer RNAs are constantly being "inserted into" (and next to, not always into) by mobile elements, a phenomenon that's been well studied not only in bacteria but in yeast and elsewhere. Over evolutionary timespans, tRNA genes are duplicated, then disrupted, over and over again, by mobile DNA elements. These elements (whether from phages, viruses, transposons, or what have you) are known to have played (and continue to play) a significant role in shaping genome diversity, across all taxa. This is not a trivial factor, in other words. Transfer RNA genes are insertion hotspots. Surely the patterns of tRNA disruption caused by gene-hopping, over the eons, cannot be unimportant in the determination of codon usage patterns.

Many examples can be found of ancient tRNA signatures inside the tail ends of protein genes, no doubt leftovers from millions of years of insertion events.

Shannon Entropy and DNA

Claude Shannon made an important finding when he realized that the information contribution of a symbol could be estimated very simply as -f(x) log(f(x)), where f(x) is the frequency of occurrence of the symbol x. For example, a series of a coin tosses can be considered a binary information stream with symbol values H and T (heads and tails). If the frequency of H is 0.5 and f(T) is 0.5, the entropy E, in bits per toss, is -0.5 times log (base 2) 0.5 for heads, and a similar value for tails. The values add up (in this case) to 1.0. The intuitive meaning of 1.0 (the Shannon entropy) is that a single coin toss conveys 1.0 bit of information. Contrast this with the situation that prevails when using a "weighted" or unfair penny that lands heads-up 70% of the time. We know intuitively that tossing such a coin will produce less information because we can predict the outcome (heads), to a degree. Something that's predictable is uninformative. Shannon's equation gives -0.7 times log(0.7) = 0.3602 for heads and -0.3 * log (0.3) = 0.5211 for tails, for an entropy of 0.8813 bits per toss. In this case we can say that a toss is 11.87% redundant.

Claude Shannon

DNA is an information stream resembling a series of four-sided-coin tosses, where the "coin" can land with values of A, T, G, or C. In some organisms, the four bases occur with equal rates (25% each), in which case the DNA has a Shannon entropy of 2.0 bits per base (which makes sense, in that a base can encode one of 2² possible values). But what about organisms in which the bases occur with unequal frequencies? For example, we know that many organisms have DNA with G+C content quite a bit less than (or in some cases more than) 50%. The information content of the DNA will be less than 2 bits per base in such cases.

As an example, let's take Clostridium botulinum (the source of "Botox" serum), a soil bacterium with unusually low G+C content, at 28%. If we go through the organism's 3,404 protein-coding genes, we find actual base contents of:

A    0.40189
T    0.30603
G    0.18255
C    0.10840

These numbers are for a single strand (the coding strand or "message" strand) of DNA, which is why A and T aren't equal. For whole DNA, of course, A =T and G = C, but that's not the case here. We're just interested in the message strand.

If we put the above base frequencies into the Shannon equation, we come up with a value of 1.8467 for the information content (in bits) of one base of C. botulinum DNA. The DNA is about 7.67% redundant on a zero-order entropy basis. The DNA may be over 70% A and T, but it's a long way from being a two-base (one bit) information stream. Each base encodes an average of 1.8467 bits of information, which is a surprising amount (surprisingly close to 2.0) for such a skewed alphabet.

Chiggers, Scrub Tyhpus, and Pseudogenes

If you've ever been bitten by tiny red bugs in the garden, you're familiar with members of the Trombiculidae, a family of mites known variously as berry bugs, harvest mites, red bugs, scrub-itch mites, aoutas, or (in the southern U.S.) "chiggers."

In the United States., the garden-variety chigger is basically harmless, but in much of the world this tiny arthropod comes with a very nasty endosymbiont known as Orientia tsutsugamushi, which is a bacterium related to the Rickettsia organisms that cause various tick-borne diseases. Throughout much of the Orient, O. tsutsugamushi infections (from chigger bites) cause scrub typhus, which begins with a rash and fever but can progress to a cough, intestinal distress, swelling of the spleen, abnormal liver chemistry, and ultimately pneumonitis, encephalitis, and/or myocarditis and even death. Treatment with doxycycline, azithromycin, or chloramphenicol is usually successful.

The "harvest mite" (chigger) can carry scrub
typhus, although U.S varieties are typically harmless.

The sequenced genome for O. tsutsugamushi is available, and if you go to this link and click on "Click for features" at the bottom of the Dataset Information box you should be able to open up a table that shows the organism as having 1,182 protein-coding genes (quite a small number), plus an additional 1,994 pseudogenes (quite a huge number, by comparison). The "DNA Seqs" links in the table will let you download the DNA sequences of all the organism's genes and pseudogenes.

This is an extremely unusual situation, in that we're dealing with a bacterium that has more pseudogenes (switched-off, defunct, damaged genes) than regular genes, something that can be said of no other bacterium of which I'm aware. The leprosy bacterium (Mycobacterium leprae) is famed for having approximately 1100 pseudogenes and 1604 "normal" genes. Astonishingly, Orientia tsutsugamushi reverses that ratio, and then some.

We don't know for sure how old Orientia tsutsugamushi's pseudogenes are. A standard rule of thumb in biology is that microbial genomes experience one spontaneous mutation per chromosome per 300 generations. But this doesn't really help us decide how old Orientia's pseudogenes are, since the pseudogenes probably didn't arise one by one, indepedently, through accumulation of random mutations. More than likely, a massive pseudogenization event caused the simultaneous deactivation of a large, unknown number of the organism's genes (of which 1100 survive today as pseudogenes), much the same as has been hypothesized for M. leprae. We have good reason to believe M. leprae's pseudogenes are at least 9 million years old. It seems likely that the pseudogenes in Orientia are also quite old, or at least not terribly new.

To get more perspective on this, I analyzed Orientia's pseudogenes from a couple of perspectives. What I found, first of all, is that the pseudogenes are shorter than their non-pseudo counterparts, averaging 700 bases in length (versus 879 for normal genes). This is similar to the case with M. leprae (where pseudogenes are 795 bases long and normal genes average 1,098). The average shorter gene length for Orientia vis-a-vis M. leprae is consistent with the fact that this is a greatly gene-reduced low-GC (30.5%) endosymbiont, whereas the Mycobacterium family is (in theory) free-living, with higher GC content (57.8% for M. leprae; 65% or more for tuberculosis species).

I've written before about the fact that in most genes, in most organisms, codons tend to begin with a purine base. Therefore I decided to look at purine usage in base one of normal-gene codons versus pseudogene codons (pseudocodons?), finding the following distribution in normal genes:

Purine usage in base one of codons in Orientia tsutsugamushi (N=346,326 codons). No pseudogenes were included in this graph. See the next graph (below) for pseudogenes.

This graph leaves little doubt that most codons begin with a purine (A or G). The median AG1 value is 63.8%. Very few proteins lie to the left of x=0.50, and frankly some of those are probably misannotated as to reading frame.

The situation with pseudogenes is quite a bit different:

Purines in codon base one (AG1) of pseudogenes (N=462,933 codons) in Orientia.

Here we see that purine usage in codon base one is not as strong (median 58.4%), although clearly, plenty of codons still show AG1 above 60%, implying that many pseudogenes are still "in frame" (not frameshifted).

Interestingly, AG1 is not only higher in normal-gene codons than in pseudogene codons, it's also higher in codons associated with proteins of known function than for "hypothetical protein" genes. Only 41.3% of pseudogene codons have AG1 greater than 60%, whereas 66.7% of "hypothetical protein" genes have AG1 > 60% and 84.3% of genes with functional assignments have codon AG1 greater than 60%. This implies that some genes annotated as hypothetical proteins may, in reality, be pseudogenes that are incorrectly annotated. I'll return to that topic some other time.

Looking for LUCA

In 1964, Emile Zuckerkandl and Linus Pauling wrote a paper (published the following year) for the Journal of Theoretical Biology suggesting the use of amino-acid and nucleic-acid sequences for deducing phylogenetic relationships. Ever since then, biologists have been trying to use sequence data to get to the root of the tree of life. Darwinian logic says that at some point, all cells had to have diverged from a Last Universal Common Ancestor (LUCA). Unfortunately, as pointed out by Doolittle and others, the quest for LUCA is greatly complicated by mutational saturation effects, reductive genome loss in important members of the most ancient taxa, convergent evolution, and non-negligible (yet difficult to estimate) amounts of horizontal gene transfer, among other serious problems.

An evolutionary tree of life based on analysis of N=420 genomes of free-living organisms. Proteomes are taxa and protein fold superfamilies are character data. Adapted from Kim and Caetano-Anollés, BMC Evolutionary Biology (2011), 11:140. Click to enlarge. See text for discussion.

The difficulty (I won't say folly) of trying to construct a well-rooted tree of life is made evident in various failed attempts to trace common descent via protein sequences. In January 2010, a few months after the sequencing of the 1000th bacterial genome, Karin Lagesen, Dave W. Ussery, and Trudy M. Wassenaar published a paper in which they expressed surprise over the fact that when they looked all 1,000 then-existing genomes (the number is now more than 10 times that), they could not find a single protein that was conserved across all bacteria. (Here, "conserved" means >50% amino-acid sequence identity.) Harris et al. took a slightly different approach, using the Clusters of Orthologous Groups (COG) database to search for universally conserved genes that follow the same phylogenetic patterns as ribosomal RNA (and therefore might constitute the ancestral genetic core of today's cells). The upshot:

Of the roughly 3100 COGs analyzed, only 80 were found to occur in all organisms. Fifty of these universally present genes showed the same phylogenetic relationships as rRNA.

Harris et al. found that the majority of universally conserved three-domain COG genes (37 of 50) are physically associated with the ribosome. Surprisingly, they found that "relatively few genes encoding proteins involved in DNA replication or transcription from DNA to RNA proved to be three-domain." In particular, RNA polymerases (except for certain subunits) did not follow rRNA distribution patterns and are not conserved across the three domains of life (archaea, bacteria, and eukaryotes). Moreover, the only component of the replicative DNA polymerase in modern cells that was found to be conserved across domains was DnaN (COG0592), the gene for the “sliding clamp.”

These disappointing results are understandable and perhaps expected, given the huge amount of deck-reshuffling that's happened in three billion years. It might well be that genome sequence data, with its constant churn, represents the wrong level of granularity for deep-phylogenetic studies. What matters for organisms, after all, is function, and function is an outcome of protein tertiary structure, not just primary structure.

With that in mind, Kyung Mo Kim and Gustavo Caetano-Anollés in 2011 published a brilliant study in BMC Evolutionary Biology called "The proteomic complexity and rise of the primordial ancestor of diversified life," relying on major structural motifs as the unit of phylogenetic discrimination. Defining protein domains at the highly conserved fold superfamily (FSF) level of structure, Kim and Caetano-Anollés used an iterative, parsimony-based phylogenomic approach to reconstructing FSF repertoires as upper and lower bounds of a presumed urancestral proteome ("ur" here meaning universal). Their conclusion:

The minimum urancestral FSF set reveals the urancestor had advanced metabolic capabilities, was especially rich in nucleotide metabolism enzymes, had pathways for the biosynthesis of membrane sn1,2 glycerol ester and ether lipids, and had crucial elements of translation, including a primordial ribosome with protein synthesis capabilities. It lacked however fundamental functions, including transcription, processes for extracellular communication, and enzymes for deoxyribonucleotide synthesis. Proteomic history reveals the urancestor is closer to a simple progenote organism but harbors a rather complex set of modern molecular functions.

The paper is quite long (14,700 words) and often relentlessly technical, but convincingly restores the quest for LUCA to the firm empirical grounding that such a quest seemed (for a while) to have been robbed of after Doolittle's "Uprooting the Tree of Life" and Dagan and Martin's "The Tree of One Percent." 

While parasitic microorganisms were found to occupy some of the most ancient branches of the superkingdom tree, Kim and Caetano-Anollés nevertheless decided to omit such organisms from their study since reductive evolution (wholesale loss of entire families of enzymes and their control systems) might otherwise queer the results. The final set of free-living organisms included 48 archaeal, 239 bacterial, and 133 eukaryotic members. To avoid potential problems with long-branch attraction, the researchers wisely sampled (at random) equal numbers of proteomes per superkingdom and replicated trees of proteomes, so that bacterial data (which of course predominated) wouldn't swamp archaea or eukaryota.

Among the many fascinating findings in the study:

• The earliest start of organismal diversification occurred sometime between 2.91 and 2.03 billion years ago.
• Translation had metabolic origins. It appeared only after the emergence "of a large number of metabolic functions, but before enzymes necessary for the synthesis of DNA."
• Proteomic analysis of extant fold superfamilies (FSFs) showed that "over 200 additional FSFs are necessary in urancestral FSF sets to account for the complexity of the simplest organism in existence today."
• None of the domains present in ribonucleotide reductase (RDR) enzymes was present in the min_set (representing the LUCA lower bound of complexity). Further, "We note that the reduction of ribonucleotides to deoxyribonucleotides involves the production of an active site thiyl radical that requires contacts with cysteines in all protein domains of the catalytic subunit of the oligomeric enzymatic complex, suggesting modern ribonucleotide reductase functions is [sic] indeed derived."
• Commenting on the known active-site domain homology between class III ribonucleotide reductase and pyruvate formate lyase (a link proposed to have mediated the RNA-to-DNA biological transition), Kim and Caetano-Anollés point out that phylogenomic analysis at the fold-family level suggests the pyruvate formate-lyase domain emerged later than its ribonucleotide reductase counterpart. Therefore it's likely that the urancestor stored genetic information as RNA and not DNA.

Kim and Caetano-Anollés note: "The urancestor had an advanced metabolic network, especially rich in nucleotide metabolism enzymes, had primordial pathways for the biosynthesis of membrane glycerol ether and ester lipids, crucial elements of translation, including amino-acyl tRNA synthases, regulatory factors, and a primordial ribosome with protein synthesis capabilities. It lacked however transcription and in advanced evolutionary stages stored genetic information in RNA (not DNA) molecules."

The authors have many interesting things to say about the evolution of archaeal and bacterial membrane-lipid chemistry (and much else). If you're a biologist and you haven't yet read the Kim and Caetano-Anollés paper, do yourself  favor and take a look at it now. It's a fascinating read, no matter what side of the LUCA fence you're on.

Which Direction Does the Gene Point?

A maddening problem in genome annotation is determining the "sense" strand for a gene, especially when the gene is short and/or the genome has a high GC content (and thus contains few or no stop codons in reverse translation). To convince yourself this is a very real and serious problem, all you have to do is browse a few genomes to see the ridiculously high number of "hypothetical proteins" (over 40% in some genomes), bogus overlaps, genes that score BLAST hits in reverse-complement mode but not frame zero, and other artifacts that are a direct result of the aforesaid problem.

I've presented examples of this problem before, but just so there's no confusion, I want to show you a particularly maddening example so you can see what I'm talking about. (Then I'll suggest a solution.) The following graphic shows a region of E. coli UTI89 in which several genes are shown as overlapping (that is to say, existing on opposite strands of DNA in the same coverage area). Small overlaps sometimes happen between genes, but whole genes rarely, if ever, overlap, and never in clusters. The situation shown below is bogus, but you see it all the time in public genomes. In fact, some of the genes shown below also show up as overlaps in Mycobacterium abscessus M93 (see gene OUW_18941), Citrobacter koseri strain ATCC BAA-895 (gene CKO_00072), and quite a few others. Glimmer has choked on this exact situation many times, in many genomes.  

A region with overlapping genes in the genome of E. coli UTI89.

The big gene on the top strand, middle, is UTI89_C4288 (DNA sequence here). It's annotated as (what else?) a "hypothetical protein." The M. abscessus version of the gene (here) is marked as a "cellobiose phosphorylase," and you can find many BLAST hits (the Rothia version gives an E-value of 6.0×10^-101) for similar "cellobiose phosphorylase" genes in other organisms at UniProt.org and elsewhere. Of course, they're all bogus and represent Glimmer choke points, but the question is how one can determine that, and be sure about it.

E. coli's hypothetical protein (UTI89_C4288) has a wobble-base GC percentage of 64.3%, whereas the gene on the opposite strand (just below it, pointing left), namely UTI89_C4287 (marked as "membrane-bound ATP synthase F1 sector alpha-subunit"), has GC3 = 54.5%. In a much higher-GC organism like Mycobacterium or Pseudomonas, you would find out which gene has the higher GC3 percentage and crown it the winner (and most of the time, you'd be right). In this case, it's not so simple. The gene with the higher GC3 value isn't necessarily the winner.

Of course, in this particular example, you can cheat and look at the identities of the genes in the immediate vicinity of the hypothetical protein, on the bottom strand, and if you do, you'll find that all of the bottom-strand genes are ATPase subunits. Mystery solved, right? Sure, in this particular case. But what about situations where overlapping genes are all shown as "hypothetical protein"? (You can find many such cases in the genome for Burkholderia pseudomallei strain 1710b, for example.) When a hypothetical overlaps a hypothetical in a low-GC genome, then what?

One of my favorite cheats (but this isn't the final solution!) is to check the gene's AG1 percentage (adenine plus guanine, codon base one). This percentage averages ~60% in something like 90% of protein-coding genes. The problem is, AG1 is often 60% whether you read the gene forward, or backward (off the antisense strand). The reverse complement of a gene usually has high AG1, because the forward AG3 is usually under 50%.

Almost any trick you can dream up will fail under edge cases. GC3 is helpful, but only in high-GC genomes. AG1 is helpful, but only sometimes. Shine Dalgarno signals are not universally used by all organisms, and even in those that do use them, they're usually reserved for highly conserved genes encoding things like ribosomal proteins. Gene context is helpful in some cases but not others.

It turns out, the best clue for positively identifying the correct strand and correct reading frame is codon usage frequency patterns. If you know what the codon frequencies are, genome-wide, for a given organism, you can use this information to good advantage, even if the genome (and therefore the codon table) contains inaccuracies. As long as the codon frequencies are approximately correct, you can use them to verify the reading frame of a protein-coding gene.

The algorithm I came up with is very simple, yet effective. For a given gene, read each triplet of bases sequentially, and score each triplet twice: keep two scores going. First, score it according to its frequency in the codon table for the organism. Then score it according to a second table developed for reverse-complement codons.

The following table shows codon frequencies in Caulobacter crescentus NA1000. If you were to encounter a "hypothetical protein" gene in Caulobacter, and you couldn't decide whether the strand assignment was correct or not, first develop a score for the gene by reading its triplets and adding the frequency value of the corresponding codon to the running total. For example, if you encounter the triplet "CTG," add 6.84 to the score (see table). For every occurrence of CTT, add 0.60, for CTC add 1.70, and so on, using the values in the table.

Codon frequencies for Caulobacter crescentus NA1000.

But you also have to create an anticodon frequency table as follows: For every codon in the original table, apply the same score to the corresponding reverse-complement codon in the "antcodon table." E.g., for CTG, the first table would contain 6.84 (as above), but the value 6.84 would apply to CAG (the reverse complement of CTG) in the second table. I call the first table the "forward" table and the second table the "back" table. One represents the frequencies of codons encountered in protein genes in the forward reading direction. The other represents those same frequencies applied to the reverse-complement of the codons (the same codons read in the reverse direction, off the opposite DNA strand).

When scoring an unknown gene, you tally a "forward table" score, and keep a separate score using the "back" table. When you're done, the gene's "forward table" score should be greater than the "back table" score. If it's not, you're reading the gene off the wrong strand.

When I scored all 3,737 C. crescentus CB15 genes using this technique, I found 136 genes that gave a "back" score higher than the "forward" score. Interestingly, when I checked the identities of those 136 putative "backwards-annotated" genes, 132 of them were listed as "hypothetical proteins." Only four genes with assigned functions gave suspect scores, and one of those (CC_0662) turns out to be a 100%-identity match for the reverse of gene CCNA_00700 in Caulobacter crescenstus NA1000.  The other three are less than 200 bases long and could well be non-coding regions.

For a more challenging test, I turned to the genome of Rothia mucilaginosa DY-18, one of the most disastrous annotation nightmares of all time. In the genome for DY-18 you will find 524 protein-coding genes (out of 1,905 total) that are involved in significant overlaps. (Some overlaps are 2-on-1, some are 1-on-1; but the genome is almost certainly overannotated by at least 260 genes.) I trained my program on the codon usage table of R. mucilaginosa M508 (which contains fewer overlaps than DY-18), then tallied codon and anticodon scores on all of DY-18's CDS genes. In the end, 276 genes gave scores indicative of a reversed reading frame. Of those, 265 were, in fact, involved in overlaps.

Codon scoring is such an effective method, I don't know why programs like Glimmer don't use it. It's quite obvious they're not using it, though, because every genome has reverse-annotated genes (by the hundreds, in some cases) that are easily detected using this simple method.

Here, for the record, are the Caulobacter crescentus CB15 genes that appear to be annotated on the wrong strand:

CC_0023
CC_0048
CC_0073
CC_0099
CC_0149
CC_0354
CC_0480
CC_0546
CC_0564
CC_0605
CC_0662
CC_0666
CC_0676
CC_0677
CC_0680
CC_0681
CC_0687
CC_0728
CC_0739
CC_0775
CC_0782
CC_0786
CC_0825
CC_0850
CC_0853
CC_0913
CC_0987
CC_0996
CC_0997
CC_1020
CC_1022
CC_1031
CC_1032
CC_1050
CC_1069
CC_1073
CC_1084
CC_1094
CC_1123
CC_1127
CC_1161
CC_1174
CC_1212
CC_1222
CC_1238
CC_1245
CC_1274
CC_1312
CC_1322
CC_1340
CC_1349
CC_1392
CC_1393
CC_1394
CC_1395
CC_1414
CC_1416
CC_1513
CC_1561
CC_1648
CC_1789
CC_1793
CC_2000
CC_2086
CC_2116
CC_2163
CC_2184
CC_2193
CC_2240
CC_2256
CC_2308
CC_2334
CC_2338
CC_2351
CC_2376
CC_2413
CC_2424
CC_2442
CC_2445
CC_2450
CC_2452
CC_2471
CC_2475
CC_2499
CC_2519
CC_2525
CC_2571
CC_2574
CC_2597
CC_2602
CC_2621
CC_2624
CC_2665
CC_2698
CC_2705
CC_2718
CC_2719
CC_2720
CC_2731
CC_2732
CC_2738
CC_2739
CC_2756
CC_2769
CC_2800
CC_2850
CC_2865
CC_2875
CC_2878
CC_2907
CC_2916
CC_2949
CC_3050
CC_3055
CC_3251
CC_3302
CC_3318
CC_3342
CC_3360
CC_3429
CC_3437
CC_3438
CC_3451
CC_3453
CC_3463
CC_3479
CC_3517
CC_3519
CC_3547
CC_3548
CC_3553
CC_3554
CC_3608
CC_3665
CC_3671
CC_3700