Molecular Analyses of the Genus Ilex (Aquifoliaceae) in Southern South America, Evidence from AFLP and ITS Sequence Data

In order to clarify the relationships among southern South American (sSA) representatives of the genus Ilex, an amplified fragment length polymorphism (AFLP) analysis was accomplished. In addition, the phylogenetic relationships of the species were studied using ribosomal internal transcribed spacer (ITS) sequence data alone and in combination with AFLP data, taking into account the possible existence of paralogous sequences and the influence of alignment parameters. To explore stability of phylogenetic hypotheses, a sensitivity analysis was performed using 15 indel-substitution models. Within each species assayed, the AFLPs allowed the recognition of several diagnostic bands. Furthermore, the AFLP analysis revealed that individuals belonging to the same morpho-species formed coherent clades. In addition, some cases of geographical association were noted. Studies on ITS sequences revealed divergence between data obtained herein and sequence data downloaded from GenBank. The sensitivity analyses yielded different interspecific hypotheses of relationships. Notwithstanding, analyses of the ITS data alone and in combination with AFLPs, rendered clades stable to variation in the analytical parameters. Topologies obtained for the AFLPs, the ITS data alone and the combined analyses, demonstrated the existence of a group formed by I. argentina, I. brasiliensis, I. brevicuspis, L integerrima, and L theezans, and that L dumosa and I. paraguariensis were distantly related to the former. Incongruence with traditional taxonomical treatments was found.

