Empty Nets? Fishing for Hard Facts on Pfiesteria

Background: Models of codon evolution have proven useful for investigating the strength and direction of natural selection. In some cases, a priori biological knowledge has been used successfully to model heterogeneous evolutionary dynamics among codon sites. These are called fixed-effect models, and they require that all codon sites are assigned to one of several partitions which are permitted to have independent parameters for selection pressure, evolutionary rate, transition to transversion ratio or codon frequencies. For single gene analysis, partitions might be defined according to protein tertiary structure, and for multiple gene analysis partitions might be defined according to a gene's functional category. Given a set of related fixed-effect models, the task of selecting the model that best fits the data is not trivial. Results: In this study, we implement a set of fixed-effect codon models which allow for different levels of heterogeneity among partitions in the substitution process. We describe strategies for selecting among these models by a backward elimination procedure, Akaike information criterion (AIC) or a corrected Akaike information criterion (AICc). We evaluate the performance of these model selection methods via a simulation study, and make several recommendations for real data analysis. Our simulation study indicates that the backward elimination procedure can provide a reliable method for model selection in this setting. We also demonstrate the utility of these models by application to a single-gene dataset partitioned according to tertiary structure (abalone sperm lysin), and a multi-gene dataset partitioned according to the functional category of the gene (flagellar-related proteins of Listeria). Conclusion: Fixed-effect models have advantages and disadvantages. Fixed-effect models are desirable when data partitions are known to exhibit significant heterogeneity or when a statistical test of such heterogeneity is desired. They have the disadvantage of requiring a priori knowledge for partitioning sites. We recommend: (i) selection of models by using backward elimination rather than AIC or AICc, (ii) use a stringent cut-off, e.g., p = 0.0001, and (iii) conduct sensitivity analysis of results. With thoughtful application, fixed-effect codon models should provide a useful tool for large scale multi-gene analyses. from First International Conference on Phylogenomics Sainte-Adèle, Québec, Canada. 15–19 March, 2006 Published: 8 February 2007 BMC Evolutionary Biology 2007, 7(Suppl 1):S5 doi:10.1186/1471-2148-7-S1-S5 <supplement> <title>
First International Conference on Phylogenomics
 </title> <editor>Hervé Philippe, Mathieu Blanchette</editor> <note>Proceedings</note> </supplement> © 2007 Bao et al; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

(page number not for citation purposes)

Background
The ratio d N /d S (ω) has proven a valuable index of the strength and direction of selection pressure. Because genetic data are typically subject to a diversity of evolutionary constraints, estimating ω as an average over many sites diminishes the effectiveness of this approach [1]. Statistical power is substantially improved, however, by accommodating variable selection pressures among sites (e.g., [2][3][4]). We follow Kosakovsky Pond and Frost [5] by placing such methods in three groups: (i) the counting methods, which estimate ω from counts of substitutions at individual sites (e.g., [3]), (ii) the random-effect models, which assume a parametric distribution of variability in the ω ratio across sites (e.g., [2]), and (iii) the fixedeffect models, which assume sites can be assigned a priori to different partitions [4]. The most generalized form of the fixed-effect models treats each site as an independent partition [5,6].
The recent growth of genome scale sequencing projects has sparked interest in using codon models to study mechanisms of innovation and functional divergence in genome-scale datasets [7]. Although the fixed-effect models were developed for analysis of multiple partitions of sites within a single gene, they are also appropriate for joint analyses of multi-gene datasets [4,8]. Fixed-effect models categorize codon sites into different classes which are allowed to have heterogeneous evolutionary dynamics, and such partitions are easily delineated on the basis of complete gene sequences. Moreover, by partitioning genes according to criteria such as their functional category, or role in a metabolic pathway, the fixed-effect models provide a statistical framework for making use of such information when analysing multi-gene datasets.
Yang and Swanson [4] introduced six fixed-effect models (Table 1) based on the codon model of Goldman and Yang [9]. The simplest model (A) assumes that the pattern of substitution is homogeneous over all sites; i.e., there are no partitions under model A. Branch lengths are included as parameters of the model. The most complex model (F) treats the different site partitions as independent datasets, having independent model parameters. As it involves a substantial increase in branch length parameters, model F is not recommended for datasets with many partitions [4]. The remaining four models (B-E in Table 1) lie between A and F in complexity. These four models scale the branch lengths of k partitions according to the parameter c k , which is a multiple of the branch lengths of the first partition; hence c 1 = 1. Models B through E differ in their treatment of parameters ω, κ (transition to transversion ratio) and π (codon frequencies) among partitions (Table 1).
We implemented 11 more fixed-effect models, which represent all the remaining combinations of heterogeneity or homogeneity among partitions for the parameters c, ω, κ and π ( Table 2). A full description of the fixed effect models and the details of our implementation are presented in the methods section. Hereafter we refer to the complete set of fixed-effect (FE) models by using the revised notation shown in table 2 (FE1-FE16). Note that a capacity to specify fixed-effects under the alternative formulation of Muse and Gaut [10] is available through the program HyPhy [11], although it has not been documented.
Given a related set of fixed-effect models (Figure 1), one is immediately faced with the non-trivial task of selecting the model that best fits the data in hand. Likelihood ratio tests (LRTs) have been shown to be a powerful and reliable means of testing site specific heterogeneity in selective pressure [8,12]. However, Figure 1 illustrates that there are 32 possible nested comparisons of models. It is not desirable to conduct 32 LRTs because computational costs are expensive for datasets with too many sequences or partitions. A popular method of model selection based on LRTs is "backward elimination" [13][14][15]. Backward elimination reduces a comparatively complex model to a simpler one in a step-by-step fashion. An alternative to "backward elimination" is the Akaike Information Criterion (AIC) [16], where the model with the smallest AIC score is chosen as the ideal model. For a small sample correction, typically when the number of observations is less than 40 times the number of parameters in the model [17], we borrowed the corrected Akaike Information Cri-
(page number not for citation purposes) terion (AICc) originally developed by Hurvich and Tsai [17] for regression settings.
Although the statistical issues surrounding model selection are well known within the field of molecular evolution [18][19][20], the established statistical techniques have not been applied to the fixed-effect setting. In this study we used computer simulation to evaluate the performance of backward elimination, AIC and AICc for selecting an optimal model from an array of models specifying different levels of heterogeneity among partitions. We then illustrated the application of these methods on two real datasets. The first was comprised of the buried and exposed sites of the abalone sperm lysin gene; this lysin partition was one of the original test cases of Yang and Swanson [4]. For the second case, we examined the evolutionary heterogeneity of a multi-gene dataset; the region of the genome encoding all the components of the flagellar system of Listeria species and several proteins of unknown function.

Simulated data
A simulation study was used to measure the accuracy of fixed-effect model selection. We simulated under the 16 different scenarios for heterogeneous codon evolution among data partitions shown in table 2 (see methods for a detailed description of the simulation study), and measured the number of cases where each procedure identified the correct generating model. The backward elimination procedure uses the likelihood ratio test (LRT) to simplify a complex model one parameter at a time; in this case we start at the top of Figure 1 (FE1) and use the LRT to remove non-significant parameters in a step-wise fashion. A more detailed description is provided in the methods. When we applied the LRT under a cut-off probability of 0.05, the backward elimination procedure provided more accurate model specification than either AIC or AICc in all the cases except for model 2 (Table 3). Among all 336 datasets, the accuracy of backward elimination was 78% whereas the accuracy of AIC and AICc was 63% and 64% respectively ( Table 3). Note that each model can be related to all other models by the number of connections, or "steps", between them in Figure 1. For all models that are wrongly specified by backward elimination, most were just one step away from the true model (85%). Among these 1-step wrong models, there was a bias in the direction of greater complexity for one of ω, κ or c; replicates heterogeneous for these parameters were never misclassified as homogenous. Taken over all replicates homogenous for ω, κ or c, this bias was generally low, with 1-step error rates of 13%, 9% and 9% respectively. Similar results were observed for AIC and AICc. Most misclassifications were 1-step errors (77% and 78%), with a bias in the direction of greater complexity for parameter ω, κ or c. Again, heterogeneous replicates were not misclassified as homogenous for these parameters. The 1-step error rates across replicates homogenous for ω, κ or c were 28%, 18% and 26% for AIC, and 26%, 17% and 22% for AICc.
Number of parameters is for the F3 × 4 method of estimating codon frequencies. b and g denote the number of branches and the number of site classes, respectively. Letters in parentheses indicate the model codes formerly used by Yang and Swanson [4].
As the number of misclassification errors ≥2-steps was much smaller than the number of 1-step errors, we examined these as an average over backward elimination, AIC and AICc. In 90% of the cases these errors resulted in too simple a model. The involved parameters were ω, c and π; the κ parameter was rarely misclassified.
We simulated under two models of codon frequencies: (i) unbiased (π i = 1/61) and (ii) biased frequencies taken from empirical frequencies of the lysin gene. In composite datasets with a 90:10 partition the number of codons in the smaller partition is too low for reliable empirical estimation of 61 different codon frequencies ("F61" method). Hence, in only those cases we used the "F3 × 4" method, which computes codon frequencies from nucleotide frequencies at the three positions of the codon [9]. In the 50:50 and 70:30 datasets we used the empirical estimates of codon frequencies (F61) in each partition. We note that such empirical estimates do not satisfy the requirements of LRT [21] and, hence, the backward elimination procedure. For backward elimination the 1-step error rate for incorrectly specifying heterogeneous π was Relationships among fixed-effect codon models Figure 1 Relationships among fixed-effect codon models. The most complex model (FE1) is located at the top and a completely homogenous model (FE16) is located at the bottom. Parameters heterogeneous among partitions for a given model are shown after the model name. Lines between models indicate "1-step" differences in complexity among the models. 6%, and for incorrectly specifying homogenous π was 14%, indicating a greater tendency towards too simple a model. For AIC and AICc, the misspecification of π was almost entirely for too simple a model. Note that most of these errors were made in the 90:10 datasets, suggesting that misspecification of codon frequency heterogeneity is mainly due to large empirical estimation-errors of codon frequencies due to the insufficient information of small partitions. Thus power is lowest to identify heterogeneity in codon frequencies when a partition consists of a small number of codon sites. We anticipate that power also will be low in larger partitions of real datasets where the difference among partitions is not as great as in our simulations.

Model 2 c
Next we investigated the possibility of tuning the cut-off p-value of the backward elimination procedure to improve the accuracy of model specification. We evaluated accuracy for cut-off p-values of 0.01, 0.001 and 0.0001. Substantial improvements were obtained, with average accuracy increasing from 78% (under the original cut-off p-value of 0.05) to 83%, 87% and 88% (Table 3) respectively. Under a cut-off value of 0.0001, 39 models were misspecified, with 33 being too simple with respect to codon frequencies in the 90:10 datasets. Among the 6 remaining misspecified models, 4 were too complex for ω and 2 were too complex for π. All but one of the misspecified models under cut-off p = 0.0001 were one-step errors. Based on these findings we used a cut-off p = 0.0001 in our application of these models to real data.

