HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) Supplementary Data Submit a response Alert me when this article is cited Alert me when eLetters are posted Alert me if a correction is posted Services Email this article to a friend Similar articles in this journal Similar articles in ISI Web of Science Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (4) Citing Articles via Google Scholar Google Scholar Articles by Gottlieb, A. M. Articles by Poggio, L. Search for Related Content PubMed Articles by Gottlieb, A. M. Articles by Poggio, L. Agricola Articles by Gottlieb, A. M. Articles by Poggio, L.

Molecular analyses of the genus Ilex (Aquifoliaceae) in southern South America, evidence from AFLP and ITS sequence data¹

¹Universidad de Buenos Aires, Facultad de Ciencias Exactas y Naturales, Ciudad Universitaria, Pabellón II, 4to piso, Departamento de Ecología, Genética y Evolución, Lab. 62, C1428EHA, Buenos Aires, Argentina; ²Instituto de Química y Metabolismo del Fármaco (IQUIMEFA), Junín 956, 1113-Buenos Aires, Argentina

Received for publication December 11, 2003. Accepted for publication September 17, 2004.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

 
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, I. integerrima, and I. theezans, and that I. dumosa and I. paraguariensis were distantly related to the former. Incongruence with traditional taxonomical treatments was found.

Key Words: AFLP • direct optimization • Ilex • ITS • sensitivity analysis

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

 
The plant genus Ilex L., which belongs to the holly family, Aquifoliaceae Bartl., comprises deciduous and evergreen bushes or trees with economic importance as crops and ornamentals. These plants are functionally dioecious; rudimentary pistils in staminate flowers and sterile stamens in pistillate flowers are borne on different plants. A remarkable morphological diversity has been described in the genus, including 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 ; Cuénoud 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, I. brasiliensis (Sprengel) Loes., I. brevicuspis Reissek, I. chamaedryfolia Reissek, I. dumosa Reissek, I. integerrima (Vell.) Reissek, I. microdonta Reissek, I. paraguariensis A. St. Hil., I. pseudobuxus Reissek, I. taubertiana Loes., and I. theezans C. Martius ex Reissek. Of those mentioned, I. paraguariensis is the most relevant from an economic perspective. In Argentina, southern Brazil, Paraguay, and Uruguay, the aerial parts of this plant are used to prepare a very popular beverage called mate, which is highly appreciated for its flavor and stimulating properties due to its caffeine and theobromine contents (Filip et al., 2001 ).

About a hundred years after Linnaeus (1753) published the diagnosis for the genus Ilex, Gray (1856) attempted to establish an infrageneric classification. In 1861, Reissek published the greatest species inventory for the time, followed by Hooker (1862) who reported the occurrence of about 145 species worldwide. Loesener (1891 , 1901 , 1908 , 1942 ) published a series of seminal studies recording almost 280 species of Ilex, dispersed globally, and recognizing several subgeneric ranks although under an unconventional hierarchical arrangement. Since then, several authors have recognized some inaccuracies in Loesener's system and proposed emendations in the systematics of the genus (Rehder, 1927 ; Hu, 1949 , 1950 ; Krüssmann, 1962 ; Baas, 1975 ; Lobreau-Callen, 1975 ; Martin, 1977 ). Species delimitation in Ilex, based on overall morphological similarity, is complicated due to the various concepts used by different authors to define the taxa and to the lack of complete information, particularly for evergreen species (Galle, 1997 ).

