Molecular Biology and Evolution 18:223-234 (2001)
© 2001 Society for Molecular Biology and Evolution
ARTICLE |
Evaluating Hypotheses on the Origin and Evolution of the New Zealand Alpine Cicadas (Maoricicada) Using Multiple-Comparison Tests of Tree Topology
*Institute for Molecular Systematics, School of Biological Sciences, Victoria University of Wellington, Wellington, New Zealand;
Department of Biology, Duke University;
Department of Ecology and Evolutionary Biology, University of Connecticut; and
§Institute of Statistical Mathematics, Tokyo, Japan
 |
Abstract |
|---|
 TOP Abstract IntroductionMaterials and MethodsResultsDiscussionAcknowledgementsliterature cited |
|---|
The statistical testing of alternative phylogenetic trees is central to evaluating competing evolutionary hypotheses. Fleming proposed that the New Zealand cicada species Maoricicada iolanthe is the sister species to the major radiation of both low-altitude and montane Maoricicada species. However, using 1,520 bp of mitochondrial DNA sequence data from the cytochrome oxidase subunit I, tRNA aspartic acid, and the ATPase subunit 6 and 8 genes, we inferred that both M. iolanthe and another low-altitude species, Maoricicada campbelli, are nested within the montane Maoricicada radiation. Therefore, we examined the stability of the inferred phylogenetic placement of these two species using the newly developed Shimodaira-Hasegawa test (SH test) implemented in a maximum-likelihood framework. The SH test has two advantages over the more commonly used Kishino-Hasegawa (KH) and Templeton tests. First, the SH test simultaneously compares multiple topologies and corrects the corresponding P values to accommodate the multiplicity of testing. Second, the SH test is correct when applied to a posteriori hypotheses, unlike the KH test, because it readjusts the expectation of the null hypothesis (that two trees are not different) accordingly. The comparison of P values estimated under the assumptions of both the KH test and the SH test clearly demonstrate that the KH test has the potential to be misleading when the issue of comparing of a posteriori hypotheses is ignored and when multiple comparisons are not taken into account. The SH test, in combination with a variety of character-weighting schemes applied to our data, reveals a surprising amount of ambiguity in the phylogenetic placement of M. iolanthe and M. campbelli.
|
Introduction |
|---|
|
TOPAbstractIntroductionMaterials and MethodsResultsDiscussionAcknowledgementsliterature cited |
|---|
The origin and evolution of the New Zealand alpine biota has long intrigued biologists (e.g., Fleming 1963, 1979
; Raven 1973
; McGlone 1985
; Given and Gray 1986
; Morgan-Richards and Gibbs 1996
; Trewick, Wallis, and Morgan-Richards 2000
). It contains many endemic and highly specialized taxa despite the fact that the alpine habitat formed less than 5 MYA (Batt et al. 2000
; Chamberlain and Poage 2000
). One notable New Zealand alpine taxon is the endemic cicada genus Maoricicada (Dugdale 1972
) (Tibicinidae: Cicadettini). The genus is characteristically montane, although five species are found in low altitude habitats such as riverbeds (M. campbelli and M. hamiltoni), eroded surfaces (M. lindsayi), coastal rock fans (M. myersi), and low-altitude forest and scrub (M. iolanthe; fig. 1
). Fleming (1971)
and Dugdale and Fleming (1978)
published a variety of hypotheses to explain the origin and radiation of the genus. In particular, Fleming (1971)
suggested that the ancestral Maoricicada species, which he called proto-iolanthe, was similar to extant M. iolanthe. This supposition was based on the fact that present-day M. iolanthe are forest dwelling and restricted to the North Island. Fleming's (1971)
hypothesis makes the implicit prediction that M. iolanthe will be a sister lineage to the remaining Maoricicada species. This hypothesis is now readily testable using DNA sequence data coupled with newly developed statistical methodology.
|
Here, we implemented the Shimodaira-Hasegawa (SH; Shimodaira and Hasegawa 1999
) test of tree topology to statistically evaluate alternative phylogenetic hypotheses regarding the evolution of Maoricicada. Although the Kishino-Hasegawa (KH; Kishino and Hasegawa 1989
) and Templeton (1983)
tests are correct when two topologies selected a priori are tested, these tests are commonly used to compare a priori and a posteriori hypotheses. In this situation, both of these tests will calculate significance levels based on the incorrect assumption that two trees selected a priori are being compared. The SH test does two things to compensate for the incorrect application of other tree comparison tests: (1) it adjusts the expected difference in log likelihoods, and (2) it samples multiple alternative topologies. In an a posteriori test that compares the most likely tree to an alternative hypothesis, the expectation of the difference in log likelihoods is not 0 as it should be, even when all the trees are equally good in expectation. Rather, it will always be >0 because the most likely tree will always have the highest likelihood score, and subtracting any other likelihood from it will give a value higher than 0 (Shimodaira and Hasegawa 1999
; Goldman, Anderson, and Rodrigo 2001
). The SH test readjusts the expectation of the null hypothesis that the two trees are not different (a process known as "centering") and is correctly one-tailed. In addition, the SH test requires sampling of all reasonable alternative hypotheses so that the true topology is always available for comparison against the maximum-likelihood (ML) topology (Goldman, Anderson, and Rodrigo 2001
). The Templeton (1983)
test for the parsimony method may be similarly modified to correct the selection bias. We used 1,520 bp of mitochondrial DNA sequence data from the cytochrome oxidase subunit I (COI), ATPase subunits 6 (A6) and 8 (A8), and the transfer RNA aspartic acid (tRNAAsp) genes to reconstruct phylogenetic relationships among species of the genus Maoricicada. Using the SH test and various character weighting schemes, we critically evaluated competing hypotheses relating to the phylogenetic position of two key species; M. iolanthe and M. campbelli, in the Maoricicada radiation.
|
Materials and Methods |
|---|
|
TOPAbstractIntroductionMaterials and MethodsResultsDiscussionAcknowledgementsliterature cited |
|---|
Species Sampling
All described species and subspecies of Maoricicada were sampled, with the exception M. otagoensis maceweni (table 1 ). At least two individuals were sequenced for each of the 20 taxa to guard against problems arising from mistaken identity, contamination, or sample mix-up. Where multiple populations were available, individuals were analyzed from geographically separate areas to control for intraspecific variation and in order to detect possible cryptic species which may be present in Maoricicada (unpublished data). Those taxa that displayed high levels of geographic variation were represented by at least two populations in the phylogenetic analyses (M. campbelli, M. cassiope, and M. tenuis). We were unable to obtain live specimens of M. iolanthe, because this species has a high-frequency call that is inaudible to adult human ears and is therefore very difficult to collect. Extractions were made from two museum specimens of M. iolanthe collected from one population in 1971 (table 1 ). The New Zealand cicada species Rhodopsalta leptomera and Kikihia scutellaris were used as outgroups in the phylogenetic analyses. These two taxa were selected because of their close evolutionary relationship with Maoricicada based on molecular phylogenetic comparisons with a range of New Zealand, Australian, and New Caledonian species (unpublished data).
|
Laboratory Protocols
Cicadas were frozen on dry ice in the field or preserved in 95% ethanol and stored in an ultrafreezer. DNA was extracted from thorax muscle for males or from ovarian tissue for females using the "salting-out" method described by Sunnucks and Hales (1996)
. Briefly, tissue was homogenized in 20 µl of proteinase K solution (20 mg/ml) and incubated in 300 µl of TNES buffer (50 mM Tris buffer, 0.5% SDS, 20 mM EDTA, and 400 nM NaCl) at 55°C for 3 h. Subsequently, 85 µl of cold 5 M NaCl was added to the homogenate and vortexed. The homogenate was then centrifuged for 10 min at top speed, and the supernatant was removed. Centrifugation was repeated to ensure removal of all debris. An equal volume of cold 100% ethanol was added to the supernatant and then centrifuged for 30 min at top speed to pellet the DNA. The resulting pellet was washed with 70% ethanol, air-dried, and resuspended in 200 µl of ddH2O.
We used the polymerase chain reaction (PCR) to amplify an 819-bp target from the COI gene with the primers C1-J-2195 and TL2-N-3014 (Simon et al. 1994
). We also amplified a second region of 771 bp from the A6, A8, transfer RNA lysine (tRNALys), and transfer RNA aspartic acid (tRNAAsp) genes with primers TK-J-3799 and A6-N-4570 (this study). All primer sequences are given in table 2 . PCR reactions contained 1.0 µl of purified DNA solution, 1.5 µl MgCl2 (50 mM), 2.5 µl dNTPs (20 mM), 5.0 µl 10 x GibcoBRL Taq buffer, 2.5 µl of each primer (10 µM), 2.5 units GibcoBRL Taq polymerase (Life Technologies), and 34.5 µl ddH2O overlaid with 35 µl of mineral oil. Thermal cycling conditions using a Perkin Elmer Cetus DNA Thermal Cycler 480 were denaturation at 94°C for 45 s, annealing at 56–57°C for 45 s, and extension at 72°C for 75 s for 30 cycles. Negative controls containing no DNA were used in all sets of PCR reactions.
|
For one of the M. iolanthe museum specimens, we extracted DNA as above but included blank extractions with no tissue in order to detect any possible contaminants introduced during the extraction procedure. The DNA pellets were redissolved in 50 µl of ddH2O. For the second M. iolanthe specimen, we extracted total genomic DNA using a CHELEX protocol (Singer-Sam, Tanguay, and Riggs 1989
). Tissue was macerated in 20 µl of 5% CHELEX solution. The extract was heated to 95°C for 10 min and centrifuged at top speed for 10 min. One microliter of the supernatant was used as a template in the PCR reactions. For both M. iolanthe individuals, we were unable to obtain amplifications using either of the above primer combinations; thus, based on conserved Maoricicada mtDNA sequences, we designed a set of internal primers and amplified homologous mitochondrial gene regions from M. iolanthe in several partially overlapping fragments (table 2
). The following primer pairs were used: TK-J-3799 + A6-N-4102, A6-J-4097 + A6-N-4325, A6-J-4268 + A6-N-4521, A6-J-4470 + A6-N-4570, C1-J-2195 + C1-N-2476, C1-J-2430 + C1-N-2789, and C1-J-2621 + TL2-N-3014. The PCR conditions used for amplification of the overlapping gene targets were the same as those listed above, except 5 µl of BSA (10 mg/ml) was added to each reaction, and primer annealing was performed at 45°C for 35 cycles. To test for possible contamination in the DNA extraction and amplification process, PCR reactions containing the blank extraction were also included. Initial amplicons were gel-purified and reamplified using the same conditions, except no BSA was added and only 30 cycles were necessary to generate sufficient template for cycle sequencing.
PCR amplicons were purified by excision from a 1% low-melting-point agarose gel and centrifuged at top speed for 10 min, and the supernatant was used directly in cycle sequencing reactions. Alternatively, PCR products were purified using the PrepAGene (Biorad) kit following the manufacturer's instructions. Purified PCR products were cycle-sequenced using the Perkin Elmer Big Dye Terminator Cycle Sequencing Ready Reaction kit following the manufacturer's instructions. Cycle sequencing products were purified by ethanol precipitation and separated by electrophoresis on an ABI Prism 377 DNA sequencer. Sequences were checked for accuracy using the ABI sequence analysis software and were manually aligned using ESEE3.2 (Cabot and Beckenbach 1989
), facilitated by the conserved amino acid sequences and the absence of indels.
Patterns of Sequence Variation
Unless otherwise noted, all analyses were conducted using PAUP*, version 4.0b2a (Swofford 1998
). We calculated the numbers of varied sites, parsimony-informative sites, and base frequencies for each gene and codon position. The null hypothesis of base frequency stationarity among sequences was evaluated using the 
2 heterogeneity test as implemented in PAUP*, version 4.0b2a (Swofford 1998
). We examined all sites and parsimony sites only in order to assess the potentially confounding effects of unvaried sites, which, by definition, have stationary base frequencies (Waddell et al. 1999
). We note that the 2 heterogeneity test ignores the lack of independence among sequences due to a shared phylogenetic history; therefore, the results may be biased.
Phylogenetic Analyses
To detect the presence of significant heterogeneity in signal among data sets, we implemented the partition homogeneity test (PHT; Farris et al. 1994
) by dividing the characters according to two criteria. First, characters were partitioned into third positions versus all other sites combined. This was done in order to detect the presence of misleading signal in the third positions, which, as we have shown elsewhere (Buckley, Simon, and Chambers 2001
), have a higher substitution rate than the first, second, and tRNAAsp sites and are therefore more likely to suffer from homoplastic substitutions (e.g., Simon et al. 1994
). Second, we partitioned the data into A6 and A8 genes combined versus the COI gene. The tRNAAsp sites were excluded and the A8 sites were combined with the A6 sites, because the tRNAAsp gene contains only a small number of parsimony-informative sites, and both the A6 and the A8 genes are part of the same enzyme complex. Following Cunningham (1997)
, invariant sites were removed before the PHT was performed, under a heuristic search with 500 replications. The PHT was repeated under six-parameter parsimony (Williams and Fitch 1989, 1990
; see below).
Phylogenetic analyses were conducted under the maximum-parsimony (MP; Fitch 1971
) and ML (Felsenstein 1981
) optimality criteria. We followed this approach because for both of these two methods there is a theoretical expectation of the behavior of the method given various sets of relative branch lengths (relative rates of evolution among taxa). Thus, comparing topologies and nodal support under different optimality criteria can provide clues to confounding artifacts present in the data. For the MP analyses, we constructed trees under equal weights, six-parameter parsimony (Williams and Fitch 1989, 1990
), and downweighting of third-position sites. To obtain appropriate step-matrices, R-matrices (underlying relative substitution rates among the four nucleotides) were estimated under the general time reversible model (e.g., Yang 1994a
) using ML optimization and transformed by taking the natural log of each relative frequency parameter (Cunningham 1997
; Stanger-Hall and Cunningham 1998
). Because they violated the assumption of triangle inequality, some of the step matrices were adjusted using PAUP*, version 4.0. Heuristic searches were performed on an initial stepwise addition tree followed by tree bisection-reconnection (TBR) branch swapping. The third positions were also progressively downweighted relative to the first, second, and tRNAAsp sites by ratios ranging from 1:2 to 1:15. We used a range of weighting ratios, because the comparison of tree lengths across weighting schemes is not valid, thus making the selection of such values arbitrary to an extent (Sullivan, Markert, and Kilpatrick 1997
). We have shown elsewhere (Buckley, Simon, and Chambers 2001
) that the third codon positions in the Maoricicada mitochondrial genes examined here evolve at a faster rate, making them potential candidates for downweighting.
For the ML analyses, we estimated substitution model parameters using a best-fit model approach as described by Frati et al. (1997)
and Sullivan, Markert, and Kilpatrick (1997)
. Substitution models tested included those of Jukes and Cantor (1969
; JC69), Kimura (1980
; K80), Hasegawa, Kishino, and Yano (1985
; HKY85), and Yang (1994a
; GTR). These substitution models were also coupled with a series of parameters describing the distribution of among-site rate variation, including invariable sites (I sites; e.g., Hasegawa, Kishino, and Yano 1985
), gamma distributed rates (
rates; Yang 1994b
), a mixed model of gamma distributed rates and invariable sites (I+; Gu, Fu and Li 1995
), and a site-specific rates model (SSR; Swofford et al. 1996
; Buckley, Simon, and Chambers 2001
). For the SSR model, each codon position for each of the three protein-coding genes and all the tRNAAsp sites together were assigned their own relative substitution rate (SSR10 model).
To calculate genetic distances, we used the estimator described in equation (4) of Waddell and Steel (1997)
as implemented in PAUP*, version 4.0. This estimator produces what are usually referred to as ML distances (e.g., Swofford et al. 1996
), in which the relative-rate matrix is optimized using ML from all of the data and fixed for each pairwise comparison. We used this method because Waddell and Steel (1997)
found through simulation studies that homogenizing R for all distance comparisons tends to yield estimates with a lower variance than the more commonly used distance estimators. Nonparametric bootstrapping (Felsenstein 1985
) was performed with 500 replicates for the MP analyses and 100 replicates for the ML analyses. To detect increases in substitution rates among various Maoricicada lineages, we implemented the ML relative-rate test from HY-PHY, version 0.7b (Kosakovsky and Muse 2000
).
We also used marginal ancestral-state reconstruction under ML (Yang, Kumar, and Nei 1995
), implemented in PAUP*, version 4.0, to estimate the spatial location and substitution type of synapomorphies that unite the species M. iolanthe with M. campbelli. This was done in order to assess the nature of the synapomorphies supporting the grouping of these two species (see below).
The Shimodaira-Hasegawa Test
We used the particular implementation of the SH test referred to as the MS method by Shimodaira and Hasegawa (1999)
, which corresponds to the posNPncd test described by Goldman, Anderson, and Rodrigo (2001)
, except that our implementation is a weighted version. We emphasize to readers that there are other possible implementations of the SH test as described by Shimodaira and Hasegawa (1999)
. In this particular implementation of the SH test, the maximum of the standardized difference of log likelihoods is used as the test statistic, and the nonparametric bootstrap with reestimated log likelihoods (RELL) approximation (Kishino, Miyata, and Hasegawa 1990
) is used for resampling the log likelihoods. The RELL approximation is used to avoid reestimation of the parameters in the nonparametric bootstrap replicates. There is room to improve the SH test while controlling the coverage probability approximately, and the P values will be between those of the KH test and those of the SH test. We also note that it is possible to construct a parametric test to perform a similar a posteriori test using parametric bootstrap techniques (Goldman 1993
; Huelsenbeck and Rannala 1997
; Goldman, Anderson, and Rodrigo 2001
).
Using the SH approach, we examined the phylogenetic position of M. iolanthe and M. campbelli by taking the ML tree estimated under the GTR+I+ model and evaluating it against a series of topologies where the M. iolanthe and M. campbelli clade was shifted to a range of alternative phylogenetic positions. In constructing the alternative topologies, we avoided breaking up monophyletic groups that were consistently recovered in all phylogenetic analyses. In all, we tested 13 alternative phylogenetic topologies against the optimal ML GTR+I+ tree. The set of topologies that we selected were therefore a mixture of a priori and a posteriori hypotheses, exactly the situation in which the KH test is inappropriate. Sitewise log likelihoods were calculated under the GTR+SSR10 and GTR+I+ models for each topology tested using PAUP*, version 4.0, and exported into programs written by H.S. (under development) for calculation of SH test P values. Substitution model parameters were reoptimized for each topology in order to maximize the likelihood score of each tree and minimize the selection bias. For comparative purposes, we also included results from pairwise KH tests performed under the GTR+I+ model and pairwise Templeton (1983)
tests performed under six-parameter parsimony. The P values calculated from the KH and Templeton tests were halved in order to convert the test to a one-sided test, because one of the topologies being tested was known to be the optimal topology (Goldman, Anderson, and Rodrigo 2001
). The significance levels of both the KH and the Templeton tests were adjusted using a Bonferroni correction.
|
Results |
|---|
|
TOPAbstractIntroductionMaterials and MethodsResultsDiscussionAcknowledgementsliterature cited |
|---|
Patterns of Sequence Variation and Substitution Model Selection
The complete data set consisted of 1,520 sites for each of the 25 individuals sequenced. The target regions contained 64 tRNAAsp sites, 753 sites from the COI gene, 156 sites from the A8 gene, and 547 sites from the A6 gene. The number of varied and parsimony-informative sites within each of the four genes are given in table 3 . The sequences are available from GenBank under the accession numbers given in table 1 . Using the
2 heterogeneity tests, we were unable to detect any significant heterogeneity in base frequencies among taxa for the parsimony-informative sites (P = 0.628). The two best fitting models were the GTR+I+ and the GTR+SSR10 models. It is not valid to quantitatively test these two models against one another using the 2 approach, because they are not nested hypotheses.
|
The results of the PHT indicated that there was no significant difference in the number of steps required by the individual and combined-gene analyses under equal weights (P = 0.482) and six-parameter parsimony (P = 0.950). Similarly, we were unable to detect any significant incongruence between the COI gene and the A6 and A8 genes combined under equal weights (P = 0.980) or six-parameter parsimony (P = 0.9940).
Phylogenetic Relationships Among Maoricicada Species
All corrected genetic distances given in this section were estimated using ML under the GTR+I+ model and represent the expected number of substitutions per site. We do not report the GTR+SSR10 distances, because we have previously shown that this model underestimates branch lengths relative to gamma and invariable-sites models (Buckley, Simon, and Chambers 2001
). The following description of phylogenetic relationships among Maoricicada species begins with the inferred root of the tree and proceeds up the tree to the more derived taxa.
Both ML and MP supported the monophyletic status of the genus Maoricicada, with all bootstrap support values 
85% (figs. 2 and 3 ). The two outgroup species, K. scutellaris and R. leptomera, were differentiated from the ingroup taxa by corrected genetic distances ranging from 0.19 to 0.28. Starting at the base of the Maoricicada phylogeny, we observed a split between M. hamiltoni, M. myersi, and M. lindsayi on the one hand and the remaining Maoricicada radiation on the other hand. The three former species are all restricted to low-altitude habitats such as coastal rock fans and associated stream channels (M. myersi), clay banks (M. lindsayi), and riverbeds (M. hamiltoni). Among these three low-altitude species, M. myersi and M. lindsayi are sister species, with all methods giving 100% estimates of bootstrap support (figs. 2 and 3
).
|
The two remaining low-altitude species, M. iolanthe and M. campbelli, did not group with M. myersi, M. lindsayi and M. hamiltoni, but were instead nested within a clade of montane species. However, when we progressively downweighted the third positions, a radically different picture emerged regarding the phylogenetic placement of M. iolanthe, and, additionally, support for a monophyletic M. campbelli increased (fig. 4 ). When the third positions were downweighted by a factor of 1:11 or greater, the three most-parsimonious trees that were recovered all placed M. iolanthe as the basal Maoricicada species, although estimates of bootstrap support were low for many nodes (fig. 5 ). Although there is no objective means for selecting among different parsimony weighting schemes, the results shown in figures 4 and 5 illustrate the sensitivity of the parsimony analyses to the weight accorded to the third positions.
|
|
Under ML and MP with equal weights and six-parameter parsimony, the sister group to M. myersi, M. lindsayi, and M. hamiltoni was supported by bootstrap values ranging from less than 50% (fig. 3B ) to 76% (fig. 2 ). All species within this clade are montane, with the exception of M. iolanthe and M. campbelli, which are low-altitude.
|
Within this radiation, the species M. mangu had a poorly supported phylogenetic position. With the exception of the Maoricicada mangu mangu population from Awakino Ski Field, all populations and subspecies of M. mangu formed a monophyletic group, supported by bootstrap values ranging from less than 50% (fig. 3B ) to 96% (fig. 2 ). The M. mangu mangu individual from Awakino was differentiated from the Porters Pass M. mangu mangu by a corrected genetic distance of 0.06. The Awakino M. mangu mangu sequence grouped with Maoricicada alticola under most phylogenetic reconstruction methods, although this was poorly supported by the bootstrap analyses. The M. mangu mangu population from Porters Pass and Maoricicada mangu multicostata were monophyletic, with estimates of bootstrap support for this grouping ranging from 85% (fig. 3 ) to 92% (fig. 2 ). Based on the data presented here, we believe that the Awakino population of M. mangu mangu represents a new species of Maoricicada. Further morphological, behavioral, and molecular studies are in progress to test this hypothesis.
Most of the phylogenetic methods placed M. cassiope and M. tenuis as sister species; however, all estimates of bootstrap support were less than 50%, with the exception of the two ML trees (fig. 3 ). Relationships among M. cassiope, M. tenuis, and M. mangu were poorly resolved under both optimality criteria (figs. 2 and 3 ).
The Maoricicada species M. campbelli, M. iolanthe, M. alticola, M. otagoensis, M. nigra, M. clamitans, and M. oromelaena, and M. phaeoptera and the M. mangu mangu from Awakino formed a monophyletic clade in the ML, equal-weighted, and six-parameter MP analyses. However, bootstrap support for this clade was >50% in the ML analyses (fig. 3 ). Within this clade, most relationships were poorly resolved. However, we consistently observed a relationship between M. nigra and M. otagoensis, with bootstrap values ranging from 58% (fig. 2 ) to 65% (fig. 3A ). Another well-supported sister species relationship existed between M. oromelaena and M. clamitans (figs. 2 and 3 ). Estimates of bootstrap support for this grouping were all >95%. A consensus tree showing well-supported relationships among Maoricicada species with habitat preferences is given in figure 6 .
|
Examining Support for the Phylogenetic Position of M. iolanthe
In table 4 , we present data that clearly demonstrate that synapomorphies that support the grouping of M. iolanthe and M. campbelli together contain large amounts of homoplasy. Using marginal ancestral-state reconstruction under ML (Yang, Kumar, and Nei 1995
) on the ML GTR+I+ tree (fig. 3A
), we determined the spatial location and substitution types of these synapomorphies. We identified 10 sites that had changed at the node uniting all of the M. iolanthe and M. campbelli sequences and that were unvaried among these four sequences (table 4
). All of these substitutions were synonymous, and all but one was located at the third position. Eight of the 10 changes were transitions, and all of the sites had experienced at least two substitutions over the entire phylogeny and were thus convergent.
|
Because the predominantly low-altitude M. iolanthe and M. campbelli were nested within a clade of montane species and the codon-weighted MP analyses indicated that their position within the tree was somewhat ambiguous, we evaluated support for alternative phylogenetic scenarios using the SH test (table 5 ). When the SH test was implemented under the GTR+I+
and GTR+SSR10 models, we were unable to reject any of the alternative topologies regarding the phylogenetic placement of the M. iolanthe and M. campbelli clade, including a basal position in the tree. Using the Templeton tests, we were able to reject a basal position for M. iolanthe; however, we have already demonstrated that the parsimony analyses were highly sensitive to the weighting scheme imposed on the data. Genetic diversity within M. campbelli ranges from 0.021 to 0.056, and the species is differentiated from M. iolanthe by a relatively large genetic distance of 0.070–0.077. We used ML relative-rate tests to assess the significance of the apparent rate acceleration in the M. iolanthe lineage. Using Maoricicada oromaelana as an outgroup, M. iolanthe was revealed to be evolving at a significantly higher rate than the Otago (P = 0.033) and South Island (P = 0.002) M. campbelli haplotypes, but not the North Island haplotype (P = 0.072).
|
|
Discussion |
|---|
|
TOPAbstractIntroductionMaterials and MethodsResultsDiscussionAcknowledgementsliterature cited |
|---|
Our analyses do not agree with the suggestions of Fleming (1971)
, who hypothesized that the lowland, forest-dwelling species M. iolanthe represents the ancestral Maoricicada lineage and that lowland species in general arose first. Fleming's (1971)
hypothesis implies that M. iolanthe should be basal in any phylogenetic analysis. Yet, with the single exception of the codon-weighted MP analyses, the topologies that we estimated indicate that M. iolanthe is nested within M. campbelli and that these two taxa are closely related to alpine species rather than to the lowland species M. myersi, M. hamiltoni, and M. lindsayi. Although the M. iolanthe sequences were obtained from pinned specimens collected in 1971, our controls suggest that the sequences presented here are genuine and are not compromised by contamination. DNA extracted from two individuals at different times produced two sequences that differed at only one site, and blank extractions indicated no contaminating DNA. If we accept the topologies obtained from our phylogenetic analyses, we are left with the conclusion that alpine species have given rise to low-altitude lineages through an evolutionary reversal in habitat preference. The close relationship between M. iolanthe and M. campbelli was not entirely unexpected—Fleming (1971)
noted similarities in song structure and morphology of the genitalia—however, our analyses have shown that the placement of these two species is ambiguous.
The result of the SH tests indicates that the precise placement of the M. iolanthe and M. campbelli clade is not statistically well supported. This clade can be shifted to a range of alternative positions in the ML tree without a significant change in likelihood score. Goldman, Anderson, and Rodrigo (2001)
have noted that the results of the SH test are highly dependent on the total number of topologies made available for simultaneous comparison. Ideally, all possible a priori topologies should be included in this set to ensure that the true tree is available for testing. However, with a data set the size of ours, this is difficult, if not impossible. If we were to include more topologies, the test would simply become more conservative, and our ultimate biological conclusion would not be altered.
The data presented here clearly demonstrate that the overconfidence of the KH test can affect the general biological conclusions of an empirical study in the absence of corrections for multiple comparisons. For example, when this bias was ignored, we were able to incorrectly reject the hypothesis that M. iolanthe and M. campbelli are sister species to the remaining Maoricicada radiation. The pairwise KH and Templeton tests seemed to offer greater resolution than the SH test, because the P values calculated using the KH and Templeton tests were much lower; however, this result is misleading, since the KH and Templeton tests give overconfidence in a topology because sampling error due to estimation of the topology is ignored (unlike in the SH test, which explicitly accounts for this problem; Shimodaira and Hasegawa 1999
). The problem of multiple testing can be compensated for using a Bonferroni correction (see also Bar-Hen and Kishino [2000] for a similar approach). For example, many of the P values obtained from both the Templeton test and the KH test were rendered nonsignificant following the appropriate adjustment. Although the Bonferroni correction is statistically valid, it is far more conservative than the SH test. Therefore, when comparing many topologies, the SH test will always be preferable. Our observations regarding the relative power of the SH and KH tests agrees with the prediction of Goldman, Anderson, and Rodrigo (2001)
, who noted that if the adjusted (one-tailed) P value from the KH test indicates acceptance of the null hypothesis, then the SH test will imply the same result. However, this is not a reason to incorrectly apply the KH test, because, as Goldman, Anderson, and Rodrigo (2001)
point out, one cannot predict the result of the SH test if the KH test indicates rejection of the null hypothesis.
How can we explain the inferred derived phylogenetic placement of M. iolanthe and M. campbelli together apart from the other lowland species? There are two classes of explanations: (1) the tree is incorrect (consistent with the SH-test) because (a) M. iolanthe is a long branch, or (b) the substitution models that we used may not fit the data well; or (2) the tree is correct (the SH test is too conservative), but (a) due to lineage sorting, the gene tree does not match the species phylogeny, or (b) the derived position results from the fact that the montane lineages could have sequentially split off M. iolanthe and M. campbelli, both of which retained the ancestral Maoricicada phenotype, or (c) M. iolanthe and M. campbelli represent evolutionary reversals regarding habitat preferences.
Our analyses indicate that M. iolanthe forms one of the longest branches in the phylogeny (fig. 3
) and is evolving at a significantly higher rate than at least two of the M. campbelli lineages. It is well known that long branches can bias phylogenetic analysis due to the accumulation of multiple substitutions obscuring phylogenetic signal (Felsenstein 1978
; Hendy and Penny 1989
; Swofford et al. 1996
). We deliberately employed methods that are known to be less susceptible to the confounding effects of long branches, such as ML using a best-fit substitution model (Swofford et al. 1996
; Cunningham, Zhu, and Hillis 1998
). However, there are numerous artifacts of the evolutionary process that have the potential to mislead our analyses. If the model of evolution fits poorly, noisy data may cause serious problems even for ML. These observations indicate that the M. iolanthe sequence analyzed here may be in a window of nucleotide variation where sufficient noise exists at third positions to be misleading, yet not enough variation is present at more conserved first, second, and tRNAAsp sites to overwhelm this noise or to strongly support an alternative hypothesis (e.g., Halanych and Robinson 1999
). Despite the results of the PHT tests, there may not be enough variation at the more conserved (i.e., first, second, and tRNAAsp) sites to support a finding of any significant incongruence. If this is the case, then we note that even the most general of the commonly used substitution models (i.e., GTR+I+ and GTR+SSR10) are unable to overcome this misleading signal, as noted in other studies (e.g., Cao et al. 1998
). However, importantly, the complex models indicate that the difference between the optimal topology and the alternative topologies are nonsignificant. Sequence data from a nuclear locus, longer mtDNA sequences containing, presumably, more informative nonsynonymous replacements, or improvements in knowledge regarding the substitution process may ultimately be required to stabilize the phylogenetic position of M. iolanthe and M. campbelli.
If the topology is correct (and the SH test is too conservative), then we must invoke alternative explanations to account for the derived position of M. iolanthe and M. campbelli. First, retention of ancestral polymorphisms could cause the inferred gene tree to be a misleading estimate of the species phylogeny. Thus, sequence data from a nuclear gene could be used to evaluate this hypothesis. However, Moore (1995, 1997)
demonstrated that the probability that the mtDNA gene tree will track the population tree is much greater than that of a gene tree derived from a single nuclear locus. Hoelzer (1997)
has raised possible exceptions to this explanation that invoke gender-biased dispersal or extremely polygynous breeding systems. Although we cannot formally discount the possibility of gender-biased dispersal or polygyny within any species of Maoricicada, dispersal is a rare phenomenon in cicadas generally (Itô and Nagamine 1981
; Williams and Simon 1995
; deBoer and Duffels 1996
) and in the few instances when it does occur, it tends to involve female rather than male dispersal (Williams and Simon 1995
).
Second, it is possible that montane lineages could have sequentially split off M iolanthe and M. campbelli, both of which retained the ancestral Maoricicada phenotype. Evolutionary processes of this sort would produce the phylogenetic pattern that we observed whereby modern M. iolanthe and M. campbelli are nested within the Maoricicada radiation, yet still represent the ancestral phenotype. Although this hypothesis is possible, we believe that it is more likely that systematic error in the data is responsible for some of the observed phylogenetic patterns within Maoricicada. Further studies are now being initiated to test this hypothesis.
|
Acknowledgements |
|---|
|
TOPAbstractIntroductionMaterials and MethodsResultsDiscussionAcknowledgementsliterature cited |
|---|
The following people contributed to this study by aiding in fieldwork or providing specimens or information on collecting localities: P. Arensburger, B. Borger, S. Chiswell, K. Donohue, J. Dugdale, the late P. Fleming, F. Frati, G. Gibbs, D. Gleeson, D. Lane, M. Lark, and M. MacEwen. J. Dugdale kindly identified several of the specimens. W. M. Boon, D. Day, L. MacAvoy, L. Milicich, and P. Sunnucks provided advice and expertise in laboratory methods. N. Goldman, P. Lockhart, M. Hasegawa, M. Steel, and J. Sullivan provided advice on the phylogenetic analyses. C. Cunningham, G. Gibbs, D. Penny, T. Schultz, K. Crandall, and two anonymous reviewers provided comments on an earlier version of the manuscript. The New Zealand Department of Conservation (Te Papa Atawhai) provided collecting permits and advice on access to field sites. This work was supported by grants from the National Geographic Society (to C.S., George Gibbs, Maxwell Moulds, and G.K.C.), Victoria University of Wellington (T.R.B. and G.K.C.), the National Science Foundation (C.S.), the Fulbright Foundation (C.S., G.K.C., and George Gibbs), the Monbusho of the Government of Japan (H.S.), and the University of Connecticut (C.S.).
|
Footnotes |
|---|
Keith Crandall, Reviewing Editor
1 Keywords: multiple-comparison test
Shimodaira-Hasegawa test
Kishino-Hasegawa test
maximum likelihood
phylogenetics
Maoricicada 
2 Address for correspondence and reprints: Thomas Buckley, Department of Biology, Box 90338, Duke University, Durham, North Carolina 27708. E-mail: tbuckley{at}duke.edu
|
literature cited |
|---|
|
TOPAbstractIntroductionMaterials and MethodsResultsDiscussionAcknowledgementsliterature cited |
|---|
Bar-Hen, A., and H. Kishino. 2000. Comparing the likelihood functions of phylogenetic trees. Ann. I. Stat. Math. 52:43–56.
Batt, G. E., J. Braun, B. P. Kohn, and I. McDougall. 2000. Thermochronological analysis of the dynamics of the Southern Alps, New Zealand. Geol. Soc. Am. Bull. 112:250–266.
Buckley, T. R., C. Simon, and G. K. Chambers. 2001. Exploring among-site rate variation models in a maximum likelihood framework: the effects of model assumptions on estimates of topology, branch lengths and bootstrap support. Syst. Biol. (in press).
Cabot, E. L., and A. T. Beckenbach. 1989. Simultaneous editing of multiple nucleic acid and protein sequences with ESEE. Comput. Appl. Biosci. 5:233–234.
Cao, Y., P. J. Waddell, N. Okada, and M. Hasegawa. 1998. The complete mitochondrial DNA sequence of the shark Mustelus manazo: evaluating rooting contradictions with living bony vertebrates. Mol. Biol. Evol. 15:1637–1646.[Abstract]
Chamberlain, C. P., and M. A. Poage. 2000. Reconstructing the paleotopography of mountain belts from the isotopic composition of authigenic minerals. Geology 28:115–118.
Cunningham, C. W. 1997. Is congruence between data partitions a reliable predictor of phylogenetic accuracy? Empirically testing an iterative procedure for choosing among phylogenetic methods. Syst. Biol. 46:464–478.[Web of Science][Medline]
Cunningham, C. W., H. Zhu, and D. M. Hillis. 1998. Best-fit maximum likelihood models for phylogenetic inference: empirical tests with known phylogenies. Evolution 52:978–987.
DeBoer, A. J., and J. P. Duffels. 1996. Historical biogeography of the cicadas of Wallacea, New Guinea and the West Pacific: a geotectonic explanation. Paleogeogr. Paleoclimatol. Paleoecol. 124:153–177.
Dugdale, J. S. 1972. Genera of New Zealand cicadidae (Homoptera). N. Z. J. Sci. 14:856–882.
Dugdale, J. S., and C. A. Fleming. 1978. New Zealand cicadas of the genus Maoricicada (Homoptera: Tibicinidae). N. Z. J. Zool. 5:295–340.
Farris, J. S., M. Källersjö, A. G. Kluge, and C. Bult. 1994. Testing significance of incongruence. Cladistics 10:315–319.
Felsenstein, J. 1978. Cases in which parsimony or compatibility methods will be positively misleading. Syst. Zool. 27:401–410.[Web of Science]
———. 1981. Evolutionary trees from DNA sequences: a maximum likelihood approach. J. Mol. Evol. 17:368–376.[Web of Science][Medline]
———. 1985. Confidence limits on phylogenies: an approach using the bootstrap. Evolution 39:783–791.
Fitch, W. M. 1971. Toward defining the course of evolution: minimum change for a specific tree topology. Syst. Zool. 20:406–416.[Web of Science]
Fleming, C. A. 1963. Age of the alpine biota. Proc. N. Z. Ecol. Soc. 10:15–18.
———. 1971. A new species of cicada from rock fans in Southern Wellington, with a review of three species with similar songs and habitat. N. Z. J. Sci. 14:443–479.
———. 1979. The geological history of New Zealand and its life. Auckland University Press, Auckland, New Zealand.
Frati, F., C. Simon, J. Sullivan, and D. L. Swofford. 1997. Evolution of the mitochondrial cytochrome oxidase II gene in collembola. J. Mol. Evol. 44:145–158.[Web of Science][Medline]
Given, D. R., and M. Gray. 1986. Celmisia (Compositae-Astereae) in Australia and New Zealand. Pp. 451–470 in B. A. Barlow, ed. Flora and fauna of alpine Australasia—ages and origins. CSIRO, Australia.
Goldman, N. 1993. Statistical tests of models of DNA substitution. J. Mol. Evol. 36:182–198.[Web of Science][Medline]
Goldman, N., J. P. Anderson, and A. G. Rodrigo. 2001. Likelihood-based tests of topologies in phylogenetics. Syst. Biol. 49:1–19.
Gu, X., Y.-X. Fu, and W.-H. Li. 1995. Maximum likelihood estimation of the heterogeneity of substitution rate among nucleotide sites. Mol. Biol. Evol. 12:546–557.[Abstract]
Halanych, K. M., and T. J. Robinson. 1999. Multiple substitutions affect the phylogenetic utility of cytochrome b and 12S rRNA data: examining a rapid radiation in leporid (Lagomorpha) evolution. J. Mol. Evol. 48:369–379.[Web of Science][Medline]
Hasegawa, M., H. Kishino, and T. Yano. 1985. Dating of the human-ape split by a molecular clock by mitochondrial DNA. J. Mol. Evol. 22:160–174.[Web of Science][Medline]
Hendy, M. D., and D. Penny. 1989. A framework for the quantitative study of evolutionary trees. Syst. Zool. 38:297–309.[Web of Science]
Hoelzer, G. A. 1997. Inferring phylogenies from mtDNA variation: mitochondrial-gene trees versus nuclear-gene trees revisited. Evolution 51:622–626.
Huelsenbeck, J. P., and B. Rannala. 1997. Phylogenetic methods come of age: testing hypotheses in an evolutionary context. Science 276:227–232.
Itô, Y., and M. Nagamine. 1981. Why a cicada, Mogannia minuta Matsumura, became a pest of sugarcane: an hypothesis based on the theory of "escape." Ecol. Entomol. 6:273–283.
Jukes, T. H., and C. R. Cantor. 1969. Evolution of protein molecules. Pp. 21–123 in N. H. Munro, ed. Mammalian protein metabolism. Academic Press, New York.
Kimura, M. 1980. A simple method for estimating evolutionary rate of base substitutions through comparative studies of nucleotide substitutions. J. Mol. Evol. 16:111–120.[Web of Science][Medline]
Kishino, H., and M. Hasegawa. 1989. Evaluation of the maximum likelihood estimate of the evolutionary tree topologies from DNA sequences, and the branching order of Hominoidea. J. Mol. Evol. 29:170–179.[Web of Science][Medline]
Kishino, H., T. Miyata, and M. Hasegawa. 1990. Maximum likelihood inference of protein phylogeny and the origin of chloroplasts. J. Mol. Evol. 30:151–160.
Kosakovsky, S. L., and S. V. Muse. 2000. HY-PHY: hypothesis testing using phylogenies. Version 0.71b. North Carolina State University, Raleigh.
McGlone, M. S. 1985. Plant biogeography and the late Cenozoic history of New Zealand. N. Z. J. Bot. 23:723–749.
Moore, W. S. 1995. Inferring phylogenies from mtDNA variation: mitochondrial-gene trees versus nuclear gene trees. Evolution 49:718–726.
———. 1997. Mitochondrial-gene trees versus nuclear-gene trees, a reply to Hoelzer. Evolution 51:627–629.
Morgan-Richards, M., and G. W. Gibbs. 1996. Colour, allozyme and karyotype variation show little concordance in the New Zealand giant scree weta Deinacrida connectens (Orthoptera: Stenopelmatidae). Hereditas 125:265–276.
Raven, P. H. 1973. Evolution of alpine and subalpine plant groups in New Zealand. N. Z. J. Bot. 11:177–200.
Shimodaira, H., and M. Hasegawa. 1999. Multiple comparisons of log-likelihoods with applications to phylogenetic inference. Mol. Biol. Evol. 16:1114–1116.[Web of Science]
Simon, C., F. Frati, A. Beckenbach, B. Crespi, H. Liu, and P. Flook. 1994. Evolution, weighting and phylogenetic utility of mitochondrial gene sequences and a compilation of conserved polymerase chain reaction primers. Ann. Entomol. Soc. Am. 87:651–701.
Singer-Sam, J., R. C. Tanguay, and A. D. Riggs. 1989. Use of Chelex to improve the PCR signal from a small number of cells. Amplifications 3:11.
Stanger-Hall, K., and C. W. Cunningham. 1998. Support for a monophyletic Lemuriformes: overcoming incongruence between data partitions. Mol. Biol. Evol. 15:1572–1577.
Sullivan, J., J. A. Markert, and C. W. Kilpatrick. 1997. Phylogeography and molecular systematics of the Peromyscus aztecus species group (Rodentia: Muridae) inferred using parsimony and likelihood. Syst. Biol. 46:426–440.[Web of Science][Medline]
Sunnucks, P., and D. F. Hales. 1996. Numerous transposed sequences of mitochondrial cytochrome oxidase I-II in aphids of the genus Sitobion (Hemiptera: Aphididae). Mol. Biol. Evol. 13:510–524.[Abstract]
Swofford, D. L. 1998. PAUP*. Phylogenetic analysis using parsimony (*and other methods). Version 4. Sinauer, Sunderland, Mass.
Swofford, D. L., G. J. Olsen, P. J. Waddell, and D. M. Hillis. 1996. Phylogenetic inference. Pp. 407–514 in D. M. Hillis, C. Motitz, and B. K. Mable, eds. Molecular systematics. 2nd edition. Sinauer, Sunderland, Mass.
Templeton, A. R. 1983. Phylogenetic inference from restriction site endonuclease cleavage sites maps with particular reference to humans and apes. Evolution 37:221–244.
Trewick, S. A., G. P. Wallis, and M. Morgan-Richards. 2000. Phylogeographic pattern correlates with Pliocene mountain building in the alpine scree weta (Orthoptera, Anostostomatidae). Mol. Ecol. 9:657–666.[Medline]
Waddell, P. J., Y. Cao, J. Hauf, and M. Hasegawa. 1999. Using novel phylogenetic methods to evaluate mammalian mtDNA, including amino acid-invariant sites-LogDet plus site stripping, to detect internal conflicts in the data, with special reference to the positions of hedgehog, armadillo, and elephant. Syst. Biol. 48:31–53.[Web of Science][Medline]
Waddell, P. J., and M. A. Steel. 1997. General time-reversible distances with unequal rates across sites: mixing and inverse Gaussian distributions with invariant sites. Mol. Phylogenet. Evol. 8:398–414.[Web of Science][Medline]
Williams, K. S., and C. Simon. 1995. The ecology, behavior and evolution of periodical cicadas. Annu. Rev. Entomol. 40:269–295.[Web of Science]
Williams, P. L., and W. M. Fitch. 1989. Finding the weighted minimal change in a given tree. Pages 453–470 in B. Fernholme, K. Bremer, and H. Jornval, eds. Nobel symposium on the hierarchy of life. Elsevier, Cambridge.
———. 1990. Phylogeny determination using the dynamically weighted parsimony method. Methods Enzymol. 183:615–626.[Web of Science][Medline]
Yang, Z. 1994a. Estimating the pattern of nucleotide substitution. J. Mol. Evol. 39:105–111.
———. 1994b. Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods. J. Mol. Evol. 39:306–314.
Yang, Z., S. Kumar, and M. Nei. 1995. A new method of inference of ancestral nucleotide and amino acid sequences. Genetics 141:1641–1650.
CiteULike 
Connotea 
Del.icio.us What's this?
This article has been cited by other articles:
|
|
 |