Abalone sperm lysin gene
Abalone sperm lysin is a reproductive protein well known for rapid evolution under strong diversifying selection [22]. We partitioned the lysin dataset into the same set of 46 buried sites and 88 solvent-exposed sites as in Yang and Swanson [4] and applied backward elimination (cutoff p = 0.0001), AIC and AICc to select among the full set of fixed-effect models. Under backward elimination, we used the likelihood ratio test to compare FE1, which assumes different κ, ω, c, and π's for buried and exposed sites, with those nested models at the next level in Figure  1 (FE2, FE3, FE5 and FE9). Each model at the next level assumes one of these four parameters (κ, ω, c or π) is homogeneous among site classes. FE9 assumes homogenous κ for both buried and solvent-exposed sites, and the likelihood ratio statistic comparing FE1 against FE9 is 2 × (4474.75-4473.88) = 1.74, which is not significant (d.f. = 1; p = 0.1871). As all other LRTs at this level are significant, we simplified our model according to κ and compared FE9 to those models nested at the next level in Figure 1 (FE10, FE11 and FE13). As subsequent LRTs involving FE9 were significant, model FE9 was selected by backward elimination. We note that even when we use a cut-off p = 0.05, we still select FE9 by backward elimination. Table 4 illustrates that FE9 is also selected by using AIC or AICc.
Yang and Swanson [4] conducted LRTs of the subset of models shown in Table 1 and found that Model E (FE1) provided the best fit to the lysin data. FE1 and FE9 are qualitatively similar in suggesting heterogeneity in ω, c and π's among the buried (b) and solvent exposed (e) sites. Moreover, these models provide similar quantitative estimates of the strength of selection and rate of evolution in these two partitions (FE1: ω b = 0.45, ω e = 1.28, c b = 1, c e = 2.54; FE9: ω b = 0.43, ω e = 1.31, c b = 1, c e = 2.56). Both models suggest that buried sites are evolving under strong purifying selection and exposed sites under diversifying selection. Note that Yang and Swanson [4] used a model that specified heterogeneous κ's (FE1), as it was not possible to test for heterogeneity in κ's independently of ω's ( Table 1). As estimates of ω were very similar under FE1 and FE9, and κ's for partitions under FE1 were very similar (κ b = 1.9 and κ e = 1.5), the use of FE1 was not problematic in this case.

Components of the Listeria flagellar system
Listeria are gram positive rod shaped bacteria which are motile between 4°C and 30°C, and grow in a wide range of pHs, temperatures, and osmotic pressures. The natural habitat of Listeria is thought to be soil rich in decaying matter; however, Listeria monocytogenes is an important food-borne pathogen of humans and animals capable of both a free-living and intracellular lifestyle. Interestingly, the motility of Listeria monocytogenes is thermoregulated, being reduced above 30°C, and completely shut down above 37°C [23], temperatures which correspond to their host intracellular environment. The consensus opinion is that the shut down of expression of flagellar related proteins, thereby shutting down motility, is an adaptation to avoid recognition by the hosts innate immune system [24]. Specifically, recognition of the flagellin protein, a product of the flaA gene, activates the host inflammatory responses through Toll-like receptor 5 (TLR5) [25].
Twenty-eight genes encoding putative flagellar related proteins, including flaA, are located together in the genome of Listeria. Several proteins having functions unrelated to flagellar machinery, or unknown functions, also are encoded in this region of the genome. We ana- Parameter estimates Parameters in parentheses are fixed. Partition 1 contained the buried sites and partition 2 contained the solvent exposed sites. ΔAIC i = AIC i -min AIC, and ΔAIC Ci = AIC Ci -min AIC C .
lysed the genes from this region with three issues in mind: (i) to test for heterogeneous evolutionary dynamics among genes, (ii) to examine the evolutionary dynamics of proteins with unknown function and determine if they have any similarities to flagellar machinery proteins or proteins having unrelated functions, and (iii) to compare selection pressures on flaA with other genes known to encode flagellar related proteins. We note that thermoregulation of motility is not always perfect, thereby raising the possibility that the host's innate immune system is occasionally able to recognize flagellin [24]. This would set up selection pressure for a "co-evolutionary chase" between host and pathogen, leading to an elevated rate of nonsynonymous evolution in flaA.
Our dataset was comprised of 37 of the 43 genes (lmo0675 -lmo0717) located contiguously within the genomes of 5 lineages of Listeria. Two genes (lmo0684 and lmo0711) were excluded because their gene trees were incompatible with the genome-tree. Four genes (lmo0677, lmo0698, lmo0709 and lmo0712) were excluded because they were less than 100 codons long. Next we partitioned the genes according to functional category: (i) flagellar machinery (7973 codons), (ii) nonflagellar functions (1427 codons) and (iii) unknown functions (2308 codons). We then applied backward elimination (cut-off p = 0.0001), AIC and AICc to select among the full set of fixed-effect models (Table 5) We used BLAST to compare the present set of unknown proteins to other bacilli and found one case (lmo0715) that was similar to a known flagellar protein (FliH). We note that the KEGG database [26] has annotated lmo0715 as a putative flagellar assembly protein. Based on genome location and levels of selection pressure, we suggest that the "unknown" genes in this dataset represent the best candidates for the unaccounted components of the Listeria flagellar machinery. Genes that evolve at high rates can be difficult to identify [27]; however this is not the case here, as estimates of ω indicate a relatively slow rate of nonsynonymous evolution. If these genes indeed encode the missing components of the Listeria flagellar machinery, we speculate an ancestor of Listeria might have acquired them via an LGT event.
The rapidly-evolving non-flagellar genes encoded (i) a protein involved in regulating chemotaxis (lmo0691), (ii) a chemotaxis-related sensory protein (lmo0692), (iii) a cell surface protein (lmo0701), and (iv) a phage-related protein similar to transglycosylase (lmo0717). To further investigate the evolutionary dynamics of these genes we applied the full set of fixed-effect models to them, with each having a unique partition. Again, the model selected by backward elimination (FE9) differed from the model selected by AIC and AICc (FE10). Results under FE9 (Table 6) and FE10 are consistent in suggesting heterogeneous ω among them, with one gene, the cell surface protein, having a substantially higher relative rate of nonsynonymous evolution. Interestingly, a genome wide survey of Listeria genes reveals that, in general, cell-surface proteins exhibit accelerated evolutionary rates as compared to housekeeping genes (unpublished data).
Lastly, we investigated the evolutionary dynamics of flaA as compared to those genes known to encode flagellar related proteins. We applied the full set of models to this subset of proteins, with flaA having a unique partition and the remaining 23 proteins placed in a second partition. In this case backward elimination, AIC and AICc selected model FE13. This indicates that, despite heterogeneity in both c and π, selection pressure on flaA does not differ significantly from the average for genes encoding a flagellar related protein. This result supports the hypothesis that thermoregulation of motility remains an effective adaptation to avoid recognition by the host's innate immune system [24], despite sometimes less than perfect control over gene expression.

Discussion
The simulation results show that under a cut-off p-value = 0.05 the likelihood ratio test is more accurate than AIC and AICc. With the exception of the π parameters, AIC and AICc chose overly complex models more frequently than did the backward elimination procedure. For the π parameters, AIC and AICc chose overly simple models more frequently than did backward elimination. The difference lies in the different cut-off values that are used to penalize the more complex model. Take the heterogeneity test of κ as an example (df = 1), the LRT statistic is defined as twice the difference in log likelihood values between a pair of nested models: Λ = 2 × (ln L( | x) -ln L(θ 0 | x)). Based on the LRT under a significance level of 0.05, we reduce the complexity of the model if Λ is smaller than the critical value 3.84. Under AIC we choose a simpler model only when Λ < 2; hence, AIC tends toward more complex models. However when we reduce a model by more than 7 parameters (e.g. homogeneous π's versus heterogeneous π's, with d.f. = 9), the critical value becomes 16.92 for the LRT, which is less than the critical value of Λ under AIC, 18. In this case, AIC will select the same or simpler model than the LRT. Note that AICc compares Λ with 2 × k × n/ (n-2) which is always greater than 2 × k used by AIC, hence AICc will always choose the same or simpler model than AIC. This property of AICc had only a small effect on the results of model selection, as AICc performed only θ  slightly better than AIC but substantially worse than backward elimination.
LRT statistics involving parameters such as ω appear to be asymptotically χ 2 distributed for random-effect codon models; such models employ a parametric distribution (e.g., the β distribution of Models M7 and M8 [2]) to accommodate among site variability in the ω ratio. However, Aagaard and Phillips [8] reported that for comparisons of Yang and Swanson's [4] models C and E (FE13 and FE1 in Figure 1), the empirical distribution of LRT statistics deviated from the expected χ 2 distribution, leading to a type I error rate in excess of that specified by the level of the test. The results of our simulation study suggest this bias might affect all tests in Figure 1 involving parameters κ, ω and c. Several authors have noted that LRT statistics derived from models that employ empirical estimates of nucleotide or codon frequencies might not be well approximated by a χ 2 distribution [8,28]. Moreover, when Aagaard and Phillips [8] repeated their simulation study under equal codon frequencies and computed LRT statistics by using models with frequency parameters fixed to the true values (π i = 1/61), they found that the LRT statistics matched the χ 2 expectation. Aagaard and Phillips [8] suggested that the approximation of codon frequencies is the source of the observed bias in the LRT.
Indeed, empirical estimates do not satisfy the requirements of LRT [21], and consequently the backward elimination procedure. To further investigate the impact of empirical estimates on model selection, we reanalysed all the simulated datasets under the true codon frequencies; i.e., those used to generate the data. We note that for a given dataset the empirical codon frequencies yielded higher likelihood scores than did the true frequencies.
This was expected, as empirical estimation will "pick up" some of the sampling errors in each simulation replicate. We found that the accuracy and bias of backward elimination, AIC and AICc under the true codon frequencies were identical to those when empirical codon frequencies were used. This suggests that bias in the model selection procedures used here did not arise from empirical estimation of π's alone.
There are several possible explanations for the bias of all three model selection methods in the direction of greater complexity for one of ω, κ or c. One possibility is that  Codon frequencies π i1 = 1/61 π i1 = lysin π i2 = 1/61 π i2 = lysin π i1 = 1/61 π i1 = lysin π i2 = lysin π i2 = 1/61 The lysin gene was used to obtain empirical estimates of codon frequencies, this is indicated in the table by π i = lysin.  potential non-independence among the values for different parameters means that AIC might not be a good approximation to the Kullback-Leibler divergence, and that the requirements for the χ 2 approximation might not be met for the LRT, thus the degree of freedom is not accurate. Also, backward elimination may find a local optimum solution. Under backward elimination the cut-off pvalue is subjectively decided before the tests, leading to the possibility that the procedure will stop too early and suggest an overly complex model. This phenomenon is sometimes seen in the regression context. Clearly, these issues require further attention; in the mean time we explored the possibility of tuning the cut-off p-value in order to improve the performance of backward elimination (see also [8]). After evaluating several cut-off values for p, we found that a substantial improvement in performance was obtained by using p = 0.0001. Moreover, we found that by using this cut-off, nearly all the tendency to select an overly complex model for ω, κ or c was avoided, and that errors for selecting overly-simplistic models happen mostly for datasets where one of the partitions was comprised of a very small number of codon sites.
Our application of fixed-effect models to real data was encouraging, having uncovered previously unrecognized heterogeneity among Listeria genes and among sites within the abalone sperm lysin gene. However, if the objective is only to identify individual positive selection sites within a gene, the a priori structural information is not likely to serve as a perfect proxy for those partitions most relevant to differences in selection pressures. For example, Yang and Swanson [4] showed that the exposed sites of lysin likely include both conserved and positively selected codon positions. Hence, averaging ω over all sites in the exposed partition yields a reduced estimate of positive selection pressure. We note this effect was the same under the best-fit model, FE9, as under FE1 (Model E) used by Yang and Swanson [4]. We agree with Yang and Swanson [4] in anticipating that the power of randomeffect models to test the strength and direction of selection pressure at sites within genes will be greater than fixedeffect models in most cases.
If the objective is to investigate heterogeneous evolution among genes, as in genome-scale analyses, then fixedeffect models are useful. The present set of models represents only a small step towards genome-scale evolutionary models. For example, decoupling synonymous and nonsynonymous rates, as in the random-effect model of Kosakovsky Pond and Muse [29], would allow users freedom to model gradients in synonymous substitution rates along a genome while allowing independent variability in nonsynonymous rates among genes. Yang and Swanson [4] made several suggestions, including the intriguing idea of enforcing a molecular clock for synonymous changes and leaving nonsynonymous changes unconstrained. We predict that joining fixed-effect codon models and datamining methods to obtain new methods analogous to model based clustering [30,31] could provide extremely useful tools for genome-scale data analysis. Lastly, there is growing interest in both the performance of codon models [32,33] and the impact of heterogeneity among genes [34,35] in multi-gene phylogenetic analysis; improved ability to model among-gene heterogeneity at the codon level could improve their utility for comparing alternative phylogenomic hypotheses.