In order to investigate the relationships among the species of Ilex found in sSA, we implemented the amplified fragment length polymorphism (AFLP) technique (Vos et al., 1995 ). In addition, we studied the phylogenetic relationships of the species 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.

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 Roselló, 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 alignments. 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 , 2003 ). By doing this, indels appear as transformations (not as states) linking ancestral and descendent nucleotide sequences (Giribet and Wheeler, 1999 ).

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

 
Specimens studied
The list of sSA Ilex taxa employed, together with their accession numbers, geographical origin and GenBank accession numbers of ITS sequences obtained in this study are presented in Appendix 1 (see Supplemental Data accompanying the online version of this article). Plant material was principally obtained from the Estación Experimental INTA Cerro Azul (EEINTACA) Germplasm Bank (Misiones, Argentina). Leaf samples were also collected from natural populations from Northern Argentina; voucher specimens of these samples were deposited at the BACP Herbarium (Centro de Estudios Farmacológicos y Botánicos, CEFyBO, Buenos Aires, Argentina). Furthermore, leaf samples of commercial lines of mate were obtained from the Establecimiento Las Marías (Corrientes, Argentina). Young leaf samples were collected from healthy plants (without evidence of fungal or insect infection) and were silica-gel preserved. Additional 42 ITS sequences were downloaded from GenBank (http://www.ncbi.nlm.nih.gov; for their geographical origin and accession numbers see Appendix 2 (Data Supplement accompanying the online version of this article). Helwingia japonica (C. P. Thunberg ex A. Murray) Dietrich and Helwingia chinensis Batalin were chosen as the outgroup. According to several molecular studies (Cuénoud et al., 2000 ; Powell et al., 2000 ; Savolainen et al., 2000 ; Albach et al., 2001 ; Bremer et al., 2002 ), Helwingia and Ilex are closely related members of Aquifoliales. Loeseners' (1891 , 1901 , 1908 , 1942 ) original dispositions of subgeneric ranks were used in this work (Fig. 19). The species Ilex microdonta, I. pseudobuxus, and I. taubertiana, were only included in the sequence analyses because accessions of these taxa were unavailable in sufficient number for AFLP analysis. Only two sSA taxa remain to be studied, namely Ilex affinis and I. chamaedryfolia, because these have not been found in the field since the 1980s (G. Giberti, Instituto de Química y Metabolismo del Fármaco, IQUIMEFA, personal communication). Attempts to obtain good quality DNA from herbarium/air-dried specimens were unsuccessful.

View larger version (38K): [in this window] [in a new window]   Fig. 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

AFLP technique
The AFLP methodology was carried out using the AFLP Analysis System I Starter Primer kit (Gibco BRL, Life Technologies, Carlsbad, California, USA) as described in the instruction manual. Eight selective primers were combined as follows: E + AGG, M + CTT; E + ACT, M + CAC; E + AGC, M + CAT; and E + ACC, M + CAG. Other primers combinations assayed were discarded because they generated profiles in which the amplification products were too dense to allow reliable scoring, or alternatively, generated too few amplification products. Selective amplification products were mixed with an equal volume of dye reagent (98% [v/v] formamide, 10 mM EDTA, 0.025% [w/v] bromophenol blue and 0.025% [w/v] xylene cyanol). Five microlitres were separated by electrophoresis in a Model S2 apparatus (Gibco BRL Sequencing System, Life Technologies) through 6% (w/v) polyacrylamide gels containing 5 M urea, in 1 x TBE buffer (89 mM Tris, 89 mM boric acid, 2 mM EDTA, pH 8). A 30–330 bp AFLP DNA Ladder (Gibco BRL, Life Technologies) size marker was included twice or three times in each electrophoresis run. Thus, the size of AFLP bands scored ranged from 50 to 330 bp. Gels were stained with silver nitrate (Bassam et al., 1991 ).

ITS amplification
Amplification of ITS 1, 5.8S gene, and ITS 2 was carried out using primers ILEXFP (5'-AACAAGGTTTCCGTAGGTGA-3', Powell et al., 2000 ) and ITS-4 (White et al., 1990 ). Amplifications were performed in 50 µL using 0.225 mM of each dNTP (Promega Corp., Madison, Wisconsin, USA), 20–25 pmoles of each primer, 5 µL of Thermophilic DNA Polymerase 10 x Buffer (Promega Corp.), 1.5 mM MgCl₂, 1.25 units of Taq DNA polymerase (Promega Corp.), 1 µL of genomic DNA template, and sterile double-distilled water. Conditions for PCR amplification were those of Vilgalys and Hester (1990) . Amplifications were carried out in an Eppendorf Mastercycler (Perkin-Elmer Corp., Foster City, California, USA). Control samples without DNA template were included in each PCR run. PCR products were checked by electrophoresis of 1-µL samples in 1% agarose gels in 1 x TAE buffer (0.04 M Tris, 0.114% glacial acetic acid, 1 mM EDTA, pH 8.0). Quantification of PCR products was estimated by comparison with known amounts of a DNA molecular-size marker (Lambda EcoRI/HindIII, Promega Corp.) included in duplicate in all gels. Gels were stained with ethidium bromide for 15–20 min and photographed under UV light. In some cases, best amplification results were achieved by adding 6% bovine serum albumin (BSA, Promega Corp.) to the PCR reaction mix. PCR products were purified using a QIAGEN Gel Extraction kit (QIAGEN Inc.). Both strands of the complete ITS region (ITS 1–5.8S-ITS 2) were sequenced in Perkin Elmer/ Applied Biosystems automatic sequencing units (ABI Prism 310 or 377) at the Molecular Biology Facility of the Universidad de Alcalá (Madrid, Spain). Electropherograms were proofread with the software BioEdit Sequence Alignment Editor (Hall, 1999 ). Boundaries of the coding and spacer regions were determined by comparison with published sequences of Ilex. In general, when two or more orthologous ITS sequences were available per sSA taxon (see Appendix 1; Supplemental Data accompanying the online version of this article), these were used to produce consensus sequences in BioEdit (Hall, 1999 ) to speed-up computational times of phylogenetic analyses.

Data analyses
AFLP
Air-dried gels were scanned and banding patterns were registered with the Gel-Pro 4 Analyzer Program (Media Cybernetics, Maryland, USA). Each AFLP band, regardless of its relative intensity, was considered as a dominant allele at a unique locus. In some cases, fragments were scored as missing data because character states could not be interpreted unambiguously. The data were extracted as a table as either present (1) or absent (0). Monomorphic bands (bands present in all individuals of a species) and diagnostic bands (monomorphic bands exclusively present in one species) were discriminated within each species and across the entire data set. Box plots of mean values; SE and SD were calculated and plotted with the program STATISTICA '99 edition (Kernel release 5.5 A, StatSoft, 1999 ). For these calculations, accessions with missing data in one or more primer combinations were excluded. Pairwise genetic distance (D) between all individuals was estimated according to the complementary value of Nei and Li's (1979) similarity coefficient, implemented in PAUP* (Swofford, 1998 ): D = 1 – (2n_ij/ [n_i + n_j]), where n_ij is the number of shared fragments between two individuals i and j, and n_i and n_j are the total number of fragments in individual i and j, respectively. From these values, mean intra- and interspecific distance values were obtained. Furthermore, the AFLP data set for 120 accessions belonging to eight taxa was analyzed using maximum parsimony in PAUP*. Initial trees in the heuristic search by tree-bisection-reconnection (TBR) branch-swapping algorithm were built stepwise with 500 random addition sequence replicates. Characters were equally weighted and unordered, and all other settings were defaults. The reliability of the clades was evaluated by bootstrap analysis (Felsenstein, 1985 ) in PAUP* with 1000 replications and all other settings as described earlier. Trees were inspected with the program TreeView (Page, 1998 ). The unrooted resulting strict consensus tree was rooted using I. paraguariensis to be consistent with ITS trees.

ITS sequence studies
We characterized and identified ITS paralogues in sSA Ilex species by examining their G + C content, secondary structure stabilities, and nucleotide substitution patterns and using the tree-based approach advocated by Bailey et al. (2003) . To discriminate the presence of presumed active sequences from paralogues and putative pseudogenes, we first followed the procedure proposed by Mayol and Roselló (2001) . Therefore, pairwise nucleotide divergences were calculated between the 21 ITS sequences obtained herein and the 10 sequences retrieved from GenBank that belong to the same sSA Ilex taxa. The values were calculated using the Kimura two-parameter model (Kimura, 1980 ) with DNAdist 3.5c program implemented in BioEdit (Hall, 1999 ). Additionally, ITS sequences were screened using BioEdit (Hall, 1999 ) for length variation, G + C content, and the presence of a structural motif in ITS 1 (i.e., GGCRY(4 to 7 n)GYGYCAAGGAA; Liu and Shardl, 1994 ). This conserved motif is often associated with part of the hairpin and loop structures that may serve as a critical signal in the enzymatic processing of the ribosomal RNA (Liu and Shardl, 1994 ). Furthermore, optimal and suboptimal folding structures, and associated free energy values (G) in Kcal/mol, of ITS 1 and ITS 2 were recorded. Fold predictions were made at the Quickfold web server (Zuker, 2003 ) of the MFOLD program, version 3.1. In all cases, foldings were conducted at 37°C by use of a search within 5% of the thermodynamic optimality set. Because each sequence could adopt more than one folding structure, maximum and minimum G values were recorded for each sequence, together with the number of all possible (optimal and suboptimal foldings) secondary structures. However, no attempts were made to study the RNA structures per se. These screenings were extended to the remaining representatives of Ilex and outgroup.

Phylogenetic and sensitivity analyses
Parsimony was chosen as the optimality criterion, as implemented in the program for Direct Optimization POY, version 3.0.11 (Gladstein and Wheeler, 1997 ; Wheeler and Gladstein, 2000 ; De Laet and Wheeler, 2003 ). To save computation time, as suggested in the POY user's manual and in Giribet (2001) , the sequence data set was partitioned. The first partition comprises 257 to 265 bp, and the second partition has ca. 340 bp. Approximately 45 positions of the 5.8S gene, representing missing data in several sequences downloaded from the GenBank, were excluded from the analyses. A parameter space of two analytical variables was examined, i.e., gap cost and transversions/transitions (Tv/Ts) ratio cost. Precisely, we used gap costs = 1, 2, 4, 8 and 10, and Tv/Ts costs = 1 (equal weights), 2 (transversions receive twice the weight than transitions) and 4. All combinations of transformation costs used were in accordance with the triangle inequality (Wheeler, 1993 , 1995 ). Each particular tree search was performed as follows: 50 independent heuristic searches starting from the best of five independent Wagner trees were conducted via subtree-pruning and regrafting (SPR) followed by TBR, retaining two trees per step, as specified in the following commands: "–norandomizeoutgroup –buildsperreplicate 5 –replicates 50 –maxtrees 2 –slop 2 –checkslop 5". Commands "–slop" and "–checkslop" were employed to reduce error derived from the heuristic operations (Wheeler, 2002 ). Step matrices (i.e., transformation matrices) were specified in POY using the command "–molecularmatrix" with an argument for the given step matrix. Strict consensus trees were calculated using the command: "–poystrictconsensustreefile". An implied alignment was generated in each case by the command: "–impliedalignment". Approximate Bremer support values (Bremer, 1988 ) were calculated with a heuristic approximation procedure implemented in POY (command: "–bremer"). Finally, strict consensus trees and single trees were used to obtain 50% majority rule consensus tree. Implied alignments and trees were deposited at TreeBase (S996). We specifically used nodal stability as defined by Giribet (2003) , i.e., the degree to which nodes are affected by variation in the analytical conditions. Robustness was used herein as a measure of stability of nodes to fluctuation in parameters values across an analytical landscape (Ogden and Whiting, 2003 ).

Two sets of analyses were performed. First, we undertook a phylogenetic analysis of all ITS sequences, including those tentatively determined as paralogues and pseudogenes. To explore the influence of alignment parameters on the relative placement of these sequences, five different combinations of the two variables were employed. In particular, these combinations were: 111 (i.e., gap cost = 1, Tv = 1 and Ts = 1), 241, 421, 811 and 1021. Then, we performed a broader sensitivity analysis using only presumably active ITS sequences under 15 analytical conditions (namely, 111, 121, 141, 211, 221, 241, 411, 421, 441, 811, 821, 841, 1011, 1021 and 1041). Each of the 15 tree searches was performed as indicated earlier.

Combined analysis
The combined analysis involved AFLP data and ITS sequences. The sensitivity analysis was performed in POY as indicated previously, with the exception that a third file consisting of AFLP data per species was included in all calculations. The new data matrix for species was extracted from the original AFLP matrix for individuals and consisted of eight taxa and 578 characters. A band was considered present in a taxon (1) if it was recorded in more than 80% of the individuals, whereas it was considered absent in a taxon (0) if it was recorded in less than 20% of the individuals. Otherwise, bands were considered as polymorphic (0,1). The weights of AFLP changes were set so that nucleotide characters were never upweighted relative to AFLP data. The computer program POY allows nonsequence character changes to be simultaneously optimized with molecular data on multiple trees.

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

 
AFLP analysis
The four primer combinations generated 836 scoreable bands, of which 22% were monomorphic (Table 1). All 120 individuals tested shared a single constant band. Primer combinations E + ACT/M + CCA and E + ACC/M + CAG yielded a nearly identical number of monomorphic and diagnostic bands. Among the seven species included in this study, 94 diagnostic bands were distinguished, which accounted for more than 50% of the monomorphic fragments recorded. The proportion of monomorphic and diagnostic bands detected varied among the species (Table 2). Ilex integerrima had the highest number of monomorphic bands, 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).