|
  S.-T. Kim and M. J. Donoghue Incongruence between cpDNA and nrITS trees indicates extensive hybridization within Eupersicaria (Polygonaceae) Am. J. Botany, September 1, 2008; 95(9): 1122 - 1135. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
A. G. Murdock Phylogeny of marattioid ferns (Marattiaceae): inferring a root in the absence of a closely related outgroup Am. J. Botany, May 1, 2008; 95(5): 626 - 641. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 T. R. Buckley, M. Cordeiro, D. C. Marshall, and C. Simon Differentiating between Hypotheses of Lineage Sorting and Introgression in New Zealand Alpine Cicadas (Maoricicada Dugdale) Syst Biol, June 1, 2006; 55(3): 411 - 425. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
D. W. Weisrock, L. J. Harmon, and A. Larson Resolving Deep Phylogenetic Relationships in Salamanders: Analyses of Mitochondrial and Nuclear Genomic Data Syst Biol, October 1, 2005; 54(5): 758 - 777. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 P. Lemey, S. Van Dooren, and A.-M. Vandamme Evolutionary Dynamics of Human Retroviruses Investigated Through Full-Genome Scanning Mol. Biol. Evol., April 1, 2005; 22(4): 942 - 951. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
T. M. Townsend, A. Larson, E. Louis, and J. R. Macey Molecular Phylogenetics of Squamata: The Position of Snakes, Amphisbaenians, and Dibamids, and the Root of the Squamate Tree Syst Biol, October 1, 2004; 53(5): 735 - 757. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
A. L. Hipp, J. C. Hall, and K. J. Sytsma Congruence Versus Phylogenetic Accuracy: Revisiting the Incongruence Length Difference Test Syst Biol, February 1, 2004; 53(1): 81 - 89. [Full Text] [PDF]
|
|
|
|
|
|
L. Hufford, M. M. McMahon, A. M. Sherwood, G. Reeves, and M. W. Chase The major clades of Loasaceae: phylogenetic analysis using the plastid matK and trnL-trnF regions Am. J. Botany, August 1, 2003; 90(8): 1215 - 1228. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
S. Guindon and O. Gascuel Efficient Biased Estimation of Evolutionary Distances When Substitution Rates Vary Across Sites Mol. Biol. Evol., April 1, 2002; 19(4): 534 - 543. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||