Conclusion
Random-and fixed-effect codon models have unique advantages and disadvantages. Random-effect models are desirable when there is no a priori knowledge by which sites might be partitioned, or when only a few sites are expected to comprise a partition of interest (but see [5]). Their disadvantage is that models for heterogeneity among sites in features such as the transition to transversion ratio (κ) and equilibrium codon frequencies (π) are unavailable. Fixed-effect models are desirable when data partitions are known to exhibit significant heterogeneity in parameters such as c, κ or π, or when a statistical test of such heterogeneity is desired. The disadvantage here is that any uncertainty in the site partition is not accommodated.
The growing importance of phylogenomics and metagenomics (e.g., [36,37]) will lead to a greater need for models suitable for testing hypotheses, and estimating rates and patterns of evolution, in large multi-gene datasets. Although considerable development remains to be done, we believe the present set of models will find many useful applications provided that results are interpreted with the inherent limitations of the methods in mind. In particular we note: (i) power can be low (see also [8]), particularly when partitions are small, (ii) the accuracy of the partitions may influence the results of model specification and (iii) the tree topology is assumed to be known without error. For the time being we make the following recommendations: (i) select among models by using backward elimination rather than AIC or AICc, (ii) use a stringent cut-off p-value; p = 0.0001 seems appropriate, and (iii) sensitivity analysis should be included in an investigation. Sensitivity of results should be investigated for robustness to tree topology and model of codon frequencies. We note that by using Akaike weights [16] to quantify the evidence in favour of a model, estimates of parameters could be obtained that accommodate model uncertainties (e.g., [19]). Where practical, we recommend that sensitivity to alternative data partitions also should be explored. Lastly, any complex model can have convergence problems or implementation errors; one must always inspect the parameter estimates for atypical results. With thoughtful application, fixed-effect codon models should provide a useful tool for large scale multi-gene analyses.

Fixed-effect models of codon evolution
The basic codon model [9,10] assumes that the process of substitutions from one codon to another is a Markov process, where the next codon state only depends on the present state, and not on any past state. The element P ij (t) in the transition matrix P(t) gives the probability of going from codon i to codon j during time t. Because they do not occur within a functional protein-coding gene, the three stop codons, UAA, UAG and UGA, are excluded. The transition matrix P(t) can be calculated by P(t) = e Qt , where Q = {q ij } is a 61 × 61 rate matrix. The element q ij denotes the instantaneous substitution rate from codon i to codon j as follows: When a change among codons involves a transition, the rate is multiplied by κ, the transition/transversion rate ratio. In the same way, if the substitution is nonsynonymous, the rate is multiplied by ω, the nonsynonymous/ synonymous rate ratio. Usage of codons within a gene can also be highly biased, and consequently the rate is multiplied by the equilibrium frequency of the targeted codon π j , which is assumed to remain unchanged between generations of the evolutionary process. Given prior information by which sites can be partitioned into classes, parameters such as ω, κ and π can be allowed to differ among partitions, with different Q matrices used for different partitions [4]. Independence among all codon sites is assumed; hence the log likelihood of the complete dataset is the sum of the log likelihood of each site [4].
For all fixed-effect codon models described in this paper ( Table 2) the tree topology is fixed a priori and, with the exception of the codon frequencies, the parameter values are estimated by numerical maximization of the likelihood function. Codon frequencies are empirically estimated from the data. For parameters heterogeneous among partitions, we use the maximum likelihood estimation implemented by Yang and Swanson [4]. For the models with identical substitution rate, we simply fix the branch length ratio c at 1 across partitions. For models with homogeneous κ or ω, we use an algorithm similar to the Expectation-Maximization (EM) algorithm [38]. Let's take a single parameter, say ω, as an example of a homogenous parameter. We first independently estimate the parameters in each partition by maximum likelihood. At the "E-step" of the algorithm we obtain the weighted average of ω over all partitions, where the weight is given by the proportion of codon sites in the corresponding partition. At the "M-step" we re-estimate parameters heterogeneous over partitions after fixing ω to its weighted average value. Following this we run the "E-step" by fixing the parameters assumed to be heterogeneous, and re-estimating ω from each partition; an updated estimate of ω is again obtained as the weighted average. The E-and Msteps are run iteratively until successive estimates of ω converge.

Model selection among a set of related fixed-effect models
Backward elimination reduces the most complex codon model, shown at the top of Figure 1, to a simpler one in a stepwise fashion. We begin from model FE1, which assumes different κ, ω, c, and π's for different site classes, and then compare it with simpler models which assume one of these four parameters to be homogeneous (see next level in Figure 1). For example, if the hypothesis of homogeneous κ is not rejected, we will go to FE9 at the next level in Figure 1, which assumes different ω, c, π's but the same κ for different site classes. We then compare FE9 with its nested models at the next level (see FE10, FE11 and FE13 in Figure 1). If more than one homogenous model cannot be rejected by the LRTs at a given level, the backward elimination procedure will choose the model with the largest p-value in the LRT.
Akaike Information Criterion (AIC) is based on minimizing the expected Kullback-Leibler divergence [39,40]; AIC = -2 × ln L( | x) + 2 × k, where k denotes the number of free parameters in the candidate model. For small samples, the AIC is corrected by a second order bias adjustment in the regression and time series settings in order to avoid over-fitting, let n denote the number of observations [17]: . This adjustment places a heavier penalty on the number of parameters when the number of observations is not much larger than the number of parameters; thus AICc tends to choose a simpler model than AIC. We note that the basis for this correction is not expected to hold outside of the regression and time-series settings; however, because we desired a correction for small samples we evaluated the performance of AICc by numerical simulations. In our analysis k denoted the number of parameters in a given model ( Table 2) and n denoted the number of codon sites. For both AIC and AICc, the model with the smallest score is chosen as the ideal model. The full set of 16 fixed-effect models ( Figure 1) is implemented in a modified version of codeml that we make freely available [41]. The original version of codeml (3.14b) is one of several programs in the PAML package [42] distributed by Ziheng Yang [43].

Simulated and real sequence data
We simulated codon evolution by using the "evolver" program of the PAML package [42]. There were 16 different scenarios, based on the heterogeneity or homogeneity of four parameters ( Table 2) among two data partitions. Parameter values are shown in Table 7. We also varied the ratio of the number of codon sites contained in each partition (50:50, 70:30 or 90:10). Data were simulated independently in two partitions and then concatenated to obtain a composite dataset. Each composite dataset contained 2000 codons for 16 species. A total of 432 composite datasets were possible; i.e., the product of 2 options on κ, 3 options on c, 4 options on π, 6 options on ω and 3 partition proportions. However, as we had two options for heterogeneous rates (1:3 or 1:9) and one possibility for the homogeneous rates (1:1) we made an adjustment to obtain similar number of datasets for each of the 16 scenarios. Specifically, we simulated under only 4 of the options on ω (first and second rows for selection pressure in Table 7) when rates were assumed to be heterogeneous among partitions (c 2 = 3 and c 2 = 9). Thus 2(κ) × 2(c) × 4(π) × 2(ω) × 3(sites proportion) = 96 possibilities were not included in our simulation. This strategy provided a grand total of 336 composite datasets for evaluating the performance of the different model selection strategies.
The first real dataset was comprised of sequences for the sperm lysin gene from 25 species of abalone. This is one of the original test cases of Yang and Swanson's paper [4], and it is distributed online by Ziheng Yang as part of the PAML package [43]. Note that this lysin dataset and phylogeny is the same as those from [44], except that a single site containing an alignment gap was excluded. The second real dataset was comprised of 37 genes located within the genomic region of Listeria bacteria that encode the components of the flagellar system (lmo0675 -lmo0717). We note that this region includes several proteins of unknown function, and for comparative purposes we included them in our dataset. Genes were partitioned according to functional categories of the ListiList database [45]. The sequence alignments, phylogenetic trees and gene ontologies for this multi-gene dataset is available online [46]. Although multi-gene datasets can be much larger than ours, it represents both a real biological problem and serves as an illustration of the types of data upon which these techniques can be used.


extfieldsofstudy: Medicine Biology

text:
The Histone Deacetylase Genes HDA1 andRPD3 Play Distinct Roles in Regulation of High-Frequency Phenotypic Switching in Candida albicans

