BACKGROUND
the cytochrome p <dig> monooxygenases  constitute a conserved gene superfamily of heme-thiolate proteins ubiquitously distributed in different life forms, including prokaryotes , unicellular eukaryotes  and multicellular eukaryotes   <cit> .
they play key roles in the metabolism of a wide variety of substrates, including endogenous chemicals such as steroids and other important small molecules, and also xenobiotic compounds including drugs, pesticides and environmental contaminants  <cit> .
substrate and functional diversity are characteristic features of p <dig> enzymes and are considered to be the consequence of evolutionary adaptation driven by different metabolic or environmental demands in different organisms.
evolution and expansion of major p <dig> branches has been suggested to be linked with the historical occurrence of important evolutionary events.
one particular example is the divergence of p450s of the common plant-animal ancestor either to synthesize biochemicals/metabolites  or to detoxify xenobiotics , followed by p <dig> gene expansions, especially in the plants  <cit> .
these may well reflect different survival strategies adopted between the two kingdoms, i.e.
plants evolved sessile systems with p <dig> enzymes with more diverse and essential roles, while animals developed higher order sensory and locomotor systems, and comparatively fewer p450s  <cit> .
although over  <dig>  p450s have been named and about  <dig>  p <dig> genes have been characterized to date, none has been systematically characterized in the phylum of ciliophora which, together with dinoflagellates and the exclusively parasitic apicomplexa, constitute the three major evolutionary lineages that make up the alveolates  <cit> .
the unicellular ciliate tetrahymena thermophila is a free-living protozoan widely distributed in freshwater and estuarine environments, elaborating typical eukaryotic components  into a highly organized cell whose structural and functional complexity is comparable to, or exceeds that, of human and other metazoan cells  <cit> .
the physiology, biochemistry and molecular biology of tetrahymena are well characterized  <cit> , and it is an excellent model organism for toxicological and ecotoxicological studies in aquatic toxicity test systems  <cit> .
results from the est project  <cit>  and the macronuclear genome sequencing project  <cit>  have shown that, although single-celled, tetrahymena possesses core processes conserved across a wide diversity of eukaryotes  that are not found in other unicellular model species such as the yeasts saccharomyces cerevisiae and schizosaccharomyces pombe.
it also contains a large number of gene families that are involved in processes associated with sensing and responding to environmental cues.
in the case of the cytochrome p <dig> gene family, s. cerevisiae and s. pombe have only three and two p450s, respectively  <cit> , while in t. thermophila the number is more than  <dig> , which is close to the typical number  in a vertebrate genome  <cit> .
in humans and other mammals, extensive studies have focused on aspects of p <dig> gene structure and biochemical properties.
important biological functions of p450s and the associated high degree of complexity in the gene polymorphism and expression patterns have been demonstrated  <cit> .
as genomic databases became available and functional genomics techniques such as dna-microarrays have been applied, investigations on p <dig> isoforms have also been extended to other organisms  .
these developments have, in turn, led to the emergence of a new field, ecotoxicogenomics  <cit> , which aims to develop effective tools for identification of possible toxic environmental pollutants by characterizing their effects on terrestrial and aquatic model organisms, such as the soil-dwelling nematode caenorhabditis elegans  <cit> , the aquatic crustacean daphnia magna  <cit> , and the small fish fathead minnow   <cit> .
we previously identified a series of differentially expressed ests of t. thermophila that respond sensitively to treatment with the organochlorine insecticide ddt  <cit> .
one est  was identified to be a p <dig> gene by homology searches, and its expression levels under different concentrations of ddt treatment were further assessed.
recently, the first genome-wide microarray platform containing the predicted coding sequences  has been established and validated in t. thermophila and was used to study gene expression during three major stages of the organism's life cycle: vegetative growth, nutrient starvation and conjugation  <cit> .
substantial progress has also been made in closure and reannotation of the mac genome sequence of this eukaryotic model organism  <cit> .
all these provided us the opportunity to investigate both the functional and evolutionary characteristics of the cytochrome p <dig> genes in t. thermophila at the genomic level.
in this study, the putative t. thermophila p <dig> genes that were previous identified both by the international nomenclature committee  and by the tigr genome annotation team were further checked by est data and through cdna sequence cloning experiments for improvement of annotation.
the expression patterns of the t. thermophila p <dig> genes during three important cell physiological/developmental conditions  were analyzed based on est and microarray data.
these results are discussed in the context of understanding the characteristics of the t. thermophila p <dig> monooxygenase isoforms and their functional and evolutionary diversity.
RESULTS
t. thermophila putative p <dig> gene sequences
the putative cytochrome p <dig> gene sequences in t. thermophila were initially identified in  <dig> by dr. nelson based on an early tigr release of its macronuclear genome assemblies, and  <dig> p450-like genes  were posted on his p <dig> website .
when the draft whole genome of t. thermophila was sequenced and reported in  <dig>   <dig> p <dig> genes were predicted by the tigr genome annotation team  <cit> .
in this study, we retrieved the putative t. thermophila p <dig> genes from above two locations and further searched against the tgd  utilizing the updated annotation .
the obtained p450-like gene sequences were thereafter compared.
then, we used the data from two est database resources  and our own rt-pcr investigations to verify the accuracy of the predictions.
a total of  <dig> putative cytochrome p <dig> gene sequences with full-length orfs were identified in t. thermophila and were assigned names as suggested by the "p <dig> nomenclature committee" .
our results showed that the gene predictions made by nelson are more accurate than those presented by the tigr annotation.
the tigr predictions have some errors in identifying the start codons, especially the signal peptide, and intron-exon boundaries of several p450s are inaccurate.
in one case, it merged  <dig> adjacent p <dig> genes  into one "monster" gene , which also includes two intergenic regions plus  <dig>  kb upstream of cyp5013a <dig> .
for many other p <dig> genes, the tigr annotations are systematically short, presumably due to its failure to recognize the signal peptide at the 5' end.
this problem may arise if the hydrophobic region is coded by at-rich codons, which are not easily distinguished statistically from the at-rich intergenic regions .
for all the  <dig> p <dig> genes, only three sequences  were modified slightly in our annotation compared with the predictions made by dr. nelson.
gene id, genbank accession numbers, scaffold genbank accession numbers, est accession numbers and predictions of sub-cellular localizations were listed.
s: secretory pathway.
*: found only in tbestdb.
there were four p <dig> pseudogenes  in nelson's predictions .
it was revealed that one  was a probable pseudogene due to the presence of in-frame stop codons, and two partial fragments  with typical p <dig> sequence features appeared to be missing the n-terminus.
these sequences were not included in the following analysis of the t. thermophila putative p <dig> genes.
however, when checking the cdna sequences of putative t. thermophila p <dig> genes, we identified one transcription product that is most identical to the predicted pseudogene "cyp5005a7p" which has an in-frame tga codon  within its orf.
compared with the "cyp5005a7p" sequence, the est of the cdna transcript has several site mutations including one nucleotide transition  within the in-frame stop codon found in the genomic sequence, and thus can be read through.
however, a blast search of the t. thermophila whole genome database failed to detect any sequence that completely matched this cdna sequence.
thus we designed a pair of primers  located in either the 5' or 3' flanking region of the "cyp5005a7p" sequence, and obtained a  <dig> kb pcr amplification product that was sequenced.
due to the possibility that different p <dig> isoforms may exist in different t. thermophila strains, we checked the pcr products both from the genomic dna of strain sb <dig>  and cu <dig>  the results showed that the pcr products from the two strains were 100% identical and were also consistent with the sequence of the unexpected cdna transcript.
we thus assigned the name cyp5005a <dig> to this newly observed p <dig> isoform.
further, blast searches in the tbestdb and ncbi est database either found one sequence in each that perfectly matched the cyp5005a <dig> sequence .
besides, the microarray data demonstrated that the corresponding cdna transcript was constantly expressed during the different cellular conditions.
all the above observations serve as evidence that cyp5005a <dig> is a functional p <dig> gene in both the t. thermophila cu and sb <dig> strains.
the erroneous prediction of the pseudogene "cyp5005a7p" by tigr was probably due to false assembly of genome sequence contigs, caused by the high similarity between the cyp5005a5p and cyp5005a <dig> sequences .
structural features and intron-exon organization
p <dig> genes are classified and annotated on the basis of amino acid sequence identity, phylogenetic homology and gene organization, and a four-digit naming system has been set up to meet the need of the increasing number of newly discovered p <dig> sequences.
all the t. thermophila p <dig> sequences are distributed into  <dig> families and  <dig> sub-families according to the nomenclature criterion of sequence similarity   <cit> , indicating p450s exist as a large gene family in this organism.
the phylogenetic analysis of t. thermophila p <dig> genes obtained using either a neighbor-joining tree or a maximum-likelihood tree gave similar topological nodes, and bootstrap testing showed good reliability of the phylogenetic tree .
the largest expanded p <dig> family in t. thermophila is the cyp <dig> family that contains  <dig> likely functional p <dig> genes  and the two predicted pseudogenes .
the cyp <dig> family  and the cyp <dig> family  were "medium"-expanded p <dig> families in t. thermophila that possessed  <dig> and  <dig> genes, respectively.
moreover,  <dig> cyp <dig> family members are located on the same scaffold , which contains a gene cluster in the form of tandem repeats with four genes  and one pseudogene  in the same orientation.
three cyp <dig> family members  form a tandem triplicate repeat on the scaffold ch <dig>  besides these, there are three other p <dig> gene pairs  organized in tandem duplication.
such a feature is consistent with the gene expansion strategy used by some other lower eukaryotic genomes, particularly those who have evolved to meet an extensive demand for generation of a broad range of metabolites in the secondary metabolism when interacting with the environmental niches.
for example, it was revealed that the majority of the multigene p <dig> families in the model white rot fungus phanerochaete chrysosporium appear to have expanded locally in its genome as a result of extensive gene duplications and rearrangements, indicating a strong need for functional divergence in response to environmental stimuli  <cit> .
the intron-exon organization of p <dig> genes exhibits a diversity of gene structure with indicating that multiple gains and losses of introns have occurred during the evolution of p <dig> genes in diverse species, with little conservation of intron positions among divergent p <dig> families  <cit> .
intron positions were mapped and characterized for all  <dig> putative t. thermophila cyp genes .
of the  <dig> introns that were identified in the p <dig> gene sequences,  <dig> were phase zero ,  <dig> were phase i  and  <dig> were phase ii .
meanwhile, based on the definition of the uip   <cit> ,  <dig> %  of the total introns were unique and the remaining  <dig> introns were present in  <dig> consensus locations  among three different family clades.
sixteen p <dig> genes lack introns entirely.
five genes  have the maximum number of four introns, respectively.
from figure  <dig>  it can be observed that there is a good correlation between the conservation of intron position and phylogenetic relationships of t. thermophila p <dig> subfamily members.
the evolution of introns in alveolates was recently studied by nguyen et al.
<cit> , who concluded that the rates of intron gain and loss were more or less constant in the last ~ <dig> myr after tetrahymena branched off.
therefore, it can be inferred that the enlargement of several p <dig> families in t. thermophila resulted from recent separate small duplication events, which is also apparent within many other tetrahymena protein families containing paralogous genes  <cit> .
our analyses of the est data to determine the intron boundaries also revealed that the cyp5013a <dig> gene might exhibit alternative intron splicing.
alternative splicing was suggested to be an uncommon phenomenon in t. thermophila, at least under the several growth conditions that have examined  <cit> .
cyp5013a <dig> had the greatest number of ests of all t. thermophila p <dig> isoforms.
among the  <dig> retrievable est sequences,  <dig> represented the correct transcriptional product of cyp5013a <dig>  while one  retained the first intron.
the different transcripts were further investigated by rt-pcr.
an intron-containing transcript was observed in pcr products both from the cdna templates of starvation and conjugation stages, but was barely detectable from vegetative growing cells .
the retained intron in cyp5013a <dig> is  <dig> bp in length and belongs to phase i type.
thus, the putative intron-containing transcript should contain a frame-shift mutation.
rna transcripts carrying such premature stop codons can be prevented by the nmd pathway  <cit> .
recently, this pathway was found to play an essential role in another ciliate, p. tetraurelia, by knocking down the two homologous genes encoding upf <dig>  a protein which is crucial in nmd, thus indicating a universal translational control of intron splicing in eukaryotes based on nmd surveillance  <cit> .
a blast search of the t. thermophila genome using p. tetraurelia upf <dig> genes as queries revealed the existence of one homolog with 70% identity .
thus, intron retention in the cyp5013a <dig> transcript might be caused by inefficient nmd activity during specific physiological/developmental stages.
however, it would lead to a mature "non-sense" transcript and could only serve to down-regulate expression of the cyp5013a <dig> gene.
whether such rarity of alternative splicing in this species may simply be tolerated or directed by other mechanisms, is currently unknown.
all known p450s appear to take on a similar folded structure, yet frequently show less than 30% sequence identity and have very different substrates  <cit> .
in t. thermophila, the average length of the p <dig> amino acid sequences in the alignment is around  <dig> amino acid residues , and the genes share different degrees of amino acid sequence conservation between different families.
to better understand the structure and functional relationship of p <dig> enzymes, mammalian p <dig> cyp3a <dig>  seven t. thermophila p <dig> sequences  and one paramecium tetraurelia p <dig> sequence  were selected and aligned .
secondary elements were then assigned using cyp3a <dig> as the template for secondary structure, based on its known crystal structures and six putative srss were identified in the aligned sequences according to gotoh's predicted models  <cit> .
srss participate in the contact with ligands and many of them locate at regions with variable structural elements.
all the t. thermophila p <dig> sequences contain the typical conserved p <dig> domains, including the heme-binding region  near the c-terminal, the perf domain  and the k-helix .
multiple sequence alignment also revealed that some residues are identical across t. thermophila p450s; most of these are distributed in the five regions that constitute the conserved p <dig> structural core, including: the c helix , the c-terminal part of the e helix , the i helix , the j helix , and the heme-binding loop between helix k  and l  .
t. thermophila p <dig> member composition
prediction of sub-cellular localization of the cyp proteins showed that no predicted mitochondrial localized p <dig> was found in t. thermophila , and all of the genes except cyp5001a <dig> encode a putative typical microsomal signal peptide of about  <dig> hydrophobic residues that likely serves as a membrane anchor in the endoplasmic reticulum.
the secretory pathway of eukaryotic cells consists of a series of discrete, membrane-bound compartments, including the er, golgi, and vacuoles  <cit> , all of which are present in tetrahymena  <cit> .
for cyp5001a <dig>  the predicted signal sequence was ambiguous and also lacked the positively charged residues in this n-terminus critical for endocellular transportation targeting.
so far, mitochondrial p450s have been described only in the animal kingdom  <cit> .
as with t. thermophila, none was found in fungi and plants.
these observations support the suggestion that mitochondrial p450s did not originate from the ancestral mitochondrial endosymbiont, but evolved later, possibly by mistargeting of microsomal p450s  <cit> .
cyp <dig> is the only p <dig> family that has orthologs in multiple phyla of the animal, plant, fungal and bacterial kingdoms, and it has been postulated to be the ancestral p <dig> isoform  <cit> .
this family is functionally conserved and has a very limited range of substrates related to the biosynthesis of sterols and their derivatives, and a high conservation within the srss is a specific feature of cyp <dig>  it was believed that two motifs, "yxxf/lxxpxfgxxvxf/yd/a" and "gq/hht/ss", present within the regions of srs <dig> and srs <dig>  respectively, can be considered as the cyp <dig> signature  <cit> .
in this study, no cyp51-like gene has been identified with such a signature in the t. thermophila p <dig> family.
so far, all studied cyp <dig> family members were found to exhibit catalytic function in the oxidative removal of the 14α-methyl group from sterol precursors formed downstream from cyclization of squalene  <dig> -epoxide, and the 14α-demethylated products are intermediates in the sterol biosynthetic pathways leading to formation of cholesterol in animals, ergosterol in fungi and a variety of functional sterols in plants and algae  <cit> .
cyp <dig> appears to be absent in insects and nematodes due to the fact that they don't make cholesterol and lack the post-squalene sterol biosynthesis pathway, and absence of the cyp <dig> gene was thought to be the consequence of gene loss during evolution  <cit> .
thus, a similar cyp <dig> gene loss may be inferred in t. thermophila.
instead of the sterols in most other eukaryotes, tetrahymena produces the pentacyclic triterpenoid alcohol 'tetrahymanol' and/or hopanoids de novo as functionally equivalent structural components of cell membranes, and its ability to synthesize and use this "primitive" substance can be considered as a metabolic relic  <cit> .
besides, although t. thermophila is able to incorporate exogenous sterols into cell membranes and convert those to various derivates, it utilize enzymes other than p450s in the downstream metabolism of certain steroid compounds, such as the cytochrome b <dig> in the direct desaturation reaction of cholesterol  <cit> .
the ciliate p. tetraurelia is the species closest to t. thermophila that has an available genome database and phylogenetic analyses between these two species can help to understand the evolutionary relationship of t. thermophila p <dig> families with other ciliate members.
a phylogenetic network was constructed of both the t. thermophila and p. tetraurelia p <dig> protein sequences .
in the branch that constitutes three p. tetraurelia p <dig> genes , two genes from t. thermophila  also present and cluster together.
however, according to the deep branching of these sequences and the relatively low bootstrap values on the branch nodes , there is yet no strong evidence supporting the possibility that they represent orthologs shared by the two species.
surprisingly, the other major clades of p <dig> genes seem to be highly specific to t. thermophila or to p. tetraurelia.
moreover, the number of p <dig> genes in p. tetraurelia is less than half of that in t. thermophila, and blast searches in different est databases revealed that only about 1/ <dig> of the p. tetraurelia p450s have retrievable est data .
it has been pointed out that while tetrahymena and paramecium are related species within the diverse clade of ciliates, their genomes have obviously evolved by significantly different mechanisms  <cit> .
the paramecium genome has gone through a series of whole genome duplications and may exhibit the most completely preserved whole-genome duplication described to date, where ~95% of the genome is located in duplicated blocks  <cit> .
by contrast, although up to 40% of the t. thermophila genome was also indicated to be located within duplicated blocks, the most of the large expansions of its gene families may have arisen through tandem duplication events  <cit> .
while differences in the p <dig> genes in the two species may reflect the difference in evolutionary strategies between the two organisms since their divergence, further investigations are needed to address the adaptive significance of these diversities relative to the well known functions of cytochrome p450s in sensing of and responding to environmental changes, particularly environmental stress.
p <dig> gene expression analysis
spatial and temporal gene expression patterns are important aspects of gene regulation and expression pattern analysis has played an important role in the study of the cytochrome p <dig> gene superfamily.
in the present study, the gene expression patterns of all t. thermophila p450s during three important cell physiological stages were analyzed based on est and microarray data.
to examine the expression profiles of the p <dig> genes, we initially retrieved the est data of t. thermophila from both in the ncbi est database and the taxonomically broad est database   that had been derived from cdna libraries made from cells in different physiological conditions.
the ests matching with a corresponding predicted p <dig> gene are listed in table  <dig>  which indicates that the gene of origin is expressed and the numbers of ests found could be considered as a first indication of the relative expression abundance of that gene.
for example, ests derived from a few p <dig> genes  were found a number of times and their relatively high expression levels were also checked by semi-quantitative rt-pcr analysis, while ests from many other p450s were not identified at any of the examined cell physiological/developmental stages.
while est analyses enabled the localization of introns and demonstration of a likely case of alternative splicing, due to random fluctuation in est numbers and differences between the sizes of each library, the total est counts in different cdna libraries cannot be rigorously compared for quantifying the relative amounts of the expressed genes in different physiological states.
therefore, a more comprehensive analysis of the gene-expression patterns of the each t. thermophila p <dig> family member was done by whole-genome microarray analyses of rnas from three different cell physiological/developmental conditions .
cluster analysis of the heat map of the p <dig> expression profiles among the different conditions showed all the p <dig> genes could be clustered into three major clades on the basis of their basal expression levels and dynamic fluctuation patterns during the three different stages.
generally, three main categories of expression patterns could be postulated.
the first category contains p <dig> genes that are silent or have relatively low expression levels during all the experimental conditions and time scales .
one surprising phenomenon is that most of the p <dig> isoforms within this group come from branches having many paralogs in the phylogenetic clade, such as the cyp <dig> and cyp <dig> family.
although these genes are possibly becoming pseudogenes after relatively recent duplication events, the transcriptional activities of these genes at a specific physiological phase other than the three major ones examined here, or under some xenobiotic stimulus cannot yet be excluded.
secondly, nearly 50% of all the p <dig> isoforms, including the mostly highly expressed isoforms  are constantly expressed at all life cycle stages, although their expression levels may vary to different extents .
these genes probably take part in constitutive, endogenous physiological processes.
it should be pointed out that, due to an erroneous prediction of the tigr genome sequence that merged three adjacent p <dig> genes  into one "monster" gene when the microarray was designed; the separate microarray data for these three genes were unavailable.
however, by manually checking the signal of  <dig> of the  <dig> oligonucleotide probes whose position correctly matched the putative cyp5013a <dig> gene, its high expression level could be inferred.
the est information of cyp5013a <dig> also supports this conclusion.
the last category consists of five p <dig> genes whose expression levels varied markedly at specific time points or during one of the three physiological conditions that have been examined, including cyp5005a <dig>  cyp5005a <dig>  cyp5005a <dig>  cyp5005a <dig> and cyp5010c <dig> .
this probably indicates their involvement in some important, stage-specific endogenous cellular process.
starvation is not only a distinct physiological state that tetrahymena likely encounters in its freshwater environment, but it also induces numerous phenotypic and behavioral changes resulting in the acquisition of competence for mating  <cit> ; the nuclear events that occur in conjugating tetrahymena have clear parallels in multicellular eukaryotes and include meiosis, pronuclear formation, pronuclear fusion, postzygotic divisions, and cytoplasmic determination of nuclear fate  <cit> .
thus, further studies are necessary to correlate the expression level of specific p <dig> isoforms with stage-specific cellular processes.
meanwhile, expression studies of t. thermophila cytochrome p <dig> genes when exposed to specific chemical substances are currently under investigation.
among the p <dig> genes that are not expressed in the three major physiological/developmental states, at least three of these isoforms  are transcriptionally active under different xenobiotic stresses , suggesting that they may play a role in the metabolism and biotransformation of some chemical compounds.
evolutionary analysis of t. thermophila p <dig> genes
codon-usage analysis
codon usage bias in genes is an important evolutionary phenomenon and has been widely examined among prokaryotic and eukaryotic organisms.
firstly, since the codon bias can be shaped by preferences at the level of nucleotide, specifically the gc content levels of the coding regions, we carried out a relative neutrality plot analysis  of the t. thermophila p <dig> genes .
unlike gc3s, gc <dig> and gc <dig> are subject to functional constraints against change because a mutation at these positions usually leads to an amino acid change, except between some codons of arginine, leucine, or serine  <cit> .
the neutrality plot of the t. thermophila p <dig> genes showed that there is a significantly positive correlation  and the correlation coefficient was  <dig> , indicating that the effect on the gc contents by the intragenic gc mutation bias was similar at all three codon positions.
since both the neutral mutation and selective constraint play roles in shaping codon usage pattern of gene sequences  <cit> , this suggests that there were relatively high mutation biases or low conservation of gc content levels among the t. thermophila p <dig> genes.
secondly, if translational selection pressure influences the shaping of codon usage, the bias would show significant positive correlation to expression levels and some translation-preferred codons should be used more frequently than others.
thus, we calculated the two often used measures of codon usage bias, the cai  <cit>  and the enc  <cit>  of each t. thermophila p <dig> gene,.
among the most highly expressed p <dig> genes , three isoforms  have a cai value above  <dig> , which is significantly above that of the average  following a student's t test .
meanwhile, four p <dig> genes , although exhibiting a similar high basal expression level, showed a cai value close to, or even less than, that of the average.
thus, no simple trend exists between the expression levels of p <dig> genes and the codons they use.
meanwhile, in the enc-plot analysis, if a given gene is only subject to g+c composition mutation constraint, it will lie above or just below the standard curve.
from figure 6b, it can be seen that while the most t. thermophila p <dig> genes had enc values slightly smaller than expected ones, the three genes with the significant high cai values displayed a more biased codon usage according to the respective gc3s, indicating that they are probably under pressure from direct expression selection.
the correspondence analysis on the rscu was also performed to avoid identification of trends in codon usage due to biased amino acid usage among the genes.
rscu of t. thermophila cyp genes detected one major trend on the first axis of inertia which accounted for  <dig> % of the total variation and was approximately three times as much of the variation as the second axis , and no other axes accounted for more than  <dig> % of the total variation, respectively.
the codons of t. thermophila cyp genes that make the major contributors to this pattern showed preference for ending in c, including ccc , ucc , cuc , guc , cac  and auc , and this might indicate that codons ending in c could lead to better translational accuracy/efficiency.
besides, axis  <dig> was significantly correlated to gc and gc3s , but not to enc .
while this result implies that the nucleotide composition mutation could be considered as the relatively major factor in shaping the codon usage in p <dig> genes, it should be taken into account that the robustness of our current analysis would be limited by the relative small quantities of calculated genes.
a recent study reported that there was a very small but significant correlation between the enc and expression level  in t. thermophila  <cit> .
however, in miao et al.
<cit>  the authors found that while most codon-biased genes  are expressed constitutively and at high levels, only a fraction of highly expressed genes  are codon-biased.
this may well reflect that the genome-wide microarray analyses covered a wide range of physiological/developmental stages and are subject to less bias than non-saturated, random analyses of cdnas.
thus, although the overall codon usage biased is determined by nucleotide composition mutation in t. thermophila, the presence of a preferred gene subset is under pressure from direct expression selection which probably results from the large effective population size of this species  <cit> .
it has been reported that in the bacteria bacillus subtilis  <cit>  and chlamydia trachomatis  <cit> , the synonymous codon usage appeared to be determined by both the neutral mutational biases and translation selection.
therefore, while genes within one species often share a single codon usage pattern, exceptions also exist, especially among microorganisms.
site-specific selection of duplicate p <dig> genes
the two largest p <dig> families in t. thermophila  were chosen for investigating whether different evolutionary pressures exist on particular sites of p <dig> isoforms.
here the rates of substitutions that occurred within the coding regions of the family members were calculated for all types of selection .
different rates of substitutions were estimated to occur within the coding regions among the two family members.
the obtained ka/ks ratios were between low  to median  along the whole sequences, and were especially low among some residues that represent some core p <dig> structure regions, indicating strong negative selection on these structural regions.
this suggests they are functionally important sites in catalytic reactions and may experience restriction of deleterious mutations.
meanwhile, no positive selection was indicated at any specific site along the rest of gene coding regions, suggesting a divergent or neutral evolution process among both the cyp <dig> and cyp <dig> family members .
in ciliate species, fast protein evolution has been reported for several genes, and was suggested to be the consequence of either relaxed functional constraint on the nucleosomes of amitotic macronuclei or of adaptive evolution through gene duplication coupled with the ciliates' highly processed macronuclear genomes  <cit> .
the low expression levels of many isoforms within the cyp <dig> family, suggested they might have undergone pseudogene formation.
however, the possibility cannot be ruled out that some members of these groups serve as inducible p <dig> isoforms that function when the cell contacts specific xenobiotic compounds from the environment.
in a recent comparative genome analysis of vertebrate cytochrome p <dig> genes  <cit> , it was suggested that all of the cyp genes that encode enzymes with known endogenous substrates are phylogenetically stable, characterized by few or no gene duplications or losses.
in contrast, the unstable p <dig> genes that are characterized by frequent gene duplications and losses were found to constitute most of those that encode enzymes that function as xenobiotic detoxifiers.
besides, adaptive evolution has also been estimated to occur restrictively on those duplicated cyp <dig> genes at the amino acid sites within the srss regions in the protein structure.
this may be due to their functional association with unstable environmental interactions such as toxin and pathogen exposure.
such specific xenobiotic-driven p <dig> gene expansion events were also observed in d. melanogaster  <cit>  and c. elegans  <cit> .
the relatively high heterogeneity among the major clades within the t. thermophila p <dig> family members, along with their diverse expression patterns suggested that a similar situation might be occurring in this organism.
expression divergence of duplicate p <dig> genes
in evolving genomes, change in gene expression is one important mechanism that can lead to retention of duplicate genes, and studying the gene expression patterns has been suggested to be an important measure of gene functions that can facilitate our understanding of the genetic basis of evolutionary change  <cit> .
high-throughput approaches, such as microarray techniques, provide an opportunity to investigate gene expression of whole genomes simultaneously, allowing studies of how different genes respond to a certain environmental stimuli and the general gene expression patterns among various gene families that were categorized into different cellular functions on genome-wide scales  <cit> .
we first investigated the gene-expression evolution patterns of the p <dig> duplicate pairs in the two largest t. thermophila p <dig> gene families, cyp <dig> and cyp <dig>  by analyzing the gene-expression data from multiple microarray experiments.
for the three physiological/developmental stages, the t. thermophila cells examined, the nonphylogenetic model was the best supported model both for the two families .
for the cyp <dig> family, the best "nonphylogenetic-free" model suggested that more closely related duplicate genes are no more likely than more distantly related genes to share similar expression patterns.
this may indicate either the gene family has little influence over physiological functions or the rapid rates of gene-expression evolution  <cit> .
meanwhile, the "nonphylogenetic-distance" model best fit the cyp <dig> family, indicating that genetic distances since last gene duplication predict change in expression, consistent with an initial coupling during evolution of expression and coding sequences, i.e.
a correlation between the genetic distance in the coding region of cyp <dig> family members and the change in gene expression level still can be detected .
this was similar to a report on duplicate genes in the yeast s. cerevisiae showing that gene expression patterns remain similar shortly after gene duplication, but the evolution of expression occurs quickly so that the patterns become distinct from each other in a relatively short period of time  <cit> .
due to the observation that many cyp <dig> family members are within the "expression-silent" category, the assumed gene pseudonization events probably obscured the possible correlation between the genetic distance and gene-expression divergence in the cyp <dig> family.
further, we estimated the relative duplication time of the two gene families.
since silent substitution rates were often used as an approximation to the neutral mutation rates  <cit> , we calculated the 4-fold substitution rates  of the synonymous sites for all the duplicate pairs for each of the two gene families.
it was indicated that except for three genes  that raised by very recent duplication events , most p <dig> genes in the cyp <dig> and cyp <dig> families have uniform mutation rates .
then we calculated the amino acid distance between duplicates for using as a proxy for evolutionary time.
the results showed that the overall distance values are highly identical , thus indicating that the duplication time of the two gene families was close to each other.
therefore, given that the enlargement of these p <dig> families in t. thermophila was probably caused by recent separate small duplication events, the divergence rate of gene expression thus may also vary between the two p <dig> gene families.
since a rapid divergence of expression among the t. thermophila duplicate p <dig> genes could be inferred, we also tried to examine the 5' flanking regions of the t. thermophila duplicate p <dig> genes for identification of likely regulatory elements.
it is well known that promoters and enhancers located upstream of the coding regions usually have critical roles in regulating gene expression levels, and that major p <dig> genes are selectively induced by different nuclear receptors in response to endogenous substances or diverse xenochemicals in multicellular organisms  <cit> .
however, the fact that the tetrahymena genome has a 78% a/t content makes it difficult to identify regulatory elements by sequence homology to well-characterized promoters from other model organisms using transcriptional regulation databases such as transfac  <cit> .
currently, except for a few regulatory elements that have been identified by deletion and mutational analysis  <cit> , little information is available on the cis-elements of specific genes in tetrahymena.
we further checked possibility whether the t. thermophila p <dig> duplicate pairs or their adjacent genes in the same chromosome  have a tendency to show similar physiological/stage-specific patterns of expression.
the results showed that most of these genes exhibited wide variations in their levels and stage-specificity of expression, i.e.
no chromosomal or sub-chromosomal level of gene regulation was indicated .
therefore, although the rapid divergence of upstream non-coding sequences of the relatively well-conserved p <dig> ofrs may contribute to the remarkably varied expression patterns among the t. thermophila duplicate p <dig> genes, further experimental investigations are needed to identify specific regulatory elements or trans-acting factors involved in the transcriptional induction of p450s in this organism.
CONCLUSIONS
in the current study, we identified  <dig> putative functional p <dig> monooxygenase genes in the model ciliate organism t. thermophila, analyzed their evolutionary relationships and characterized their expression based on both est and microarray data, using bioinformatics tools.
the current microarray data provide background information of t. thermophila p <dig> gene expression in normal physiological states.
our analyses provide information on the characteristics and diversity of the p <dig> genes in the ciliophora, and will facilitate further functional analysis to understand the roles of individual p <dig> isoforms either in cellular physiological metabolism or the possible oxidative detoxification catalysis under environmental toxic exposures in this model ciliate species.

 the cytochrome p monooxygenases constitute a conserved gene superfamily of heme-thiolate proteins ubiquitously distributed in different life forms, including prokaryotes, unicellular eukaryotes and multicellular eukaryotes. they play key roles in the metabolism of a wide variety of substrates, including endogenous chemicals such as steroids and other important small molecules, and also xenobiotic compounds including drugs, pesticides and environmental contaminants. substrate and functional diversity are characteristic features of p enzymes and are considered to be the consequence of evolutionary adaptation driven by different metabolic or environmental demands in different organisms. evolution and expansion of major p branches has been suggested to be linked with the historical occurrence of important evolutionary events. one particular example is the divergence of p450s of the common plant-animal ancestor either to synthesize biochemicals/metabolites or to detoxify xenobiotics, followed by p gene expansions, especially in the plants. these may well reflect different survival strategies adopted between the two kingdoms, i.e. plants evolved sessile systems with p enzymes with more diverse and essential roles, while animals developed higher order sensory and locomotor systems, and comparatively fewer p450s. although over p450s have been named and about p genes have been characterized to date, none has been systematically characterized in the phylum of ciliophora which, together with dinoflagellates and the exclusively parasitic apicomplexa, constitute the three major evolutionary lineages that make up the alveolates. the unicellular ciliate tetrahymena thermophila is a free-living protozoan widely distributed in freshwater and estuarine environments, elaborating typical eukaryotic components into a highly organized cell whose structural and functional complexity is comparable to, or exceeds that, of human and other metazoan cells. the physiology, biochemistry and molecular biology of tetrahymena are well characterized, and it is an excellent model organism for toxicological and ecotoxicological studies in aquatic toxicity test systems. results from the est project and the macronuclear genome sequencing project have shown that, although single-celled, tetrahymena possesses core processes conserved across a wide diversity of eukaryotes that are not found in other unicellular model species such as the yeasts saccharomyces cerevisiae and schizosaccharomyces pombe. it also contains a large number of gene families that are involved in processes associated with sensing and responding to environmental cues. in the case of the cytochrome p gene family, s. cerevisiae and s. pombe have only three and two p450s, respectively, while in t. thermophila the number is more than, which is close to the typical number in a vertebrate genome. in humans and other mammals, extensive studies have focused on aspects of p gene structure and biochemical properties. important biological functions of p450s and the associated high degree of complexity in the gene polymorphism and expression patterns have been demonstrated. as genomic databases became available and functional genomics techniques such as dna-microarrays have been applied, investigations on p isoforms have also been extended to other organisms. these developments have, in turn, led to the emergence of a new field, ecotoxicogenomics, which aims to develop effective tools for identification of possible toxic environmental pollutants by characterizing their effects on terrestrial and aquatic model organisms, such as the soil-dwelling nematode caenorhabditis elegans, the aquatic crustacean daphnia magna, and the small fish fathead minnow. we previously identified a series of differentially expressed ests of t. thermophila that respond sensitively to treatment with the organochlorine insecticide ddt. one est was identified to be a p gene by homology searches, and its expression levels under different concentrations of ddt treatment were further assessed. recently, the first genome-wide microarray platform containing the predicted coding sequences has been established and validated in t. thermophila and was used to study gene expression during three major stages of the organism's life cycle: vegetative growth, nutrient starvation and conjugation. substantial progress has also been made in closure and reannotation of the mac genome sequence of this eukaryotic model organism. all these provided us the opportunity to investigate both the functional and evolutionary characteristics of the cytochrome p genes in t. thermophila at the genomic level. in this study, the putative t. thermophila p genes that were previous identified both by the international nomenclature committee and by the tigr genome annotation team were further checked by est data and through cdna sequence cloning experiments for improvement of annotation. the expression patterns of the t. thermophila p genes during three important cell physiological/developmental conditions were analyzed based on est and microarray data. these results are discussed in the context of understanding the characteristics of the t. thermophila p monooxygenase isoforms and their functional and evolutionary diversity.  t. thermophila putative p gene sequences the putative cytochrome p gene sequences in t. thermophila were initially identified in by dr. nelson based on an early tigr release of its macronuclear genome assemblies, and p450-like genes were posted on his p website. when the draft whole genome of t. thermophila was sequenced and reported in p genes were predicted by the tigr genome annotation team. in this study, we retrieved the putative t. thermophila p genes from above two locations and further searched against the tgd utilizing the updated annotation. the obtained p450-like gene sequences were thereafter compared. then, we used the data from two est database resources and our own rt-pcr investigations to verify the accuracy of the predictions. a total of putative cytochrome p gene sequences with full-length orfs were identified in t. thermophila and were assigned names as suggested by the "p nomenclature committee". our results showed that the gene predictions made by nelson are more accurate than those presented by the tigr annotation. the tigr predictions have some errors in identifying the start codons, especially the signal peptide, and intron-exon boundaries of several p450s are inaccurate. in one case, it merged adjacent p genes into one "monster" gene, which also includes two intergenic regions plus kb upstream of cyp5013a. for many other p genes, the tigr annotations are systematically short, presumably due to its failure to recognize the signal peptide at the 5' end. this problem may arise if the hydrophobic region is coded by at-rich codons, which are not easily distinguished statistically from the at-rich intergenic regions. for all the p genes, only three sequences were modified slightly in our annotation compared with the predictions made by dr. nelson. gene id, genbank accession numbers, scaffold genbank accession numbers, est accession numbers and predictions of sub-cellular localizations were listed. s: secretory pathway. *: found only in tbestdb. there were four p pseudogenes in nelson's predictions. it was revealed that one was a probable pseudogene due to the presence of in-frame stop codons, and two partial fragments with typical p sequence features appeared to be missing the n-terminus. these sequences were not included in the following analysis of the t. thermophila putative p genes. however, when checking the cdna sequences of putative t. thermophila p genes, we identified one transcription product that is most identical to the predicted pseudogene "cyp5005a7p" which has an in-frame tga codon within its orf. compared with the "cyp5005a7p" sequence, the est of the cdna transcript has several site mutations including one nucleotide transition within the in-frame stop codon found in the genomic sequence, and thus can be read through. however, a blast search of the t. thermophila whole genome database failed to detect any sequence that completely matched this cdna sequence. thus we designed a pair of primers located in either the 5' or 3' flanking region of the "cyp5005a7p" sequence, and obtained a kb pcr amplification product that was sequenced. due to the possibility that different p isoforms may exist in different t. thermophila strains, we checked the pcr products both from the genomic dna of strain sb and cu the results showed that the pcr products from the two strains were 100% identical and were also consistent with the sequence of the unexpected cdna transcript. we thus assigned the name cyp5005a to this newly observed p isoform. further, blast searches in the tbestdb and ncbi est database either found one sequence in each that perfectly matched the cyp5005a sequence. besides, the microarray data demonstrated that the corresponding cdna transcript was constantly expressed during the different cellular conditions. all the above observations serve as evidence that cyp5005a is a functional p gene in both the t. thermophila cu and sb strains. the erroneous prediction of the pseudogene "cyp5005a7p" by tigr was probably due to false assembly of genome sequence contigs, caused by the high similarity between the cyp5005a5p and cyp5005a sequences. structural features and intron-exon organization p genes are classified and annotated on the basis of amino acid sequence identity, phylogenetic homology and gene organization, and a four-digit naming system has been set up to meet the need of the increasing number of newly discovered p sequences. all the t. thermophila p sequences are distributed into families and sub-families according to the nomenclature criterion of sequence similarity, indicating p450s exist as a large gene family in this organism. the phylogenetic analysis of t. thermophila p genes obtained using either a neighbor-joining tree or a maximum-likelihood tree gave similar topological nodes, and bootstrap testing showed good reliability of the phylogenetic tree. the largest expanded p family in t. thermophila is the cyp family that contains likely functional p genes and the two predicted pseudogenes. the cyp family and the cyp family were "medium"-expanded p families in t. thermophila that possessed and genes, respectively. moreover, cyp family members are located on the same scaffold, which contains a gene cluster in the form of tandem repeats with four genes and one pseudogene in the same orientation. three cyp family members form a tandem triplicate repeat on the scaffold ch besides these, there are three other p gene pairs organized in tandem duplication. such a feature is consistent with the gene expansion strategy used by some other lower eukaryotic genomes, particularly those who have evolved to meet an extensive demand for generation of a broad range of metabolites in the secondary metabolism when interacting with the environmental niches. for example, it was revealed that the majority of the multigene p families in the model white rot fungus phanerochaete chrysosporium appear to have expanded locally in its genome as a result of extensive gene duplications and rearrangements, indicating a strong need for functional divergence in response to environmental stimuli. the intron-exon organization of p genes exhibits a diversity of gene structure with indicating that multiple gains and losses of introns have occurred during the evolution of p genes in diverse species, with little conservation of intron positions among divergent p families. intron positions were mapped and characterized for all putative t. thermophila cyp genes. of the introns that were identified in the p gene sequences, were phase zero, were phase i and were phase ii. meanwhile, based on the definition of the uip, % of the total introns were unique and the remaining introns were present in consensus locations among three different family clades. sixteen p genes lack introns entirely. five genes have the maximum number of four introns, respectively. from figure it can be observed that there is a good correlation between the conservation of intron position and phylogenetic relationships of t. thermophila p subfamily members. the evolution of introns in alveolates was recently studied by nguyen et al., who concluded that the rates of intron gain and loss were more or less constant in the last ~ myr after tetrahymena branched off. therefore, it can be inferred that the enlargement of several p families in t. thermophila resulted from recent separate small duplication events, which is also apparent within many other tetrahymena protein families containing paralogous genes. our analyses of the est data to determine the intron boundaries also revealed that the cyp5013a gene might exhibit alternative intron splicing. alternative splicing was suggested to be an uncommon phenomenon in t. thermophila, at least under the several growth conditions that have examined. cyp5013a had the greatest number of ests of all t. thermophila p isoforms. among the retrievable est sequences, represented the correct transcriptional product of cyp5013a while one retained the first intron. the different transcripts were further investigated by rt-pcr. an intron-containing transcript was observed in pcr products both from the cdna templates of starvation and conjugation stages, but was barely detectable from vegetative growing cells. the retained intron in cyp5013a is bp in length and belongs to phase i type. thus, the putative intron-containing transcript should contain a frame-shift mutation. rna transcripts carrying such premature stop codons can be prevented by the nmd pathway. recently, this pathway was found to play an essential role in another ciliate, p. tetraurelia, by knocking down the two homologous genes encoding upf a protein which is crucial in nmd, thus indicating a universal translational control of intron splicing in eukaryotes based on nmd surveillance. a blast search of the t. thermophila genome using p. tetraurelia upf genes as queries revealed the existence of one homolog with 70% identity. thus, intron retention in the cyp5013a transcript might be caused by inefficient nmd activity during specific physiological/developmental stages. however, it would lead to a mature "non-sense" transcript and could only serve to down-regulate expression of the cyp5013a gene. whether such rarity of alternative splicing in this species may simply be tolerated or directed by other mechanisms, is currently unknown. all known p450s appear to take on a similar folded structure, yet frequently show less than 30% sequence identity and have very different substrates. in t. thermophila, the average length of the p amino acid sequences in the alignment is around amino acid residues, and the genes share different degrees of amino acid sequence conservation between different families. to better understand the structure and functional relationship of p enzymes, mammalian p cyp3a seven t. thermophila p sequences and one paramecium tetraurelia p sequence were selected and aligned. secondary elements were then assigned using cyp3a as the template for secondary structure, based on its known crystal structures and six putative srss were identified in the aligned sequences according to gotoh's predicted models. srss participate in the contact with ligands and many of them locate at regions with variable structural elements. all the t. thermophila p sequences contain the typical conserved p domains, including the heme-binding region near the c-terminal, the perf domain and the k-helix. multiple sequence alignment also revealed that some residues are identical across t. thermophila p450s; most of these are distributed in the five regions that constitute the conserved p structural core, including: the c helix, the c-terminal part of the e helix, the i helix, the j helix, and the heme-binding loop between helix k and l. t. thermophila p member composition prediction of sub-cellular localization of the cyp proteins showed that no predicted mitochondrial localized p was found in t. thermophila, and all of the genes except cyp5001a encode a putative typical microsomal signal peptide of about hydrophobic residues that likely serves as a membrane anchor in the endoplasmic reticulum. the secretory pathway of eukaryotic cells consists of a series of discrete, membrane-bound compartments, including the er, golgi, and vacuoles, all of which are present in tetrahymena. for cyp5001a the predicted signal sequence was ambiguous and also lacked the positively charged residues in this n-terminus critical for endocellular transportation targeting. so far, mitochondrial p450s have been described only in the animal kingdom. as with t. thermophila, none was found in fungi and plants. these observations support the suggestion that mitochondrial p450s did not originate from the ancestral mitochondrial endosymbiont, but evolved later, possibly by mistargeting of microsomal p450s. cyp is the only p family that has orthologs in multiple phyla of the animal, plant, fungal and bacterial kingdoms, and it has been postulated to be the ancestral p isoform. this family is functionally conserved and has a very limited range of substrates related to the biosynthesis of sterols and their derivatives, and a high conservation within the srss is a specific feature of cyp it was believed that two motifs, "yxxf/lxxpxfgxxvxf/yd/a" and "gq/hht/ss", present within the regions of srs and srs respectively, can be considered as the cyp signature. in this study, no cyp51-like gene has been identified with such a signature in the t. thermophila p family. so far, all studied cyp family members were found to exhibit catalytic function in the oxidative removal of the 14α-methyl group from sterol precursors formed downstream from cyclization of squalene -epoxide, and the 14α-demethylated products are intermediates in the sterol biosynthetic pathways leading to formation of cholesterol in animals, ergosterol in fungi and a variety of functional sterols in plants and algae. cyp appears to be absent in insects and nematodes due to the fact that they don't make cholesterol and lack the post-squalene sterol biosynthesis pathway, and absence of the cyp gene was thought to be the consequence of gene loss during evolution. thus, a similar cyp gene loss may be inferred in t. thermophila. instead of the sterols in most other eukaryotes, tetrahymena produces the pentacyclic triterpenoid alcohol'tetrahymanol' and/or hopanoids de novo as functionally equivalent structural components of cell membranes, and its ability to synthesize and use this "primitive" substance can be considered as a metabolic relic. besides, although t. thermophila is able to incorporate exogenous sterols into cell membranes and convert those to various derivates, it utilize enzymes other than p450s in the downstream metabolism of certain steroid compounds, such as the cytochrome b in the direct desaturation reaction of cholesterol. the ciliate p. tetraurelia is the species closest to t. thermophila that has an available genome database and phylogenetic analyses between these two species can help to understand the evolutionary relationship of t. thermophila p families with other ciliate members. a phylogenetic network was constructed of both the t. thermophila and p. tetraurelia p protein sequences. in the branch that constitutes three p. tetraurelia p genes, two genes from t. thermophila also present and cluster together. however, according to the deep branching of these sequences and the relatively low bootstrap values on the branch nodes, there is yet no strong evidence supporting the possibility that they represent orthologs shared by the two species. surprisingly, the other major clades of p genes seem to be highly specific to t. thermophila or to p. tetraurelia. moreover, the number of p genes in p. tetraurelia is less than half of that in t. thermophila, and blast searches in different est databases revealed that only about 1/ of the p. tetraurelia p450s have retrievable est data. it has been pointed out that while tetrahymena and paramecium are related species within the diverse clade of ciliates, their genomes have obviously evolved by significantly different mechanisms. the paramecium genome has gone through a series of whole genome duplications and may exhibit the most completely preserved whole-genome duplication described to date, where ~95% of the genome is located in duplicated blocks. by contrast, although up to 40% of the t. thermophila genome was also indicated to be located within duplicated blocks, the most of the large expansions of its gene families may have arisen through tandem duplication events. while differences in the p genes in the two species may reflect the difference in evolutionary strategies between the two organisms since their divergence, further investigations are needed to address the adaptive significance of these diversities relative to the well known functions of cytochrome p450s in sensing of and responding to environmental changes, particularly environmental stress. p gene expression analysis spatial and temporal gene expression patterns are important aspects of gene regulation and expression pattern analysis has played an important role in the study of the cytochrome p gene superfamily. in the present study, the gene expression patterns of all t. thermophila p450s during three important cell physiological stages were analyzed based on est and microarray data. to examine the expression profiles of the p genes, we initially retrieved the est data of t. thermophila from both in the ncbi est database and the taxonomically broad est database that had been derived from cdna libraries made from cells in different physiological conditions. the ests matching with a corresponding predicted p gene are listed in table which indicates that the gene of origin is expressed and the numbers of ests found could be considered as a first indication of the relative expression abundance of that gene. for example, ests derived from a few p genes were found a number of times and their relatively high expression levels were also checked by semi-quantitative rt-pcr analysis, while ests from many other p450s were not identified at any of the examined cell physiological/developmental stages. while est analyses enabled the localization of introns and demonstration of a likely case of alternative splicing, due to random fluctuation in est numbers and differences between the sizes of each library, the total est counts in different cdna libraries cannot be rigorously compared for quantifying the relative amounts of the expressed genes in different physiological states. therefore, a more comprehensive analysis of the gene-expression patterns of the each t. thermophila p family member was done by whole-genome microarray analyses of rnas from three different cell physiological/developmental conditions. cluster analysis of the heat map of the p expression profiles among the different conditions showed all the p genes could be clustered into three major clades on the basis of their basal expression levels and dynamic fluctuation patterns during the three different stages. generally, three main categories of expression patterns could be postulated. the first category contains p genes that are silent or have relatively low expression levels during all the experimental conditions and time scales. one surprising phenomenon is that most of the p isoforms within this group come from branches having many paralogs in the phylogenetic clade, such as the cyp and cyp family. although these genes are possibly becoming pseudogenes after relatively recent duplication events, the transcriptional activities of these genes at a specific physiological phase other than the three major ones examined here, or under some xenobiotic stimulus cannot yet be excluded. secondly, nearly 50% of all the p isoforms, including the mostly highly expressed isoforms are constantly expressed at all life cycle stages, although their expression levels may vary to different extents. these genes probably take part in constitutive, endogenous physiological processes. it should be pointed out that, due to an erroneous prediction of the tigr genome sequence that merged three adjacent p genes into one "monster" gene when the microarray was designed; the separate microarray data for these three genes were unavailable. however, by manually checking the signal of of the oligonucleotide probes whose position correctly matched the putative cyp5013a gene, its high expression level could be inferred. the est information of cyp5013a also supports this conclusion. the last category consists of five p genes whose expression levels varied markedly at specific time points or during one of the three physiological conditions that have been examined, including cyp5005a cyp5005a cyp5005a cyp5005a and cyp5010c. this probably indicates their involvement in some important, stage-specific endogenous cellular process. starvation is not only a distinct physiological state that tetrahymena likely encounters in its freshwater environment, but it also induces numerous phenotypic and behavioral changes resulting in the acquisition of competence for mating; the nuclear events that occur in conjugating tetrahymena have clear parallels in multicellular eukaryotes and include meiosis, pronuclear formation, pronuclear fusion, postzygotic divisions, and cytoplasmic determination of nuclear fate. thus, further studies are necessary to correlate the expression level of specific p isoforms with stage-specific cellular processes. meanwhile, expression studies of t. thermophila cytochrome p genes when exposed to specific chemical substances are currently under investigation. among the p genes that are not expressed in the three major physiological/developmental states, at least three of these isoforms are transcriptionally active under different xenobiotic stresses, suggesting that they may play a role in the metabolism and biotransformation of some chemical compounds. evolutionary analysis of t. thermophila p genes codon-usage analysis codon usage bias in genes is an important evolutionary phenomenon and has been widely examined among prokaryotic and eukaryotic organisms. firstly, since the codon bias can be shaped by preferences at the level of nucleotide, specifically the gc content levels of the coding regions, we carried out a relative neutrality plot analysis of the t. thermophila p genes. unlike gc3s, gc and gc are subject to functional constraints against change because a mutation at these positions usually leads to an amino acid change, except between some codons of arginine, leucine, or serine. the neutrality plot of the t. thermophila p genes showed that there is a significantly positive correlation and the correlation coefficient was, indicating that the effect on the gc contents by the intragenic gc mutation bias was similar at all three codon positions. since both the neutral mutation and selective constraint play roles in shaping codon usage pattern of gene sequences, this suggests that there were relatively high mutation biases or low conservation of gc content levels among the t. thermophila p genes. secondly, if translational selection pressure influences the shaping of codon usage, the bias would show significant positive correlation to expression levels and some translation-preferred codons should be used more frequently than others. thus, we calculated the two often used measures of codon usage bias, the cai and the enc of each t. thermophila p gene,. among the most highly expressed p genes, three isoforms have a cai value above, which is significantly above that of the average following a student's t test. meanwhile, four p genes, although exhibiting a similar high basal expression level, showed a cai value close to, or even less than, that of the average. thus, no simple trend exists between the expression levels of p genes and the codons they use. meanwhile, in the enc-plot analysis, if a given gene is only subject to g+c composition mutation constraint, it will lie above or just below the standard curve. from figure 6b, it can be seen that while the most t. thermophila p genes had enc values slightly smaller than expected ones, the three genes with the significant high cai values displayed a more biased codon usage according to the respective gc3s, indicating that they are probably under pressure from direct expression selection. the correspondence analysis on the rscu was also performed to avoid identification of trends in codon usage due to biased amino acid usage among the genes. rscu of t. thermophila cyp genes detected one major trend on the first axis of inertia which accounted for % of the total variation and was approximately three times as much of the variation as the second axis, and no other axes accounted for more than % of the total variation, respectively. the codons of t. thermophila cyp genes that make the major contributors to this pattern showed preference for ending in c, including ccc, ucc, cuc, guc, cac and auc, and this might indicate that codons ending in c could lead to better translational accuracy/efficiency. besides, axis was significantly correlated to gc and gc3s, but not to enc. while this result implies that the nucleotide composition mutation could be considered as the relatively major factor in shaping the codon usage in p genes, it should be taken into account that the robustness of our current analysis would be limited by the relative small quantities of calculated genes. a recent study reported that there was a very small but significant correlation between the enc and expression level in t. thermophila. however, in miao et al. the authors found that while most codon-biased genes are expressed constitutively and at high levels, only a fraction of highly expressed genes are codon-biased. this may well reflect that the genome-wide microarray analyses covered a wide range of physiological/developmental stages and are subject to less bias than non-saturated, random analyses of cdnas. thus, although the overall codon usage biased is determined by nucleotide composition mutation in t. thermophila, the presence of a preferred gene subset is under pressure from direct expression selection which probably results from the large effective population size of this species. it has been reported that in the bacteria bacillus subtilis and chlamydia trachomatis, the synonymous codon usage appeared to be determined by both the neutral mutational biases and translation selection. therefore, while genes within one species often share a single codon usage pattern, exceptions also exist, especially among microorganisms. site-specific selection of duplicate p genes the two largest p families in t. thermophila were chosen for investigating whether different evolutionary pressures exist on particular sites of p isoforms. here the rates of substitutions that occurred within the coding regions of the family members were calculated for all types of selection. different rates of substitutions were estimated to occur within the coding regions among the two family members. the obtained ka/ks ratios were between low to median along the whole sequences, and were especially low among some residues that represent some core p structure regions, indicating strong negative selection on these structural regions. this suggests they are functionally important sites in catalytic reactions and may experience restriction of deleterious mutations. meanwhile, no positive selection was indicated at any specific site along the rest of gene coding regions, suggesting a divergent or neutral evolution process among both the cyp and cyp family members. in ciliate species, fast protein evolution has been reported for several genes, and was suggested to be the consequence of either relaxed functional constraint on the nucleosomes of amitotic macronuclei or of adaptive evolution through gene duplication coupled with the ciliates' highly processed macronuclear genomes. the low expression levels of many isoforms within the cyp family, suggested they might have undergone pseudogene formation. however, the possibility cannot be ruled out that some members of these groups serve as inducible p isoforms that function when the cell contacts specific xenobiotic compounds from the environment. in a recent comparative genome analysis of vertebrate cytochrome p genes, it was suggested that all of the cyp genes that encode enzymes with known endogenous substrates are phylogenetically stable, characterized by few or no gene duplications or losses. in contrast, the unstable p genes that are characterized by frequent gene duplications and losses were found to constitute most of those that encode enzymes that function as xenobiotic detoxifiers. besides, adaptive evolution has also been estimated to occur restrictively on those duplicated cyp genes at the amino acid sites within the srss regions in the protein structure. this may be due to their functional association with unstable environmental interactions such as toxin and pathogen exposure. such specific xenobiotic-driven p gene expansion events were also observed in d. melanogaster and c. elegans. the relatively high heterogeneity among the major clades within the t. thermophila p family members, along with their diverse expression patterns suggested that a similar situation might be occurring in this organism. expression divergence of duplicate p genes in evolving genomes, change in gene expression is one important mechanism that can lead to retention of duplicate genes, and studying the gene expression patterns has been suggested to be an important measure of gene functions that can facilitate our understanding of the genetic basis of evolutionary change. high-throughput approaches, such as microarray techniques, provide an opportunity to investigate gene expression of whole genomes simultaneously, allowing studies of how different genes respond to a certain environmental stimuli and the general gene expression patterns among various gene families that were categorized into different cellular functions on genome-wide scales. we first investigated the gene-expression evolution patterns of the p duplicate pairs in the two largest t. thermophila p gene families, cyp and cyp by analyzing the gene-expression data from multiple microarray experiments. for the three physiological/developmental stages, the t. thermophila cells examined, the nonphylogenetic model was the best supported model both for the two families. for the cyp family, the best "nonphylogenetic-free" model suggested that more closely related duplicate genes are no more likely than more distantly related genes to share similar expression patterns. this may indicate either the gene family has little influence over physiological functions or the rapid rates of gene-expression evolution. meanwhile, the "nonphylogenetic-distance" model best fit the cyp family, indicating that genetic distances since last gene duplication predict change in expression, consistent with an initial coupling during evolution of expression and coding sequences, i.e. a correlation between the genetic distance in the coding region of cyp family members and the change in gene expression level still can be detected. this was similar to a report on duplicate genes in the yeast s. cerevisiae showing that gene expression patterns remain similar shortly after gene duplication, but the evolution of expression occurs quickly so that the patterns become distinct from each other in a relatively short period of time. due to the observation that many cyp family members are within the "expression-silent" category, the assumed gene pseudonization events probably obscured the possible correlation between the genetic distance and gene-expression divergence in the cyp family. further, we estimated the relative duplication time of the two gene families. since silent substitution rates were often used as an approximation to the neutral mutation rates, we calculated the 4-fold substitution rates of the synonymous sites for all the duplicate pairs for each of the two gene families. it was indicated that except for three genes that raised by very recent duplication events, most p genes in the cyp and cyp families have uniform mutation rates. then we calculated the amino acid distance between duplicates for using as a proxy for evolutionary time. the results showed that the overall distance values are highly identical, thus indicating that the duplication time of the two gene families was close to each other. therefore, given that the enlargement of these p families in t. thermophila was probably caused by recent separate small duplication events, the divergence rate of gene expression thus may also vary between the two p gene families. since a rapid divergence of expression among the t. thermophila duplicate p genes could be inferred, we also tried to examine the 5' flanking regions of the t. thermophila duplicate p genes for identification of likely regulatory elements. it is well known that promoters and enhancers located upstream of the coding regions usually have critical roles in regulating gene expression levels, and that major p genes are selectively induced by different nuclear receptors in response to endogenous substances or diverse xenochemicals in multicellular organisms. however, the fact that the tetrahymena genome has a 78% a/t content makes it difficult to identify regulatory elements by sequence homology to well-characterized promoters from other model organisms using transcriptional regulation databases such as transfac. currently, except for a few regulatory elements that have been identified by deletion and mutational analysis, little information is available on the cis-elements of specific genes in tetrahymena. we further checked possibility whether the t. thermophila p duplicate pairs or their adjacent genes in the same chromosome have a tendency to show similar physiological/stage-specific patterns of expression. the results showed that most of these genes exhibited wide variations in their levels and stage-specificity of expression, i.e. no chromosomal or sub-chromosomal level of gene regulation was indicated. therefore, although the rapid divergence of upstream non-coding sequences of the relatively well-conserved p ofrs may contribute to the remarkably varied expression patterns among the t. thermophila duplicate p genes, further experimental investigations are needed to identify specific regulatory elements or trans-acting factors involved in the transcriptional induction of p450s in this organism.  in the current study, we identified putative functional p monooxygenase genes in the model ciliate organism t. thermophila, analyzed their evolutionary relationships and characterized their expression based on both est and microarray data, using bioinformatics tools. the current microarray data provide background information of t. thermophila p gene expression in normal physiological states. our analyses provide information on the characteristics and diversity of the p genes in the ciliophora, and will facilitate further functional analysis to understand the roles of individual p isoforms either in cellular physiological metabolism or the possible oxidative detoxification catalysis under environmental toxic exposures in this model ciliate species. 