trees to creeping prostrate plants or epiphytic shrubs. It has more than 400 species, dispersed mainly in America and Eurasia, but also in Oceania and Africa (Giberti, 1994b).
South America is considered one of the main areas of diversification of Ilex, together with East Asia (Loesener, 1942;Cuenoud et al., 2000). In southern South America (sSA) the species are mostly found in areas with tropical or subtropical climates, mainly in northeastern Argentina, southeastern Brazil and eastern Paraguay. Twelve species have been described for sSA (Loesener, 1901;Lillo, 1911;Edwin and Reitz, 1967;Giberti, 1998), namely, Ilex affinis Gardner, I. argentina Lillo, L brasiliensis (Sprengel) Loes., I. brevicuspis Reissek, I. chamaedryfolia Reissek, L dumosa Reissek, L integerrima (Vell.) Reissek, L microdonta Reissek, L paraguariensis A. St. H I. pseudobuxus Reissek, I. taubertiana Loes., and I theez C. Martius ex Reissek. Of those mentioned, I. paraguariensi is the most relevant from an economic perspective. In Arge tina, southern Brazil, Paraguay, and Uruguay, the aerial pa of this plant are used to prepare a very popular beverage ca mate, which is highly appreciated for its flavor and stimulat properties due to its caffeine and theobromine contents (F et al., 2001).
About a hundred years after Linnaeus (1753) published th diagnosis for the genus Ilex, Gray (1856) attempted to esta lish an infrageneric classification. In 1861, Reissek publishe the greatest species inventory for the time, followed by Hoo (1862) who reported the occurrence of about 145 speci worldwide. Loesener (1891Loesener ( , 1901Loesener ( , 1908Loesener ( , 1942 published series of seminal studies recording almost 280 species of Il dispersed globally, and recognizing several subgeneric rank although under an unconventional hierarchical arrangemen Since then, several authors have recognized some inaccurac in Loesener's system and proposed emendations in the system atics of the genus (Rehder, 1927;Hu, 1949Hu, , 1950Krtissman 1962;Baas, 1975;Lobreau-Callen, 1975;Martin, 1977). Sp cies delimitation in Ilex, based on overall morphological sim ilarity, is complicated due to the various concepts used different authors to define the taxa and to the lack of comp information, particularly for evergreen species (Galle,199 In order to investigate the relationships among the speci of Ilex found in sSA, we implemented the amplified fragme length polymorphism (AFLP) technique (Vos et al., 1995). In addition, we studied the phylogenetic relationships of the sp cies using ribosomal internal transcribed spacer (ITS) quence data alone and in combination with AFLP data, takin into account the possible existence of paralogous sequen and the influence of alignment parameters. ' Manuscript received 11 December 2003;revision accepted 17 September 2004. Beside the general use of AFLP markers in the estimation of genetic diversity and differentiation among individuals, populations, and plant species, they have also been successfully used in phylogenetic analysis between closely related species (Cervera et al., 1998;Kardolus et al., 1998;Mace et al., 1999;Negi et al., 2000;Ren and Timko, 2000;Koopman et al., 2001). Because of the rapidity with which reliable high resolution markers can be generated, the AFLP technique is considered a powerful tool for molecular studies (Mueller and Wolfenbarger, 1999). In addition, the high ratio of polymorphism generated per PCR experiment (the multiplex ratio) (Powell et al., 1996) and the random distribution of most AFLP bands across the genome (Zhu et al., 1998) are considered key features of this technique. Several drawbacks have been detected for AFLPs (Robinson and Harris, 1999), which introduce bias into the estimates, whether analyzed by phenetics or cladistics.
A fundamental requirement in phylogenetic studies based on nucleic acid or protein sequences is that the genes compared are orthologous (originating from organismal cladogenesis) (Alvarez and Wendel, 2003). Plant rDNA consists of thousands of copies or paralogues (derived from gene duplication events) spanning several loci, which evolve through gene conversion, unequal crossing over, and perhaps repeat amplification (Baldwin et al., 1995). ITS sequences have proven useful in species-level phylogenetic studies of a wide range of taxa (Alice and Campbell, 1999;Downie et al., 2000;Lia et al., 2001, among others). Within-individual rDNA polymorphisms may occur when concerted evolution is not fast enough to homogenize repeats in the face of high rates of mutation and/or recent interspecific hybridization (Campbell et al., 1997). Intraspecific variation concerning ITS sequences has been detected in several plant groups (Baldwin et al., 1995 and references therein;Buckler and Holtsford, 1996a, b;Buckler et al., 1997;Campbell et al., 1997;Hartmann et al., 2001;Mayol and Rosell6, 2001). Active ITS regions have functional constraints that can help discriminate between functional and nonfunctional paralogous (i.e., pseudogene) sequences (Buckler et al., 1997). For instance, nonfunctional ITS paralogues could accumulate random substitutions, which will destabilize hairpins and reduce secondary structure stability. When nonfunctional paralogous ITS sequences are unknowingly included in the phylogenetic analyses (to the exclusion of appropriate orthologous comparisons), the resulting gene tree may confound organismal divergence events with a tracking of the history of gene duplication (Alvarez and Wendel, 2003). Erroneous assessment of orthology and paralogy will lead to phylogenetic incongruence. Hence, discrimination of inactive paralogues and pseudogenes becomes crucial.
Another critical step in molecular phylogenetic studies is the alignment of sequence data, in which primary homology statements are being established. However, it is typically performed without considering the effects of analytical parameters, such as indel and base transformation values, on the phylogenetic inference. In the case of multiple sequence alignments, different parameters may yield alternative alignments, and these primary homology hypotheses may support distinct topologies (Doyle and Davies, 1998). Indeed, Morrison and Ellis (1997) showed that different multiple sequence alignment methods were responsible for the variation encountered among resultant phylogenetic trees topologies, and this variation was more important than the one found when different methods of phylogenetic reconstruction were applied to the same align-ments. Because parameter choice is arbitrary but unavoidable when doing algorithmic DNA sequence comparisons, its exploration becomes essential (Giribet and Wheeler, 1999;Giribet and Ribera, 2000). Herein, we used ITS sequences to study the phylogenetic relationships of sSA taxa and to explore the influence of applying different insertion/deletion (indel or gap) and substitution models (sensitivity analysis sensu Wheeler, 1995) on phylogenetic inference. Phylogenetic analyses were accomplished through direct optimization of DNA sequences (Wheeler, 1996). The sensitivity analysis is a way to explore the data and to discern between stable relationships (those supported throughout the parameters range) and unstable relationships (those appearing only under particular conditions), allowing the formulation of more robust hypotheses (Giribet and Ribera, 2000) than in customary phylogenetic analyses. Direct optimization is a method that does not disconnect the alignment step from the tree-building step and that generalizes phylogenetic character analysis to include indel events (Giribet, 2001(Giribet, , 2003. By doing this, indels appear as transformations (not as states) linking ancestral and descendent nucleotide sequences (Giribet and Wheeler, 1999).

MATERIALS AND METHODS
Specimens studied-The list of sSA flex taxa employed, together accession numbers, geographical origin and GenBank accession nu ITS sequences obtained in this study are presented in Appendix 1 plemental Data accompanying the online version of this article). Plan was principally obtained from the Estaci6n Experimental INTA C (EEINTACA) Germplasm Bank (Misiones, Argentina). Leaf sam also collected from natural populations from Northern Argentina specimens of these samples were deposited at the BACP Herbariu de Estudios Farmacol6gicos y Botfinicos, CEFyBO, Buenos Aires, A Furthermore, leaf samples of commercial lines of mate were obta the Establecimiento Las Marias (Corrientes, Argentina). Young le were collected from healthy plants (without evidence of fungal or fection) and were silica-gel preserved. Additional 42 ITS sequen downloaded from GenBank (http://www.ncbi.nlm.nih.gov); for t graphical origin and accession numbers see Appendix 2 (Data Su accompanying the online version of this article). Helwingia japon Thunberg ex A. Murray) Dietrich and Helwingia chinensis Batalin sen as the outgroup. According to several molecular studies (Cu6n 2000; Powell et al., 2000;Savolainen et al., 2000;Albach et al., 200 et al., 2002), Helwingia and Ilex are closely related members of Aq Loeseners' (1891,1901,1908,1942) original dispositions of subgen were used in this work (Fig. 19). The species flex microdonta, I. pseu and I. taubertiana, were only included in the sequence analyses b cessions of these taxa were unavailable in sufficient number for AFLP analysis. Only two sSA taxa remain to be studied, namely flex affinis and I. chamaedryfolia, because these have not been found in the field since the the amplification products were too dense to allow reliable scor natively, generated too few amplification products. Selective a products were mixed with an equal volume of dye reagent (98% amide, 10 mM EDTA, 0.025% [w/v] bromophenol blue and 0 xylene cyanol). Five microlitres were separated by electrophoresis S2 apparatus (Gibco BRL Sequencing System, Life Technolog 6% (w/v) polyacrylamide gels containing 5 M urea, in 1 X TBE mM Tris, 89 mM boric acid, 2 mM EDTA, pH 8). A 30-330 b Ladder (Gibco BRL, Life Technologies) size marker was inclu three times in each electrophoresis run. Thus, the size of AFLP ranged from 50 to 330 bp. Gels were stained with silver nitrat al., 1991).
ITS amplification-Amplification of ITS 1, 5.8S gene, and ITS ried out using primers ILEXFP (5'-AACAAGGTTTCCGTAG Powell et al., 2000) and ITS-4 (White et al., 1990) (Hall, 1999). Boundaries of the coding and spacer r determined by comparison with published sequences of Ilex. In ge two or more orthologous ITS sequences were available per sS Appendix 1; Supplemental Data accompanying the online versio ticle), these were used to produce consensus sequences in BioEdi to speed-up computational times of phylogenetic analyses.
Data analyses-AFLP-Air-dried gels were scanned and bandin were registered with the Gel-Pro 4 Analyzer Program (Media Maryland, USA). Each AFLP band, regardless of its relative int considered as a dominant allele at a unique locus. In some cases were scored as missing data because character states could not be unambiguously. The data were extracted as a table as either pr absent (0). Monomorphic bands (bands present in all individuals and diagnostic bands (monomorphic bands exclusively present in were discriminated within each species and across the entire d plots of mean values; SE and SD were calculated and plotted with t STATISTICA '99 edition (Kernel release 5.5 A, StatSoft, 1999) calculations, accessions with missing data in one or more primer c were excluded. Pairwise genetic distance (D) between all ind estimated according to the complementary value of Nei and Li' ilarity coefficient, implemented in PAUP* (Swofford, 1998): D [n; + ni]), where n,i is the number of shared fragments between uals i and j, and n, and n, are the total number of fragments in and j, respectively. From these values, mean intra-and interspeci values were obtained. Furthermore, the AFLP data set for 120 belonging to eight taxa was analyzed using maximum parsimon Initial trees in the heuristic search by tree-bisection-reconn branch-swapping algorithm were built stepwise with 500 rand sequence replicates. Characters we other settings were defaults. The bootstrap analysis (Felsenstein,1 all other settings as described earl TreeView (Page, 1998). The unroote ed using I. paraguariensis to be c ITS sequence studies-We charac sSA Ilex species by examining t stabilities, and nucleotide substit proach advocated by Bailey et a presumed active sequences from first followed the procedure pro fore, pairwise nucleotide diverg sequences obtained herein and the belong to the same sSA Ilex taxa mura two-parameter model (Kim plemented in BioEdit (Hall, 1999) using BioEdit (Hall, 1999) for leng ence of a structural motif in ITS Liu and Shardl, 1994). This conser the hairpin and loop structures t zymatic processing of the riboso more, optimal and suboptimal fold values (AG) in Kcal/mol, of ITS 1 were made at the Quickfold web gram, version 3.1. In all cases, fo search within 5% of the thermo quence could adopt more than one AG values were recorded for each possible (optimal and suboptimal f attempts were made to study th were extended to the remaining r Phylogenetic and sensitivity an timality criterion, as implement POY, version 3.0.11 (Gladstein an 2000;De Laet and Wheeler, 2003) in the POY user's manual and in partitioned. The first partition partition has ca. 340 bp. Approxi resenting missing data in several were excluded from the analyses iables was examined, i.e., gap cost tio cost. Precisely, we used gap = 1 (equal weights), 2 (transversi tions) and 4. All combinations of dance with the triangle inequal tree search was performed as f starting from the best of five ind subtree-pruning and regrafting ( per step, as specified in the foll -buildsperreplicate 5 -replicate Commands "-slop" and "-checks rived from the heuristic opera transformation matrices) were s lecularmatrix" with an argument f trees were calculated using the An implied alignment was gen "-impliedalignment". Approxima were calculated with a heuristic POY (command: "-bremer"). Fina were used to obtain 50% majorit and trees were deposited at Tr stability as defined by Giribet (2 affected by variation in the analy whereas the highest number of diagnostic bands was obtained for I. argentina and I. brevicuspis. In the box plot in Fig. 1, I. argentina gave the highest mean number of AFLP bands (mean = 175, SD = 29.9), whereas the lowest number of bands was found in I. brevicuspis (mean = 124, SD = 6.1).
Intra-and interspecific distance values derived from AFLP data are presented in Table 3. Excluding I. dumosa, for which the highest intraspecific distance value was obtained (0.1223, SD = 0.0408), remaining taxa had similar values. When interspecific mean distances were compared, I. integerrima and I. theezans appeared to be the least distant pair of taxa (0.1014, Table 3).   The strict consensus tree derived from AFLP data showed that all individuals of the same taxon grouped accordingly (Fig. 2). Accessions of L theezans yielded a coherent group (BV, bootstrap value; BV = 93%), closely related to the clade L integerrima (BV = 95%). Within L brasiliensis (BV = 66%), two subclades were distinguished, one comprised accessions from Argentina and Paraguay (BV = 83%), and the other subclade involved individuals from Brazil (BV = 82%). In contrast, specimens of L argentina formed a highly supported clade (BV -97%). The clade L brevicuspis (BV = 100%) showed several internal nodes unresolved and poorly supported. Within the I. dumosa clade (BV = 99%), some geographical structuring was detected. The first subclade (BV < 50%) involved individuals from Argentina and Paraguay representing both varieties of I. dumosa (I. dumosa var. dumosa and L dumosa var. guaranina, the latter indicated with an arrow in Fig. 2), while the other subclade comprised solely Brazilian accessions (BV 68%). Regarding L paraguariensis (BV = 100%), the majority of Argentine, Brazilian, and Paraguayan accessions of mate appeared intermingled (BV = 96%); still, most internal relationships were poorly supported. Commercial lines of mate (SI-16, SI-19, SI-167 and G-18) formed a clade with a specimen from Paraguay (BV = 59%). Another subclade involved lines SI-49, Y-383, and two accessions from Paraguay (BV 66%). The relationship between L brasiliensis, I. integerrima, and L theezans was highly supported by the data (BV -100%). Likewise, the relationship of L argentina and I. brevicuspis to the aforementioned taxa was also well-supported by the data (BV = 91%).

ITS sequence studies-Nucleotide divergence calculations
showed that in several cases sSA ITS sequences retrieved from GenBank were highly divergent from those obtained in this report (Table 4). For instance, similar divergence values (ca. 14-15%) were found for ITS sequences of L brasiliensis, I. integerrima, and I. pseudobuxus. While I. brevicuspis AJ492663 diverged with our ITS data in ca. 11%, accession AJ492656 of I. dumosa, was ca. 21% divergent from L dumosa AJ492657 and with the five sequences of L dumosa obtained herein. On the other hand, for L theezans, L microdonta, I. brevicuspis AJ492662, and I. dumosa AJ492657, the range of nucleotide variation detected was approximately under 5%.
In all cases, divergence values calculated among co-specific sequences obtained in this study were lower than 2%. The length, G + C content, free energy values of the RNA transcripts, and the presence of a conserved structural motif were also evaluated in ITS 1 and ITS 2 of sSA and worldwide distributed Ilex taxa (Table 5). For I. microdonta and I. thee-zans, no distinct differences were observed between our sequences and those from GenBank. Instead, major differences were found between our nucleotide data and sequences I brasiliensis AJ492661, I. brevicuspis AJ492663, I. dumosa AJ492656, L integerrima AJ492664, and I. pseudobuxus AJ492660. These downloaded data showed a markedly low G + C content and a considerable reduction in stability of the most probable secondary structure, both for ITS 1 and ITS 2. Ilex brevicuspis AJ492662 and I. dumosa AJ492657 were not completely comparable to other co-specific data due to a difference in length. Unfortunately, these sequences were incomplete lacking at least, 60 bp from the 5' end of ITS 1. Notwithstanding, their ITS 2 had comparable values to our sequences in all factors considered. Excluding these incomplete sequences, no significant differences in length were encountered. Among North American and East Asian sequences, some variations in the screened parameters were found (Table  5). In particular, the values obtained for ITS 2 of L kiusiana make us question its functionality. The lack of additional sequences impeded us from making a final decision.
The universal ITS 1 core motif of Liu and Schardl (1994) was screened in all 55 sequences, and some variants were detected. Particularly in Ilex, 13 variants of this sequence were recognized and arbitrarily named ( Table 6). The types B and D were the most frequent ones (21/53 and 11/53, respectively).
This scrutiny revealed that all six I. dumosa sequences (AY183479, AY183480, AY183481, AY183488, AY183489, and AJ492657) had a different motif that lacks the final adenine. As to the GenBank data (L brasiliensis AJ492661, L brevicuspis AJ492663, L dumosa AJ492656, I. integerrima AJ492664, and I. pseudobuxus AJ492660), these exhibited a motif that substantially departs from the core sequence, contributing to a less stable secondary structure (Tables 5 and 6).
Phylogenetic and sensitivity analyses-The input data set for the phylogenetic analysis comprised 53 ingroup sequences and two outgroup, also including putative divergent paralogues and pseudogenes. Analytical conditions 241, 421, and 1021 yielded a single MPT of length 2437, 2796, and 4707, respectively (not shown), whereas conditions 811 and 111 rendered two and three MPTs of length 3617 and 1346, respectively (not shown). The 50% majority rule consensus tree of Fig. 3, condenses the information display in all MPT encountered. Regarding sSA taxa, sequences of L argentina, I. brasiliensis, L brevicuspis, I. integerrima, I. microdonta, I. pseudobuxus, L taubertiana, L theezans, I. microdonta AJ492665, and L theezans AJ492666 appeared related in all the conditions assayed. However, under conditions 111 and 241, se-  AJ492660, could be considered as pseudogenes: the high nucleotide divergence with co-specific data, the low G + C content, the departure from the conserved motif in ITS 1, the less thermodynamic stability in the RNA structure, and their basal position on phylogenetic trees with respect to sSA counterparts. Therefore, only L theezans AJ492666, I. microdonta AJ492665, and L argentina AJ492655, retrieved from Gen-Bank, were further used as alternative functional ITS types.
The phylogenetic analysis restricted to active ITS sequences, comprised 46 ingroup sequences and two outgroup. Strict consensus trees and single MPTs obtained under each combination of gap cost and Tv/Ts ratio cost are shown in Figs. 4-18. The 50% majority rule consensus tree (Fig. 19, on the left), condenses the information display in all MPT encountered. The overall topology of the trees was not altered by the exclusion of the pseudogenes, although its resolution increased. Regarding sSA taxa, the clade composed by I. argentina, I. brasiliensis, L brevicuspis, L integerrima, L microdonta, L pseudobuxus, L taubertiana, and L theezans, including sequences I. microdonta AJ492665 and L theezans AJ492666 (hereafter, sSA clade) appeared in 100% of the conditions assayed. This clade received Bremer support values that ranged from 16 to 58, depending on the analytical condition. These values represent the third to fifth highest values compared to all values obtained in each condition. As before, the pairs (L brasiliensis + I. theezans) and (I. integerrima + I. theezans AJ492666) were found throughout the explored weights, although both subclades received Bremer supports from 7 to 24, which represent the seventh to less than tenth highest values.
Similarly, the subclade (L microdonta + I. pseudobuxus + I. microdonta AJ492665) appeared in 100% of the analyses, but its Bremer support was either the tenth highest value or lower. The formation and placement of subclade (I. brevicuspis + I. taubertiana) was dependent on analytical conditions; it was found in 80% of the conditions, and it formed a monophyletic group with (I. microdonta + I. pseudobuxus + L microdonta AJ492665) in only 53% of the conditions. The sister sequence to the sSA clade, in 67% of the parameters sets assayed, was represented solely by Asian L pedunculosa, but with low Bremer support values. Solely under analytical condition 811 (Fig.   13), I. paraguariensis, both varieties of I. dumosa, and East Asian I. crenata and L mutchagara, appeared as sister group to the sSA clade, with a Bremer support of 29 (seventh highest    Liu and Schradl (1994) found in 55 ITS 1 Ilex sequences and outgroup sequences (for detailed distribution of each type, see Table 5). value). In the analytical space explored, L dumosa varieties consistently yielded a clade with the second highest Bremer value (18-99). The varieties appeared more closely related to East Asian I. crenata and L mutchagara (73% of the conditions) than to other sSA taxa, excluding I. paraguariensis. The latter taxa appeared basal in relation to remaining ingroup taxa, in eight of the 15 parameters sets assayed, whereas under six analytical conditions, it was closely allied to L dumosa, I. crenata and L mutchagara.
Concerning North American taxa, the ITS sequences of L collina, I. decidua, and L vomitoria were related across all the analytical space, although with variable Bremer support (from 3 to 54), i.e., values ranging from the fourth to less than 10th highest. The pair of taxa I. cassine and I. opaca formed a clade with I. argentina AJ492655 in 100% of the conditions, with Bremer support values within the third and seventh highest (13-49). Ilex amelanchier and L mucronata appeared monophyletic in 80% of the conditions. The location of L glabra varied considerably under different analytical conditions. Likewise, the placement of Hawaiian L anomala relies on analytical conditions, appearing basal to several East Asian taxa in 53% of the conditions. The Asian taxa formed two stable groups in all of the analyses performed. One was composed of sequences of L beecheyii, L integra, L kiusiana, I. mertensii, and L percoriacea, and had variable Bremer support values ranging from 5 to 87 (second to ninth higher values). Instead, the other stable clade formed by I. cornuta, I. dimorphophylla, and L maximowicziana had lower Bremer support values. Whereas, Asian taxa I. buergeri, I. latifolia, I. liukiuensis, and I. warburgii appeared monophyletic in 93% of the conditions. The position of the remaining Asian taxa (i.e., I. macropoda, I. micrococca, I. rotunda, L serrata, and I. yunnanensis) varied with alignment parameters. ITS sequences of African L mitis and Macaronesian L canariensis formed a clade in 60% of the parameter space.
Combined analysis-Additional analyses were carried out combining AFLP data for eight taxa and nucleotide sequences of active ITS alleles. A 50% majority rule consensus tree is shown as a summary of the sensitivity analysis (Fig. 19, on the right). The lengths of strict consensus trees and single MPT             Figure 17. 1 e,,,,, 19. Fifty percent majority rule consensus tree derived from the sensitivity analyses of ITS sequences alone (on the left) and from the combined analysis (on the right). Numbers below branches indicate percentage of formation of the node. Numbers above branches represent ranges of Bremer support values. Clades found throughout the parameter space explored are marked with thick lines. Gbk indicates sSA sequences downloaded from the GenBank and discussed in the text. Loesener's subgeneric disposition is presented. obtained, ranged from 1356 to 7479. When compared with the ITS alone analysis, we observed that in 10/15 of the explored parameter space, the number of MPTs found remained unchanged; in 4/10 situations, the number of MPT was lower than in the ITS analysis alone; and only in one case, was the number of MPT higher (condition 411 rendered 10 MPT when combining and 4 in ITS alone). In general, the combined to-pology was congruent with topology of ITS sequence alone, as were Bremer support values for stable nodes. Main differences were found within the sSA clade, though some Asian and North American taxa showed alternative placement. Within sSA, the relationship between I. brasiliensis, I. integerrima, and I. theezans was resolved in the combined analysis as with AFLP data alone. In contrast, the combination of data did not resolve the relationship between both I. microdo (which appeared monophyletic in 47% of the condit I. pseudobuxus. The placement of I. argentina appear slightly less dependent on parameters. Similarly, th of AFLP data slightly stabilized the relationship of with I. crenata and I. mutchagara. Notwithstanding, combination did not help in stabilizing the paramet dent position of I. paraguariensis. Neither the analysis of ITS data alone nor the co analysis, yielded trees that completely agreed with taxonomic classification. Although the trees derived analysis under parameter sets 221 (Fig. 8), 441 ( Fig.  821 (Fig. 14) had the highest number of clades in co with current classification, none of the supraspecific peared monophyletic. DISCUSSION AFLP analyses-The AFLP technique is highly efficient for detecting polymorphisms of restriction fragments by PCR amplification. The difficulty in identifying homologous markers, the potential non-independence of bands and the dominant nature of AFLP markers, are the main disadvantages of the technique (Mueller and Wolfenbarger, 1999). However, AFLP allows the simultaneous screening of many different high-resolution markers randomly distributed throughout the genome.
The combinations of selective AFLP primers employed in this study allowed the recognition of several diagnostic bands within each of the seven species assayed. These would be useful for molecular identification of the plants, particularly at early developmental stages. A relationship between the number of AFLP bands scored for each species and its ploidy level was detected, even though it was not statistically tested ( Fig.  1). For instance, L argentina had the most complex banding pattern, and it is the sole tetraploid species (2n = 80) detected in the region (Barral et al., 1995).
Our AFLP analysis revealed that individuals belonging to the same morpho-species formed coherent clades (with moderate to strong bootstrap support). Therefore, seven clades were distinguished corresponding to each of the species surveyed. The species I. brasiliensis, I. integerrima, and I. theezans had a strongly supported relationship. The branching order in which I. argentina and I. brevicuspis associate with I. brasiliensis, I. integerrima, and I. theezans was not undoubtedly defined by the data; nevertheless, the monophyly of the five species was very well-supported. Instead, the relationship to I. dumosa and I. paraguariensis of the aforementioned species was not clearly established.
An important aspect of our analysis is the observation that AFLP data do not confirm the traditional distinction, established by Loesener (1901), between I. dumosa var. dumosa and I. dumosa var. guaranina, which was based on the inflorescence type. It seems possible that morphological features used to distinguish and define the varieties (Loesener, 1901;Giberti, 1994a) are not associated with the differences observed herein. Instead, our results suggest that within I. dumosa, Argentine and Paraguayan populations were derived from a common ancestral population a posteriori of the divergence from Brazilian population. A similar association was observed within I. brasiliensis. On the contrary, no geographical pattern was noted within I. paraguariensis, for which internal nodes appeared, in general, poorly supported and resolved. For this species, calculated mean intraspecific distance was among the lowest values obtained. Similarly, Gaue Cavalli-Molina (2000), using RAPD markers, found vergence among natural populations of L paraguariens Brazil. It seems reasonable that the arguments used by authors (the occurrence of gene flow favored by frui birds or recent fragmentation caused by deforestation reduction of natural forests) are responsible for the re both studies.
ITS sequence studies-The nucleotide divergence comparison revealed that in several cases, divergence between sequence data obtained in this study and sequence data downloaded from GenBank, far exceeded the values reported in the literature, usually below 5% (Baldwin et al., 1995). Results presented here, as a rule, do not support the view that the ITS sequences retrieved from GenBank (those of Manen et al., 2002), and our ITS sequences reflect intraspecific polymorphism. In general, the sSA ITS sequences of Manen et al. (2002) were highly divergent from those obtained herein, have lower G + C content, and were less thermodynamically stable. Previously, Manen et al. (2002) reported two of their South American ITS sequences as paralogous (L dumosa AJ492656 and L brevicuspis AJ492663). In this study, we gathered solid evidence indicating that those sequences might be treated as pseudogene rDNA alleles, or as highly divergent paralogues alleles (sensu Baldwin, 1995), or as class II paralogues (sensu Buckler et al., 1997) or deep paralogues (sensu Bailey et al., 2003). Furthermore, we report three additional ITS sequences from Manens' original data set as pseudogene rDNA alleles, i.e., I. brasiliensis AJ492661, L integerrima AJ492664, and L pseudobuxus AJ492660. In contrast, L microdonta AJ492665 and L theezans AJ492666 might represent alleles within different active ribosomal arrays. Our phylogenetic analysis including divergent ITS alleles, confirmed our results based on examining factors correlated to functionality. Different laboratory methods were employed in Manen's work and ours, concerning mainly the amplification primers (universal primer plus specific primer used herein, and both universal primers in Manen's) and the sequencing (direct sequencing of double stranded DNA was employed herein vs. cloning by Manen et al., 2002). These might account for the differential recovery of pseudogenes.
Our results support the idea that phylogeneticists using nuclear ITS sequences should be aware of the potential risks of downloading sequences from GenBank and using them for phylogenetic analysis without prior corroboration of their orthologous nature. As stated by Mayol and Rosell6 (2001), without a detailed inspection of some basic features, errors in evolutionary inferences could be easily overlooked. Incorporating ITS paralogues can distort the phylogenetic signal (Buckler et al., 1997) yielding spurious clades.
Phylogenetic inferences and sensitivity analyses--Direct optimization of DNA sequence data yields improved testing of phylogenetic hypotheses because provisional homologies are optimized for each cladogram individually, not a priori and universally as with the static homologies of multiple alignments (Wheeler, 2003). The sensitivity analysis as proposed by Wheeler (1995), results in the knowledge of nodal stability, showing how an alignment parameter affects the outcome tree topology (Giribet, 2003). Stability, logically different from support, can only be defined by reference to some factors (Goloboff et al., 2003). A group found in all the parameter space may nonetheless be poorly supported, and at the same time, a group with high "support" may be unstable to changes in some parameters (Goloboff et al., 2003). A drawback of sensitivity analysis, is that the number of parameters to be explored is unbounded, and thus new or more gap : change cost ratios could always be added (Giribet and Wheeler, 1999). Trying to generalize an evolutionary model for all taxa, all loci, and all regions may be simplistic; however, evaluating two general parameters (gap/change ratio and Tv/Ts) is a starting point (Maxmen et al., 2003).
In this group of plants, analytical conditions rendered different hypotheses of relationships between species, indicating that the alignment parameters did influence the phylogenetic inference. Notwithstanding, for the ITS sequence data alone and in combination with AFLP data, the sensitivity analysis showed that some clades were stable to the variation in analytical parameters and that those clades mainly agreed with geographical origin of the taxa involved. For instance, the clade composed solely by southern South American taxa (the sSA clade) namely, L argentina, I. brasiliensis, L brevicuspis, L integerrima, L microdonta, L pseudobuxus, L taubertiana, and L theezans was stable to parameter fluctuations (i.e., robust) in both analyses of ITS data alone and combined (Fig.   19). It is noteworthy that neither the sequences of L microdonta nor those of I. theezans appeared monophyletic, which could be caused by ancestral polymorphism, introgression, or sampling errors. Another stable clade consisted of East Asian taxa L beecheyi, L integra, L kiusiana, L mertensii, and L percoriacea. In addition, both stable clades were highly supported by the data in all conditions. On the other hand, East Asian clade formed by L cornuta, L dimorphophylla, and L maximowicziana and the clade concerning North American taxa I. collina, L decidua, and I. vomitoria were robust to fluctuations but showed variable support values. Many examples exist in the literature showing cases like the ones presented here, decoupling support from stability in both directions (Giribet, 2003).
Results presented here contradict those of Manen et al. (2002). A possible explanation for the incongruence of both studies would be based on the use of several highly divergent ITS sequences by Manen's team and on the analytical methodology employed. In this last respect, in the work of Manen et al. (2002), sequences were aligned manually, gaps were treated as missing data and indels were coded as present or absent, but without an explicit criterion for doing so. Manual alignments have been criticized for the lack of objectivity, lack of repeatability, and lack of criterion for exclusion data (e.g., Gatesy et al., 1993;Giribet and Wheeler, 1999). In our study, the analysis proceeds directly from sequence data to phylogenetic reconstruction treating gaps as relevant evolutionary transformations (Giribet and Wheeler, 1999).
Herein, majority rule consensus trees were presented because the lack of an appropriate external criterion (i.e., character congruence, as suggested by Wheeler, 1995 andWheeler, 1999) prevents us from choosing among competing phylogenetic hypotheses. Taxonomic congruence is probably not a good external criterion when groups are defined arbitrarily (Giribet and Wheeler, 1999).
The AFLP topology (Fig. 2) and the majority rule consensus cladograms of ITS data alone and combined data (Fig. 19) were remarkably congruent in that they demonstrated the existence of a group formed by L argentina, I. brasiliensis, I.
brevicuspis, L integerrima, and L theezans and that I. dumosa and I. paraguariensis were distantly related to the former species. These groupings were found to agree with the phenological data obtained by Coelho and Mariath (1996). Two main strategies of leaf longevity were distinguished by these authors, although some taxa were not included: subdeciduous (for I. brevicuspis, I. integerrima, I. theezans, I. microdonta, and L taubertiana) and perennifoliar (for I. dumosa and I. paraguariensis).
Apart from performing separate analyses of ITS sequence data and AFLPs, we carried out a combined phylogenetic analysis. We found that robust clades obtained with ITS data alone were not affected by the addition of AFLP data. Only sSA clade exhibited any changes, such as a better resolution of the relationship between I. brasiliensis, I. integerrima, and I theezans. Loesener (1901Loesener ( , 1908Loesener ( , 1942 and Baas (1975), based on shared morphological characters (for instance, inflorescences in fasciculate corymbs, flower and fruit of large sizes, coriaceous limbs with entire or pauci-aserrate leaf margins), proposed a close relationship among these taxa. Other clades showed in the combined analysis an increment in their percentage of formation throughout the analytical space. However, no major differences in support values were observed between sequence data alone and the combined analysis.
In spite of their taxonomic classification, I. pseudobuxus and I. taubertiana failed to form a clade, suggesting that morphological features shared by these taxa (the occurrence of long and slender flower peduncles; cfr. Edwin and Reitz, 1967) might have been developed by convergence. Similarly, I. microdonta and I. brevicuspis did not group in accordance with their taxonomic disposition, and thus morphological characters that distinguish these taxa from the rest could have arisen by convergence.
Furthermore, our analyses revealed that changes in alignment parameters resulted in alternative taxa being placed as sister to the sSA clade and that Asian taxa were involved in the majority of analytical conditions. The placement of I. paraguariensis was likewise highly sensitive to alignment parameter variation, appearing either related to the I. dumosa-I. crenata-I. mutchagara group, or as a third independent sSA lineage (Figs. 7,8,10,11,14,15,17,and 18). With regard to dumosa, our results conflict with the work of Manen et al. (2002) in which L dumosa appeared closely related to I. argentina, with both nuclear and plastid sequence data. We found no evidence of such a relationship, neither with AFLP data nor with ITS data. Instead, our results indicate a close relationship of I. argentina specimens to the remaining taxa of the sSA clade. It is worth mentioning that Manen's sequence of I. argentina AJ492655 appeared herein, robustly related and highly supported by the data to North American I. cassine and I. opaca. The fact that both ITS sequences of I. argentina failed to form a monophyletic clade could be explained by the coexistence within this tetraploid species of many types of active ITS copies due to an incomplete homogenization or to a recent formation of the polyploid. Alternatively, the possibility of mislabeled samples cannot be discarded. It would be interesting to reanalyze all the data set used by Manen et al. (2002), including plastid sequence data, through optimization and sensitivity analyses, in order to evaluate the stability of the phylogenetic hypotheses proposed. Present results support the idea that current classification of I. argentina, closely allied to I. paraguariensis (Baas, 1975), may not be justified and that a re-classification should be accomplished.
Our phylogenetic analysis provides no examples plete agreement between stable clades and the ta treatment of Loesener (1891Loesener ( , 1901Loesener ( , 1908Loesener ( , 1942. Oth such as those of Cu6noud et al. (2000) and Ma (2002), indicated an inconsistency of Loesener's infr system with sequence data. This inconsistency was a with Galle's system (not shown), which basically sum previous taxonomic works, but using contradictory The phylogenetic study of Setoguchi and Watanabe comprises solely Asian Ilex taxa followed the taxono ment of Yamazaki (1989). In our opinion, a discrepan between Yamazaki's treatment and sequence data pre those authors. Contradictions between morphologica lecular data are somehow expected because, as exp Systma (1990), reliance on overall morphological sim an approximation of phylogenetic relationships rathe shared derived character states can generate mislead mates.
Our results draw attention to the need of a thorough taxonomic revision of Ilex that would be based on a simultaneous cladistic analysis of all morphological and molecular characters available, without disregarding the influence of analytical parameters on phylogenetic reconstruction.