ABSTRACT Five histone deacetylase genes (HDA1, RPD3, HOS1, HOS2, and HOS3) have been cloned fromCandida albicans and characterized. Sequence analysis and comparison with 17 additional deacetylases resulted in a phylogenetic tree composed of three major groups. Transcription of the deacetylases HDA1 and RPD3 is down-regulated in the opaque phase of the white-opaque transition in strain WO-1. HOS3 is selectively transcribed as a 2.5-kb transcript in the white phase and as a less-abundant 2.3-kb transcript in the opaque phase. HDA1 andRPD3 were independently deleted in strain WO-1, and both switching between the white and opaque phases and the downstream regulation of phase-specific genes were analyzed. Deletion ofHDA1 resulted in an increase in the frequency of switching from the white phase to the opaque phase, but had no effect on the frequency of switching from the opaque phase to the white phase. Deletion of RPD3 resulted in an increase in the frequency of switching in both directions. Deletion ofHDA1 resulted in reduced white-phase-specific expression of the EFG1 3.2-kb transcript, but had no significant effect on white-phase-specific expression of WH11 or opaque-phase-specific expression of OP4, SAP1, andSAP3. Deletion of RPD3 resulted in reduced opaque-phase-specific expression of OP4, SAP1, and SAP3 and a slight reduction of white-phase-specific expression of WH11 and 3.2-kbEFG1. Deletion of neither HDA1 norRPD3 affected the high level of white-phase expression and the low level of opaque-phase expression of the MADS box protein gene MCM1, which has been implicated in the regulation of opaque-phase-specific gene expression. In addition, there was no effect on the phase-regulated levels of expression of the other deacetylase genes. These results demonstrate that the two deacetylase genes HDA1 and RPD3 play distinct roles in the suppression of switching, that the two play distinct and selective roles in the regulation of phase-specific genes, and that the deacetylases are in turn regulated by switching.