View larger version (26K): [in this window] [in a new window]   Fig. 1. Comparison of mean number of AFLP bands obtained for selected Ilex species.

View larger version (22K): [in this window] [in a new window]   Fig. 2. Maximum parsimony strict consensus tree of selected Ilex species derived from AFLP data. Numbers above branches indicate bootstrap values (1000 replicates). Accessions, as presented in Appendix 1 (see Supplemental Data accompanying the online version of this article), are indicated between brackets for clearness. F, stands for female; M, stands for male. The arrow indicates terminals belonging to I. dumosa var. guaranina. *, indicates specimens used for sequencing

View larger version (20K): [in this window] [in a new window]   Fig. 3. Fifty percent majority rule consensus tree derived from the sensitivity analysis of Ilex ITS sequences, including divergent copies. 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

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, I. brevicuspis, I. integerrima, I. microdonta, I. pseudobuxus, I. taubertiana, and I. theezans, including sequences I. microdonta AJ492665 and I. 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 (I. 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 (I. 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 + I. 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 I. 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 I. mutchagara, appeared as sister group to the sSA clade, with a Bremer support of 29 (seventh highest value). In the analytical space explored, I. 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 I. 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 I. dumosa, I. crenata and I. mutchagara.

View larger version (45K): [in this window] [in a new window]   Figs. 4–7. Results from the sensitivity analyses of Ilex taxa via direct optimization of ITS sequences. 4. Strict consensus tree of four most parsimonious trees (MPT) of length 1013, under parameters 111. 5. Strict consensus tree of eight MPT of length 1174, under parameters 121. 6. Strict consensus tree of 11 MPT of length 1297, under parameters 141. 7. Strict consensus tree of three MPT of length 1301, under parameters 211. Figure abbreviations: The transformation matrix stands for the step matrix used in each search (i.e., 241: gap cost = 2, Tv cost = 4, Ts cost = 1). Gbk, indicates sSA sequences retrieved from GenBank and discussed in the text

View larger version (44K): [in this window] [in a new window]   Figs. 12–15. Results from the sensitivity analyses of Ilex taxa via direct optimization of ITS sequences. 12. Strict consensus tree of two MPT of length 2540, under parameters 441. 13. Single MPT of length 2655, obtained under parameter 811. 14. Single MPT of length 3047, obtained under parameter 821. 15. Strict consensus tree of three MPT of length 3629, under parameters 841

The Asian taxa formed two stable groups in all of the analyses performed. One was composed of sequences of I. beecheyii, I. integra, I. kiusiana, I. mertensii, and I. 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 I. 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, I. serrata, and I. yunnanensis) varied with alignment parameters. ITS sequences of African I. mitis and Macaronesian I. 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 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 topology 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. microdonta alleles (which appeared monophyletic in 47% of the conditions) and I. pseudobuxus. The placement of I. argentina appeared here slightly less dependent on parameters. Similarly, the addition of AFLP data slightly stabilized the relationship of I. dumosa with I. crenata and I. mutchagara. Notwithstanding, the data combination did not help in stabilizing the parameter-dependent position of I. paraguariensis.

Neither the analysis of ITS data alone nor the combined analysis, yielded trees that completely agreed with Loesener's taxonomic classification. Although the trees derived from ITS analysis under parameter sets 221 (Fig. 8), 441 (Fig. 12), and 821 (Fig. 14) had the highest number of clades in concordance with current classification, none of the supraspecific ranks appeared monophyletic.

View larger version (45K): [in this window] [in a new window]   Figs. 8–11. Results from the sensitivity analyses of Ilex taxa via direct optimization of ITS sequences. 8. Single MPT of length 1533, obtained under parameter 221. 9. Strict consensus tree of two MPT of length 1842, under parameters 241. 10. Strict consensus tree of four MPT of length 1797, under parameters 411. 11. Strict consensus tree of six MPT of length 2083, under parameters 421

	DISCUSSION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

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, I. 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, Gauer and Cavalli-Molina (2000) , using RAPD markers, found low divergence among natural populations of I. paraguariensis from Brazil. It seems reasonable that the arguments used by those authors (the occurrence of gene flow favored by fruit-eating birds or recent fragmentation caused by deforestation and/or reduction of natural forests) are responsible for the results of 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 (I. dumosa AJ492656 and I. 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, I. integerrima AJ492664, and I. pseudobuxus AJ492660. In contrast, I. microdonta AJ492665 and I. 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 Roselló (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, I. argentina, I. brasiliensis, I. brevicuspis, I. integerrima, I. microdonta, I. pseudobuxus, I. taubertiana, and I. 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 I. 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 I. beecheyi, I. integra, I. kiusiana, I. mertensii, and I. percoriacea. In addition, both stable clades were highly supported by the data in all conditions. On the other hand, East Asian clade formed by I. cornuta, I. dimorphophylla, and I. maximowicziana and the clade concerning North American taxa I. collina, I. 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 and Giribet and Wheeler, 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 I. argentina, I. brasiliensis, I. brevicuspis, I. integerrima, and I. 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 I. 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 (1901 , 1908 , 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 I. dumosa, our results conflict with the work of Manen et al. (2002) in which I. 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.