extfieldsofstudy: Chemistry Medicine

text:
Physico-Chemical Studies on the Evaluation of the Antioxidant Activity of Herbal Extracts and Active Principles of Some Indian Medicinal Plants

Understanding of the efficacy and mechanism for the reaction of the biologically important radicals with natural and/or synthetic antioxidants is the first step towards the development of future therapeutic agents. The kinetic parameters e.g., formation and decay rate constants predict the efficacy of an antioxidant and its fate after reaction. These parameters also dictate the ease with which competing reactions would occur in a bio-environment. The spectroscopic parameters provide the clue to the site of free radical attack to these antioxidants. Here, in this article an attempt has been made to show the use of physico-chemical methods in the evaluation of antioxidant activity of some important medicinal plants commonly used in India and the subcontinent. The systems chosen here for discussions are herbal extracts as such, curcumin from turmeric, methoxy phenols from Indian spices, dehydrogingerdione from ginger and bakuchiol from Psoralea corylifolia. All the examples shown in this article illustrate the potential of the pulse radiolysis coupled with kinetic spectroscopy and other physicochemical techniques for the study of antioxidants either in the form of mixture as in herbal extract or as an isolated compound.

Introduction
Oxidation of biological molecules has been postulated to induce a variety of pathological events such as atherogenesis, carcinogenesis and ageing [1]. It has been proven beyond doubt that these damaging events are caused by free radicals [2]. The most important free radicals identified to induce such oxidative damage are commonly termed as reactive oxygen species (ROS) such as ˙OH, HO2, O2 −˙, H2O2 etc and reactive nitrogen species (RNS) such as NO˙, NO2˙, ONOO − etc. Therefore the concept of pharmacological supplementation of a defense against ROS/RNS with antioxidants has become an intensive area of research and here lies the interest in the chemistry of antioxidants. The use of herbal products for improving one's health has augmented in modern time. Indian medicinal plants provide a rich source of such products. Ayurveda and other systems of Indian medicine have identified a number of plants and their components that have curative values.
The measure of the antioxidant properties of the herbal products and their active components is important. Different biochemical as well as chemical methods can be employed *To whom correspondence should be addressed. Tel: 91-  Fax: 91-  for these studies. In vitro animal model systems such as rat brain and liver mitochondria, various cell culture systems are some of the subjects where these studies are carried out. The methodologies that have been used are immunocytochemistry, confocal microscopy and electron microscopy as morphological techniques to name a few. Several standard biochemical and molecular biological techniques are used to assess the damage caused by free radicals and its repair by the antioxidants. All these studies provide gross pictures about the events that take place involving free radicals and antioxidants in in vitro systems. Under this scenario, molecular level studies become indispensable. The physico-chemical parameters like redox potential, reactivity, half-life, chemical transformation and stable end products determine the positive and adverse effects of the test molecules in use. These studies might provide a clear view of the electron transfer processes that occur both in free radical induced damage as well as repair by antioxidants. Advanced techniques such as pulse radiolysis and cyclic voltammetry (CV) have been used extensively to determine redox properties of molecules in solution. The details of pulse radiolysis are given below. CV was extensively used by Chevion et al. [3] for the evaluation of antioxidant capacity of blood plasma, tissue homogenates and extracts of edible plants. Also some cyclic voltammograms were reported for beer samples as an example of a food product containing polyphenolic compounds [4]. It is also imperative to understand the mechanism involved in the antioxidant activity of such compounds. Kinetic spectroscopy in combination with pulse radiolysis is generally being used for such studies. In this review we will focus on some recent results obtained in our laboratory. Mostly kinetic, spectroscopic and electrochemical methods were used to determine the antioxidant efficacy and mechanism of selected herbal extracts and their active ingredients. It must be noted that no attempt at a comprehensive survey of literature is made, nor is any claim made that the literature citations given are fully representative.

Pulse Radiolysis Setup and Measurement
In the pulse radiolysis setup used in this group at Bhabha Atomic Research Center [5], 7 MeV electrons come out from the waveguide through a titanium window with mean diameter of 2 mm. The machine is capable of producing electron pulses of widths between 50 ns to 2 µs with peak current ranging from 400 mA to 70 mA. The block diagram of the machine is shown in Fig. 1. In this technique free radicals are generated by employing electron pulse and the fate of these radicals are followed with time in absence and presence of scavengers. Pulse radiolysis is the most convenient method to generate detectable instantaneous concentrations of free radicals. Ever since its discovery, this technique has been used to study reactions of biological molecules. The simple reason is that one can mimic to some extent biological situations by studying the free radical induced reactions of bio-molecules in solutions. The intermediates, which are formed in vivo can be mimicked by studying redox reactions of biological molecule in vitro and can be assigned by the method like transient optical absorption technique in the present case. By knowing the intermediates, mechanism of the reaction can be delineated.
Accurate measurements of kinetic parameters require knowledge of the concentrations of the reactants and the rate constants for the reactions so that importance of the potential competing reactions can be understood. The rate constants are determined by two methods, one is direct determination and the other is by competition kinetic method. In direct determination the principle of pseudo first order kinetics is generally followed. The substrate is used in large excess to scavenge all the free radicals (R˙). The equations follow:

Studies with Herbal Extracts
Now a days plant extracts and herbal formulations, are gaining importance as less toxic antioxidants and radioprotectors [6,7]. While there are several attempts to study the bio-chemical and pharmacological activities, there is a need to develop methods for evaluating the free radical scavenging ability of such extracts. Extracts contain mixture of compounds, with unknown compositions and also there is no direct method to estimate their concentration. In such case, a general practice is to express total concentration of the extract as µg/ml. Because of the multi-constituent mixture, it is not possible to determine the absolute rate constants for the reactions of extracts with free radicals. However, it is possible to estimate their relative reactivity in comparison with a standard, which can indirectly be used to indicate the total scavenging ability of the extract. Reactivity is a product of rate constant and the concentration and is a quantity, which can be directly compared when mixtures participate in a chemical or biological reaction. Some of the methods employed for studying free radical reactions are given below.
For the determination of hydroxyl radical (˙OH) scavenging ability, potassium thiocyanate (KSCN or SCN − ) is used as standard using pulse radiolysis. In this method ˙OH is allowed to react with known concentration of KSCN at pH 7, in the absence and in the presence of various concentrations of extracts. In the absence of extract, the ˙OH reacts completely with SCN − to produce (SCN)2˙−. The formation of (SCN)2˙− can be observed by monitoring its transient absorbance in microsecond time scales at 500 nm [8]. In the presence of extract, capable of reacting with ˙OH, the absorbance due to (SCN)2˙− at 500 nm decreases. From the extent of decrease in the absorbance at 500 nm, rate of scavenging of ˙OH by extract (E) can be calculated. For these studies with ˙OH radicals, it is necessary to use advanced techniques like pulse radiolysis, however it is possible to estimate the free radical scavenging ability using other radicals like 2,2'-diphenyl-1-picrylhydrazyl (DPPH) and 2,2'-azinobis(3-thylbenzthiazoline-6-sulfonic acid (ABTS˙−). DPPH is purple in colour (λmax 517 nm) and in presence of a substance capable of donating an electron or a hydrogen atom, it loses its radical character and becomes colourless [9]. The total reactivity of the extract towards DPPH radical can be studied by following the kinetics of DPPH radical reaction in presence of the extracts. In general such reactions take place in time scales of fraction of a second to several hundred seconds and stopped-flow kinetic spectrometer is best suited for this. ABTS˙− radical is also a stable radical like DPPH, shows absorption in the wavelength region 400 nm to 750 nm, with absorption maxima at 417 nm (ε = 3.47 × 10 4 M −1 cm −1 ), 645 nm (ε = 1.35 × 10 4 M −1 cm −1 ) and at 728 nm (ε = 1.5 × 10 4 M −1 cm −1 ) [10]. The reaction of ABTS˙− radical with extracts can be followed by steady state or by time resolved methods. However, unlike DPPH, ABTS˙− radical reactions involve electron transfer and take place at a much faster rate and hence depending on the time scale of reaction, it can be monitored by either pulse radiolysis or stopped-flow spectrometer.
Employing the above techniques, we were able to study a number of extracts such as Phyllanthus emblica (T1), Terminalia chebula (T2), Terminalia belerica (T3) and triphala (T4) an equi-proportional mixture of fruits of the three herbs [11][12][13][14][15]. Table 1 gives the ascorbate equivalents estimated by this method in different extracts studied in our laboratory. Further we were able to compare these results with their antioxidant activity like inhibition of lipid peroxidation in liposomes induced by γ-radiation. The IC50 value, i.e. the concentration of the extract required to inhibit lipid peroxidation by 50% was estimated and these are listed in Table 1. The values correlate well within experimental limits except in case of T2, which is being addressed in future. Cyclic votammetric studies were used to obtain the redox potential of the extracts. The reactions of the above extracts with different radicals are possible only when the reduction potential of constituents present in the extract is lower as compared to the reduction potential of the radicals. To determine the reduction potential, cyclic voltammetric traces of the extracts were recorded. None of them showed reversibility but broad peaks were observed at the peak potential values of 0.25 V, 0.3 V, and 0.24 V for T1, T2, and T3, respectively. The broad peaks have been attributed to the presence of large amount of oxidisable substrates present in these extracts. The low redox potential (0.2 to 0.3 V) of T1, T2, and T3 suggest that the extracts contain easily oxidisable substrates that help in neutralizing the oxidizing free radicals via electron transfer reaction. T4 being a mixture showed similar features, with major peak at 0.3 V.

Studies with Curcumin
Turmeric is the root of the perennial herb Curcuma longa, belonging to the ginger family is cultivated in tropical regions of south and southeast Asia. It is commonly used as a spice and coloring agent in Indian cooking. Its dietary consumption is estimated to be ~1.5 g per person per day. Its utility in traditional Indian medicine is referred as old as second millennium BC. It is a household medicine widely used to cure biliary disorders, anorexia, cough, diabetic wounds, hepatic disorders, rheumatism and sinusitis [16][17][18][19]. The most active component of turmeric is curcumin, which makes upto 2-5% of the spice. Several studies have shown curcumin to be a potent antioxidant with its activity 10 times more than that of the well-known antioxidant vitamin E [20]. Recent phase I clinical trials indicate that human beings can tolerate a dose of curcumin as high as 8 g/day with no side effects [16,21]. Curcumin has been shown to be an effective antioxidant under several conditions, where free radicals are generated and since the antioxidant activity is always associated with free radical scavenging activity, curcumin has been examined for free radical scavenging studies.
Using pulse radiolysis, reactions of superoxide, halocarbon peroxyl, lipidperoxyl, methyl and methyl peroxyl radicals, glutathione radicals and azide radicals have been reported [20,22,23]. The rate constants for the reactions of these radicals were determined in aqueous solutions by kinetic spectrophotometry, where the formation of transient at suitable wavelength is followed as a function of time. A transient having a lifetime of a few hundred microseconds absorbing from 250 to 600 nm with maximum absorption at 500 nm was characterized in all the reactions, with extinction coefficient of ~7300 M −1 cm −1 . In aqueous solutions the radicals of curcumin decay by second order kinetics with 2 k values of 7.5 × 10 8 M −1 s −1 . The reaction takes place by the radical-radical recombinations leading to the formation of products, one of them being a two electron oxidized species, which undergoes subsequent cleavage to produce species like, vanillin, ferulic acid etc. The second order radical decay suggests that the decay and lifetime of curcumin radical depend on the initial concentration of radicals (t1/2 = 0.693/2 k [radical]). Under cellular conditions, where very low concentration of radicals is formed (less than micromolar concentration), such reactions between two curcumin radicals do not take place, rather the radicals react with other cellular antioxidants and ultimately may be regenerated. In fact increase in lifetime by an order of magnitude was observed in TX-100 micellar systems.
Using stopped-flow spectrometry the reactions of DPPH radicals and superoxide radicals were studied and the rate constants reported [22,24]. In all these cases decrease in the absorption of either the radical species or the parent curcumin were followed as a function of time. Using transient luminescence detection, the rate constant for the reaction of curcumin with singlet molecular oxygen was determined by following the decay of singlet oxygen emission at 1270 nm. All these rate constants are listed in Table 2. The radicals of curcumin absorbing at 500 nm have been attributed to phenoxyl radicals, as the decay of the radicals was not affected by the presence of oxygen and the radical absorption spectrum is not influenced by the change in pH from 3 to 10 and such radicals acquire resonance stability over the extended conjugation system.
Phenoxyl radicals are formed by one of the two mechanisms such as direct hydrogen abstraction or by initial electron transfer to the free radical, where by a radical cation is formed, which subsequently undergoes a proton loss to produce a phenoxyl radical [25]. Jovanovic et al. also studied the free radical reaction mechanism of curcumin by pulse radiolysis and flashphotolysis and concluded that H atom transfer is the preferred antioxidant mechanism of curcumin [26]. They proposed that due to the availability of two methylenic hydrogens, the preferred site of free radical reaction is the CH2 group in the heptadienone link and not the phenolic OH group, which accounted for only 15% of the reaction. This was sharply contradicted later by many others. Berclay et al. [27], determined the inhibition of controlled initiation of styrene oxidation by curcumin and other synthetic non-phenolic curcuminoids having basic curcumin skeleton and do not possess phenolic groups and confirmed that curcumin is a classical phenolic chainbreaking antioxidant, donating H atoms from the phenolic OH groups and not the CH2 group. Later Jovanovic et al. [28] did some more supporting experiments and corrected their claims. Based on laser flash photolysis studies of curcumin and its mono methylated and trimethylated derivatives, they proposed that the radical formed is that of a phenoxyl radical, which is formed by direct one-electron oxidation of a methoxyphenol moiety or by that the initially formed β-oxo-alkyl radical undergoes an intramolecular Hatom transfer from the methoxyphenol. Later a conclusive evidence has been given by our group by the systematic comparison of the antioxidant activity and the free radical reactions of curcumin and dimethoxy curcumin, in which the two phenolic OH groups are blocked by the alkylation and the diketone moiety is intact, which makes an excellent comparison for this [22]. Three independent experiments were performed and compared for both the compounds. In these experiments, curcumin showed three times more efficacy in inhibiting lipid peroxidation than the dimethoxy derivative. Rate constants for hydrogen abstraction reaction with DPPH radical was 1800 times more than that with the dimethoxy derivative. Pulse radiolysis induced one-electron oxidation of dimethoxy derivative did not produce the characteristic 500 nm absorbing transient instead a week band at 310 nm was observed, that decayed in less than few microseconds. All these studies confirmed that phenolic OH makes major contribution in both antioxidant activity and hydrogen abstraction reactivity, while the methylenic group makes minor contribution for this and the 500 nm band observed in microseconds during one-electron oxidation corresponds to the phenoxyl radical. This was further confirmed by density functional theory (DFT), studies. The bond dissociation energies of the phenolic OH and the methylenic CH2 group in diketo form and that of CH group in the enolic form were estimated and found that the O-H bond dissociation energy of the phenolic group is 5.04 kcal/ mol lower than the C-H bond dissociation energy (BDE) of the central CH2 of the diketo curcumin, suggesting that the phenolic OH is the most vulnerable target for free radicals in curcumin. The site for the free radical reaction of curcumin and the resonant structures of phenoxyl radicals are given in Fig. 2.
From the above studies, it was confirmed that the phenolic site is essential for the antioxidant activity. The diketone moiety of curcumin possesses chelating abilities towards transitions metals. This chelation of curcumin towards redox-active metals like iron and copper has been found to be useful in the treatment of Alzheimer's disease and it is interesting to note that among the Indian populations, regular usage of turmeric is one of the reasons responsible for the reduced risk of Alzheimer's disease. There is a similarity between drugs used for the treatment of Alzheimer's disease and their usefulness as super oxide dismutase (SOD) mimics. This prompted us to develop new copper-curcumin complexes, which showed promising SOD activity and are found to be better antioxidants than curcumin at cellular level [29].

Studies with Other Methoxy Phenols: Correlating In-vitro Activity with Redox Potentials
In another study, four different o-methoxy phenols from Indian spices and plant products, viz., dehydrozingerone, brompentenone, eugenol and isoeugenol were examined for free radical scavenging reactions by pulse radiolysis [30]. In all these cases, one-electron oxidation produced phenoxyl radicals, having absorption spectrum in the range from 300 to 400 nm with absorption maximum at~350 nm. Further, the redox potentials of the compounds to form phenoxyl radicals from their parent phenols at pH 7 were estimated by pulse radiolysis. The reduction potentials are direct measure for their ability to undergo oxidation on reaction with oxidizing free radicals. This in turn is dependent on the O-H bond dissociation energy. This data was compared with the compound in vitro antioxidant ability like inhibition of free radicals induced lipid peroxidation in model systems. The results confirmed that a phenol, which can undergo easy oxidation, is more potent antioxidant. The comparative data is given in Table 3.

Studies with Dehydrogingerdione
Ginger (Zingiber officinale Roscoe) is one of the world's best-known spices, and it is universally used throughout history for its health benefits. The dried extract of ginger contains monoterpenes and sesquiterpenes. The main antioxidant principles in ginger are the gingerols, shogaols and some related phenolic ketone derivatives. Ginger extract possesses antioxidative characteristics, since it can scavenge superoxide anion and hydroxyl radicals. Ginger alone contains up to 12 compounds providing as much as 40-times higher antioxidant activity than vitamin E. The National Cancer Institute has identified ginger as one of the top-10 anticancer foods. As with other botanicals and herbs, ginger also offers numerous other benefits and could easily be placed in the cardiovascular category. Dehydrogingerdione (1), one of the constituents of ginger, possesses a partially methylated catechol moiety, a 1,3-diketo functionality and an extended conjugation (Fig. 3).
In principle, compounds possessing multiple double bonds and especially with active hydrogen can act as radical scavengers via addition to double bonds and/or abstraction of hydrogen atom from the allylic position as is the case with lipids. The positive role of the β-diketone moiety in the antioxidant action of four ginger-derived phenolic diketones (1-4) (Fig. 3) against lipid peroxidation were reported [31] using biologically relevant model systems and pulse radiolysis. In iron-independent peroxidation, however, 1 showed activity comparable to that of curcumin. This was ascribed due to its higher affinity for the lipid peroxyl radical due to higher hydrophobicity. The antioxidant/prooxidant activity of a test compound in iron-mediated lipid peroxidation in a biological system is critically governed by factors such as the concentration of endogenous reducing systems and/or redox couples. Thus, the true antioxidant activity of 1 was assessed by use of a biologically relevant liposomal system, devoid of reductants such as ascorbate, reduced form of  , and by carrying out the Fe 2+ -mediated peroxidation at pH 5.0, at which value Fe 2+ remains unoxidized, for at least 3 h, during which considerable peroxidation took place. This study not only ensured minimization of the prooxidant effect of ascorbate, but also assisted in the quantitative comparison of the antioxidant activities of 1 and ascorbate, individually and in combination. The protection offered by vitamin C (15 µM) alone was considerably less than that offered by 1 even at lower concentrations (5 and 10 µM). More interestingly, the antioxidant activity of 1 increased drastically in the presence of vitamin C at increasing concentrations. These results strongly suggest that 1 scavenges the LOO˙ radical more quickly than vitamin C does, but the radical generated from 1 then reacts with vitamin C to regenerate 1. This established a synergistic effect of vitamin C on the antioxidant activity of 1. The above results were also confirmed by pulse radiolytic studies of the reaction between 1 and the trichloromethylperoxyl radical (CCl3O2˙). In spite of its higher reduction potential value (1.5 V), compared to the physiologically relevant peroxides, the CCl3O2˙ radical is extensively used as a representative peroxyl radical because of the inherent simplicity in performing the experiments, and it was indeed used earlier to study a very important biophysical phenomenon, namely free radical interaction between vitamin E and vitamin C [32]. The radical can be generated in aerated water/isopropanol/acetone mixtures containing carbon tetrachloride. A very high rate constant (2 × 10 9 M −1 s −1 ) for the reaction between 1 and CCl3O2˙ established the superior antioxidant efficacy of 1 in comparison to vitamins C and E. This study also established a synergistic effect of vitamin C on the antioxidant activity of 1. On the basis of pulse radiolysis study, liposomal peroxidation study and HPLC analysis of the products, a mechanism for the antioxidant action of 1 has been proposed, suggesting the contribution of the phenolic group as well as the active methylene group of the 1,3diketones [31]. The pulse radiolysis studies revealed its high reactivity with the hydroxyl and glutathiyl radicals. With the hydroxyl radicals, it initially formed three species, an adduct, a phenoxyl and a methylenic radical. The bimolecular rate constants for the formation of these species were the same (6.8 × 10 9 M −1 s −1 ). At a later stage a phenoxyl radical is produced through an intramolecular transformation from the initial radicals [33]. In the case of glutathiyl-radical-induced oxidation, a carbon-centered radical was inferred from the result of the oxygen effect on the decay of the radical. A suitable mechanism for all the reactions mentioned for the oxidation of 1 has been proposed (Fig. 4).

Studies with Bacuchiol
Besides phenolics, antioxidant activities of other classes of compounds including terpenoids are also reported [34,35] in literature. In this section the mechanism for the antioxidant activity of the meroterpene, bakuchiol, an important constituent of the plant Psoralea corylifolia (P. corylifolia) is discussed. P. corylifolia is an ancient medicinal plant found in India and China and is still being cultivated for medicinal purposes in China today. It is abundant in various parts of India, and its seeds have enjoyed an honoured place in Ayurvedic medicine, particularly for the treatment of inflammatory and skin diseases. Bakuchiol, one of the important constituents of P. corylifolia seeds, is also reported to possess anti-inflammatory and lipid protective activities. The protective activity of bakuchiol against oxidative damages to lipids and proteins has been investigated [36] and rationalized based on the scavenging activity against various oxidizing radicals (CCl3OO˙, LOO˙, hydroxyl and glutathiyl radicals). The rate constants of the scavenging reactions (Table 4), transients formed in these reactions and their mechanistic pathways had been probed using optical pulse radiolysis (PR) technique.
Bakuchiol inhibits microsomal lipid peroxidation in a concentration-dependent manner showing 74.7% protection at a concentration of 10 µM. In view of its solubility in both lipid and water (at higher pH), it is expected to be distributed in both these phases. This may account for its low IC50 = 6.1 ± 0.2 µM value against lipid peroxidation. Thus, the protective effect of 5 against lipid peroxidation was explained on the basis of radical scavenging only. Surprisingly, the methyl ether of bakuchiol (6) also prevented Fenton-mediated peroxidation of rat brain homogenate (RBH), with an IC50 value of 50.2 ± 5.3 µM. The activity of 6 was explained by abstraction of its allylic hydrogen by the lipid peroxyl radical (LOO˙). This was more evident from the subsequent studies on the reactions of different radicals with 5.
Bakuchiol reacted with CCl3O2˙ producing radical with absorption in the visible region. Based on the absorption, formation and decay constants it was proposed that only the phenoxyl radical of bakuchiol was generated. Similarly the azidyl and hydroxyl radicals also produce the phenoxyl radicals.
Studying the dynamics of the reaction between 5 and glutathiyl radical (GS˙) revealed that two species are formed initially, one being the phenoxyl radical other is an allylic radical. Further confirmation of the generation of the allylic radical in 5 was confirmed by studying the reaction of GSẇ ith the protected phenol derivative 6. As expected, 6 produced only one transient, the C-centered radical at the C-5 position of its terpenoid chain. This was further corroborated based on the reaction between glutathiyl radicals and two other congeners, 7 and 8 (Fig. 5). On the basis of strong optical absorption signals with 5-7 as compared to that in the case of 8, formation of the allylic radical adjacent to the trisubstituted olefin function in 5-7 was envisaged. In the case of 5, the allylic radical was found to be transformed into the phenoxyl radical at a later stage. It revealed that 5 not only scavenged GS˙ radicals but might also repair it by abstracting a hydrogen atom from its terpenoid chain. This repair phenomenon is very attractive especially with respect to its biological activity. These results showed the importance of the terpenoid moiety of bakuchiol in controlling its antioxidant action via radical scavenging [36].

Conclusions
Measurements of kinetic parameters and physico-chemical properties can be useful in understanding the radical scavenging efficiency of herbal extracts. This in turn provides an idea about the antioxidant activity of the extracts as such. Even though the extracts are generally mixture of many compounds, rate constant for the reaction of hydroxyl radical with the extract can be determined by competition kinetic method. Such studies help in corroborating the results obtained from same sample in biochemical and chemical assays for antioxidant activities.  kf and kd are the bimolecular rate constant for the formation of the scavengers radical and the decay constant of this radical, respectively.
In the case of known antioxidant like curcumin, methoxy phenols, gingeroids, and bakuchiol a proper mechanism can be delineated based on the kinetic spectroscopic results. Certainly parallel studies on the redox potential, products formed during the reaction and theoretical modeling should substantiate the proposed mechanism. Towards the development of future therapeutic agents as antioxidants and radioprotectors, the methodology described in the article may be very important but preliminary steps to achieve the goal.


extfieldsofstudy: Medicine Biology

text:
Tandem gene arrays in Trypanosoma brucei: Comparative phylogenomic analysis of duplicate sequence variation

Background The genome sequence of the protistan parasite Trypanosoma brucei contains many tandem gene arrays. Gene duplicates are created through tandem duplication and are expressed through polycistronic transcription, suggesting that the primary purpose of long, tandem arrays is to increase gene dosage in an environment where individual gene promoters are absent. This report presents the first account of the tandem gene arrays in the T. brucei genome, employing several related genome sequences to establish how variation is created and removed. Results A systematic survey of tandem gene arrays showed that substantial sequence variation existed across the genome; variation from different regions of an array often produced inconsistent phylogenetic affinities. Phylogenetic relationships of gene duplicates were consistent with concerted evolution being a widespread homogenising force. However, tandem duplicates were not usually identical; therefore, any homogenising effect was coincident with divergence among duplicates. Allelic gene conversion was detected using various criteria and was apparently able to both remove and introduce sequence variation. Tandem arrays containing structural heterogeneity demonstrated how sequence homogenisation and differentiation can occur within a single locus. Conclusion The use of multiple genome sequences in a comparative analysis of tandem gene arrays identified substantial sequence variation among gene duplicates. The distribution of sequence variation is determined by a dynamic balance of conservative and innovative evolutionary forces. Gene trees from various species showed that intraspecific duplicates evolve in concert, perhaps through frequent gene conversion, although this does not prevent sequence divergence, especially where structural heterogeneity physically separates a duplicate from its neighbours. In describing dynamics of sequence variation that have consequences beyond gene dosage, this survey provides a basis for uncovering the hidden functionality within tandem gene arrays in trypanosomatids.

Background
Evolutionary biology has begun to utilise the abundance of genome sequence data, applying a comparative approach to the evolutionary dynamics of genome structure and sequence [1][2][3]. Within this approach genomic position is a novel source of contextual information, and an additional criterion for assessing homology, independent of gene sequences. The debate has recently concentrated on the relative contributions of evolutionary processes affecting sequence divergence, such as gene conversion, positive and negative selection and concerted evolution, and processes affecting gene quantity and order [4,5]. This study concerns the evolution of tandem gene arrays in the kinetoplastid parasite Trypanosoma brucei, the cause of sleeping sickness in humans and various diseases in animals, for which a genome sequence was recently completed [6]. A comparative method, using genome sequence from multiple trypanosomatid species, was applied to tandem-duplicated genes throughout the T. brucei genome to test for various evolutionary forces. The results reveal the variety of dynamics among tandem arrays and clarify the contribution of concerted evolution and gene conversion in regulating sequence diversity among repetitive genes.
Long before the appearance of genome sequences, tandem gene arrays were known to affect genome structure and gene function. Ribosomal rDNA is commonly tandemly-duplicated across the taxonomic spectrum, and is typically thought to be homogenised as a result [7]. In contrast, metazoan hox genes are also tandemly arrayed, but contain substantial variation that demonstrates the functional consequences of the array organisation [8]. Tandem arrays are thought to evolve through a process of replication slippage during DNA synthesis or unequal crossing over during meiosis [9], a phenomenon conceptually consistent with the rapid evolution of mini-satellites and other highly repetitive regions in many genomes through mis-alignment [10]. These processes often result in contiguous gene copies evolving in concert toward a consensus sequence, such that duplicate sequences within tandem arrays may fail to diverge as expected under a neutral model [9]. In these circumstances gene duplicates within a species look more alike than orthologous loci in different species, which is indicative of concerted evolution [11]. Concerted evolution of duplicated genes was first recognised among the rDNA of Xenopus [12] and subsequently shown to be widespread among the globins of primates [13,14] and heat-shock proteins in Drosophila [15]. Furthermore, unequal crossing-over (UCO) and gene conversion (GC) were suggested as mechanisms for this homogenisation within arrays, although their relative contributions were unknown [16], and remain so.
UCO results from mis-alignment of repetitive DNA and occurs between sister chromatids during mitosis and between homologous chromosomes during meiosis [9]; it principally affects gene numbers and can cause both the addition and subtraction of repeat units, as shown in Additional file 1. Consequently, tandem arrays affected by UCO can show intraspecific copy number variation; but while copy number can vary, gene sequences are relatively invariant because UCO adds or subtracts whole repeat units. UCO was implicated in the formal description of concerted evolution among tandem α-globins in primates [18] and is also known to occur in trans, for example, among the five unlinked rDNA loci in humans [19]. The removal of variation by UCO was demonstrated in Streptococci where the functional polymorphism of rDNA tandem copies, which facilitates antibiotic resistance, was abolished when diversifying selection was removed from the bacteria [17]. GC describes the reciprocal or directed exchange of material between alleles on homologous chromosomes or between gene duplicates; these may be at the same locus (allelic gene conversion, AGC) or at different loci (ectopic gene conversion, EGC). Reciprocal exchange between alleles on homologous chromosomes is commonly understood as 'recombination', while non-reciprocal exchanges within or between loci are termed 'biased gene conversion'. These processes are also illustrated in Additional file 1. This study primarily concerns exchanges among gene duplicates, within or between loci, regardless of their reciprocity (although most are biased events). First observed in fungi as a mechanism for homogenisation of alleles [20], GC can homogenise sequences, for instance among mammalian [21] and nematode [22] heat-shock proteins, but also promote polymorphism, for example, in the diversification of both immunoglobulin [23,24] and multi-histocompatibility complex alleles [25,26].
Detecting the action of concerted evolution is essentially about identifying similarities where none are expected; previously, concerted evolution has been identified by eye from obvious, anomalous similarities [27][28][29][30]. Genomic position is central to showing that a lack of divergence was not caused by purifying selection or recent gene duplications. The contextual information provided by genome sequences show that, first, there is a substantial excess of mutations between duplicates in different species, compared to duplicates in the same genome, and second, because gene duplicates display conserved synteny with surrounding loci, they were established long ago. Quantitative methods for the treatment of GC and UCO have been developed in the last three decades (for a review see [31]) and were applied in this study under four criteria: • Cladistic criterion: homologous genes from several organisms are arranged in a phylogenetic tree to estimate their cladistic relationships; concerted evolution is inferred where homologous genes cluster by species, rather than by locus, i.e., paralogs appear more similar than expected given that their closest relatives are found in other genomes.
• Phylogenetic criterion: a 'sliding window' analysis [32,33] is applied to a multiple alignment of homologous genes, with phylogenetic trees estimated for each window. A significant departure in phylogenetic signal between windows identifies two regions with different evolutionary histories, i.e., a recombination 'breakpoint'.
• Distribution criterion [34]: the distribution of silent polymorphisms along a sequence alignment is scored for areas where the distribution is 'condensed', i.e., matches between two otherwise discordant sequences. Where the length of a corresponding region is significantly larger than expected after permutation, a gene conversion event is inferred.
• Compatibility criterion [35]: two homologous genes are compared on a site-by-site basis, breakpoints are inferred where neighbouring, informative sites have incompatible phylogenetic signals.
The importance of concerted evolution has been demonstrated through whole genome comparisons of gene duplicates, across diverse taxa [36][37][38][39], but in most detail for yeast [31,40,41]. Most of these studies have compared duplicated sequences within genomes, but have not provided a comparative context for sequence evolution. In this study, tandem gene arrays were chosen because they offer a structural paradigm to infer evolutionary events, and one peculiar to trypanosomatids since they typically lack individual gene promoters and transcribe polycistronically [42,43,6]. Therefore, the argument that tandem duplication serves to increase gene dosage in the absence of promoters [44] could be assessed. A comparative approach was also specifically chosen, in which variation among T. brucei duplicate sequences was analysed in relation to homologous sequences in five other kinetoplastids, ranging in relation from a subspecies (T. b. gambiense) to another genus (L. major). This provided the historical perspective necessary to identify the origins and causes of sequence variation. Hence, the specific aims were to (i) catalogue all tandem gene arrays in the genome sequence, establish their phylogenetic distribution in other trypanosomatid species and quantify the variation found among the duplicate sequences of each; (ii) assess the relative contributions of concerted evolution, UCO and GC, through the application of multiple tests; and (iii) use this information to characterise the causes of variation (or lack of) in each array.

Results
After surveying the Trypanosoma brucei genome sequence, 47 tandem arrays with at least four gene copies were identified and characterised. This is unlikely to be an exhaustive list of tandem-duplicated loci in the T. brucei genome since some lengthy tandem arrays may be collapsed to tandem pairs or triplets if gene number could not be resolved during genome assembly. A full list of tandemduplicates, regardless of copy number, is available from the author. Duplicate sequences were combined with homoeologous sequences (i.e., homologous sequences occurring at the same genomic position) from five other trypanosomatid species and several criteria were applied to the multiple alignments to identify the molecular signatures of concerted evolution and gene conversion. Each of the methods scrutinised multiple alignments for unexpected similarity between otherwise dissimilar gene sequences, which could not be explained through occasional conservation or derivation. The evidence suggested that tandem arrays display a variety of evolutionary dynamics and often contain substantial sequence variation; among them examples of conserved arrays, evolving in copy number through unequal crossing-over, moderately variable duplicates evolving in concert, those affected by allelic gene conversion, and also structurally segregated arrays, with different duplicates exposed to distinct evolutionary forces.

Variation within tandem arrays
Values for D n and D s in Table 1 show that for many cases there was substantial variation among the duplicate genes within tandem arrays. Indeed, in a few cases D s exceeded 1 (6 cases), probably due to highly divergent copies in these arrays that showed little or no sequence homology over part of the alignment. Typically, tandem arrays showed moderate variation among copies, where D s always exceeded D n , but D n was non-zero (31 cases). However, a sizable minority showed no variation at all (10 cases).
Turning first to the invariant arrays, these included histones (Tb4, 39, 57, 107), tubulin (Tb3), cysteine peptidase (Tb43), glycerol kinase (Tb89), paraflagellar rod proteins (Tb18, 73) and S-adenosylmethionine synthetase (Tb47). The coding sequences of these arrays were identical, except for very occasional silent changes. The non-coding parts were typically identical also, although there was some variation in copy number of repetitive tracts, for instance in the 3' UTR of Tb4 (although the ability to precisely determine tract length is limited). Generally, the terminal UTRs of the arrays were identical to internal IGSs but curtailed; this shortening was assumed to reflect the proportion of the IGS properly assigned as UTR.
Among moderately variable cases, variation was typically continuous among copies and greater in coding regions. Notable examples of this included Tb93 (BARPs), where the corrected genetic distance between coding sequences (CDS) was 0.335 sub/site, while equivalent values for 5' and 3' UTRs were 0.073 and 0.038 respectively. In other arrays, such as Tb24 (S-T phosphatase) and Tb62a, the non-coding regions were invariant, while their associated CDS showed significant divergence. The only instances to show the opposite pattern, where non-coding regions were more variable than CDS were the histone arrays (Tb4, 39 and 107), albeit because the CDSs were invariant.
Tb927.6.1300 X 0.0220 S-adenosylmethionine synthetase X X X 0.0000 0.0000 500 Presence of a homologous locus within conserved gene order is shown by a cross X. † Sequence divergence is expressed using the rate of synonymous (D s ) and non-synonymous (D n ) substitutions per site. # Monophyly of array copies was assessed through non-parametric bootstraps of a ML phylogeny (500 replicates) and posterior probabilities in a Bayesian phylogeny ('na' denotes that monophyly was not recovered by the Bayesian tree). SH tests compared the likelihood of an optimal topology with 'orthologous' and 'concerted' scenarios; significant p-values are shown in bold and indicate that the scenario concerned was significantly less likely that the optimal topology. The cladistic criterion was only applied to arrays with widespread distributions.
Where variation existed, the phylogenetic signal among coding sequences was typically inconsistent with 5' and 3' UTRs. Additional file 2 shows phylogenetic networks for three components of Tb55 (retrotransposon hot-spot protein). In frame A, copies 2-4 are robustly divergent from other duplicates, but this was not repeated with either 5' or 3' UTRs (frames C and D respectively). Not only are copies 2-4 dispersed in these latter networks, but nowhere do the UTRs show the same pattern at CDSs. Hence, these regions of the repeat unit showed different histories and evolutionary rates. Similarly, the relationships between Tb61 CDSs, shown in Additional file 3, are not reflected in the affinity of UTRs, of which there are three distinct and unrelated types. The exceptions to this inconsistency were those arrays showing discontinuous and mixed patterns of variation (see below).
Some tandem arrays containing highly variable duplicates were idiosyncratic. Figure 1 describes Tb109c, which contained four copies with extensive length differences that were genuinely divergent from one another. Other arrays with high average D s obtained this value due to discontinuous variation, that is, they all included an atypical duplicate that did not represent the patterns shown by other copies. Instances where divergence was significantly greater between some duplicates than others are referred to as having 'mixed dynamics'. Tb24 (see Additional file 4) included a divergent copy at its 3' end that lacked the C-and N-termini of other repeats. Tb32 (75 kDa invariant surface glycoprotein) included a divergent copy at its 5' end that differed from other duplicates through frequent indels. Tb113 (cation transporter) and Tb82 also possessed divergent copies at their 5' ends. In these cases, the divergent gene was the exception and other copies varied continuously. The divergent copies were always accompanied by radically remodelled IGSs, which failed to align with other non-coding repeats within the array. For instance, Tb105a is shown in Figure 2; its last copy had a distinctive C-terminus and both 5' and 3' UTRs bore no resemblance to others in the array. This often resulted in obvious structural change within the array, with the divergent copy being physically segregated from the remainder. Consequently, CDS and IGS relationships usually coincided for these cases.
In summary, duplicate genes within tandem arrays typically contained variation and not all elements of the array provided an identical phylogenetic signal. Some duplicates were physically segregated from the remaining array by divergent non-coding sequence; where this had occurred, these genes diverged substantially and variation was distributed discontinuously.

Evidence for concerted evolution
A cladistic criterion was applied to each sequence alignment 'C', to infer concerted evolution, i.e., where paralogous sequences within a species looked more alike than their putative homoeologs in other species. It was trivial to confirm this for invariant arrays, since all copies were identical, but the likelihood of concerted evolution was reinforced by the large genetic distances between sequences evolving in concert in T. brucei and homologs in other species. Additional file 5 shows the ML phylogeny for Tb43 (cysteine peptidase), which is highly conserved within T. brucei and other species. Conspecific sequences obviously cluster together, but there is considerable sequence divergence between clades. In fact, the ratio of interspecific to intraspecific divergence was consistently two orders of magnitude greater for these invariant arrays than more variable cases. Tb43 had a ratio of 217.6, while values for Tb26 and Tb75, both moderately variable, were 3.22 and 1.81 respectively. For all of the 'invariant' arrays, the lack of change observed within clades was not matched by the length of branches relating them; these genes had diverged, but not while contained within a single genome.
Variable arrays also displayed the phylogenetic signal of concerted evolution. Almost all cases for which homoeologous arrays existed in other species produced the pattern of clustering by species. However, this did not extend to T. brucei and T. brucei gambiense; sequences from these organisms were usually interspersed. The pattern presented by Tb17 in Figure 3 is representative of the general trend; high bootstrap values support species-specific clades and this is confirmed by SH tests, described in Table 1. A topology in which Tb17 sequences group by species is not significantly worse than the most optimal tree (ΔlnL = -10.4, p = 0.448), whereas a topology in which presumed orthologs (based on position) from different species are forced to cluster together is significantly less likely (ΔlnL = -3002.3, p < 0.0001). 16 out of 47 cases failed to reject the concerted evolution scenario and Tb80a confirmed this to the extent that T. brucei and T. b. gambiense sequences were reciprocally monophyletic also, (the only exception to the previous statement). In 14 out of 47 instances, the cladistic criterion could not be applied because the array was found only in T. brucei and T. b. gambiense; for these cases it was not possible to infer concerted evolution.
For the remaining tandem arrays, the SH test either rejected the concerted evolution scenario in favour of a pattern of orthology, as for Tb64a; or, as in all other instances, both scenarios were rejected. The optimal tree topology for Tb82 sequences was significantly different to topologies describing concerted evolution and orthology respectively; this occurred because most copies had their Mixed patterns of variation within array Tb109c (hypothetical protein, Tb10.389.0740) Figure 1 Mixed patterns of variation within array Tb109c (hypothetical protein, Tb10.389.0740). A. GENECONV output after analysis of an alignment of gene copies 1-3. Similarity was inferred from the non-random distribution of shared polymorphisms among any two sequences; the likelihood of each phylogenetic grouping is represented by coloured rectangles: copies 1 and 2 as closest relatives (green); 1 and 3 (red); 2 and 3 (blue). B. Diagram of the chromosomal position of Tb109c. DNA strands are shown in grey, scale is in base-pairs. The tandem array is coloured red, other annotated genes are green and hypothetical genes are coloured blue. The identity of UTRs around each gene duplicate is denoted by geometric symbols; identical UTRs possess the same symbol. C. Phylogenetic tree for Tb109c homologs from four species; orthology is retained by the divergent fourth duplicate, producing a distinct clade (shaded). Values are bootstrap proportions out of 500.
Ectopic gene conversion affecting array Tb105a (expression-site associated gene, Tb10.6k15.0910), identified by distribution criterion closest relative in a different species (i.e., still displayed orthology), while some copies in T. vivax, T. cruzi and L. major were homogenised or independently derived. Equally, each copy within Tb106 had its closest relative at a different unlinked locus; Tb106 was only represented in T. brucei, T. b. gambiense and T. congolense and there was no clear indication of the relationships between homoeologous copies. Thus, it is reasonable that this data set should reflect neither concerted evolution nor orthology. A concerted evolution scenario was also rejected by five arrays showing mixed dynamics (Tb24, 64a, 80, 105a, 109c). In these cases the divergent copy retained orthology with homoeologous divergent copies in other organisms (see Figure 1), while the remaining genes in each array were monophyletic. This was sufficient for concert to be rejected, and, conversely, for orthology to be accepted for Tb64a.
In summary, most phylogenies of homologs from various species (alignment 'C', see methods) showed some evidence of concerted evolution, although many arrays could not be tested because they were only present in T. brucei.
In those that could, monophyly of conspecific gene duplicates was usually significantly more likely than orthology with other species. Even in cases where orthology was retained by certain duplicates, others displayed the signature of concerted evolution.

Evidence for allelic gene conversion
Three different criteria were applied to examine tandem gene duplicates from T. brucei only for allelic gene conversion (alignment 'B', see methods). The phylogenetic criterion, applied through the PDM and HMM methods in TOPALi, is illustrated in Figure 3 (Tb17). The distribution criterion, applied through the GENECONV (GCV) method in RDP, is illustrated in Figures 1 and 2 (Tb109c and Tb105a). The compatibility criterion, in the form of RETICULATE, executed in RDP, is illustrated in Figure 4 (Tb62). The consensus of all this information, critically examined by eye, produced a total of 76 GC events, in 21 arrays. These are recorded by sequence in Additional file 6 (Table S1) and by event in Additional file 7 (Table S2).
Like many arrays, Tb17 displayed variation among its duplicate sequences but, despite this, the duplicates clustered together when compared to homoeologs in other species. In addition, Figure 3(A) shows that the pattern of variation was indicative of a gene conversion event at the C-terminus. From here onwards, gene copies will be defined by their position from the 5' end; hence, the third copy of the Tb17 array is written Tb17.3. It follows from  GCV was successful in identifying three regions in the Tb61 multiple alignment where similarities in the distribution of variation were significant. Additional file 3 (A) shows the output of an SSC analysis, which was executed to confirm the GCV results. When only copies 1, 2 and 8 were compared, three distinct breakpoints were visible from a sliding window analysis that otherwise suggested Tb61.1 and Tb61.2 were always closest relatives. In the left-most anomaly, Tb61.2 and Tb61.8 share significant similarity (p = 2.14 e-9 ), which was due to the conversion of Tb61. The majority of putative genetic exchanges were biased, rather than reciprocal. Of 76 observations, only 7 were best explained through the equitable exchange of homologous sequence between two sequences (shown in Additional file 7). All other events involved a given sequence being uncharacteristically similar to unrelated sequences for a given region, and thereby altering its phylogenetic position for the length of that region. These events were interpreted as biased gene conversion.
In many instances, it was possible to delimit the CDS and UTRs within an array, either because EST data was availa-Allelic gene conversion between within array Tb62 (adenylate cyclase, Tb927.6080), identified by compatibility criterion Figure 4 Allelic gene conversion between within array Tb62 (adenylate cyclase, Tb927.7.6080), identified by compatibility criterion. A. SISCAN output, showing the changes in phylogenetic signal along a multiple alignment of gene copies 2, 3 and 4. Genetic distances were calculated to identify the closest relatives from three sequences using a sliding window approach. Coloured lines represent the three phylogenetic resolutions: copies 2 and 3 as closest relatives (red), 2 and 4 (blue) and 3 and 4 (green). The significance of observed similarity is assessed through randomisation of the sequences within a given window and expressed as a Z-score; if the similarity persists over many replicates, outweighing the contribution of any compositional bias, the observation is significant (where the Z-score exceeds 1.96). B. Diagram of the chromosomal position of Tb62. Dashed lines link gene copies to putative breakpoints. C. RETICULATE output. The multiple alignment of all copies is represented by both axes and shaded areas identify where there is incompatibility between the phylogenetic signal of two different regions.
ble, or through comparison with the terminal UTRs of the array. Therefore, it was shown that GC frequently affected the UTRs. Tb62 had three putative breakpoints, identified in Figure 4; the event shown at left spanned the start codon, the second event spanned the stop codon, while the peak shown at right corresponded to a region exclusively within the 3' UTR. Similarly, of the 10 GC events identified for Tb93, only three occurred within the coding sequence while another spanned the stop codon.

Evidence for ectopic gene conversion
The application of GC tests to the 'A' alignments (i.e., tandem gene duplicates from T. brucei plus any other gene family members) showed that ectopic gene conversion was very rare among these tandem gene arrays. Only two gene arrays were implicated in EGC by several methods. The first instance concerned Tb105a, where there were significant regions of similarity between the third gene copy and Tb10.70.2860; these loci are separated by 1.17 mb on chromosome 10. Analysis by SSC identified a region between 519-628 bp that showed significant identity (p = 9.88 e-04 ), which is indicated by the GCV output in Figure  2(A). These elements of similarity stand out despite Tb105a.3 otherwise being closely related to the first and second gene copies in the array, as shown in Figure 2(C). The second instance of EGC involved Tb32, where the tandem array of 75 KDa invariant surface glycoproteins was preceded upstream by an additional tandem pair of gene duplicates, on the opposite strand. Phylogenetic analysis showed that the tandem pair was most likely to be a segmental duplication of Tb32.3 and Tb32.4. However, while the pair still bore these resemblances, one copy showed regions of similarity to Tb32.2. One of these regions was significant (408-603 bp, p = 1.21 e-11 ). Additionally, the tandem pair seemed to have concerted at the N-terminus as a significant GC event was identified here also (136-164 bp, p = 9.61 e-04 ).

Categorisation by evolutionary dynamic and base composition
The tandem arrays in T. brucei were categorized according to the evidence presented above, in order to make sense of the various patterns observed and to test the relationship between double strand-break repair mechanisms and base composition (see Additional file 8). Most arrays fell into one of four categories (see methods for explanation): CE-AGC, CE+AGC, AGC and UCO. Exceptions included Tb105a (Figure 2), which was affected by EGC alone, and Tb82, which failed the cladistic criterion (i.e., did not reflect concerted evolution) and retained orthology for all copies. Five arrays fell into the 'mixed dynamics' category, where not all duplicates were monophyletic but a proportion had evolved in concert. Having grouped arrays according to their evolutionary dynamic, the data set was used to evaluate the theory that gene conversion, through double-strand break repair mechanisms, results in G-C enrichment of genomes. The proportion of G-C at third codon positions (G-C3) was calculated for each category and for 100 random, singleton genes from across the genome. This showed that the average G-C3 for singleton genes was 0.416, while the value for UCO arrays was 0.676; this was a highly significant elevation in third position G-C content (df = 10, t = -11.89, p = 1.57 e-07 ). A significant elevation was also observed among CE-AGC arrays, where pG-C3 was 0.454 on average (df = 8, t = -2.511, p = 0.018). However, neither the CE+AGC (G-C3 = 0.426) or AGC (G-C3 = 0.435) categories produced significantly different G-C3 to single genes.

Discussion
Homoeology, the conservation of genomic position, was crucial to establishing the literal equivalence between loci in different organisms, given that these gene sequences were repetitive by nature, and subject to gene conversion that could generate misleading phylogenies. Using this comparative approach, different levels of variation among duplicates were identified among the 47 tandem arrays studied, reflecting various evolutionary dynamics. Where gene duplicates were invariant, UCO was suggested as the mechanism for concerted evolution, but concert was also widespread among variable tandem arrays. The combination of sequence variation and concerted evolution indicated that these arrays were exposed to GC while still diverging. Finally, arrays with copies displaying different patterns in variation were exceptions that proved the rule: copies within homogeneous repeat structures conformed to a pattern indicative of concerted evolution, while others physically segregated from the repetitive region retained orthology with loci in other species.

Patterns of sequence variation indicative of concerted evolution
For most arrays that could be tested the cladistic criterion firmly rejected orthology between homoeologous gene copies in different organisms. For tandem arrays present only in T. brucei it was not possible to find evidence for concerted evolution, given that lack of variation itself is not evidence. If one accepts that where there is conserved gene order around a homologous array in multiple genomes, this proves that the loci are orthologous (i.e., directly inherited from an ancestor with an equivalent character), the failure of these orthologs to appear monophyletic in a phylogeny is evidence for the alteration of conspecific sequences, i.e., concerted evolution within individual genomes. If one does not accept this as the cause of the unexpected homogeneity, convergence due to positive selection or, more plausibly, the removal of substitutions due to strong structural conservation (i.e., purifying selection) are possible alternatives.
It has been suggested that a mixture of recent gene duplication and purifying selection can explain why tandem gene duplicates appear similar [5]. In this study, purifying selection is not a sufficient explanation since the gene copies have evolved; including those 'invariant' arrays showing large genetic distances between species. Whatever small amount of variation exists should still reflect the orthologous relationships indicated by shared genomic position. Therefore, purifying selection would only be adequate if posited after the isolation of a genome, coinciding with the origin of an array by gene duplication de novo. Rapid gene duplication would then result in homogenous gene copies, clustering by species in a phylogeny. Yet, many arrays in this study were clearly not exposed to purifying selection and were evolving neutrally; even invariant arrays displayed large interspecific genetic distances and so had diverged within each lineage. And it is implausible that an array present in five different species was inherited by each as a single gene and was independently expanded five times. It is more parsimonious that the array present in all species was an array in their ancestor.
Among the 'invariant' tandem arrays, concerted evolution was the result of unequal crossing-over between homologous chromosomes at meiosis. The absence of any variation within the non-coding parts of the arrays, near-zero D n and D s values in the CDSs and, conversely, substantial divergence from homologs in other species, suggest that variation has been purged from these sequences irrespective of codon position, rather than constrained by purifying selection. In other words, while UCO results in dispersal of polymorphisms among gene duplicates and thus 'invariant' arrays within individual genomes, divergence between homologs in different genomes continues.
In agreement with previous observations, the G-C3 content is significantly elevated among these arrays, supporting the theory that double-strand break repair mechanisms cause G-C enrichment [45,46]. In T. brucei, homogeneity has previously been observed among the highly-conserved heat-shock proteins [29,47]. These and other arrays in Trypanosoma spp., such as histones [48,49] and cysteine protease [50,51], emulate the homogeneity among important structural or 'housekeeping' genes in other organisms [9].
The idea that purifying selection could explain concert among arrays with greater levels of variation is even more unlikely, given that non-coding regions were once again often more conserved, and many genes were evolving neutrally. The coincidence of variation and concerted evolution demands a different mechanism: gene conversion. These arrays suggest that repeated gene conversion can cause tandem duplicates to evolve in concert, in the presence of neutral divergence (or even positive selection).
Although no GC was observed in 7 arrays where concerted evolution was inferred under the cladistic criterion, the incidence of GC suggests that both homogenising and diversifying processes operate simultaneously within these genomic structures. It is not clear which would occur more frequently, or where the balance of forces would lie; however, increased divergence would inevitably reduce the probability of conversion between two sequences and less conservative regions of a gene may rapidly escape conversion after duplication, allowing unfettered divergence.
In summary, the action of concerted evolution in tandem arrays is betrayed by the loss of orthology between homoeologous genes in different organisms, yet this does not seem to prevent divergence between and within genomes. These results show that paralogous sequences can become more similar over time, even in the presence of divergence, and that such concerted evolution is only one conservative process affecting duplicate sequences, which may be simultaneously affected by other conservative or innovative forces.
Allelic gene conversion as a mechanism of sequence evolution 21 out of 47 tandem arrays contained chimaeric gene copies, which generated putative cases of partial gene conversion. In the majority of cases this gene conversion was biased and presumably occurred with the effect of propagating one sequence motif at the expense of another. On a few occasions, recombination between two sequences was the best explanation for incompatible motifs in an alignment. The extent of AGC may be difficult to determine since it cannot be known how many exchanges have been 'overwritten' by subsequent events. And as no distinction was made between homologous chromosomes in the T. brucei genome sequence, it is not known what proportion of events occurred between alleles during meiosis. The authenticity of these putative GC events was checked by collecting voucher sequences from those read libraries that were available (chromosomes 9, 10 and 11). Of all GC events originating on these chromosomes, only four could not be substantiated by single read sequences for three gene copies containing incompatible motifs (see Additional file 7). For these cases, the contradiction between one part of the molecule and another could be due to the assembly of individual reads from different copies. However, single reads were retrieved for most putative events that covered the breakpoint and represented both incompatible motifs, eliminating the possibility of such mis-assembly.
AGC is a plausible mechanism for homogenisation of many of these tandem arrays. However, in certain circumstances, where sequences diverge rapidly, AGC may have a further diversifying effect by helping to disperse sequence motifs among duplicates. Tb93 and Tb17, both coding for surface-expressed proteins, may be examples of where the rate of AGC is outstripped by sequence divergence, preventing a net homogenising effect of GC (although, by sharing around sequence motifs, array copies would still cluster by species in a phylogeny). Contingency gene systems are a dramatic illustration of the power of GC to generate allelic diversity through the creation of chimaeric sequences. Contingency genes have evolved independently in various pathogens and parasites to provide a rapid solution to immunity raised against an existing epitope. Typically they comprise a large number of whole and partial surface antigens, only a fraction of which are expressed. A novel gene is generated by transposing an inactive motif into an expression locus, replacing the existing sequence, for example, in the generation of novel pilin in Neisseria [52], translocation of VSG pseudogenes into expression sites in T. brucei [53,54] and similarly in Babesia bovis [55] and the cestode Echinococcus [56]. These systems employ EGC in the generation of chimaeras; this, and the concerted evolution of conserved gene families, such as β-tubulin in Leishmania [57], suggested a priori that EGC played a prominent role in the evolution of gene duplicate sequences across the genome.
However, in contrast to AGC between tandem duplicates physically associated within an array, exchange between gene family members in trans, or separated from the array, was seldom observed. It is reported that conversion frequency is lower between distant and dissimilar gene copies [37,38,58,59], which suggests that gene conversion is unlikely to affect the sequence homogeneity of gene families in general. Only two arrays, Tb32 and Tb105a, produced a substantiated incidence of EGC, despite initially promising patterns in Tb12 and Tb26. The latter two demonstrated the difficulty in confirming gene conversion events; the presence of recently duplicated gene copies near to the tandem array gave the impression of gene conversion with array copies. However, any similarity between genes in and out of the array was interpreted as common inheritance from a duplication event, and any incompatibilities were introduced subsequently by allelic gene conversion between newly duplicated copies rather than through ectopic exchange with the array.

Functional differentiation and the spatial context of concerted evolution
UCO and GC are recognised mechanisms of concerted evolution [9]. Both processes occur due to the aberrant alignment of repetitive sequences in cis or trans. Since these processes involve the interaction of repetitive sequences, the frequencies of UCO and GC decline with physical distance and sequence dissimilarity [60,61], which is corroborated by the paucity of EGC here. The limitations on concerted evolution occur because it is an explicitly spatial phenomenon [62], for instance, it is largely limited to tandem arrays in T. brucei, and does not extend to gene families in general. In this study, the structural limitations on concerted evolution are very well demonstrated by arrays with mixed dynamics. While gene duplicates in arrays such as Tb109c and Tb24 show physical proximity and sequence homology typical of other arrays, the repeat structure of the array is disrupted through extensive divergence in flanking and, sometimes, coding sequences. In Additional file 4, Tb24.4 and 5 have evolved unique UTRs and are physically segregated from Tb24.1-3 by the insertion of two unrelated genes; perhaps as a consequence, they retained orthology with homoeologs in other species. Similarly, Tb109c.4 ( Figure 1) is physically separated from the array by an unknown CDS and is shown to retain orthology with homoeologs in T. cruzi, T. congolense and L. major. And likewise, Tb105a.4 ( Figure 2) is segregated from the other, concerted duplicates by its UTRs, which show no affinity with others in the array, and is more closely related to other genes around the genome. These cases link the disruption of the repeat structure, by rearrangements to IGS or CDS boundaries, with divergent gene copies that retain orthology across species, rather than converging with conspecific paralogs.
Previous studies have shown that tandem arrays have the capacity to harbour important functional variation. The phosphoglycerate kinase (PGK) array in T. brucei comprises three distinct sequence types, where the second copy is a chimaera of the first and third, created through gene conversion [63,28]; this variation is responsible for the expression of specific glycosomal and cytoplasmic isoforms [64,65]. In a similar instance, the hexose transporter array in T. brucei reflects the spatial and functional partition of two different isoforms, each with several identical duplicates, which may again have been affected by gene conversion between the isoforms [30]. In T. cruzi, duplicates within the pyr gene array controlling pyrimidine biosynthesis are also functionally distinct and expressed in particular life-cycle stages [66]. Sequence variation in the 3' UTRs of β-tubulin in trypanosomatids has been directed related to differential regulation of gene expression in different life stages [67,68]. Hence, given the abundance of sequence variation evident here, functional differentiation of tandem duplicates may be commonplace, suggesting a subtle functionality within tandem arrays not yet appreciated in T. brucei. The presence of both gene duplicates evolving in concert and others retaining orthology in mixed dynamic arrays demonstrates how gene conversion is a spatial phenomenon. Rearrangements of the primary or secondary structures of divergent gene copies make them permanently distinct and reduce the probability of recombination with a nearby duplicate. This offers a model for how functional differentiation can be reconciled with the ubiquity of concerted evolution observed here.

Conclusion
Contrary to many previous descriptions from trypanosomatids, tandem gene arrays in T. brucei contain substantial genetic variation among gene duplicates. Despite this, concerted evolution is pervasive, and evidence from various loci indicates that both unequal-crossing over and allelic gene conversion contribute to homogenisation. Therefore, gene duplicates can both evolve in concert, relative to homoeologs in other species, and diverge within an array. This indicates that multiple evolutionary forces are apparent, with sequence variation determined by a balance of conservative and innovative pressures. Structurally heterogeneous tandem arrays demonstrated that the influence of concerted evolution is mediated by the spatial environment of gene duplicates; it can be prevented by disruption of the repeat structure in the array. Functional differentiation of duplicates can and does occur through such rearrangements in intergenic sequences. Together, the dynamics of sequence variation among tandem duplicates show how various forces affecting gene sequences must be integrated with typical 'birthand-death' ideas of gene family evolution, to produce an adequate model. Abundant variation among gene duplicates indicates that tandem duplication has greater consequence than simply increasing gene dosage; the systematic account given here provides a basis for uncovering the hidden functionality within tandem gene arrays in trypanosomatids.

Methods
The primary resource for all analyses was the completed genome sequence for Trypanosoma brucei [6], which was available from the Wellcome Trust Sanger Institute's interface 'GeneDB' [69]. When evaluating concerted evolution, the completed sequences of Trypanosoma cruzi [70] and Leishmania major [44], and the draft genome sequences for other trypanosomatids (T. brucei gambiense, T. congolense, T. vivax), were also utilised through the same interface. The T. brucei genome consists of 11 chromosomes ranging between 1 and 5 mb. Each chromosome was scanned by eye for duplicate genes positioned in tandem; where the similarity was not obvious from annotation suspected duplicates were tested by reciprocal BLAST and sequence alignment. A catalogue of all tandem gene pairs and gene arrays is available as supplementary material (see Additional file 9).
Since this study concerned the variation among gene copies and the approach taken was inherently phylogenetic, only those tandem arrays with >3 copies in the annotated genome sequence were selected for analysis; these included tandem segmental duplications where the repeat units comprised multiple dissimilar, contiguous genes. The component coding sequences (typically, trypanosomatid genes do not contain introns, although they do occasionally occur [71]) for each duplicated locus were then extracted. The protocol for analysis of sequence variation is described in Additional file 8. Essentially, after preparation of multiple alignments, the variation among sequences was quantified and characterised using phylogenetic estimation. The cladistic, phylogenetic, distribution and compatibility criteria described previously were applied to detect concerted evolution, AGC and EGC in turn. Each array was then categorised according to the results of these analyses to examine the relationship between evolutionary dynamic and base composition.

Data preparation
Duplicate sequences were aligned using the ClustalW tool [72] for each tandem array, within the BioEdit platform [73]. Sequences were translated prior to alignment to preserve codon structure and then returned to nucleotide form for analysis. Each alignment was checked by eye. Three types of alignment were created: • A: duplicate sequences from the array were combined with all other gene family members. Untranscribed regions (UTRs) were not included because they usually could not be aligned. These alignments were intended to test for monophyly of array copies and to test for ectopic exchange between unlinked loci in the same gene family. These alignments were not possible in 15/47 cases because no further loci existed beyond the tandem array.
• B: the coding sequence and UTRs of each copy from the array were aligned. UTRs were defined using either expressed sequence tag (EST) data, where this existed, or on the basis of identity between internal and external intergenic sequences (IGS). The first and last copies of the array usually had curtailed (sometimes entirely different) 5' and 3' non-coding regions respectively. The length over which these regions displayed identity with the internal IGS was assumed to reflect the UTR. These alignments were created to estimate variation and phylogeny among array copies and to evaluate the evidence for allelic gene conversion.
• C: the coding sequence for each array copy in T. brucei was combined with homologs from up to five other species, where available. UTRs were not included because these often do not align over large phylogenetic distances. Homologs were found by searching each genome sequence with the gene sequence of surrounding genes in T. brucei. For a homolog to be accepted, it had to display 'homoeology' with T. brucei, i.e., conservation of surrounding gene order.

Analysis of concerted evolution
Variation among duplicate sequences within arrays was quantified with the numbers of synonymous (D s ) and non-synonymous (D n ) substitutions per site, using DNAsp [74]; these measures were favoured because they reflect purifying selection, an alternative to UCO in explaining low divergence. Concerted evolution should result in a clade of tandem duplicates, excluding homologs in other species. Hence, the cladistic criterion was applied by estimating a phylogenetic network for Alignment C. A network has advantages as an initial analysis because ambiguity in a phylogenetic hypothesis remains explicit, but it is also realistic where sequences have mixed histories (as many genes may have among these data sets). The networks were estimated using the neighbour-net method, a clustering algorithm executed in Splitstree v4.0 [75]. To correct for potential base composition bias, a logdet genetic distance matrix was used after removing parsimony-uninformative sites.
Where all array copies were monophyletic, this was confirmed by estimating maximum likelihood phylogenetic trees using PHYML [76], with 500 non-parametric bootstraps. The trees were unrooted, a GTR + I + Γ model was applied and parameters were estimated from the data [77]. Statistical evaluation of the concerted evolution hypothesis was obtained using a Shimodaria-Hasegawa test (SH, [78]), applied to alignment C; the test was executed using DNAml in the PHYLIP 3.65 package [79]. The SH test compared the most optimal topology without constraints with two tree topologies in which duplicate sequences were monophyletic by species (i.e., concerted) and by position in the array (i.e., retained orthology) respectively. The 'concerted' and 'orthologous' scenarios could be rejected where there was a significant departure in the likelihood between optimal and constrained topologies. Failure to reject the null hypothesis that the constrained topology was sub-optimal showed that any apparent differences were explicable due to sampling error.

Analysis of allelic gene conversion
AGC is defined as recombination between gene duplicates within an array, either between copies on the same chromosome or between homologous chromosomes. Phylogenetic, compatibility and distribution criteria were applied in detecting AGC events. All putative GC events were checked by eye and many were rejected when related back to the multiple alignment.
The phylogenetic criterion detected GC by identifying a significant difference in tree topology between neighbouring 'windows' in a sliding window analysis. This was conducted using the TOPALi platform [80], with the probabilistic divergence measure (PDM, [81]) and hidden Markov model (HMM, [82]) methods. The PDM method scans the entire alignment, estimating phylogenies for fixed sequence windows with a Markov Chain Monte Carlo process. A significant change in the posterior probabilities of two neighbouring windows is interpreted as a breakpoint. A step size of 1 and window size of 10 were applied, with otherwise default settings. Possible breakpoints were further analysed using the HMM method, which detects recombination breakpoints among four sequences by modelling the probabilities of transitions between the three tree topologies. This determines which topology, at any point along the alignment, is most probable. Both PDM and HMM models are designed to compare tree topologies only and minimise the bias introduced by rate heterogeneity among and within sequences [80].
The compatibility criterion detected breakpoints by identifying neighbouring sites in the multiple alignment that gave incompatible phylogenetic signals. These sites were plotted on to a 2-dimensional space using the program RETICULATE [35], within the RDP v2.0 platform [83].
Only binary sites were used and 100 randomised matrices were estimated. These breakpoints were again used to inform subsequent methods of GC detection.
The distribution criterion detected regions of significant similarity between sequence pairs by comparing the distribution of silent polymorphisms along the alignment, i.e., those changes unaffected by purifying selection. Significant clustering of identical polymorphisms was identified by GENECONV (GCV, [84]), executed within RDP v2.0. The analysis compared sequence triplets, using default settings. P-values for GC events were calculated by permutation of the data set, with 1000 randomisations. GENECONV produced a large number of candidate gene conversions, but many of these were created by rate differences within and among sequences. Each candidate event was checked to remove the obvious mistakes caused by regions of high conservation or derivation. Furthermore, each event was checked using a second method in the RDP package, SISCAN (SSC, [85]). This method conducts a sliding window analysis with a phylogenetic criterion, calculating a p-value through sequence permutation. Default settings were applied.
Putative AGC events required validation. Tandem gene arrays are among the most challenging parts of genome sequences to assemble, for much the same reason that non-coding repetitive regions are difficult to resolve; assembling sequence reads into contigs relies on unique motifs to define position. It is especially difficult to determine the precise number and order of copies within an array, where the copies differ by less than 5%. Consequently, in tandem arrays with low variation, gene copies may be mis-assembled and appear chimaeric because parts of different copies have been conjoined. To consider this, the original sequence reads were collected for every three sequences involved in a putative GC event on chromosomes 9, 10 and 11, to act as vouchers. If sequence X and Y were closest relatives, except for one region where X was more akin to Z, the appropriate regions of all three sequences had to be represented in sequence reads (each covering 0.5-1 Kb) to confirm that a phylogenetic conflict truly existed.

Analysis of ectopic gene conversion
Where alignment A was available, tests for GC were repeated as described above for each array. These comparisons were designed to identify significant similarity between gene duplicates at different loci that may be caused by EGC. Putative EGC events were validated by collecting voucher sequences as before, although the potential for mis-assembly is reduced for single-copy genes.

Relationship between gene conversion and base composition
T. brucei tandem arrays were classified according to the results of the preceding tests. The categories shown in Additional file 8 classify arrays by the presence or absence of AGC in conjunction with concerted evolution (CE+AGC and CE-AGC arrays); the presence or absence of EGC in conjunction with concerted evolution (CE+EGC and CE-EGC arrays), those with some duplicates retaining orthology, and others concerted ('mixed dynamics'); those unique to T. brucei and showing AGC; those that have retained orthology; and those with identical copies within T. brucei but divergence between species (i.e., UCO or 'invariant' arrays). To test the hypothesis that double strand-break repair mechanisms are biased towards guanine-cytosine (G-C) and result in G-C enrichment of regions undergoing recombination or biased gene conversion [45,46], the average G-C content at the third codon position (G-C3) of each class of tandem array was calculated using the base composition function in GCUA [86]. This was then compared to the average G-C3 value for 100 singleton genes chosen at random from all chromosomes.

Additional material
Additional File 1 Figure S1. Models