Most strains of Candida albicans and related species switch, or can be induced to switch, between two or more general phenotypes distinguishable by colony morphology (14,41,42,44,48). Switching regulates a number of phenotypic characteristics, including putative virulence traits, and for that reason, it has been considered a higher-order virulence factor (42). Switching has been demonstrated to occur at sites of commensalism (21) and infection (48)(49)(50) and to occur at higher frequencies in strains causing infections than in commensal strains (21,24,61). Therefore, it has been proposed that switching provides colonizing populations with phenotypic variants for rapid adaptation in response to environmental challenges, such as drug therapy and the immune response, as well as changes in host physiology, the competing microflora, and anatomical locale (34,45).
By using the simple white-opaque transition in C. albicans strain WO-1 (42) as a model system, it has been demonstrated that switching involves the transcription of a number of phasespecific genes (3,13,22,32,33,46,(53)(54)(55) and that the differential expression of these genes is accomplished through phase-regulated transcription factors (30,46,47,51,56), leading to a model in which a reversible, spontaneous change at a master switch locus leads to the downstream activation and deactivation of unlinked phase-regulated genes (46,47). Since the cis-acting regulatory sequences of coordinately expressed genes differ, models have been developed that accommodate different DNA binding proteins (46). However, in spite of recent progress in functionally characterizing the promoters of phase-regulated genes, little is known about the actual switch event or how genes are coordinately regulated in the whiteopaque transition. Because switching occurs spontaneously and reversibly at relatively high frequencies (10 Ϫ4 to Ͼ10 Ϫ2 ), it had been suggested (37,44) that it either represents a spontaneous and reversible reorganization of DNA at a master switch locus, which occurs in bacterial switching systems (15), or a reversible alteration in the state of chromatin at a master switch locus, which occurs in yeast when genes are placed at sites adjacent to subtelomeric regions (17,36). Because the "silent information regulator" gene SIR2 plays a role in gene repression in the latter mechanism, it was suggested that a C. albicans homolog of SIR2 may play a role in the regulation of switching at a master switch locus in C. albicans (37,44). Perez-Martin et al. (35) subsequently demonstrated that deletion of a C. albicans SIR2 homolog in strain CAI4 did in fact result in up-regulation of the switching system in that strain, which is analogous to the first switching system described in C. albicans strain 3153A (41,54). In addition to the "silent information regulator" (Sir) proteins, other classes of proteins are involved in the repres-sion or silencing of developmentally regulated genes in eukaryotes. One such class of genes, the histone deacetylases, has been demonstrated to regulate chromatin structure through selective histone deacetylation, which in turn affects chromatin folding and interactions between DNA and DNA-binding proteins (2,63). The role of the deacetylases in gene regulation has been demonstrated by generating mutants with loss of function or with dominant-negative effects in both Saccharomyces cerevisiae (25,39,58) and human cell lines (18,20).
Recently, we tested whether the specific deacetylase inhibitor trichostatin A (10, 65) affected the white-opaque transition (26). We found that the inhibitor caused a selective increase in the frequency of switching in the white-to-opaque transition, but had no effect on the frequency of switching in the opaqueto-white transition, suggesting that deacetylation through a trichostatin-sensitive deacetylase selectively suppresses switching in one direction (26). Since the deacetylase Hda1p is highly sensitive to trichostatin A (7), we deleted the HDA1 gene and obtained a mutant phenotype similar to that of trichostatin A-treated cells, supporting the conclusion that deacetylation through Hda1p suppresses switching selectively in the whiteto-opaque direction (26). In addition to playing a role in switching, deacetylation may also play a role in the regulation of phase-specific gene transcription. We therefore cloned five of the major C. albicans histone deacetylase genes with homology to the genes that encode the five known histone deacetylases in S. cerevisiae (HDA1, RPD3, HOS1, HOS2, and HOS3) and analyzed their deduced protein sequences and expression patterns in the white-opaque transition. We have also generated, in addition to the HDA1 deletion mutant, a deletion mutant of RPD3. We present evidence that the expression of the deacetylase genes is affected by switching and that, while HDA1 plays a selective role in suppressing the basic switch from white to opaque, but not from opaque to white, RPD3 plays a role in suppressing the basic switch events in both directions. In addition, both HDA1 and RPD3 play indirect but distinct roles in regulating the levels of expression of select phase-specific genes.

MATERIALS AND METHODS
Maintenance of stock cultures. C. albicans strain WO-1 (42) was periodically cloned from a frozen stock culture on agar containing modified Lee's medium (4). Strains TS3.3, a ura3 auxotroph (55), and TUI7, a URA3 prototrophic derivative of TS3.3, were maintained on agar containing modified Lee's medium with or without 0.01 mM uridine, respectively. Both the HDA1 and RPD3 homozygous mutant strains were maintained on agar containing modified Lee's medium. The genotypes of all strains used in this study are described in Table 1.
Cloning the deacetylase genes HDA1 and RPD3. Based on short DNA sequences reported in the Stanford C. albicans genome database (www-sequence .stanford.edu), two sets of primers, DeACT5Ј/DeACT3Ј and RPD5Ј/RPD3Ј (Table  2), were designed to amplify by PCR DNA fragments containing HDA1 and RPD3 sequences, respectively, with genomic DNA of strain WO-1 as a template and with Taq polymerase. The PCR-derived fragments were then used to screen a EMBL3A genomic library of C. albicans strain WO-1 (53). Of approximately 100,000 plaques screened for the HDA1 fragment or for the RPD3 fragment, 30 to 40 putative clones were identified. Ten putative clones of each gene were purified by a secondary screen and analyzed by Southern blot hybridization with the respective probe to assess purity. Based on these preliminary results, DA3.1 and RP2.1 were chosen to represent CaHDA1 and CaRPD3, respectively.
To characterize HDA1, DA3.1 was digested with SalI, and a 2.1-kb fragment that hybridized with the respective probe was subcloned into pGEM3Z to generate pDA14.1. The SalI DNA fragment was sequenced in both directions with an ABI model 373A automatic sequencing system and fluorescent Big Dye termination chemistry (PE-ABI Inc., Foster City, Calif.). In addition, approximately 1 kb of the 5Ј upstream sequence and 0.5 kb of the 3Ј downstream sequence of the SalI fragment were determined with DA3.1 as a template. To characterize RPD3, a 4-kb sequence of RP2.1 was directly sequenced with the same protocol applied to pDA14.1. Then, by using the specific primer pair FANRPD5Ј/ FANRPD3Ј (Table 2), the full-length RPD3 open reading frame (ORF) was amplified by PCR and subcloned into pGEM5Z between end-repaired SacI sites to generate pGRP16/5.
Analysis of protein sequences and construction of a phylogenetic tree. Alignment of multiple protein sequences was performed with the Clustal W/Jalview multiple alignment editor, version 4 (60). Pairwise comparisons between protein sequences were performed with the PROTDIST program of the PHYLIP package, version 3.57C (http://evolution.genetics. washington.edu/Phyl.p.html). The unrooted dendrograms were generated by using the FITCH program of the PHYLIP package (the Fitch-Margoliash least-squares distance method) (11). The data set was subjected to a bootstrap analysis (1,000 replicates) by sequential use of the SECBOOT, PROTDIST, and Consensus programs. Genetic distances were derived by the Day'hoff PAM matrix algorithm (8).
Construction of homozygous HDA1 deletion mutants. Two different deletion cassettes were constructed, each spanning the essential deacetylation motifs. To generate the heterozygote, a chloramphenicol acetyltransferase (CAT)-URA3-CAT-based cassette (55) was constructed. The pDA14.1 plasmid containing the 2.1-kb SalI fragment of HDA1 (Fig. 1A) was digested with BclI, deleting 1,147 bp of the HDA1 ORF, end repaired with T4 DNA polymerase, and dephosphorylated with shrimp alkaline phosphatase (SAP). The linearized vector was purified by Tris-borate-EDTA (TBE)-agarose gel electrophoresis and ligated with the 3.5-kb CAT-URA3-CAT cassette from the plasmid pCUC (55). The CAT-URA3-CAT cassette was isolated by digesting pCUC with BamHI, followed by end repair with T4 DNA polymerase. The resulting plasmid, pA48.1, contained the deleted version of HDA1 (Fig. 1A). pA48.1 was digested with SalI, and 25 g was  (Table 1). Recovered transformants were tested for heterozygosity by digestion of genomic DNA with XbaI and analyzed by Southern blot hybridization with the HDA1 ORF. One confirmed heterozygous clone was subjected to 5-fluoroorotic acid (5-FOA) treatment in order to induce "pop-outs" of the URA3 gene. To generate the homozygote, a hisG-URA3-hisG-based cassette (12) was constructed. The BclI-BclI deletion fragment of plasmid pDA14.1 was isolated, end repaired with T4 DNA polymerase, dephosphorylated with SAP, and subcloned at the EcoRV site of pGEM5Z to generate pC25.1 (Fig. 1A). pC25.1 was digested with EcoRV and BglII to delete 702 bp of DNA, end repaired with T4 DNA polymerase, dephosphorylated with SAP, and ligated with the hisG-URA3-hisG cassette from plasmid pMB7 (55). The hisG-URA3-hisG cassette was prepared by digesting PMB7 with SalI and BglII and end repaired with T4 DNA polymerase. The resultant disruption cassette, pC88.10 ( Fig. 1A), was then digested with ApaI and SacI. This DNA fragment could only target the functional HDA1 allele. It was used to transform the selected heterozygous ura3 Ϫ strain. Transformants, obtained on selection plates, were tested for homozygosity by Southern analysis as described earlier.
Construction of homozygous RPD3 deletion mutants. The hisG-URA3-hisG cassette was used to delete both alleles of RPD3. To construct the disruption cassette, pGRP16/5 was digested with SacI and BglII to delete 697 bp. The 4.0-kb hisG-URA3-hisG cassette was then inserted in its place to generate the plasmid pC7a.4 (Fig. 1B). The RPD3 disruption cassette was prepared by digesting pC7a.4 with ApaI and NsiI and used in both the first and second rounds of transformation to obtain first the heterozygous strain and second the homozygous strain. Prior to the second round of transformation, the selected heterozygous clone was subjected to the 5-FOA "pop-out" protocol to derive ura3 auxotrophy. The heterozygous and homozygous strains were confirmed by probing Southern blots of PstI-digested DNA with the RPD3 ORF.
Northern and Southern blot analyses. Northern analyses were performed according to methods previously described (54,55). To ensure that the growth conditions for the compared cell preparations were similar, cells of each cell type were removed from agar cultures and grown to the stationary phase in liquid culture. White-phase cells were then diluted to 10 6 and opaque-phase cells were diluted to 5 ϫ 10 5 in fresh growth medium and grown to mid-log phase (9 ϫ 10 6 and 5 ϫ 10 6 per ml, respectively) prior to harvesting. Total RNA was extracted with the RNeasy mini kit according to manufacturer's specifications (Qiagen, Inc., Santa Clarita, Calif.). The fold differences in transcript levels were measured according to methods recently described (55). The hybridization probes for WH11, SAP1, EFG1, and OP4 contained the full-length ORFs, derived either by PCR with specific primers (Table 2) or by digesting them from the respective plasmids containing the cognate DNA inserts (32,33,53,55). The probe for SAP3 was derived by PCR with the proximal activation sequence pAS3 as a template and the specific primers SAP3-5Ј and SAP3-3Ј (Table 2). Probes for HDA1 encompassed either the full-length ORF or the 702-bp EcoRV-BglII deletion fragment from pC25.1, as described earlier. Probes for RPD3 encompassed either the full-length ORF or the 697-bp SacI-BglII deletion fragment from pGRP16/5, as described earlier. The probes for HOS1, HOS2, HOS3, MCM1, and TUP1 encompassed the full-length ORFs obtained by PCR with the primer pairs previously described ( Table 2). Southern analyses were performed according to methods previously described (33,52,55). Specific restriction enzymes and probes are described in Results.

RESULTS
Cloning and characterization of the five histone deacetylase genes. A search of the C. albicans genome database for homology at the amino acid level initially revealed five genes homologous to each of the S. cerevisiae histone deacetylase genes HDA1, RPD3, HOS1, HOS2, and HOS3. HDA1 and RPD3 were cloned by probing a lambda genomic library of strain WO-1 (42) with PCR-generated DNA fragments based on nucleotide sequences reported in the Stanford C. albicans genome database. The three HOS genes were amplified from WO-1 genomic DNA by PCR with primers based on reported nucleotide sequences. The cloned HDA1, RPD3, HOS1, HOS2, and HOS3 homologs contained ORFs encoding proteins with sizes of 653, 478, 392, 454, and 713 amino acids, respectively (Fig. 2).
To test whether the five C. albicans histone deacetylase proteins represented distinct members of the histone deacetylase superfamily, alignment of these proteins and 12 additional fungal histone deacetylases was performed with the Clustal W multiple alignment protocol described in Materials and Methods. The alignment of approximately 330 amino acids included the nine sequence blocks that identify histone deacetylase subtypes (29). The phylogenetic analysis (Fig. 3) revealed that the five C. albicans deacetylases grouped into three distinct classes: the Hda1 class (C. albicans Hda1p), the Rpd3 class (C. albicans Rpd3p, HOS1p, and HOS2p), and the Hos3 class (C. albicans HOS3p). In every class, the C. albicans deacetylases were most similar to the corresponding S. cerevisiae deacetylase (Fig. 3). Further comparison of the five C. albicans histone deacetylase proteins with S. cerevisiae Hda1p and ScRpd3p revealed that a 200-amino-acid stretch in the middle of all seven contained highly conserved amino acid residues (Fig. 2). In particular, a highly conserved functional deacetylation motif of 53 to 63 amino acid residues resided in the middle of this stretch (see box in Fig. 2). Pairwise comparison of the seven proteins revealed four highly conserved histidines (presented in boldface in Fig. 2) organized in pairs and separated by roughly the same number of amino acids (5,10,18,26,65). Comparison of the deacetylation motifs of the five C. albicans deacetylases and 17 additional deacetylases revealed 12 identical amino acid residues and overall a high level of conservation (data not shown).
Expression of the deacetylase genes in the white-opaque transition. Northern analysis revealed that expression of the deacetylases differed between white-and opaque-phase cells. The transcript levels of HDA1 and RPD3 were significantly lower in the opaque phase (Fig. 4). The transcript levels of HOS1 and HOS2 were slightly lower in the opaque phase (Fig.  4). In the case of HOS3, both the levels and the molecular size of the transcripts differed between the two phases. In the white phase, the molecular size of the transcript was 2.5 kb, and in the opaque phase, the molecular size of the less-abundant transcript was 2.3 kb. Northern analysis with two independent RNA samples revealed similar results.
Creating deletion mutants of HDA1 and RPD3. To create deletion mutants of HDA1 and RPD3, we employed the "urablast" gene knockout strategy (12) in strain WO-1 (42). To create the HDA1 deletion mutant, the ura3 Ϫ strain TS3.3 was first transformed with the CAT-URA3-CAT-based HDA1 disruption cassette from pA48.1 (Fig. 1A). Of 40 transformants, 3 proved to be heterozygous by Southern analysis with the full-  (60). The rectangular unshaded box in the center represents the highly conserved deacetylation motif. The gray shaded area, including the unshaded rectangular box, contains amino acid residues conserved to various degrees in the deacetylases. The identical residues among all deacetylases are denoted by stars, while conservative replacements of amino acid residues are denoted by either two stacked solid circles (based on similar functional groups) or one solid circle (based on similar effects on secondary structure). Conserved histidines are presented in boldface. The accession numbers for HDA1 and RPD3 in GenBank are AF377894 and AF377895, respectively. Other accession numbers are as follows: ScHda1p, Saccharomyces cerevisiae Z71297; SCRpd3p, S. cerevisiae P32561. Ca, Candida albicans. homozygous transformant HDho15 was selected from the first transformation, and HDho11 and HDho19 were selected from the second. When Southern blots were probed with the deleted HDA1 fragment, the parental strain TS3.3 exhibited the two allelic bands at 9.8 and 3.6 kb, and the homozygous strain HDho15 exhibited neither band, confirming that both alleles of HDA1 were deleted in HDho15 (Fig. 5A). A similar analysis proved that both alleles of HDA1 were deleted in HDho19 and HDho11 (data not shown).
To create the RPD3 deletion mutant, one urablast cassette based on hisG-URA3-hisG was used to knock out both alleles. TS3.3 was transformed with the disruption cassette from pC7a.4 (Fig. 1B). Of 60 transformants, 44 (73%) proved to be heterozygous for RPD3. When Southern blots of PstI-digested TS3.3 DNA were probed with the RPD3 ORF, a single band was identified at 5.8 kb (Fig. 5B). When the URA Ϫ heterozygote strain RPheF4 was retransformed with the pC7a.4 cassette, 40 transformants were obtained. Six proved to be homozygous disruptants. Two transformants, RPho19 and RPho3, were chosen for further analysis. The 5.8-kb band was missing in Southern blots of RPho19 (Fig. 5B) and RPho3 (data not shown) probed with the RPD3 ORF. Both contained the 6.2and 9.1-kb molecular size bands. When probed with the deleted RPD3 fragment, the parent exhibited a 5.8-kb band (Fig.  5B), while RPho19 (Fig. 5B) and RPho3 (data not shown) exhibited no bands, confirming that both alleles of RPD3 were deleted in the latter strains.
Effects of HDA1 deletion on switching and phase-specific gene expression. We recently demonstrated that treatment with the deacetylase inhibitor trichostatin A (TSA) or deletion of the most TSA-sensitive gene, HDA1, had a selective effect on switching (26). While both TSA-treated cells and cells of the two independent mutant strains HDho15 and HDho19 switched from the opaque to the white phenotype at frequencies comparable to that of untreated wild-type cells (ϳ10 Ϫ3 ), both TSA-treated and mutant cells switched from the white-to the opaque-phase phenotype at frequencies more than an order of magnitude greater than that of wild-type cells (ϳ3 ϫ 10 Ϫ2 ) (26). These results demonstrated that although the deletion of HDA1 selectively increased the frequency of switching in the white-to-opaque direction, it had no effect on the unique signature morphology of opaque-phase cells.
To assess the effects of deleting HDA1 on phase-specific gene expression, Northern blots of total cellular RNA of whiteand opaque-phase cells of the parent strain TU17 and the  6A) and HDho19 (data not shown) were probed with the white-phase-specific genes WH11 (53) and EFG1 (55) and the opaque-phase-specific genes OP4 (32), SAP1 (33), and SAP3 (22,62). Deletion of HDA1 had little effect on the basic developmental regulation of most of the phase-regulated genes tested (Fig. 6A). It also had no significant effect on the level of expression of WH11 in the white phase or on the level of expression of OP4, SAP1, or SAP3 in the opaque phase (Fig. 6A). Deletion did, however, reduce the level of the EFG1 transcript in the white phase fivefold. Similar results were obtained in a repeat experiment in which RNA was extracted from independent growth cultures of TU17 and both mutant strains HDho15 and HDho19. When opaque-phase cell cultures are shifted from 25°C to 42°C, they switch en masse to the white phase (32,38,42,53,55). The average cell commits to the white-phase phenotype between 3 and 4 h, concomitant with the second cell doubling and a switch in phase-specific gene expression. At the time of phenotypic commitment, the expression of WH11 is activated and the capacity to induce OP4 and SAP1 expression by shifting cells from 37°C back to 25°C is lost (32,42,53,55). To test whether HDA1 plays a role in this transition, opaque-phase cells of the control strain TU17 and the homozygous deletion mutant HDho15 were shifted from 25°C to 42°C. After 1 h (prior to phenotypic commitment) and 7 h (after phenotypic commitment), the samples were shifted back to 25°C. After 1 h at 42°C, both TU17 and HDho15 cells contained no detectable OP4 message. A shift to 25°C after 1 h at 42°C induced reexpression of OP4 in both TU17 and HDho15. After 7 h at 42°C, TU17 cells contained no OP4 message, and HDho15 cells contained a negligible level (data not shown). A subsequent shift to 25°C did not reactivate OP4 expression in either strain (data not shown), demonstrating that the normal switch in OP4 regulation had occurred. After 1 h at 42°C, both TU17 and HDho15 cells contained no detectable WH11 transcript, and a subsequent shift to 25°C did not up-regulate WH11 expression in either strain (data not shown). After 7 h at 42°C, both TU17 and HDho15 cells expressed a low level of WH11 transcript, and a shift to 25°C resulted in up-regulation of WH11 in both strains (data not shown), demonstrating that the switch in WH11 regulation had occurred in the mutant. These results demonstrate that deletion of HDA1 has no discernible effect on the basic activation and deactivation of phase-specific genes at the point of phenotypic commitment.
The effect of RPD3 deletion on switching and phase-specific gene expression. Deletion of RPD3 also had an effect on the frequency of switching, but the effect in this case was on switching in both directions. While the frequency of opaque-phase CFU in 5-day-old white-phase colonies of the parental strain TU17 was 4 ϫ 10 Ϫ4 , those of the homozygous deletion mutants  RPho19 and RPho3 were 1 ϫ 10 Ϫ2 and 2 ϫ 10 Ϫ2 , respectively, representing 25-and 50-fold increases, respectively (Table 3). And while the frequency of white-phase CFU in 5-day-old opaque-phase colonies of the parental strain TU17 was 6 ϫ 10 Ϫ4 , those of the homozygous deletion mutants RPho19 and RPho3 were 5 ϫ 10 Ϫ2 and 10 Ϫ2 , representing 83-and 17-fold increases, respectively ( Table 3). The increases in the frequencies of the transition from white to opaque and opaque to white in RPD3 deletion mutants are evident in cultures in which cells from white and opaque colonies are streaked across nutrient agar (Fig. 7). The same differences were also evident in the frequency of sectored colonies (Table 3). Scanning electron micrographs of white-and opaque-phase cells revealed that the RPD3 deletion had no effect on either the smoothsurfaced, round budding phenotype of white-phase cells or the pimpled, elongate phenotype of opaque-phase cells (data not shown).
To assess the effect of deleting RPD3 on phase-specific gene expression, Northern blots of total cellular RNA of white-and opaque-phase cells of the parent strain TU17 and the deletion mutant RPho19 were probed with the white-phase-specific genes WH11 (53) and EFG1 (55) and the opaque-phase-specific genes OP4 (32), SAP1 (33), and SAP3 (22,62). Deletion of RPD3 had no effect on the developmental regulation of any of the phase-specific genes (Fig. 6B). It had a small effect on the transcript levels of the white-phase-specific genes WH11 and EFG1 and the opaque-phase-specific gene SAP1, all of which were approximately twofold lower (Fig. 6B). Deletion of RPD3 caused 10-to 15-fold decreases in the transcript levels of the two opaque-phase-specific genes OP4 and SAP3 (Fig. 6B). Similar results were obtained in a repeat experiment in which RNA was extracted from independent growth cultures of TU17 and two deletion strains, RPho19 and RPho3 (data not shown).
The expression of MCM1 and TUP1 in HDA1 and RPD3 deletion mutants. The opaque-phase-specific genes OP4 (30) and SAP3 (46; S. Lockhart and D. R. Soll, unpublished observations) contain MADS box protein consensus binding sites associated with the cis-acting activation sequences in their promoters, suggesting that a MADS box protein is involved in the regulation of these opaque-phase-specific genes. The MADS box binding sequences exhibited greatest homology to the MCM1 binding site in S. cerevisiae (30). We therefore cloned MCM1 by PCR using a sequence identified in the Stanford C. albicans database. Homology to the S. cerevisiae MCM1 gene was confirmed by partial sequencing. The cloned sequence was then used to probe Northern blots of the parental strain TU17, the HDA1 deletion mutants HDho15 and HDho19, and the RPD3 deletion mutants RPho3 and RPho19. We first found that the MCM1 transcript was high in whitephase cells and low in opaque-phase cells in parental strain TU17 (Fig. 8A). Neither deletion of HDA1 nor deletion of RPD3 affected the levels of expression in the white or opaque phase (Fig. 8A). The same blot was probed with TUP1, a transcription factor gene (6) expressed in white-and opaquephase cells, but at slightly lower levels in the latter (R. Zhao, S. R. Lockhart, and D. R. Soll, unpublished data). Neither deletion of HDA1 nor deletion of RPD3 affected levels of expression in the white or opaque phase (Fig. 8A).
Expression of the deacetylases in the HDA1 and RPD3 deletion mutants. The deletion mutants HDA1 and RPD3 showed no hybridization bands with the respective probes, demonstrating that these strains were devoid of the respective deacetylase mRNAs (Fig. 8B). The deletion of HDA1 had no effect on the transcript levels of RPD3 and HOS3 (Fig. 8B) or HOS1 and HOS2 (data not shown) in either the white or opaque phase. Likewise, deletion of RPD3 had no effect on the transcript levels of HDA1 and HOS3 (Fig. 8B) or HOS1 and HOS2 (data not shown).

DISCUSSION
The molecular mechanisms involved in the downstream regulation of phase-specific genes in the process of high-frequency phenotypic switching have begun to emerge in recent years. In the white-opaque transition, which has been used as an experimental model system to study these mechanisms (46,47), it has been demonstrated that phase-regulated genes are activated through the interaction of phase-specific trans-acting factors with specific cis-acting sequences in the promoters of these genes (30,51,55,56). Because the cis-acting sequences in coordinately activated genes of the same phase differ, the simplest model of regulation involving a common trans-acting DNA-binding protein has been excluded (46), and alternative more complex models have been entertained (46). The emerging complexity of the circuitry that regulates the expression of a Cells from white-or opaque-phase colonies were plated and incubated at 25°C for 5 days, and then the proportions of colony phenotypes (white, opaque, or sectored) were counted. phase-specific genes in the white-opaque transition should have been expected, given the complexity of the molecular mechanisms regulating eukaryotic gene expression in general. Indeed, in the regulation of a gene, it has become clear that in addition to activators and repressors, chromatin-modulating machinery is recruited that is involved in establishing transcriptionally active or inactive states (57). Among the classes of proteins involved in chromatin modifications, the deacetylases have been demonstrated to function in the repression of gene loci through the selective deacetylation of histones H3 and H4 (2,63). To test whether deacetylases also play a role in switching, we initially performed an experiment with the deacetylase inhibitor TSA and found that it caused a dramatic and selective increase in the frequency of switching in the white-to-opaque direction, but had no effect on switching in the opaque-towhite direction (26). Since TSA preferentially inhibits the major deacetylase Hda1p (7), we deleted the gene coding for this protein in C. albicans strain WO-1 and found that the mutant phenotype was similar to that of TSA-treated cells (26). Switching was selectively upregulated in the white-to-opaque direction only, suggesting that deacetylation through Hda1p functioned either at a "master switch locus" to suppress the event or at the site of an activator of the switch event (26). Here we have extended these studies by identifying the mem-bers of the deacetylase family in C. albicans, examining the effects of switching on expression, testing whether HDA1 plays a role in phase-specific gene expression, and testing whether a second major deacetylase, RPD3, also plays a role in switching and phase-specific gene expression.
The histone deacetylases in C. albicans. We cloned five histone deacetylases from C. albicans with homology to the five deacetylases in S. cerevisiae. By generating phylogenetic trees based on homology comparisons between the amino acid sequences of the five cloned C. albicans deacetylase genes, the five S. cerevisiae deacetylase genes and deacetylase genes from additional fungi, three groups were identified. One contained C. albicans Hda1p; the second contained C. albicans Hos1p, Hos2p, and Rpd3p; and the third contained C. albicans Hos3p. In each group, the deacetylase with the highest homology to each C. albicans deacetylase was the S. cerevisiae homolog. A comparison of the deacetylation motif of 22 deacetylases that included the five C. albicans deacetylases revealed 12 identical amino acid residues and very high levels of homology overall. By using a more rigorous search of the protein database, a sixth C. albicans deacetylase was recently identified in the Stanford C. albicans genome database with homology to Rpd3p (data not shown). We are in the process of characterizing the role of Both HDA1 and RPD3 play roles in the switching process. We had previously demonstrated the selective role of HDA1 in suppressing switching in the white-to-opaque direction (26), so we considered the possibility that a second deacetylase might selectively suppress switching in the opaque-to-white direction. We were therefore surprised to find that deletion of RPD3 resulted in an increase in the frequency of switching in both the white-to-opaque and opaque-to-white directions. The observation that deletion of either deacetylase results in an increase in switching in the white-to-opaque direction, but only deletion of RPD3 results in an increase in switching in the opaque-to-white direction, suggests the mechanisms for switching in the two directions are distinct. Other observations support this conclusion. First, an increase in temperature from 25 to Ͼ37°C leads to the mass conversion of opaque to white, but has no apparent effect on switching in the white-to-opaque direction (32,38,42,53,55). A decrease in temperature also selectively stimulates switching in the opaque-to-white direction only (42). Second, Kolotila and Diamond (27) demonstrated that leukocytes and oxidants selectively stimulated switching in the white-toopaque direction, like TSA treatment and deletion of HDA1 (26). Third, misexpression of the white-phase-specific gene WH11 selectively stimulates switching in the opaque-to-white direction, but not the white-to-opaque direction (28). Finally, UV treatment stimulates switching in both the white-toopaque direction and the opaque-to-white direction (32), as we have demonstrated here for RPD3 deletion. The effects of these different treatments on the frequencies of switching in the two directions are summarized in Fig. 9. It should be noted that although the characteristics of switching in the alternate directions differ, we cannot rule out the likely possibility that the reversible switch event occurs at a single locus.
Since deacetylases play a role in the repression of gene expression (20,40,59,64), we must first entertain the possibility that Hda1p and Rpd3p suppress switching by removing acetyl groups from histone 3 and/or 4 at the actual site of the switch event. In this case, a switch event would be affected by increased or altered acetylation patterns through one of the many histone acetyltransferases in an activation complex at the  9. Diagram of the different treatments, environmental conditions, and deacetylase mutations that cause increases (ϩ) in the frequency of switching in the white-to-opaque and opaque-to-white directions. UV treatment (31) and treatment with neutrophils and oxidants (27) are shown. An up arrow shows a shift to high temperature (temp) (32,38,42,51), and a down arrow shows a shift to low temperature (42). misexp., misexpression of WH11 in the opaque phase (28).

VOL. 183, 2001
DEACETYLASES IN CANDIDA SWITCHING master switch locus. However, we must also entertain the alternate possibility that one or both histone deacetylases function to suppress one or more trans-activators of switching, which act directly upon the master switch locus. In this case, deletion of the deacetylase leads to up-regulation of the activator and, in turn, to up-regulation of the primary switch event.
The roles of HDA1 and RPD3 in the regulation of phasespecific gene expression and other deacetylases. The histone deacetylases and SIR2, which was recently demonstrated to have NAD-dependent histone deacetylase activity (23), play roles in the expression of large numbers of different genes in S. cerevisiae. Recently, Bernstein et al. (5) performed a bioinformatic analysis of genes up-regulated at least 1.5-fold in rpd3, hda1, and sir2 mutants. They found that while the gene coding for Hda1p plays a more prominent role in regulating carbon metabolite and carbohydrate transport and utilization, RPD3 plays a role in cell cycle progression, and SIR2 plays a role in amino acid biosynthesis. There were also overlapping roles of the three genes in these three functional categories. In the more limited study of gene expression that we performed, we found no instance in which deletion of a deacetylase resulted in up-regulation of a phase-specific gene (Table 4). We did, however, observe down-regulation. In the case of the HDA1 deletion, white-phase expression of EFG1 was down-regulated, but there was no significant effect on the white-phase-specific expression of WH11 or HOS3 or the opaque-phase-specific expression of OP4, SAP1, and SAP3 (Table 4). There was also no effect on the white-phase-enriched expression of MCM1 (Table  4). In the case of RPD3 deletion, there was a small decrease in the expression of the white-phase genes WH11 and EFG1, as well as a more dramatic decrease in the opaque-phase genes OP4, SAP1, and SAP3 (Table 4). As in the case of the HDA1 deletion, there was no effect on the phase-regulated expression of MCM1 (Table 4). Neither deletion of HDA1 nor deletion of RPD3 affected the phase-regulated expression of the four other deacetylases (Table 4).
Our observation that the deletion of HDA1 or RPD3 led only to down-regulation, not up-regulation, of select phase-specific genes is interesting, but not unique. Deletion of RPD3 in S. cerevisiae leads to increased repression of reporter gene expression at the mating type loci, telomeres, and ribosomal DNA (9,20,39), all loci under the regulation of Sir proteins (9,17,44). Increased silencing at these loci results presumably from up-regulation of silencing genes usually suppressed by the activity of the deacetylases. We assume that the down-regulation of select phase-specific genes in the deacetylase deletion mutants of C. albicans is effected by similar mechanisms.
Conclusion. We have presented evidence that the relationship between high-frequency phenotypic switching and the deacetylases is complex. First, both HDA1 and RPD3 play roles in the suppression of switching: the former in the selective suppression of switching in the white-to-opaque direction and the latter in both the white-to-opaque and opaque-to-white directions. The results do not distinguish in either case whether suppression is mediated by deacetylation directly at the site of the switch event or indirectly at the site of a gene that encodes an activator of the switch event. Second, switching affects the expression of deacetylases. For four of the five deacetylases, expression is higher in the white phase. Third, deletion of HDA1 or RPD3 in no case results in up-regulation of a phase-specific or phase-enriched gene, and in select cases, it results in downregulation. Our results therefore demonstrate that the deacetylases play distinct roles not only in the suppression of switching, but also in the activation of select phase-regulated genes, in the latter case presumably through the down-regulation of suppressor genes. ϩϩϩϩ ϩϩϩϩ ϩϩϩϩ ϩϩϩ ϩϩϩ ϩϩϩ a ϩϩϩϩ, full expression; ϩϩϩ, slightly reduced expression; ϩϩ, very reduced expression; ϩ, extremely reduced expression; Ϫ, no expression. A down arrow in parentheses represents significant down-regulation.


extfieldsofstudy: Biology Medicine

text:
Specific ligands of the peripheral benzodiazepine receptor induce apoptosis and cell cycle arrest in human colorectal cancer cells

The peripheral benzodiazepine receptor (PBR) has been implicated in growth control of various tumour models. Although colorectal cancers were found to overexpress PBR, the functional role of PBR in colorectal cancer growth has not been addressed to date. Using primary cell cultures of human colorectal cancers and the human colorectal carcinoma cell lines HT29, LS174T, and Colo320 DM we studied the involvement of PBR in the growth control and apoptosis of colorectal cancers. Both mRNA and protein expression of PBR were detected by RT-PCR and flow cytometry. Using confocal laser scanning microscopy and immunohistochemistry the PBR was localized in the mitochondria. The specific PBR ligands FGIN-1-27, PK 11195, or Ro5-4864 inhibited cell proliferation dose-dependently. FGIN-1-27 decreased the mitochondrial membrane potential, which indicates an early event in apoptosis. Furthermore, FGIN-1-27, PK 11195 or Ro5-4864 increased caspase-3 activity. In addition to their apoptosis-inducing effects, PBR ligands induced cell cycle arrest in the G 1/G 0-phase. Thus, our data demonstrate a functional involvement of PBR in colorectal cancer growth and qualify the PBR as a possible target for innovative therapeutic approaches in colorectal cancer. © 2001 Cancer Research Campaign http://www.bjcancer.com

Tumour cell growth is characterized by an imbalance between cell division and cell death. A common step in tumorigenesis is the loss of the capability of cells to undergo apoptosis. The localization of PBR in the mitochondrial membrane and its implication in the permeability transition pore (Zorov, 1996;Fennell et al, 2001) suggest that the PBR takes part in the regulation of the mitochondrial permeability and induction of apoptosis. However, the role of PBR in the regulation of apoptosis is not yet understood. Specific PBR ligands were shown to either induce apoptosis directly (Marchetti et al, 1996a;Tanimoto et al, 1999;Fischer et al, 2001), or to facilitate apoptosis by inhibiting the antiapoptotic effects of BCL-2 (Hirsch et al, 1998;Larochette et al, 1999;Ravagnan et al, 1999). In contrast, apoptosis-protective effects of PBR ligands have been reported as well (Bono et al, 1999). Besides their putative role in the regulation of apoptosis,

Specific ligands of the peripheral benzodiazepine receptor induce apoptosis and cell cycle arrest in human colorectal cancer cells
PBR ligands were shown to induce cell cycle arrest in the G 1 /G 0and G 2 /M-phase in breast carcinoma and melanoma cell lines (Landau et al, 1998;Carmel et al, 1999). Similarly, the new benzazepine BBL22, classified as a PBR-specific ligand, induced arrest in the G 2 /M-phase in different tumour cell lines of epithelial or haematopoietic origin followed by an induction of apoptosis (Xia et al, 2000).
As PBR was shown to be overexpressed in colorectal tumours (Katz et al, 1990b), PBR may well play a functional role in colorectal carcinogenesis. However, an involvement of PBR in colorectal cancer growth has not been investigated so far. In this study we show that specific PBR ligands inhibit proliferation of colorectal cancer cells, which is associated with an induction of apoptosis and cell cycle arrest.

Cell culture
The human colorectal adenocarcinoma cell lines HT29 and Colo320 DM (neurorendocrine-differentiated) were grown in RPMI 1640 medium supplemented with 10% fetal calf serum. The human colorectal adenocarcinoma cell line LS174T was grown in Dulbecco's minimal essential medium supplemented with 10% fetal calf serum. Cell lines were cultured in a humidified atmosphere containing 5% CO 2 at 37˚C .
Surgically resected specimens of primary colorectal carcinomas (3 rectum, 3 colon) were obtained from 4 female and 2 male patients who underwent surgery in the Department of Surgery, Benjamin Franklin University Hospital, Free University Berlin. The human tumour material was used according to the standards set by the Ethical Committee of the Benjamin Franklin University Hospital, Free University of Berlin. The age of the patients ranged from 23 to 78 years. Primary cell cultures were prepared by mechanical dissection using a Medimachine with 50 µm Medicons (Becton Dickinson, Heidelberg, Germany) according to the manufacturer's instructions. Isolated cells were cultured in RPMI 1640 medium supplemented with 10% fetal calf serum, penicillin (100 U ml -1 ), and streptomycin (100 µg ml -1 ) and cultured in a humidified atmosphere containing 5% CO 2 at 37˚C .

Reverse transcriptase chain reaction (RT-PCR)
RNA isolation, reverse transcription, and PCR reactions were carried out as described . Amplification consisted of 38 cycles with the following conditions for denaturation, annealing and extension: 94˚C for 45 s, 60˚C for 45 s and 72˚C for 2 min. The primers for amplification of cDNA were designed using the PRIMER program (Whitehead Institute for Biomedical Research, Cambridge, MA, USA) based on the cDNA sequences of human 18 kDa PBR subunit obtained from the GenBank (accession number: M36035). Primers for amplification were: Forward: 5′-CACGCTCTACTCAGCCATGG-3′ Reverse: 5′-GCAGTAGTTGAGTGTGGTCGC-3′ The expected PCR product size was 298 bp. PCR products were sequenced on an ABI 310 sequencer (Applied Biosystems, Foster City, CA, USA) and specificity of the transcripts was confirmed by the NCBI (National Center for Biotechnology Information, Bethesda, MD, USA) BLASTIN 2.0 search program (Altschul et al, 1997).

Immunofluorescence labelling and confocal microscopy
Cells were plated on glass cover slips for 24 h. Cells were stained with the mitochondrial dye CMTMRos (100 nM in medium, 30 min, 37˚C) (Molecular Probes, Eugene, OR, USA). After washing with PBS, the samples were fixed by incubation with 4% paraformaldehyde in PBS for 30 min at room temperature, washed with PBS, and permeabilized with cold methanol (20 min, -20˚C) followed by acetone (10 s, -20˚C). Cells were incubated for 1 h at room temperature with primary anti-PBR monoclonal antibody 8D7 (5 µg ml -1 in PBS) (Dussossoy et al, 1996) or isotypic control mouse IgG1 (DAKO, Hamburg, Germany), respectively. Thereafter, cells were washed with PBS twice and incubated with 4 µg ml -1 secondary Alexa TM 488-labelled goat-anti-mouse IgG antibody (Molecular Probes, Eugene, OR, USA) for 1 h at room temperature. Fluorescence and transmission images were obtained using the inverted confocal microscope LSM 410 with a 63x /1,2 W Korr objective (Zeiss, Oberkochen, Germany).
Immunohistochemistry 23 paraffin-embedded colorectal carcinomas were studied. These were 1 well differentiated, 13 moderately differentiated, and 9 poorly differentiated cancers. 2 µm sections of the tumours were deparaffinized and rehydrated in a series of alcohol solutions of decreasing concentrations (Grabowski et al, 2000(Grabowski et al, , 2001. Then sections were transferred into a robotic machine (Chemo-mate, DAKO, Heidelberg, Germany) and staining procedure was automatically performed under following standard conditions. Sections were incubated with the anti-PBR antibody 8D7 (0.5 µg ml -1 ) for 30 min at room temperature. After washing samples were incubated with the second anti-mouse IgG (1:20 dilution) (DAKO, Heidelberg, Germany) for 30 min at room temperature. The APAAP complex (DAKO, Heidelberg, Germany) was incubated for 30 min at a dilution 1:50 in RPMI medium containing 10% fetal calf serum and 1% sodium azide. Staining was detected using fast-red system (DAKO, Heidelberg, Germany) and samples were counterstained with haemalaun.

Cell proliferation assay
The ability of the PBR ligands to modulate cell proliferation was studied using the crystal violet method (Gillies et al, 1986). In brief, cells were seeded on 96-well plates at a density of 1000 cells well -1 (HT29 and Colo320 DM) or 2000 cells well -1 (LS174T). After 72 h the PBR ligands FGIN-1-27, PK 11195 (Tocris, Bristol, UK), or Ro5-4685 (Sigma, Deisenhofen, Germany) were added at concentrations of 1 nM to 100 µM. The FGIN-1-52 (Kozikowski et al, 1993) or the GABA A receptor ligand clonazepam (Sigma, Deisenhofen, Germany) were used as controls. Each concentration group consisted of 10 wells. The incubation medium was changed every day. Cell quantification was performed after 0, 24, 48, 72 and 96 h of incubation. Cells of each well were washed (200 µl PBS) and fixed (100 µl, 1% glutaraldehyde in PBS for 15 min at room temperature). After another washing step (200 µl PBS), cells were stained with 0.1% crystal violet in PBS for 30 min at room temperature. The unbound dye was removed by washing with H 2 O for 30 min. Crystal violet which had absorbed onto the cells was solubilized with 100 µl 0.2% Triton-X-100 in PBS for at least 24 h at 37˚C. The crystal violet containing solution was measured spectrometrically at 570 nm using an ELISA-Reader. In the range of 500 to 20 000 cells well -1 , the measured extinction was a linear function of the cell number .

Cell cycle analysis
Cell cycle analysis was performed by the method of Vindelov and Christensen (Vindelov and Christensen, 1990). 5 × 10 5 cells per well were cultured for 24 h and then exposed to PBR ligands for another 24 h. Cells were trypsinized, washed and the nuclei were isolated using CycleTest PLUS DNA Reagent Kit (Becton Dickinson, Heidelberg, Germany). According to the manufacturer's instructions DNA was stained with propidium iodide. The DNA content of the nuclei was detected by flow cytometry and analysed using CellFit software (Becton Dickinson, Heidelberg, Germany).

Statistical analysis
Comparison of multiple means was performed with nonparametric ANOVA. Comparison of individual drug treatments to control treatments was performed with the unpaired, 2-tailed Mann-Whitney U test for proliferation, ∆Ψ M measurements, and caspase-3 activity experiments. Data are expressed as mean percentage of control ± SEM. For cell cycle analysis the unpaired student's t-test was used. P values were considered to be significant at < 0.05.

PBR expression in colorectal cancer cells
The mRNA expression of PBR was investigated by RT-PCR. Both in HT29 cells and in primary cultures of colorectal carcinomas, PBR mRNA was detected ( Figure 1A). The specificity of the obtained PCR products was confirmed by direct sequencing of cDNAs. HT29 cells did not express any β-chain of the GABA Areceptor (data not shown), which are essential for benzodiazepines binding to the GABA A -receptor. To investigate the expression of the PBR protein, cells were stained with the monoclonal anti-PBR antibody 8D7 and fluorescence was analysed by flow cytometry. The expression of PBR protein was detected in primary cell cultures of colorectal cancer cells ( Figure 1B), as well as in the cell lines HT29, LS174T and Colo320 DM ( Figure 1C-E). In nonpermeabilized cells, no specific 8D7 fluorescence was observed indicating an intracellular localization of PBR ( Figure 1C-E).

PBR localization in mitochondria
To further characterize the subcellular localization of PBR, HT29 cells were simultaneously stained with the mitochondrial dye CMTMRos and the anti-PBR antibody 8D7. Confocal laser scanning microscopy imaged the subcellular distribution of fluorescence. Green fluorescence of 8D7-labelled PBR was detected within the cytoplasm but neither in the cell membrane nor in the nucleus (Figure 2A, E). The staining pattern was comparable to the one obtained by the mitochondrial dye CMTMRos ( Figure 2B, I). Superposition of the images of 8D7-labelled PBR and CMTMRos stained mitochondria resulted in a yellow colour ( Figure 2C) indicating a co-localization of PBR and mitochondria. Likewise, PBR was localized in the mitochondria in the colorectal cell lines LS174T and Colo320 DM as well as in primary cultures of 6 resected colorectal cancers (data not shown). To further verify the mitochondrial localization of PBR, subcellular PBR expression was immunohistochemically determined in 23 other colorectal cancers. The differentiation grade of these tumours ranged from well (n = 1), moderately (n = 13), to poorly (n = 9) differentiated.
In all 23 cancers, the specific PBR staining was observed unevenly distributed within the cytoplasm ( Figure 2K, L). No specific PBR staining was observed in the plasma membrane nor in the nuclei of cells, indicating that in colorectal cancers PBR is located in the mitochondria.

Induction of mitochondrial alterations by specific PBR ligands
PBR is located mainly in the outer mitochondrial membrane and is thought to form the permeability transition pore, which plays an important role in the regulation of ∆ψ M . Therefore, we investigated if specific PBR ligands can modulate the ∆ψ M and mitochondrial volume. FGIN-1-27 significantly depolarized the mitochondria of HT29 cells in a dose-dependent manner as assayed by the  Figure 4A). The results were confirmed by using the mitochondrial dye CMTMRos. CMTMRos stained HT29 cells displayed mean fluorescence of 83 ± 6% (10 µM), 68 ± 4% (50 µM), and 64 ± 6% (100 µM) of control, respectively, when incubated for 6 h with FGIN-1-27. Simultaneously, an increase of mitochondrial volume of HT29 cells and primary cultures was detected ( Figure 4A), indicating that FGIN-1-27 not only induced a mitochondrial membrane depolarization but also an increase in mitochondrial volume. In contrast to FGIN-1-27, PK 11195 or Ro5-4864 showed no effects on HT29 cells after 1 h to 16 h of incubation. Moreover, the non-PBR-specific compound FGIN-1-52 did not alter the ∆ψ M (data not shown). In control experiments the K + ionophore valinomycin was used to reduce the JC-1 590 nm fluorescence (Cossarizza et al, 1993) resulting in a mean 590 nm fluorescence of one third of corresponding control.

DISCUSSION
A tightly regulated balance between cell division and cell death is a prerequisite for normal tissue development. Impaired regulation of either process can lead to tumorigenesis. In this study we show that specific PBR ligands induce both an arrest of the cell cycle and an increase in apoptosis in human colorectal cancer cells. Our data indicate that PBR might be involved in both the regulation of cell division and cell death, suggesting an important role of PBR    (Kozikowski et al, 1993) in colorectal cancer growth. Moreover, we provide evidence that the cell cycle interfering and proapoptotic effects of specific PBR ligands are associated with an inhibition of proliferation of colorectal cancer cells.
The PBR is mainly located in the outer membranes of mitochondria (Gavish et al, 1989). Recently, a perinuclear or nuclear localization has been described in breast cancer and glioma cells, which was shown to be associated with an aggressive tumour type (Hardwick et al, 1999;Brown et al, 2000). Using confocal laser scanning microscopy and immunohistochemistry, we studied PBR expression in human colorectal cancers of distinct grades of differentiation, and in several colorectal cancer cell lines. In all 29 cancers and in the 3 cell lines studied, PBR was found to be localized in the mitochondria. Thus, PBR is primarily targeted to mitochondrial membranes of either well, moderately or poorly differentiated colorectal cancers.
To investigate the functional involvement of PBR in colorectal cancer growth we used specific exogenous ligands. These ligands were shown to interfere with the regulation of both the cell cycle and apoptosis. The specific PBR ligands induced cell cycle arrest. Upon ligand treatment, the proportion of colorectal cancer cells in the G 1 /G 0 -phase increased markedly, indicating that the cells Means as percentage of control ± SEM of 4 independent experiments are shown. In addition, caspase-3 activity was determined of primary tumour cells of 6 patients after incubation with 50 µM FGIN-1-27. Mean ± SEM are shown. *P < 0.05 versus control without agent stayed longer in this phase of the cell cycle. This suggests that PBR ligands act at the classical G 1 checkpoint, preventing cells from entering the S-phase. All 3 PBR ligands studied induced cell cycle arrest at this restriction point, suggesting a common signalling pathway. Cell cycle interfering effects of PBR ligands have been shown previously. In breast cancer PBR ligands induced a cell cycle arrest at both major restriction points, the G 1 /S-and the G 2 /M-junction (Carmel et al, 1999), whereas in lung and melanoma cells an accumulation in the G 2 /M-phase was observed (Camins et al, 1995;Landau et al, 1998). These differences may reflect tissue-specific PBR signal transduction.
In addition to their cell cycle-arresting effects, we here showed for the first time that specific PBR ligands can induce apoptosis in non-haematopoietic cancer cells. The PBR ligands FGIN-1-27, PK 11195 and Ro5-4864 stimulated caspase-3 activity in the 3 cell lines and in all 6 primary cultures of colorectal cancers. Recently, mitochondrial permeability transition has been described as initiating event in the process of apoptosis . Within the mitochondria, PBR is known to be located in the membrane forming the permeability transition pore (Zorov, 1996). This pore is considered to play a central role in the initiation of apoptosis (Zamzami et al, 1995). In this study, the PBR ligand FGIN-1-27 decreased ∆ψ M and caused mitochondrial swelling. Whether FGIN-1-27 induced mitochondrial alterations are a prerequisite for the induction of apoptosis in colorectal cancer cells is not yet understood. It has been reported that mitochondrial permeability transition occurs independently of the initiation of the apoptotic cascade (Lemasters et al, 1998). The actual role of the mitochondrial transition in PBR ligand-induced apoptosis needs to be further investigated.
In this study, we showed that the cell cycle-arresting and proapoptotic effects resulted in an inhibition of cell proliferation.
In line with previous studies in a variety of tumour models (Garnier et al, 1993;Neary et al, 1995;Landau et al, 1998;Carmel et al, 1999;Beinlich et al, 1999), the specific PBR ligands FGIN-1-27, PK 11195 and Ro5-4864 exhibited antiproliferative effects in colorectal cancer cells at micromolar concentrations. To clarify the PBR specificity of the observed effects, we used the benzodiazepine clonazepam or the indoleacetamide FGIN-1-52. Despite their very similar structure to Ro5-4864 or FGIN-1-27, respectively, neither clonazepam nor FGIN-1-52 affected apoptosis, cell cycle or the proliferation of human colorectal cancer cells. This indicates that the effects are specific for PBR. Nevertheless, there is a quantitative discrepancy between the micromolar ligand concentrations necessary to inhibit cell proliferation and the nanomolar-binding affinities. This discrepancy might be due to the existence of a putative 'low-affinity PBR' which has previously been suggested to exist in other tissues (Pawlikowski et al, 1988;Kunert-Radek et al, 1994). This putative low-affinity binding site has been implied in the regulation of calcium channels (Cantor et al, 1984;Taft and DeLorenzo, 1984), suggesting that PBR ligands might affect calcium influx responsible for growth inhibition. However, other factors including cellular absorption, ligand metabolism, as well as differences in the conditions used for the cell proliferation and binding assays (Wang et al, 1984b) have yet to be ruled out.
Even though the exact targets and mechanisms of the proapoptotic and antiproliferative effects of the specific PBR ligands remain to be elucidated, the ability of these ligands to modulate both apoptosis and cell cycle regulation qualify them as promising agents for innovative treatment strategies in colorectal cancer disease.