gene expression patterns of all t. thermophila p450s during three important cell physiological stages were analyzed based on est and microarray data, and three main categories of expression patterns were postulated. t. thermophila possess many advantages as a eukaryotic model organism and it exhibits rapid and sensitive responses to xenobiotics, making it an ideal model system to study the evolutionary and functional diversity of the p monooxygenase gene family. cytochrome p monooxygenases play key roles in the metabolism of a wide variety of substrates and they are closely associated with endocellular physiological processes or detoxification metabolism under environmental exposure. evolutionary analysis including codon usage preference, site-specific selection and gene-expression evolution patterns were investigated and the results indicated remarkable divergences among the t. thermophila p genes. the characteristics of both the conserved intron-exon organization and scaffold localization of tandem repeats within each p family clade suggested that the enlargement of t. thermophila p families probably resulted from recent separate small duplication events. the characterization, expression and evolutionary analysis of t. thermophila p monooxygenase genes in the current study provides useful information for understanding the characteristics and diversities of the p genes in the ciliophora, and provides the baseline for functional analyses of individual p isoforms in this model ciliate species.
 cytochrome p monooxygenases play key roles in the metabolism of a wide variety of substrates and they are closely associated with endocellular physiological processes or detoxification metabolism under environmental exposure. to date, however, none has been systematically characterized in the phylum ciliophora. t. thermophila possess many advantages as a eukaryotic model organism and it exhibits rapid and sensitive responses to xenobiotics, making it an ideal model system to study the evolutionary and functional diversity of the p monooxygenase gene family.  a total of putative functional cytochrome p genes were identified and could be classified into families and sub-families according to standard nomenclature. the characteristics of both the conserved intron-exon organization and scaffold localization of tandem repeats within each p family clade suggested that the enlargement of t. thermophila p families probably resulted from recent separate small duplication events. gene expression patterns of all t. thermophila p450s during three important cell physiological stages were analyzed based on est and microarray data, and three main categories of expression patterns were postulated. evolutionary analysis including codon usage preference, site-specific selection and gene-expression evolution patterns were investigated and the results indicated remarkable divergences among the t. thermophila p genes.  the characterization, expression and evolutionary analysis of t. thermophila p monooxygenase genes in the current study provides useful information for understanding the characteristics and diversities of the p genes in the ciliophora, and provides the baseline for functional analyses of individual p isoforms in this model ciliate species. 