Phylogeny of the Avian Family Ciconiidae (Storks) Based on CytochromebSequences and DNA–DNA Hybridization Distances

MOLECULAR PHYLOGENETICS AND EVOLUTION

Vol. 8, No. 3, December, pp. 275–300, 1997 ARTICLE NO. FY970431

Phylogeny of the Avian Family Ciconiidae (Storks) Based on Cytochrome b Sequences and DNA–DNA Hybridization Distances Beth Slikas1 Academy of Natural Sciences, 1900 Benjamin Franklin Parkway, Philadelphia, Pennsylvania 19103-1195; and Biology Department, University of Pennsylvania, Leidy Labs, Philadelphia, Pennsylvania 19104-6018 Received October 18, 1996; revised April 28, 1997

This study is a phylogenetic analysis of the avian family Ciconiidae, the storks, based on two molecular data sets: 1065 base pairs of sequence from the mitochondrial cytochrome b gene and a complete matrix of single-copy nuclear DNA–DNA hybridization distances. Sixteen of the nineteen stork species were included in the cytochrome b data matrix, and fifteen in the DNA– DNA hybridization matrix. Both matrices included outgroups from the families Cathartidae (New World vultures) and Threskiornithidae (ibises, spoonbills). Optimal trees based on the two data sets were congruent in those nodes with strong bootstrap support. In the best-fit tree based on DNA–DNA hybridization distances, nodes defining relationships among very recently diverged species had low bootstrap support, while nodes defining more distant relationships had strong bootstrap support. In the optimal trees based on the sequence data, nodes defining relationships among recently diverged species had strong bootstrap support, while nodes defining basal relationships in the family had weak support and were incongruent among analyses. A combinable-component consensus of the best-fit DNA–DNA hybridization tree and a consensus tree based on different analyses of the cytochrome b sequences provide the best estimate of relationships among stork species based on the two data sets. r 1997 Academic Press

INTRODUCTION The avian family Ciconiidae, the storks, includes 19 extant species distributed primarily in tropical and subtropical regions around the world. Recognized as a family since 1901, the storks form a well-defined group characterized by several behavioral and morphological features: long legs with half the tibiotarsus bare, relatively short toes with small webs, 12 rectrices, 12

1 Current address: Molecular Genetics Laboratory, National Zoological Park, Smithsonian Institution, 3001 Connecticut Avenue NW, Washington, DC 20008.

primaries, bare portions on the head, feathered oil gland, common pectoral musculature associated with soaring flight, stout bill, predominantly black and white plumage, young with two down coats, air sacs under the neck skin, lack of powder down, similar stereotypical social behaviors, and urohydrosis (Kahl, 1963; Vanden Berge, 1970; Sibley and Ahlquist, 1990; Hancock et al., 1992). Relationships within the family have been assessed based on a review of morphological features and behavior (Verheyen, 1959), an extensive and detailed documentation and analysis of stereotypical social behaviors (Kahl, 1971a, 1972d), and a phenetic analysis of osteological measurements and Kahl’s behavioral data (Wood, 1983, 1984). These previous studies produced classifications (Verheyen, 1959; Kahl, 1972d) and phenetic clustering diagrams (Wood, 1983, 1984), but no phylogeny. Because storks are large and conspicuous, and many species thrive in zoos, these birds have been the focus of numerous ecological, behavioral, and morphological studies (Sibley and Ahlquist, 1990; Coulter et al., 1991). In particular, the behavioral data collected by M. Philip Kahl (1966, 1971b, 1972a–c,e, 1973) are unique and impressive. Detailed data were collected for all species, most from field observations, and the behaviors are described in a consistent fashion, facilitating comparisons among species. A reliable estimate of phylogeny is needed to complement the abundance of neontological data and to provide insight into the evolution of behavior, life-history traits, and morphological features. In addition, storks have a good fossil record (Howard, 1942; Olson, 1985) which, combined with a phylogeny, permits testing of hypotheses concerning the historical biogeography of the family. The goal of this study, therefore, was to produce a reliable estimate of phylogenetic relationships among the extant stork species. In general, congruence among trees based on different data sets is the best measure of the accuracy of a phylogenetic estimate (Penny and Hendy, 1986; Swofford, 1991; Sheldon and Bledsoe, 1993; Miyamoto and Fitch, 1995). For this reason, two molecular data sets were collected: single-copy nuclear DNA–DNA hybrid-

275

1055-7903/97 $25.00 Copyright r 1997 by Academic Press All rights of reproduction in any form reserved.

276

BETH SLIKAS

TABLE 1 List of Study Species, DNA Preparations, and Sources of Blood/Tissue

Species

DNA preparation no.

Source/source ID number

Species

Source/source ID number

90* 121 117

SD WAP AO893249 SD WAP AO893249 Miami MetroZoo AO3047

C. maguari Maguari stork

18 19* 27 42 97

Audubon Zoo 141L (1193) Audubon Zoo 140 (1192) Audubon Zoo 140 (1192) ANSP 2077 Audubon Zoo 140 (1192)

C. nigra Black stork

16*

Fort Wayne Children’s Zoo 4209 Fort Wayne Children’s Zoo 2097 Fort Wayne Children’s Zoo 2098 Fort Wayne Children’s Zoo 4210

80 101

Fort Wayne Children’s Zoo 90071 Fort Wayne Children’s Zoo 3282 Fort Wayne Children’s Zoo 1878 Fort Wayne Children’s Zoo 1985 Fort Wayne Children’s Zoo 1453 Fort Wayne Children’s Zoo 93001 Miami MetroZoo 10096 SD WAP AO998320-168

Ciconia stormi Storm’s stork

115

JAW Kirsch/U. Wisconsin 2834

Ephippiorhynchus asiaticus Blacknecked stork

106 112 113*

Mus. of Victoria, Australia Miami MetroZoo 611 Miami MetroZoo 612

Mycteria americana Wood stork

29* 124 71 104

ANSP 2070 ANSP 2070 LSU 5670 LSU 5670

E. senegalis Saddlebill stork

22 26* 94 95 114

Audubon Zoo 1618 Audubon Zoo 225R Audubon Zoo 1617 Audubon Zoo 1618 Miami MetroZoo AO3706

Audubon Zoo 1662 Audubon Zoo 1661 Audubon Zoo 1661

Jabiru mycteria Jabiru stork

47* 60 61 68 103

Audubon Zoo ‘‘Henry’’ LSU B10933 LSU B13373 LSU B13373 LSU B13373

Audubon Zoo 426 Audubon Zoo 429 LSU B10200 Audubon Zoo B95R, APZ63L LSU B10200

Leptoptilos crumeniferus Marabou stork

15

Ciconia abdimii Abdim’s stork

C. boyciana Oriental white stork

C. ciconia White stork

7* 8 46 99 111 75

2* 3 4 54 56 57

M. cinerea Milky stork

M. ibis Yellowbilled stork

M. leucocephala Painted stork

Anastomus lamelligerus African openbill stork Plegadis falcinellus Glossy ibis

24* 25 96

9* 17 58 100 105

78 79* 86 88 89 102 2687* 62* 120 127

Audubon Zoo 1575R Audubon Zoo 1582L Audubon Zoo 1622 Audubon Zoo 1515 Audubon Zoo 1575R San Diego Zoo SDZ28-14L

SDWAP WAP20-182 SDWAP WAP20-395 SDWAP WAP20-273 SDWAP AO893172 SDWAP AO893180 SDWAP AO893172 C. G. Sibley 2687 LSU B5273 ANSP 3835 ANSP 3835

C. episcopus Woollynecked stork

DNA preparation no.

31 45 55

L. javanicus Lesser adjutant stork

74 108 109

Fort Wayne Children’s Zoo 1421 LSU B6751 Fort Wayne Children’s Zoo 1420 Audubon Zoo 156 Miami MetroZoo 3738 Miami MetroZoo 3737

5 59 66 91 110*

San Diego Zoo 161R LSU B10925 LSU B10925 San Diego Zoo 161R San Diego Zoo 162R

70 72*

PHYLOGENY OF STORKS

TABLE 1—Continued

Species Cathartes aura Turkey vulture

DNA preparation no.

Source/source ID number

50* 118 119

ANSP 395 ANSP 395 ANSP 395

Note. Two DNA preparations were made for some individuals. Asterisks mark DNA preparations radiolabeled and used as tracers in the DNA–DNA hybridization matrix. DNA preparations that were sequenced and included in the cytochrome b data matrix are underlined. The following institutions provided samples: Academy of Natural Sciences of Philadelphia (ANSP), PA; Audubon Zoo, LA; Fort Wayne Children’s Zoo, MI; Louisiana State University (LSU), LA; Miami MetroZoo, FL; Philadelphia Zoo, PA; San Diego Zoo, CA; San Diego Wild Animal Park (SDWAP), CA.

ization distances and mitochondrial cytochrome b sequences. Morphological and behavioral data have been collected for species in the family, so the molecular data provide an additional source of phylogenetic information, independent in the sense that the molecular, morphological, and behavioral characters differ in evolutionary rates and in patterns and probabilities of transformation among character states. In many respects, the two molecular data sets are complementary. The cytochrome b sequences yield a mitochondrial gene tree; the DNA–DNA hybridization data provide a measure of sequence divergence in the single-copy portion of the nuclear genome. While DNA–DNA hybridization data generally cannot resolve relationships among closely related species, cytochrome b sequences usually resolve such relations well. The estimates of phylogeny obtained from these molecular data sets are compared with each other and previous assessments of relationships based on morphological and behavioral data. MATERIALS AND METHODS Choice of Taxa Table 1 lists the species included in this study, sources of blood/tissue, and DNA preparation numbers. The study includes individuals of 17 of the 19 species (and all genera) in the family Ciconiidae, with Anastomus oscitans and Leptoptilos dubius missing. Fourteen species are represented by more than one individual. Most DNA preparations are from blood samples from captive (zoo) birds. About 15% of the zoo birds are wild-caught individuals; 65%, captive-born; and 30%, of unknown birth type. For the wild-caught birds, locality information is vague or lacking. Outgroups are an essential ingredient for any phylogenetic analysis. With DNA–DNA hybridization data, inclusion of an outgroup is required to assess relative rates of genetic change among ingroup species, allowing the partitioning of symplesiomorphic and synapo-

277

morphic similarity (Sheldon, 1994). For sequence data, outgroups are needed to polarize character state changes (Maddison et al., 1984; Smith, 1994). In this study, outgroups were selected from two families: Threskiornithidae, the ibises and spoonbills, and Cathartidae, the New World vultures. The monophyly of storks with respect to the cathartids and ibises is supported by numerous morphological and behavioral features (Kahl, 1972d; Hancock et al., 1992), but the pattern of relationships among storks, ibises, and cathartids remains a subject of considerable debate (Sibley and Ahlquist, 1990; Griffiths, 1994; Seibold and Helbig, 1995). Storks and cathartids share morphological and behavioral traits, including a common pectoral musculature, an absence of intrinsic syringeal muscles, and the habit of defecating on the legs to dissipate body heat (Ligon, 1967; Olson, 1979). Storks and ibises traditionally have been considered close relatives based on similarities in behavior, nesting habits, and the pattern of feather tracts on the body (Sibley and Ahlquist, 1990), although Olson (1979) argues, based on osteological features, that storks and ibises are distantly related. Given the available evidence and conflicting opinions, I included representative(s) from both the cathartid and ibis/ spoonbill families as outgroups in both the sequencing and the DNA–DNA hybridization data sets. DNA Extraction For most individuals, whole-genomic DNA was extracted from erythrocytes in 10% EDTA; a few extractions were from frozen tissue (280°C) or alcoholpreserved tissue. DNA extraction followed standard protocols (Maniatis et al., 1982). Frozen tissue (0.5 mg) was homogenized in 2–5 ml of cold STE (sucrose/Tris/ EDTA) buffer; alcohol-preserved tissue (0.5 mg) was vacuum-dried, soaked in liquid nitrogen, immediately ground to a fine powder, and mixed with 2–3 ml of STE buffer. Approximately 3–7 ml of tissue homogenate in STE or 1–2 ml of blood diluted with 2–4 ml of STE was digested overnight at 37°C with 1% SDS (sodium dodecyl sulfate) and 5–10 mg/ml of pronase. Following digestion, the extract was purified by mixing once with phenol:chloroform (1:1) and twice with chloroform. DNA was precipitated with cold ethanol. DNA extracted from tissue was also incubated with RNase. DNA–DNA Hybridization Methods The DNA hybridization data included representatives from 15 of the 19 stork species; all six genera were represented. Not included were 2 species that were formerly considered subspecies (Ciconia (ciconia) boyciana, Ciconia (episcopus) stormi) and two species for which no DNA was available (A. oscitans, L. dubius). Plegadis falcinellus (glossy ibis; family Threskiornithidae) and Cathartes aura (turkey vulture; family Cathartidae) were included as outgroups.

278

BETH SLIKAS

Biochemistry. Whole-genomic DNA was extracted as described above. General procedures for sonification of the purified extract and preparation and labeling of single-copy DNA followed Sibley and Ahlquist (1990) and Sheldon et al. (1992), with the following specifications. All purified DNA extracts were sonified according to a standard protocol refined and tested in our lab to produce a modal fragment size of about 500 base pairs; sonified extracts were not sized. Single-copy tracer DNA was prepared to C0t 1000 (Werman et al., 1990) and labeled with tritiated thymidine by random priming (Feinberg and Vogelstein, 1983). Hybridization reactions were set up with 40,000 to 150,000 dpm of labeled single-copy tracer DNA (i.e., about 0.004 µg) and 10 µg of sonified, unlabeled driver DNA (tracer: driver about 1:2500). Reactions were boiled for 10 min, incubated at 60°C in 0.48 M phosphate buffer to C0t $ 22,000, and then fractionated on 1.0-ml hydroxylapatite (HAP) columns placed in a thermal-elution device similar to those described by Sibley and Ahlquist (1981) and Kirsch et al. (1990). Fractions were collected at 2.5°C intervals beginning at 60°, ending at 95°C. Sample melting curves are reproduced in Fig. 1. Matrix design. For each of the 17 species in the study, a representative individual was labeled and hybridized to itself (homoduplex hybrid) and several different individuals of the same and different species (heteroduplex hybrids) to produce a complete matrix of replicated comparisons. A single fractionation run included 35 hybrids, all with the same labeled (tracer) DNA. Two runs were performed for each tracer, each run including 3 homoduplex hybrids and 1–5 heterodu-

FIG. 1. Sample DNA–DNA hybridization melting curves. Each curve represents a single hybrid, and the data points are percentage counts eluted at each temperature. For these curves, the tracer (labeled DNA) is a Ciconia ciconia individual (Sample 2). The curve marked by closed circles is the homoduplex, Ciconia ciconia 2 hybridized to itself. The other curves are heteroduplexes: Ciconia ciconia 2 hybridized to (1) Ciconia maguari (inverted closed triangles), (2) Mycteria americana (inverted open triangles), and (3) Plegadis falcinellus (open squares).

plex hybrids for each unlabeled (driver) species. Melting statistics are known to vary among individuals and among different DNA preparations from the same individual (Bleiweiss and Kirsch, 1993). To assess this variation, driver DNA from different individuals or DNA preparations was used for each of the replicate heteroduplexes in most matrix cells. For Anastomus lamelligerus, only a single DNA preparation was available. Calculation of indices. Tm, Tmode, T50H, and normalized percentage reassociation (NPR) were calculated for every hybrid, as described by Sheldon and Bledsoe (1989). Tmode was determined by fitting a curve to the melting data (% counts eluted versus fractionation temperature), using the Jandel Scientific computer package ‘‘Peakfit’’; two equations were used to perform the fit: a modified Fermi-Dirac and an asymmetric double sigmoid. Because both sets of fitted DTmode values yielded the same best-fit topology, trees presented here are based on the modified Fermi-Dirac modes for consistency with previous DNA–DNA hybridization studies (e.g., Sheldon et al., 1992; Sheldon and Winkler, 1993). Measures of genetic dissimilarity, delta (D) values, were calculated for every heteroduplex hybrid by subtracting the heteroduplex value for a given index from the corresponding homoduplex value, the latter averaged over the replicate homoduplexes in the fractionation run. Uncorrected average values and standard deviations for DTmode values are provided in Table 2. Certain data were excluded a priori from analysis. Entire runs were omitted if (1) major mechanical problems occurred during the run or (2) the average homoduplex Tm was lower than 82°C, because such a low Tm suggests that the labeled DNA sample is short-stranded (Springer and Kirsch, 1991). Individual hybrids with unusually low NPR were omitted, because low NPR suggests incomplete hybridization. The minimal acceptable NPR was set at 80% for ingroup hybrids (stork-to-stork) and 70% for hybrids involving an outgroup (stork-to-Cathartes or Plegadis, Cathartes-toPlegadis). If a hybrid melting curve could not be fit reasonably with either of two fitting equations (modified Fermi-Dirac or asymmetric double sigmoid) to determine its mode, then the hybrid was omitted from all analyses under the presumption that the data were poorly distributed (e.g., because of bimodality; see Sarich et al., 1989). Goodness-of-fit was measured by the sum-of-squares error; hybrids with a sum-ofsquares error greater than 2.0 were omitted. However, if the hybrid melting curve could be fit reasonably well with one of the two fitting equations, then only the fitted mode with the large sum-of-squares error (.2.0) was omitted; all other indices for that hybrid were retained. All data for certain individuals or DNA preparations which appeared to be short-stranded were excluded. The excluded preparations were Samples 25, 70, 71,

279

PHYLOGENY OF STORKS

104, 89, and 106 (Table 1). When used as drivers, these preparations regularly yielded much larger delta values than conspecific drivers, with the discrepancy often exceeding 1°C. Two of these preparations (71, 106) were used as labels and gave low homoduplex Tm’s (80.3 and 81.1°C for 71; 80.6 and 81.3°C for 106), corroborating the assumption that these preparations are shortstranded. Samples 104 and 71 are different preparations from the same individual. One apparently short-stranded sample was not excluded, A. lamelligerus 2687, because it is the only individual representing the species and the genus. As expected for a short-stranded sample, reciprocal distances involving Anastomus were asymmetric (Sheldon and Bledsoe, 1989); measurements obtained from hybrids in which Anastomus was the driver gave consistently larger distances than the reciprocals with Anastomus as the tracer (average difference in DTmode 5 2.92°C, range 2.41–4.64°C). Hybrids including a shortstranded sample melt at lower temperatures, which increases calculated distances when the short-stranded sample is the driver. When the short-stranded sample is the tracer, calculated distances are not increased, because the homoduplex also melts at lower temperatures, and distances are calculated by subtracting the melting index of each heteroduplex from that of the homoduplex. Due to the substantial difference between reciprocal measurements for Anastomus and the likely cause of this asymmetry, all distances with Anastomus as the driver were excluded from analysis. The empty matrix cells were filled with the reciprocal distance measures. For several hybrids, Tm and DTm values were excluded because the initial 60°C elution was missed during the fractionation run, making calculation of Tm impossible. Finally, one homoduplex hybrid that showed an anomalous secondary peak in its melting distribution was omitted, and one hybrid that showed unusual ‘‘bumpiness’’ at the low-temperature end of its melting curve (i.e., an unusually large number of counts coming off below the melting peak) was excluded. In the end, about 15% of the hybrids were excluded prior to phylogeny estimation. The final phylogenetic results presented here are based on analyses of the distance measure DTmode. The measures DTm and DTmode are highly correlated (correlation coefficient R 2 5 0.96) and yield virtually identical branching topologies. However, DTm values become compressed when labeled DNA samples are shortstranded. The short-stranded label effectively increases the stringency of the hybridization, and fewer divergent hybrids form than if the fragment length of the label were longer. DTmode is not affected by compression, because Tmode is determined by the position of the peak in the melting curve and is not influenced by the low-temperature tail of the melting distribution (Sheldon and Bledsoe, 1989; Springer and Kirsch, 1991). In

this study, one stork species, A. lamelligerus, is represented by a single, short-stranded sample, so phylogenetic analyses were based on DTmode to avoid the problem of compression. The distance measure DT50H is more variable than either DTm or DTmode, rendering it less useful in phylogenetic analyses, except in comparisons of very distantly related taxa (Sheldon and Bledsoe, 1989; Goodman et al., 1990; Schmid and Marks, 1990). Phylogenetic analyses were performed on matrices of uncorrected (DTmode ) and corrected (DTmodeC) distances. The correction procedure consisted of two steps. First, DTmode was converted to percentage sequence difference (d), using the empirical formula of Springer et al. (1992): d 5 1.18(DTmode /100).

(1)

The d values were then corrected for multiple substitutions at sites using the Jukes-Cantor correction (Jukes and Cantor, 1969), assuming a 60:40 AT:GC ratio (Arthur and Strauss, 1978; Epplen et al., 1978; Swofford and Olsen, 1990): DTmC 5 (100)(20.74)(ln 51 2 1.35d6).

(2)

These corrections were intended to increase the additivity of the distances, because the tree-building algorithms assume additivity (Springer and Krajewski, 1989a,b). Tree construction. Trees were constructed by the FITCH algorithm in PHYLIP (version 3.5c; Felsenstein, 1995), using unweighted least squares (CavalliSforza and Edwards, 1967). The unweighted leastsquares method is appropriate when measurement error does not increase with increasing distance values, as is the case with these data (Fig. 2) and most DNA hybridization data sets (e.g., Sheldon, 1987a; Bleiweiss et al., 1994). In any case, similar best-fit trees were obtained with weighted least squares (Fitch and Margoliash, 1967), differing only in the branching order among the three Old World Mycteria species, which are the most weakly supported nodes in the tree. The FITCH program was run with the input order of taxa randomized and the subreplicate and global optimization options enabled. Negative branch lengths were not allowed. The FITCH algorithm does not assume uniformity of rates across taxa. The DTmode and DTmodeC data matrices also were run through the neighbor-joining tree-building algorithm (Saitou and Nei, 1987) in PHYLIP. Neighbor-joining is a sequential clustering method that produces a single final tree; the neighborjoining method does not assume equal rates of change among lineages (Swofford and Olsen, 1990). Robustness of tree topology with respect to random experimental error was assessed by bootstrapping (Fel-

280

BETH SLIKAS

TABLE 2 Matrix of Average DTmode (C) Values

senstein, 1985; Lanyon, 1985). Bootstrap pseudomatrices were constructed by sampling, with replacement, from replicate distance measures within each cell of the DTmode and DTmodeC data matrices. Each pseudomatrix was then analyzed with the FITCH program, with randomized input order of OTUs and P 5 0; 1000 pseudomatrices were generated and analyzed (Krajewski and Dickerman, 1990; Dickerman, 1991). The result-

ing trees were input into the CONSENSE program in PHYLIP, which produces a majority-rule consensus tree, with bootstrap percentages given for each monophyletic group. The bootstrap percentages indicate the percentage of trees generated from the pseudoreplicate matrices in which the designated group occurs. Jackknifing, sequential removal of single OTUs and subsequent analysis of the reduced data matrices

PHYLOGENY OF STORKS

281

TABLE 2—Continued

(Lanyon, 1985), was performed to assess the effect of individual OTUs on tree stability. Jackknifing runs were performed with both the DTmode and the DTmodeC data matrices. Sequencing Methods Amplification. An approximately 1.1-kb fragment of the mitochondrial genome, including most of the cytochrome b gene, was amplified via the polymerase

chain reaction (PCR; Saiki et al., 1988). Doublestranded amplifications were performed in 50-µl reaction volumes including 13 buffer mix (Promega; 103 buffer includes 500 mM KCl, 100 mM Tris–HCl, pH 9.0 at 25°C, and 1.0% Triton X-100), 200 µM concentrations of each deoxynucleotide, 0.8 µM concentrations of each primer, 1.5–3.0 mM MgCl2, 0.25 unit of Taq polymerase (Promega), and 25 ng of whole-genomic DNA as template. All samples were amplified using the primer pair L14990 and H16065 (Fig. 3), where ‘‘L’’ and ‘‘H’’ refer to

282

BETH SLIKAS

TABLE 2—Continued

FIG. 2. Plot of standard deviation versus DTmode. Correlation coefficient R 2 5 0.11.

turation at 94°C, and the last cycle was followed by an 8-min extension step at 72°C. The double-stranded PCR product was run out on a 1% agarose gel and stained with ethidium bromide. The product band was excised from the gel and purified of

Note. Standard deviation/number of measurements are listed beneath the mean in each cell.

the light and heavy strands of the mitochondrial genome, respectively, and numbers refer to the location of the 38 end of the primer in the Gallus gallus sequence (Desjardins and Morais, 1990). The PCR mix was subjected to 35 cycles of denaturation at 93°C (50 s), annealing at 41°C (1 min 15 s), and extension at 72°C (2 min 15 s) in an automated thermocycler (M. J. Research). The first cycle was preceded by a 2-min dena-

FIG. 3. Diagram depicting location of primers used for amplification and sequencing of the cytochrome b gene. The numbers below the tick marks indicate nucleotide position in the Gallus gallus mitochondrial genome (Desjardins and Morais, 1990). Numbers above and below the arrows refer to particular primers; primer position and sequence are listed below the diagram. Light-strand primers are deignated with an ‘‘L,’’ and heavy-strand primers with an ‘‘H’’; numbers following strand designation indicate the position of the 38 end of the primer in the G. gallus mitochondrial genome (Desjardins and Morais, 1990).

283

PHYLOGENY OF STORKS

remaining nucleotides and excess primers with glassmilk (GeneClean, Bio 101). The double-stranded product was sequenced directly using dideoxy chain termination with T7 DNA polymerase (Sequenase 2.0; United States Biochemical) and 33P-labeled dATP (Sanger et al., 1977). Annealing of the sequencing primer was accomplished by boiling the annealing mix and then rapidly immersing it in ultracold ethanol (280°C). Primers used in sequencing are listed in Fig. 3. Complete sequences were obtained for the heavy strand only; partial sequences were obtained from the light strand for some samples. Between 1005 and 1065 bp of the mitochondrial cytochrome b gene were sequenced for 16 stork species and an outgroup, P. falcinellus. The sequence data matrix (Table 3) begins at position 14,894 in the G. gallus mitochondrial sequence and ends at the cytochrome b termination codon (Desjardins and Morais, 1990). For all species except Mycteria americana, C. boyciana, C. stormi, and P. falcinellus, at least a partial sequence from a second individual was obtained ($750 bp), primarily to check for possible sample misidentification or sequencing errors. The sequences of Leptoptilos crumeniferus and three outgroup species from the family Cathartidae (C. aura, turkey vulture; Sarcoramphus papa, king vulture; and Vultur gryphus, Andean condor) were taken from GenBank (Seibold and Helbig, 1995). All sequences could be aligned unambiguously with each other and the G. gallus sequence (Desjardins and Morais, 1990). Sequence entry, editing, and alignment were performed using the GCG software package (Genetics Computer Group, 1993), available through the National Cancer Institute computer facility. PCR amplification of mitochondrial genes frequently results in the amplification of supposed nuclear homologues, or pseudogenes (e.g., Smith et al., 1992; Arctander, 1995). To confirm that the sequences obtained were from functional mitochondrial cytochrome b genes, rather than nuclear copies, sequences were converted to amino acids and checked for internal stop codons and substitutions at conserved amino acid sites (i.e., sites that are invariant in other vertebrates, based on Esposti et al., 1993). In addition, the pattern of substitutions across codon positions was checked for each sequence, relative to other sequences in the data matrix. Functional cytochrome b genes accumulate substitutions more rapidly at third codon positions than at first and second positions combined, while nuclear pseudogenes are expected to show a more even distribution of substitutions across codon positions (Arctander, 1995). Phylogenetic analyses of the cytochrome b sequence data included maximum likelihood, parsimony with transversions weighted 53 transitions, and a weighted least-squares analysis of distances calculated using the maximum-likelihood formula in DNADIST (Felsenstein, 1995). In the maximum-likelihood and parsi-

mony analyses, question marks in the data matrix were treated as missing data. When distances were calculated, for each pair of taxa compared, sites with a question mark for either taxon were excluded from the calculation. RESULTS DNA–DNA Hybridization Characterization of data. A total of 1247 DNA hybrids consisting of 18,705 melting-curve fractions contributed to the final estimate of phylogenetic relationships. Comparisons included 61 DNA preparations from 50 individuals of storks, 3 preparations from one C. aura individual, and 3 preparations from 2 P. falcinellus individuals. Average values of the distance DTmode range from 0.43 to 2.89°C between ingroup taxa and from 5.67 to 7.98°C between ingroup and outgroup taxa (Table 2; excluding the distance measures obtained with A. lamelligerus as a driver, which are inflated by the degraded quality of the sample—see DNA Hybridization Methods, above). The average standard deviation, averaged across all off-diagonal cells of the data matrix, is 0.31°C (range 0.006–0.933°C) for DTm, 0.32°C for DTmode (0.002–0.930°C), and 0.53°C (0.002–2.115°C) for DT50H. These standard deviations fall near the upper end of the range of values obtained in other DNA hybridization studies of avian families (0.20–0.30°C for DTm and DTmode; Sheldon, 1987b; Krajewski, 1989; Sheldon and Bledsoe, 1989; Sheldon et al., 1992; Sheldon and Winkler, 1993; Bleiweiss et al., 1994). As is typical for DNA–DNA hybridization data, measurement precision for DT50H is less than that for either DTmode or DTm (Sheldon and Bledsoe, 1989; Sheldon and Winkler, 1993). Reciprocal measurement discrepancy was assessed using the SYMBOOT program (A. J. Dickerman, personal communication), which calculates and corrects matrix asymmetry via the method of Sarich and Cronin (1976). Discrepancy between reciprocal measurements can arise for several reasons, but for DNA hybridization data at this taxonomic level the most probable cause of discrepancy is short-strandedness of individual DNA preparations (Springer and Kirsch, 1991). Percentage nonreciprocity is defined as 1003 the reciprocal differences between matrix cells divided by the reciprocal sums. For these data, mean percentage nonreciprocity was 7.58% for DTm and 6.25% for DTmode. Because percentage nonreciprocity varies inversely with absolute distance, comparison of values across studies is problematic, but the values obtained for these data are similar to those obtained in other DNA hybridization studies using similar protocols (e.g., Springer et al., 1990; Sheldon et al., 1992; Sheldon and Winkler, 1993; Bleiweiss et al., 1994). Rates. Homogeneity of evolutionary rates among stork species was assessed by the relative rate test,

284

BETH SLIKAS

TABLE 3 Cytochrome b Sequence Data Matrix

PHYLOGENY OF STORKS

TABLE 3—Continued

285

286

BETH SLIKAS

TABLE 3—Continued

PHYLOGENY OF STORKS

TABLE 3—Continued

287

288

BETH SLIKAS

TABLE 3—Continued

Note. Sequences begin at position 14971 in the Gallus gallus mitochondrial genome (Desjardins and Morais, 1990) and end at position 16035 (termination codon). Sequences have been submitted to GenBank under Accession Numbers U72771–U72786 and U70822.

which compares measured distances from an outgroup species to each ingroup species (Sarich and Wilson, 1967). DTmode distances from labeled P. falcinellus to each stork species were compared using ANOVA and the nonparametric Kruskal-Wallis test (Sokal and Rohlf,

1981). With all species included, the ANOVA indicated significant rate differences (F14,30 5 4.33, P , .001), but with A. lamelligerus excluded from the test, no significant differences were found (F13,30 5 0.88, P . 0.5). The distance from A. lamelligerus to Plegadis is large

PHYLOGENY OF STORKS

compared to distances to other stork species, probably due simply to the short-strandedness of the sample. With the Kruskal-Wallis test, no significant difference in rates was indicated, even with all species included (H 5 15.74, v 5 14, P . 0.5). Comparisons were made to Plegadis rather than Cathartes, because the average homoduplex Tm’s of the labeled Cathartes were relatively low (81.9 and 82.4°C), so distances between Cathartes and stork species were probably compressed and therefore less likely to show rate differences (Springer and Kirsch, 1991). Phylogeny. The best-fit FITCH tree, based on an unweighted least-squares analysis of the 17 3 17 corrected DTmode matrix is depicted in Fig. 4. The same topology was obtained with either Plegadis or Cathartes declared as the outgroup or with only one of these two taxa included in the analysis. FITCH analysis of the corrected DTm matrix yielded a topology differing from the DTmode tree only in the position of A. lamelligerus, which shifted to the base of the tree, as the sister taxon to the rest of the family. This result is probably an artifact of distance compression between the short-stranded Anastomus sample and the outgroups. The effect of compression increases with increasing distance and would be most severe in comparisons to outgroup taxa, with the possible result that a shortstranded sample would be pulled toward the base of the FITCH tree, near the outgroup(s). Because DTmode distances are not affected by compression (see Materials and Methods), the tree based on DTmode is considered more reliable. The corrected DTmode matrix was subjected to bootstrap resampling of replicate distance measures and jackknife resampling of taxa, using the FITCH algorithm with unweighted least squares (i.e., P 5 0). The

289

majority-rule bootstrap consensus tree is identical in topology to the best-fit tree, except that Mycteria leucocephala pairs with M. ibis, rather than with M. cinerea. Bootstrap percentages are indicated on the best-fit tree (Fig. 4). The jackknife strict-consensus tree is congruent with the best-fit tree, but less resolved. Nodes in the best-fit tree that appear in fewer than 85% of bootstrap replicates collapse in the jackknife strict-consensus tree, except the node defining the Mycteria–Anastomus– Leptoptilos clade. The DNA–DNA hybridization data strongly support the monophyly of the Ciconiidae and of the genera defined by Kahl (1971a, 1972d). The Ciconiidae and all genera represented by more than a single species occur as monophyletic groups in 100% of the trees based on bootstrap resampling of the data. Additionally, each of the species pairs Ciconia ciconia–C. maguari and C. episcopus–C. abdimii has 100% bootstrap support. The pairings of A. lamelligerus with the Mycteria species (88%) and Jabiru mycteria with the Ephippiorhynchus species (96%) are also strongly supported. The data also support two clades composed of the genera (1) Ciconia– Ephippiorhynchus–Jabiru and (2) Leptoptilos–Anastomus–Mycteria (bootstrap of 84%). Within the genus Mycteria, relationships among the three Old World species (M. ibis, M. leucocephala, and M. cinerea) are weakly resolved; the position of the New World M. americana as the sister to its congeners is moderately supported. The position of Ciconia nigra with respect to its congeners is unresolved. Cytochrome b Sequences Characterization of data. For the stork Leptoptilos crumeniferus, a cytochrome b sequence was taken from

FIG. 4. Best-fit FITCH tree based on corrected DTmode distances. FITCH analysis was performed with the following parameters: P 5 0, negative branches not allowed, subreplicate option enabled, global branch-swapping, jumbled input order of taxa (two repititions). Branches are drawn proportionally to the amount of change along the branch. Numbers above the branches are bootstrap percentages; 1000 bootstrap pseudomatrices were generated and analyzed. The same topology results with either Plegadis or Cathartes designated as the outgroup.

290

BETH SLIKAS

GenBank for use in phylogenetic analyses (Seibold and Helbig, 1995). Two individuals of L. crumeniferus were sequenced in this study (Samples 15 and 70), but the sequences are suspect. First, the two differ from each other by 5.12%, an unusually large intraspecific divergence (see Table 4). Compared to the GenBank sequence, Sample 70 is 2.24% divergent and Sample 15, 7.34%. In addition, the sequences for Samples 70 and 15 share a unique amino acid replacement (proline = serine) at a site that is almost invariant in other animals (residue 286 in G. gallus; Desjardins and Morais, 1990); the only other reported variant is a deletion in Xenopus laevis (African Toad; Esposti et al., 1993). Despite these abnormalities, the sequences obtained for Samples 15 and 70 are unambiguous and have no insertions, deletions, or internal stop codons. In comparisons to other sequences in the data matrix, both show a substitution pattern typical for a mitochondrial coding gene, with the majority of substitutions occurring at third codon positions. The sequences for Samples 15 and 70 and the GenBank L. crumeniferus sequence all show similar percentage divergence compared to other stork sequences, in the range of 8–10%. This combination of features suggests that the sequences obtained for Samples 15 and 70 might be from a recently incorporated nuclear pseudogene, rather than from the mitochondrial cytochrome b gene (Smith et al., 1992; Arctander, 1995). If all three sequences are included in a phylogenetic analysis, the three group together as a clade. The GenBank sequence was selected as the single representative L. crumeniferus sequence in the final phylogenetic analyses, because the sequences for Samples 15 and 70 have one unusual amino acid replacement, noted above. For the cytochrome b data, sequence divergence between conspecific individuals (excluding L. crumeniferus; see above) ranged from 0.0 to 2.1%. The two largest intraspecific divergences were between individuals of Ephippiorhynchus asiaticus (2.1%) and of C. episcopus (0.91%). Both these species have widespread distributions: E. asiaticus occurs in India, southeast Asia, southern New Guinea, and northern Australia; C. episcopus occurs in tropical and subtropical Africa, India, southeast Asia, the Sunda Islands, and Celebes (Hancock et al., 1992). Of the two E. asiaticus samples sequenced, one (106) is an Australian bird, and one C. episcopus sample (117) is an Asian bird; the source localities of the other two samples are unknown. For both species, it is possible that the two individuals sequenced were from widely separated populations, possibly on different continents. Alternatively, the relatively large intraspecific differences between these sequences might indicate that one or both in each pair are erroneous, but these sequences have no other aberrant features. No multiple bands or ambiguities were observed on the sequencing gels, and the translated amino acid sequences exhibit no improbable substitu-

tions (Arctander, 1995). Due to the relatively large intraspecific divergence, both individuals sequenced were included in the data matrix. Sequence divergence spanned a wide range between species within the stork genera (Table 4): 3.0–12.1% within Ciconia (7 spp.), 0.9–6.7% within Mycteria (4 spp.), 7.9–8.4% within Ephippiorhynchus (2 spp.), and 8.2% within Leptoptilos (2 spp.). Between species in different genera, divergences spanned a narrower range (7.71–11.63%) and averaged only slightly greater than the largest intrageneric values (intrageneric average 5 8.01%, intergeneric average 5 9.75%), possibly due to saturation effects. Distances to the outgroups averaged 13.73% for P. falcinellus and 15.52% for C. aura and did not differ significantly among stork species, suggesting that rates of sequence change in cytochrome b among stork lineages are equal. However, the probable saturation of sequence divergences between stork species and the outgroups would tend to mask any existing rate differences. Patterns of variability in this data set are typical for cytochrome b (Kocher et al., 1989; Edwards et al., 1991; Irwin et al., 1991). Of the 1065 nucleotide sites in the data matrix, 311 sites were variable and 231 were phylogenetically informative among stork species. Of the 311 variable sites, 247 occurred at third codon positions, 49 at first positions, and 15 at second positions. Of the 231 phylogenetically informative sites, 190 occurred at third positions, 31 at first positions, and 10 at second positions. In terms of amino acids, of a total of 355 sequenced, 47 were variable and 30 phylogenetically informative within storks; 57 were variable and 40 phylogenetically informative among the storks and C. aura; and 52 were variable and 41 phylogenetically informative among the storks and P. falcinellus. Across codons, as well as within codons, nucleotide variability in cytochrome b appears to be constrained by functional properties. The cytochrome b molecule weaves back and forth across the mitochondrial membrane, so it can be divided into inner membrane, outer membrane, and transmembrane regions. Empirically, the transmembrane regions exhibit the greatest amount of variability and the outer regions, the least. The lower variability of regions lying outside the membrane is believed to reflect the need for proper protein–protein interactions between cytochrome b and an iron–sulfur subunit that plays a role in oxidation of ubiquinol (Esposti et al., 1993). I determined amino acid variability as a function of position along the cytochrome b gene for the data matrix including 16 stork species, P. falcinellus, and C. aura (using MEGA; Kumar et al., 1993). The observed pattern of variability was consistent with that reported for other vertebrate groups (Irwin et al., 1991; Esposti et al., 1993). The transmembrane regions of the molecule had the highest percentage of variable amino acids (57/171 5 33.3%), followed by the inner regions

PHYLOGENY OF STORKS

291

TABLE 4 Distances Calculated from Cytochrome b Sequence Data

Note. Numbers below the diagonal are uncorrected percentage divergences; sites with missing data were excluded on a pairwise basis. Numbers above the diagonal are maximum-likelihood distances calculated with the DNADIST program PHYLIP (Felsenstein, 1995) with the following parameters: transition:transversion ratio of 5:1; three categories of substitution rates, of relative magnitude 1:0.3:30; rates were assigned to sites based on the output of a maximum-likehood analysis with the same parameters.

292

BETH SLIKAS

(13/72 5 18.1%) and the outer regions (8/111 5 7.2%). Among the transmembrane regions, regions E and H were the most variable (52.2 and 43.5% of residues variable, respectively), and A, B, and C were the least variable (10, 5.3, and 8.0% of residues variable). The C-terminus end of the molecule (adjacent to the threonine tRNA gene) was the most variable portion, again in accord with data from other animal groups (Esposti et al., 1993). The nucleotide composition of the cytochrome b gene is skewed. A strong bias against guanine at third codon positions appears to be common to all vertebrates; a significant, but less drastic bias against thymine at third positions appears to be characteristic of mammals and birds, but not fish or amphibians. At second codon positions, a deficiency of guanine and surplus of thymine are common to at least birds and mammals (Kocher et al., 1989; Irwin et al., 1991; Edwards et al., 1991; Graybeal, 1993). The compositional bias in the stork data set coincides with biases observed in cytochrome b sequences of other avian groups (e.g., Edwards et al., 1991; Lanyon and Hall, 1994; Cicero and Johnson, 1995; Nunn and Cracraft, 1996). Base composition is most variable at third codon positions, as reported previously for birds and other vertebrates (Irwin et al., 1991; Lanyon and Hall, 1994; Cicero and Johnson, 1995). A plot of the number of transitions and transverions, partitioned by codon position, versus DNA–DNA hybridization distance between species pairs (Fig. 5) suggests saturation of transitional changes, particularly at third codon positions. The transition curves plateau, and even decline with increasing divergence between taxa, a typical pattern for cytochrome b sequence comparisons (e.g., Irwin et al., 1991; Moritz et al., 1992; Cicero and Johnson, 1995; Hackett, 1996; Lara et al., 1996). Phylogeny. Phylogenetic analyses of the cytochrome b sequence data included maximum likelihood, parsimony with transversions weighted 53 transitions, and a weighted least-squares analysis of distances calculated using the maximum-likelihood formula in DNADIST (Felsenstein, 1995). The maximum-likelihood analysis of the cytochrome b sequence data included 20 taxa: 16 stork species, 3 cathartids (C. aura, S. papa, and V. gryphus) and the designated outgroup, P. falcinellus. The analysis was performed with the DNAML program in PHYLIP (Felsenstein, 1995). This program requires input of the following parameters: (1) transition:transversion (ts: tv) ratio and (2) number, relative magnitude, and probability of different substitution rate categories across sites. To obtain reasonable estimates of these parameters for this data set, a preliminary tree was estimated using a FITCH weighted least-squares analysis of Kimura-corrected sequence divergence (ts/tv 5 2; Kimura, 1980). This tree was input as a user tree in the DNAML program with the cytochrome b sequence data.

FIG. 5. Plot of number of cytochrome b transitions (above) and transversions (below) against DNA–DNA hybridization distance (corrected DTmode ) for all possible species pairs. Points below 4°C represent stork–stork pairs, and points above 8°C represent stork– outgroup pairs. Third codon positions (open circles) are plotted separately from first and second codon positions (open triangles). The solid lines are regressions. Transition data were fitted with secondorder regression; R 5 0.70 for third positions, and R 5 0.79 for first and second positions combined. Transversion data were fitted with a linear regression; R 5 0.96 for third positions, and R 5 0.84 for first and second positions combined. Both transition curves suggest saturation of changes.

I then conducted several maximum-likelihood analyses varying the ts:tv ratio and rate parameters to find the combination that gave the best ln(likelihood) for the input tree. The combination of a ts:tv ratio of 5:1 and three rate categories with equal probability (0.333, 0.333, 0.333) and relative magnitudes 1:0.3:30 gave at least a locally maximum likelihood. These parameters were then input in a search for the best tree; the search was done with global rearrangements, and P. falcinellus designated as the outgroup. The maximum-likelihood tree is depicted in Fig. 6. A parsimony analysis of 16 stork species (18 individuals, Table 3) was performed using PAUP (Swofford, 1993), with P. falcinellus, C. aura, S. papa, and V. gryphus declared as outgroups. A stepmatrix weighting transversions over transitions by 5:1 was applied to all characters; this weighting reflects the expectation that

PHYLOGENY OF STORKS

293

FIG. 6. Maximum-likelihood tree based on cytochrome b sequences in Table 3. Analysis was performed using the DNAML algorithm in PHYLIP (Felsenstein, 1995), assuming a 5:1 transition/transversion ratio and three categories of substitution rate, of relative magnitude 1:0.3:30 and equal probability. Circled branches are not significantly different from zero length. All branches are drawn proportionally to the amount of change along the branch.

transitions occur 5 times more frequently than transversions, which is the ratio used in the maximumlikelihood analysis. An heuristic search, using TBR branch-swapping, with a 100-replicate, random-addition sequence of taxa and the MULPARS option in effect, produced a single most-parsimonious tree (Fig. 7). I used the Bremer support index for each node in the most-parsimonious tree to assess reliability (Bremer,

1994). In addition, the data were bootstrapped 100 times (Felsenstein, 1985), and each bootstrap pseudomatrix was analyzed with the same search options used for the original data matrix, except that the number of addition-sequence replicates was reduced from 100 to 10. On the most-parsimonious tree, 856 changes (676 transitions) occurred at third codon positions, 124 (102

FIG. 7. Single most-parsimonious tree resulting from a PAUP (Swofford, 1993) analysis of the cytochrome b sequences in Table 3. A heuristic search was performed with random addition of taxa (10 addition-sequence replicates) and TBR branch-swapping; a stepmatrix was also imposed weighting transversions five times more heavily than transitions. Numbers above the branches are Bremer support indices, and numbers below the branches are bootstrap percentages. Branches are drawn proportionally to the amount of change along the branch. Plegadis falcinellus, Cathartes aura, Vultur gryphus, and Sarcoramphus papa were declared as outgroups.

294

BETH SLIKAS

ts) at first positions, and 37 (35 ts) at second positions. A plot of transformation frequencies shows that C = T changes were by far the most frequent, followed by A = G changes. A = C, T = C, and G = A changes were the next most frequent and were about equally common. Other types of changes were relatively rare. This pattern of transformation frequencies is typical for cytochrome b in birds (Hackett, 1996; Cicero and Johnson, 1995). Although the most common changes included all possible transitions, different types of transition were not equally probable, and one transversion (A = C) occurred as frequently as T = C and G = A transitions. The relative infrequency of the latter transition types is probably due to the low proportion of thymine and guanine nucleotides at third codon positions in cytochrome b, where most changes occurred. A distance analysis of the sequence data was performed, using distances calculated by the maximumlikelihood formula in DNADIST (PHYLIP; Felsenstein, 1995), with a ts:tv ratio of 5:1 and three rate categories of relative magnitude 1:0.3:30 and equal probability (0.333, 0.333, 0.333). The distance matrix was input to FITCH, and a weighted least-squares search was performed (i.e., P 5 2; Fitch and Margoliash, 1967) with global rearrangements and a jumbled input order of taxa (2 reps); P. falcinellus was designated as the outgroup. The resulting best-fit tree includes all clades having 49% or greater bootstrap support in the single most-parsimonious tree (Fig. 7), but the ordering of relationships among these clades conflicts with the maximum-likelihood (Fig. 6) and most-parsimonious (Fig. 7) trees. The maximum likelihood, parsimony, and distance analyses yielded optimal topologies that are generally congruent in relationships among species within genera, but incongruent with respect to basal relationships in the family. To assess congruence quantitatively, I used a test suggested by Kishino and Hasegawa (1989), available in the DNAML program in PHYLIP (Felsenstein, 1995). For a given data matrix and choice of model parameters, the Kishino-Hasegawa test compares the ln(likelihood) of the maximum-likelihood tree to that of any selected input tree(s). If the difference in ln(likelihood) is smaller than 1.96 standard errors, then the trees are not significantly different. According to this test, neither the maximum-parsimony tree (Fig. 7) nor the optimal tree based on maximum-likelihood distances is significantly different from the maximumlikelihood tree (Fig. 6). That is, given the cytochrome b data and the choice of model parameters (ts:tv of 5:1 and three substitution rate categories, with relative magnitude 1:0.3:30 and equal probabilities), none of the three topologies can be excluded as a possible phylogeny for the storks. A combinable-component consensus of these three is considered the best-supported estimate of relationships among species in the Ciconi-

idae based on the cytochrome b data (Bremer, 1990; Fig. 8). DISCUSSION Comparison of DNA–DNA Hybridization and Sequence-Based Trees Although the DNA–DNA hybridization data and cytochrome b sequences produce conflicting best-fit trees (Figs. 4, 6, and 7), the conflicting nodes are weakly supported. I used the Kishino-Hasegawa (1989) test to assess the congruence of the optimal DNA–DNA hybridization tree (Fig. 4) and the maximum-likelihood tree based on the cytochrome b data (Fig. 6). According to this test, the two trees are not significantly different; i.e., the fit of the the cytochrome b data to the optimal DNA–DNA hybridization tree is not significantly worse than the fit to the maximum-likelihood tree. The optimal DNA–DNA hybridization tree and the consensus tree based on different analyses of the cytochrome b data (Figs. 4 and 8, respectively) are entirely congruent and complementary. A combinable-component consensus (Bremer, 1990) of these two trees represents the best estimate of relationships within the Ciconiidae based on the two data sets. The combinablecomponent consensus is identical in topology to the best-fit DNA–DNA hybridization tree, except that the relationship of C. nigra to its congeners is unresolved in the consensus tree. This polytomy within Ciconia probably reflects a rapid divergence of the unresolved lineages, because the loss of resolution in the consensus tree results from the collapse of short branches in the

FIG. 8. Combinable-component consensus of three optimal trees based on the cytochrome b data matrix: maximum-likelihood (Fig. 6), maximum-parsimony (Fig. 7), and FITCH analysis of maximumlikelihood distances.

PHYLOGENY OF STORKS

optimal trees (Figs. 4, 6, and 7). Short branches indicate a small amount of divergence, which translates to a short time interval, assuming that divergence is proportional to time. The two molecular data sets are complementary in some respects. For example, the cytochrome b data resolve relationships among species in the genus Mycteria that are weakly resolved in the DNA–DNA hybridization tree. In general, data from single-copy nuclear DNA–DNA hybridization are uninformative about relationships among populations within species or between very closely related species, and the relatively small cytochrome b divergences within Mycteria (Table 4) suggest that this radiation is more recent than among species in other stork genera. The cytochrome b data, on the other hand, give poor resolution of more ancient divergences that are well resolved by the DNA–DNA hybridization data. For example, the cytochrome b data fail to pair the two Leptoptilos species. These two species are separated by a large sequence divergence for congeners, but their pairing is strongly supported in the best-fit DNA–DNA hybridization tree. Also, resolution of basal relationships among storks is strongly supported in the best-fit DNA–DNA hybridization tree (Fig. 4), but weakly resolved in the cytochrome b-based trees and incongruent between analyses (Figs. 6 and 7). Comparison to Previous Estimates of Stork Relationships Several stork species are included in the Sibley and Ahlquist (1990) phylogeny of the birds of the world based on DNA–DNA hybridization data. The species included and their relationships are as follows: (J. mycteria (L. crumeniferus (C. ciconia, (E. asiaticus, M. americana)))). This pattern of relationships conflicts with the best-fit DNA–DNA hybridization tree obtained in this study (Fig. 4). In the Sibley and Ahlquist (1990) tree, however, the nodes supporting relationships within the Ciconiidae are short, particularly the node supporting the Ephippiorhynchus–Mycteria pair. Furthermore, the Sibley and Ahlquist (1990) tree is not based on a complete data matrix. M. americana was the only stork species labeled, and the relationships within storks were estimated from only 11 hybrids. Uncertainty in the topology due to experimental error must be considered large. The placement of Jabiru as the sister to a clade including all the other storks could be due to incomplete hybridization, possibly due to a shortstranded DNA sample. The average distance (DT50H ) between Jabiru and the other storks is listed in the Sibley and Ahlquist (1990) tree as 4.4°C, a distance considerably larger than any that I obtained for comparisons of Jabiru to other storks (average DT50H 5 2.38°C, range 1.09–3.35°C). Sibley and Ahlquist (1990) also include a matrix of distances (DT50H ) and a corresponding FITCH tree for

295

five cathartid genera (Cathartes, Sarcoramphus, Vultur, Gymnogyps, and Coragyps), storks, ibises, herons, and cormorants. The distances are averages of all available measurements for species in these groups. The reported distance between Cathartes and storks is 9.23°C (n 5 9), and between ibises and storks is 9.86°C (n 5 12). From the current study, the average DT50H distance between C. aura and storks is 7.28°C (range 4.47–10.86°C, n 5 104) and between P. falcinellus and storks is 9.11°C (range 7.95–11.02°C, n 5 103). The best estimates of phylogeny from this study are highly congruent with previous assessments of relationships among stork species based on behavior and morphology (Verheyen, 1959; Kahl, 1971a, 1972d; Wood, 1983, 1984). In particular, the genera defined by Kahl (1971a, 1972d) based on behavioral traits are supported as monophyletic groups. The molecular data also support the pairing of the genera Anastomus and Mycteria, which compose Kahl’s tribe Mycteriini, and the pairing of J. mycteria with the Ephippiorhynchus species. Kahl (1971a) considered Jabiru to be most closely allied to Ephippiorhynchus based on behavioral traits, but most closely allied to Leptoptilos based on morphological features. Kahl’s tribe Leptoptilini, including Jabiru, Ephippiorhynchus, and Leptoptilos, does not appear as a monophyletic group in any of the molecular-based trees. A close relationship between Jabiru and the Ephippiorhynchus species was suggested by Verheyen (1959), who recommended that Jabiru be moved to the genus Ephippiorhynchus. Wood (1984) supported Verheyen’s recommendation, based on phenetic analyses of behavior and osteological measurements. Wood (1984) concluded that the Jabiru–Ephippiorhynchus clade is most closely related to the Ciconia clade and suggested that the former be moved into the Ciconiini, and the tribe Leptoptilini limited to the genus Leptoptilos. The pairing of Jabiru–Ephippiorhynchus and Ciconia is supported in the best-fit DNA–DNA hybridization tree with a bootstrap value of 84%, but is not supported by any of the cytochrome b-based trees. In general, relationships among major clades within the family are weakly resolved by the cytochrome b data. In his classification of the Ciconiidae, Kahl (1971a) arranged the list of species in the genus Ciconia to reflect his assessment of relationships: C. nigra, abdimii, episcopus, maguari, ciconia. Consistent with Kahl’s arrangement, the molecular data support the species pairs maguari–ciconia and abdimii–episcopus; these occur with strong support in all the optimal trees (Figs. 4, 6, and 7). C. nigra pairs with maguari–ciconia– boyciana in the maximum-likelihood tree (Fig. 6) and the best-fit cytochrome b-distance tree, but with abdimii–episcopus–stormi in the maximum parsimony (Fig. 7) and best-fit DNA–DNA hybridization trees (Fig.

296

BETH SLIKAS

4). In all four trees, the relationship of C. nigra to its congeners is weakly supported and is unresolved in the consensus tree (Fig. 8). Kahl considered C. nigra to be the most primitive (basal) member of the genus. Previous classifications paired C. nigra with C. ciconia (Peters, 1931; Verheyen, 1959), a pairing which is contradicted by the molecular data. In all cytochrome b-based trees, the pairings of C. ciconia/C. boyciana and C. episopus/C. stormi are strongly supported; C. boyciana and stormi were not included in the DNA–DNA hybridization study. The members of each of these pairs were once considered conspecifics and were only recently separated into distinct species (Hancock et al., 1992). Within the genus Mycteria, Kahl was unable to resolve relationships. Peters (1931) and Verheyen (1959) placed M. americana as the most basal member of the genus; Verheyen further arranged the three Old World species as (ibis, (leucocephala, cinerea)). The latter arrangement is supported in the four molecular trees, although all resolution within the genus is weakly supported in the best-fit DNA–DNA hybridization tree (Fig. 4). Analysis Options for Sequence Data Analysis of sequence data is challenging due to the numerous options available and the dearth of welldefined criteria for choosing among these options. This study provides an opportunity to test empirically different methods of analysis for the sequence data, because the DNA–DNA hybridization tree provides an independent estimate of phylogeny for evaluating the accuracy of trees constructed from the sequence data. Such a comparison assumes that the DNA–DNA hybridization tree is accurate. This assumption cannot be directly evaluated, but the DNA–DNA hybridization tree is concordant with heuristic assessments of relationships based on behavioral data (Kahl, 1971a, 1972d), and such congruence suggests accuracy (Sheldon and Bledsoe, 1993; Miyamoto and Fitch, 1995). In addition, bootstrap analysis indicates that most resolved nodes in the best-fit tree are robust with respect to measurement error (Fig. 4). Nodes in the DNA-hybridization tree supported by bootstrap percentages of #80% will be considered unresolved in following comparisons to the sequence-based trees. Analysis options that were explored for the sequence data include parsimony with differential weighting of character-state transformations or weighting of characters and distance analyses using alternate distance estimates. I used the same search parameters in all parsimony analyses: heuristic search with TBR branch-swapping, random addition sequence of taxa (10 replicates), and the MULPARS option in effect. Basic Fitch parsimony, with all characters weighted equally and all character-state transformations requiring a single step (Fitch, 1971), yielded six most-

parsimonious trees, none of which is completely congruent with the best-fit DNA–DNA hybridization tree (Fig. 4). The strict consensus of the six trees does not conflict with the DNA–DNA hybridization topology, but many nodes in the consensus are unresolved. The term ‘‘generalized parsimony’’ refers to alternatives to Fitch parsimony that incorporate a stepmatrix assigning different weights (number of steps) to different character-state transformations (Swofford et al., 1996). The weights are inversely related to the expected probability of occurrence of the corresponding transformations. For example, the observed predominance of transitions over transversions in mitochondrial sequences (Moritz et al., 1987) can be modeled by giving transversions a greater weight than transitions to reflect the lower frequency of occurrence of the former (e.g., Lanyon and Hall, 1994; Christidis et al., 1996; Yoder et al., 1996; Lara et al., 1996). Parsimony analyses were performed imposing stepmatrices with tv:ts weightings of 3:1, 4:1, 5:1, 10:1, and 20:1. All yielded a single, identical most-parsimonious tree (Fig. 7). This tree differs from the DNA–DNA hybridization topology only in the placement of the genus Leptoptilos, and the position of this genus is rather weakly supported in the parsimony trees (Fig. 7). Thus, analyses with transversions weighted over transitions by at least 3:1 yielded trees more congruent with the DNA–DNA hybridization tree than unordered Fitch parsimony. Furthermore, the topology of the most-parsimonious tree was insensitive to the exact weighting factor over a broad range (3:1 to 20:1). An alternative form of generalized parsimony is to use the data matrix itself (Wheeler, 1990) or a tree estimated from the data matrix (Williams and Fitch, 1989) as a source of character-transformation weights. The goal of these approaches is to provide a detailed model of sequence change tailored to the particular data set, a model more precise and accurate than simply weighting transversions over transitions (Fitch, 1993). An analysis using iterative weighting of character-state transformations was performed, using the MacClade software to construct stepmatrices from trees (Maddison and Maddison, 1992). First, the sequence data were analyzed with all characters unordered. A stepmatrix was then constructed based on the frequency of character-state changes (unambiguous changes only) on one of the resulting most-parsimonious trees, randomly selected. The assigned weights were proportional to the inverse of the frequency of a given transformation; the stepmatrix was not symmetrized. The data were then reanalyzed with this stepmatrix imposed. This procedure was performed three times. The first analysis with a stepmatrix imposed yielded a single most-parsimonious tree, differing from the DNA–DNA hybridization tree in two aspects: (1) C. nigra paired with L. crumeniferus rather than positioning within the Ciconia clade; and (2) Leptoptilos javani-

PHYLOGENY OF STORKS

cus positioned as the sister to the genus Mycteria (with L. crumeniferus/C. nigra as the sister group to L. javanicus/Mycteria). The next iteration yielded two most-parsimonious trees, both differing from the DNA– DNA hybridization tree in that the genus Ciconia paired with the genus Leptoptilos, and the Ephippiorhynchus/Jabiru clade paired with the genus Mycteria. The final iteration yielded a single most-parsimonious tree, differing in several aspects from the DNA–DNA hybridization topology, including a split in the genus Ciconia. No further iterations were attempted, although neither the stepmatrix nor the tree topology had converged to stability. Judged by congruence with the DNA–DNA hybridization tree, the trees obtained in this iterative stepmatrix approach were less accurate than the trees obtained with a simpler 5:1 tv:ts weighting. Further, accuracy of trees did not improve with successive iterations. A basic problem with estimating weights from the data or trees based on the data is the accuracy of the estimate. Weights derived from the data matrix do not take into account multiple substitutions at sites. Tree-based weights depend on several arbitrary choices, including selection of an initial tree and selecting among unambiguous changes only, or minimum, maximum, or average number of changes on the tree to calculate weights (Maddison and Maddison, 1992). The relatively poor performance of the iterative stepmatrix approach is probably due to inaccuracy of the imposed weights and the sensitivity of the analysis to the precise values of the weights in a full, asymmetric stepmatrix. By imposing a specific model of sequence evolution in the phylogenetic analysis, differential weighting of nucleotide changes provides a means of accounting for the different frequency of occurrence of particular transformations. An alternative approach is to omit or reduce the weighting of characters that are expected to have a high frequency of change, usually third codon positions in coding sequences. The rationale behind the latter approach is to eliminate or reduce the effects of characters that are likely to be homoplastic due to the accumulation of multiple substitutions (Swofford et al., 1996). Characters can also be weighted iteratively; beginning with an unweighted analysis, each character is subsequently weighted inversely to some measure of its homoplasy, such as consistency or retention index, on the tree from the previous analysis (Farris, 1969, 1989; Carpenter, 1988). A parsimony analysis of the stork sequence data was performed with third codon positions excluded, and transversions weighted five times over transitions. The analysis resulted in six most-parsimonious trees. The strict consensus does not support monophyly of the genus Ciconia, the pairing of J. mycteria with the Ephippiorhynchus species, or the pairing of M. cinerea and M. leucocephala, although the latter is strongly supported in all analyses of the cytochrome b data

297

based on all codon positions. Because most phylogenetically informative sites occur at third codon positions (82%), eliminating these positions greatly reduces the size of the data matrix, which apparently removes informative data. A second analysis was performed in which each site was weighted by the maximum value of its rescaled consistency index (Farris, 1989); the consistency indices were taken from the set of mostparsimonious trees produced from an analysis with all characters equally weighted and all transformations unordered. A single most-parsimonious tree resulted, which differed in several aspects from the DNA–DNA hybridization tree, including the splitting of the genus Ciconia. Judging from these trials with the Ciconiidae sequence data, differential weighting of characters based on their expected or estimated homoplasy is relatively ineffective in recovering accurate phylogenies, probably due to the fact that a large percentage of the informative data is discounted. Another method for obtaining a phylogeny from sequence data is to calculate distances based on the sequences and to analyze the distance matrix. Several distance estimates have been proposed for sequences, based on different models of nucleotide substitution (Swofford et al., 1996). The most frequently used distances include the following: (1) percentage sequence divergence; (2) Jukes-Cantor distance, based on a model of independent change at all sites, with equal rates of change among all bases and equal nucleotide frequencies (Jukes and Cantor, 1969); (3) Kimura distance, a modification of the Jukes-Cantor, which allows for different rates of substitution for transitions and transversions (Kimura, 1980); and Felsenstein’s maximumlikelihood distance, a special case of the Kimura distance, which incorporates unequal base frequencies (Kishino and Hasegawa, 1989; Felsenstein, 1995). Distances are intended to estimate time since divergence between two taxa. As such, the model-based distances (Jukes-Cantor, Kimura, Felsenstein’s maximum-likelihood) are preferable to sequence divergence, because the former accommodate the possibility of multiple substitutions at sites. Using the DNADIST program in PHYLIP, JukesCantor, Kimura, and maximum-likelihood distances were calculated from the Ciconiidae sequence data and analyzed using the weighted least-squares option in FITCH (Felsenstein, 1995). The Jukes-Cantor distance matrix yielded a tree differing substantially from the DNA–DNA hybridization tree. For example, the genus Ciconia was split into two nonsister clades, and J. mycteria paired with L. javanicus, rather than with the Ephippiorhynchus species. Kimura and maximumlikelihood distances were calculated using a 5:1 ts:tv ratio and (for the maximum-likelihood distances) three categories of rate substitution, with relative magnitude 1:0.3:30 and equal probability; these parameter values were used to obtain the maximum-likelihood tree from

298

BETH SLIKAS

the sequence data (Fig. 6). FITCH analysis of the Kimura and maximum-likelihood distances gave the same best-fit tree. The topology was a better match to the DNA–DNA hybridization tree: the Ciconia species grouped as a clade, and J. mycteria paired with the Ephippiorhynchus species. However, the two Leptoptilos species did not pair together, and the arrangement of the Ciconia, Ephippiorhynchus–Jabiru, and Mycteria clades was not congruent with that in the DNA– DNA hybridization tree. A neighbor-joining analysis (Saitou and Nei, 1987) of the maximum-likelihood distance matrix yielded the same topology as the FITCH analysis. Distance estimates that incorporate a transition/ transversion bias (Kimura and maximum-likelihood) produced trees more congruent with the DNA–DNA hybridization tree than Jukes-Cantor distances. This outcome parallels the results of character-based analyses, where analyses weighting transversions over transitions by at least 3:1 yielded trees more congruent with the DNA–DNA hybridization tree than unordered or less biased analyses (see above). For the trees based on Kimura or maximum-likelihood distances, conflict with the best-fit DNA–DNA hybridization tree is limited to weakly supported nodes. Thus, analyses of Kimura or maximum-likelihood distances produce trees qualitatively as accurate as comparable character-based analyses. Because distance-based analyses tend to be much faster than character-based, converting sequences to Kimura or maximum-likelihood distances and analyzing the resulting distance matrix might be preferable to parsimony or maximum-likelihood analysis of characters for very large data sets, where the character-based analyses would be practically limited to superficial heuristic searches. ACKNOWLEDGMENTS I thank the following people and institutions who provided material for this study: Academy of Natural Sciences (Philadelphia, PA); Atlanta Zoo (Atlanta, GA); Audubon Zoo (New Orleans, LA); Fort Wayne Children’s Zoo (Detroit, MI); Louisiana State University (Baton Rouge, LA); Miami MetroZoo (Miami, FL); Philadelphia Zoo (Philadelphia, PA); San Diego Zoo (San Diego, CA); San Diego Wild Animal Park (San Diego, CA); L. Christidis (Museum of Victoria, Australia); J. A. W. Kirsch (University of Wisconsin, Madison, WI); C. G. Sibley (Sonoma, CA). I thank F. H. Sheldon and F. B. Gill for intellectual and financial support throughout the project. The manuscript was improved by comments from P. Arctander, J. Garcı´aMoreno, F. B. Gill, F. H. Sheldon, and an anonymous reviewer. This work was funded by grants from the Frank M. Chapman Memorial Fund (AMNH), the AOU Research Awards program, and NSF Grants BSR 9020183 and 9207991 to F. H. Sheldon and F. B. Gill.

REFERENCES Arctander, P. (1995). Comparison of a mitochondrial gene and a corresponding nuclear pseudogene. Proc. R. Soc. London B 363: 13–19.

Arthur, R. R., and Straus, N. A. (1978). DNA sequence organization in the genome of the domestic chicken (Gallus domesticus). Can. J. Biochem. 56: 257–263. Avise, J. C., Nelson, W. S., and Sibley, C. G. (1994). DNA sequence support for a close phylogenetic relationship between some storks and New World vultures. Proc. Natl. Acad. Sci. USA 91: 5173– 5177. Bleiweiss, R., and Kirsch, J. A. W. (1993). Experimental analysis of variance for DNA hybridization. II. Precision. J. Mol. Evol. 37: 514–524. Bleiweiss, R., Kirsch, J. A. W., and Matheus, J. C. (1994). DNA hybridization evidence for subfamily structure among hummingbirds. Auk 111: 8–19. Bremer, K. (1990). Combinable component consensus. Cladistics 6: 369–372. Bremer, K. (1994). Branch support and tree stability. Cladistics 10: 295–304. Carpenter, J. M. (1988). Choosing among multiple equally parsimonious cladograms. Cladistics 4: 291–296. Cavalli-Sforza, L. L., and Edwards, W. F. (1967). Phylogenetic analysis: Models and estimation procedures. Evolution 21: 550– 570. Christidis, L., Leeton, P. R., and Westerman, M. (1996). Were bowerbirds part of the New Zealand fauna? Proc. Natl. Acad. Sci. USA 93: 3898–3901. Cicero, C., and Johnson, N. K. (1995). Speciation in sapsuckers (Sphyrapicus). III. Mitochondrial-DNA sequence divergence at the cytochrome b locus. Auk 112: 547–563. Coulter, M. C., Bryan, A. L., Jr., Young, D. P., Jr., Brouwer, K., Kahl, M. P., King, C. E., Kushlan, J. A., Luthin, C. S., and van Wieringen, V. (1991). ‘‘A Bibliography of Storks, Ibises, and Spoonbills,’’ Savannah River Ecology Laboratory, Aiken, SC. Desjardins, P., and Morais, R. (1990). Sequence and gene organization of the chicken mitochondrial genome. J. Mol. Biol. 212: 599–634. Dickerman, A. W. (1991). Among-run artifacts in DNA hybridization. Syst. Zool. 40: 494–499. Edwards, S. V., Arctander, P., and Wilson, A. C. (1991). Mitochondrial resolution of a deep branch in the genealogical tree for perching birds. Proc. Roy. Soc. London B 243: 99–107. Epplen, J. T., Leipoldt, M., Engel, W., and Schmidtke, J. (1978). DNA sequence organization in avian genomes. Chromosoma 69: 307– 321. Esposti, M. D., De Vries, S., Crimi, M., Ghelli, A., Patarnello, T., and Meyer, A. (1993). Mitochondrial cytochrome b: Evolution and structure of the protein. Biochim. Biophys. Acta 1143: 243–271. Farris, J. S. (1969). A successive approximations approach to character weighting. Syst. Zool. 18: 374–385. Farris, J. S. (1989). The retention index and the rescaled consistency index. Cladistics 5: 417–419. Feinberg, A. P., and Vogelstein, B. (1983). A technique for radiolabeling DNA restriction endonuclease fragments to high specific activity. Anal. Biochem. 132: 6–13. Felsenstein, J. (1985). Confidence limits on phylogenies: An approach using the bootstrap. Evolution 39: 783–791. Felsenstein, J. (1995). ‘‘PHYLIP (Phylogenetic Inference Package), version 3.5c,’’ Department of Genetics, University of Washington, Seattle, WA. Fitch, W. M. (1971). Toward defining the course of evolution: Minimal change for a specific tree topology. Syst. Zool. 20: 406–416. Fitch, W. M. (1993). The triangle inequality and character analysis. Mol. Biol. Evol. 10: 713–714. Fitch, W. M., and Margoliash, E. (1967). Construction of phylogenetic trees. Science 155: 279–284.

PHYLOGENY OF STORKS Genetics Computer Group. (1993). ‘‘Program Manual for the GCG Package, version 8,’’ Genetics Computer Group, Madison, WI. Goodman, M. D., Tagle, A., Fitch, D. H. A., Bailey, W., Czelusniak, J., Koop, B. F., Benson, P., and Slightom, J. L. (1990). Primate evolution at the DNA level and a classification of the hominoids. J. Mol. Evol. 30: 260–266. Graybeal, A. (1993). The phylogenetic utility of cytochrome b: Lessons from bufonid frogs. Mol. Phylogenet. Evol. 2: 256–269. Griffiths, C. S. (1994). Monophyly of the Falconiformes based on syringeal morphology. Auk 111: 787–805. Hackett, S. J. (1996). Molecular phylogenetics and biogeography of tanagers in the genus Ramphocelus (Aves). Mol. Phylogenet. Evol. 5: 368–382. Hancock, J. A., Kushlan, J., and Kahl, M. P. (1992). ‘‘Storks, Ibises, and Spoonbills of the World,’’ Academic Press, San Diego. Howard, H. (1942). A review of the American fossil storks. Carnegie Inst. Wash. Publ. 530: 187–203. Irwin, D. M., Kocher, T. D., and Wilson, A. C. (1991). Evolution of the cytochrome b gene of mammals. J. Mol. Evol. 32: 128–144. Jukes, T. H., and Cantor, C. R. (1969). Evolution of protein molecules. In ‘‘Mammalian Protein Metabolism’’ (H. N. Munro, Ed.), pp. 21–123, Academic Press, New York. Kahl, M. P. (1963). Thermoregulation in the wood stork, with special reference to the role of the legs. Physiol. Zool. 36: 141–151. Kahl, M. P. (1966). Comparative ethology of the Ciconiidae. Part 1. The Marabou Stork, Leptoptilos crumeniferus (Lesson). Behaviour 27: 76–106. Kahl, M. P. (1971a). Social behavior and taxonomic relationships of the storks. Living Bird 10: 151–170. Kahl, M. P. (1971b). Observations on the Jabiru and Maguari storks in Argentina, 1969. Condor 73: 220–229. Kahl, M. P. (1972a). Comparative ethology of the Ciconiidae. Part 2. The Adjutant storks, Leptoptilos dubius (Gmelin) and L. javanicus (Horsfield). Ardea 60: 97–111. Kahl, M. P. (1972b). Comparative ethology of the Ciconiidae. The wood storks (genera Mycteria and Ibis). Ibis 114: 15–29. Kahl, M. P. (1972c). Comparative ethology of the Ciconiidae. Part 4. The ‘‘typical’’ storks (genera Ciconia, Sphenorhynchus, Dissoura, and Euxenura). Z. Tierpsychol. 30: 225–252. Kahl, M. P. (1972d). A revision of the family Ciconiidae (Aves). J. Zool. (London) 167: 451–461. Kahl, M. P. (1972e). Comparative ethology of the Ciconiidae. Part 5. The openbill storks (genus Anastomus). J. Ornithol. 113: 121–137. Kahl, M. P. (1973). Comparative ethology of the Ciconiidae. Part 6. The blacknecked, saddlebill, and Jabiru storks (genera Xenorhynchus, Ephippiorhynchus, and Jabiru). Condor 75: 17–27. Kimura, M. (1980). A simple method for estimating evolutionary rate of base substitutions through comparative studies of nucleotide sequences. J. Mol. Evol. 16: 111–120. Kirsch, J. A. W., Springer, M. S., Krajewski, C., Archer, M., Aplin, K., and Dickerman, A. W. (1990). DNA–DNA hybridization studies of the carnivorous marsupials. I. The intergeneric relationships of bandicoots (Marsupialia: Perameloidea). J. Mol. Evol. 30: 434–448. Kishino, H., and Hasegawa, M. (1989). Evaluation of the maximumlikelihood estimate of the evolutionary tree topologies from DNA sequence data, and the branching order in Hominoidea. J. Mol. Evol. 29: 170–179. Kocher, T. D., Thomas, W. K., Meyer, A., Edwards, S. V., Paabo, S., Villablanca, S. X., and Wilson, A. C. (1989). Dynamics of mitochondrial DNA evolution in animals: Amplification and sequencing with conserved primers. Proc. Natl. Acad. Sci. USA 86: 6196–6200. Kornegay, J. R., Kocher, T. D., Williams, L. A., and Wilson, A. C. (1993). Pathways of lysozyme evolution inferred from the sequences of cytochrome b in birds. J. Mol. Evol. 37: 367–379.

299

Krajewski, C. (1989). Phylogenetic relationships among cranes (Gruiformes: Gruidae) based on DNA hybridization. Auk 106: 603–618. Krajewski, C., and Dickerman, A. W. (1990). Bootstrap analysis of phylogenetic trees derived from DNA hybridization distances. Syst. Zool. 39: 383–390. Kumar, S., Tamura, K., and Nei, M. (1993). ‘‘MEGA (Molecular Evolutionary Genetics Analysis), version 1.0,’’ Pennsylvania State Univ., University Park, PA. Lanyon, S. M. (1985). Detecting internal inconsistencies in distance data. Syst. Zool. 34: 397–403. Lanyon, S. M., and Hall, J. G. (1994). Reexamination of barbet monophyly using mtDNA sequence data. Auk 111: 389–397. Lara, M. C., Patton, J. L., and da Silva, M. N. F. (1996). The simultaneous diversification of South American echimyid rodents (Hystricognathi) based on complete cytochrome b sequences. Mol. Phylogenet. Evol. 5: 403–413. Ligon, J. D. (1967). Relationships of the cathartid vultures. Occas. Pap. Mus. Zool. U. Mich. No. 651. Maddison, W. P., Donoghue, M. J., and Maddison, D. R. (1984). Outgroup analysis and parsimony. Syst. Zool. 33: 83–103. Maddison, W. P., and Maddison, D. R. (1992). ‘‘MacClade, version 3,’’ Sinauer, Sunderland, MA. Maniatis, T., Fritsch, E. F., and Sambrook, J. (1982). ‘‘Molecular Cloning: A Laboratory Manual,’’ Cold Spring Harbor Publications, Cold Spring Harbor, NY. Miyamoto, M. M., and Fitch, W. M. (1995). Testing species phylogenies and phylogenetic methods with congruence. Syst. Biol. 44: 64–76. Moritz, C., Dowling, T. E., and Brown, W. M. (1987). Evolution of animal mitochondrial DNA: Relevance for population biology and systematics. Annu. Rev. Ecol. Syst. 18: 269–292. Moritz, C., Schneider, C. J., and Wake, D. B. (1992). Evolutionary relationships within the Ensatina eschscholtzii complex confirm the ring species interpretation. Syst. Biol. 41: 273–291. Nunn, G. B., and Cracraft, J. (1996). Phylogenetic relationships among the major lineages of the birds-of-paradise (Paradisaeidae) using mitochondrial DNA gene sequences. Mol. Phylogenet. Evol. 5: 445–459. Olson, S. L. (1979). Multiple origins of the Ciconiiformes. Proc. Colonial Waterbird Group 1978: 165–170. Olson, S. L. (1985). The fossil record of birds. In ‘‘Avian Biology’’ (D. S. Farner, J. R. King, and K. C. Parkes, Eds.), Vol. 8, pp. 79–238, Academic Press, New York. Penny, D., and Hendy, M. (1986). Estimating the reliability of evolutionary trees. Mol. Biol. Evol. 3: 403–417. Peters, J. L. (1931). ‘‘Checklist of the Birds of the World I,’’ Harvard Univ. Press, Cambridge, MA. Saiki, R. K., Gelfand, D. H., Stoffel, S., Scharf, S. J., Higuchi, R., Horn, G. T., Mullis, K. B., and Erlich, H. A. (1988). Primer-directed enzymatic amplification of DNA with a thermostable DNA polymerase. Science 239: 487–491. Saitou, N., and Nei, M. (1987). The neighbor-joining method: A new method for reconstructing phylogenetic trees. Mol. Biol. Evol. 4: 406–425. Sanger, F., Nicklen, S., and Coulson, A. R. (1977). DNA sequencing with chain-terminating inhibitors. Proc. Natl. Acad. Sci. USA 74: 5463–5467. Sarich, V. M., and Cronin, J. E. (1976). Molecular systematics of the primates. In ‘‘Molecular Anthropology’’ (M. Goodman and R. E. Tashian, Eds.), pp. 141–170, Plenum, New York. Sarich, V. M., Schmid, C. W., and Marks, J. (1989). DNA hybridization as a guide to phylogenies: A critical analysis. Cladistics 5: 3–32. Sarich, V. M., and Wilson, A. C. (1967). Immunological time scale for hominoid evolution. Science 158: 1200–1203.

300

BETH SLIKAS

Schmid, C. W., and Marks, J. (1990). DNA hybridization as a guide to phylogeny: Chemical and physical limits. J. Mol. Evol. 30: 237– 246. Seibold, I., and Helbig, A. J. (1995). Evolutionary history of New and Old World vultures inferred from nucleotide sequences of the mitochondrial cytochrome b gene. Phil. Trans. R. Soc. London B 350: 163–178. Sheldon, F. H. (1987a). Phylogeny of herons estimated from DNA– DNA hybridization data. Auk 104: 97–108. Sheldon, F. H. (1987b). Rates of single-copy DNA evolution in herons. Mol. Biol. Evol. 4: 56–69. Sheldon, F. H. (1994). Advances in the theory and practice of DNA-hybridization as a systematic method. In ‘‘Molecular Approaches to Ecology and Evolution’’ (B. Schierwater, B. Streit, G. Wagner, and R. DeSalle, Eds.), pp. 284–297, Birkhauser Verlag, Basel, Switzerland. Sheldon, F. H., and Bledsoe, A. H. (1989). Indexes to the reassociation and stability of solution DNA hybrids. J. Mol. Evol. 29: 328–343. Sheldon, F. H., and Bledsoe, A. H. (1993). Avian molecular systematics, 1970’s to 1990’s. Annu. Rev. Ecol. Syst. 24: 243–278. Sheldon, F. H., Slikas, B., Kinnarney, M., Gill, F. B., Zhao, E., and Silverin, B. (1992). DNA–DNA hybridization evidence of phylogenetic relationships among major lineages of Parus. Auk 109: 173–185. Sheldon, F. H., and Winkler, D. W. (1993). Intergeneric phylogenetic relationships of swallows measured by DNA–DNA hybridization. Auk 110: 798–824. Sibley, C. G., and Ahlquist, J. E. (1981). The phylogeny and relationships of the ratite birds as indicated by DNA–DNA hybridization. In ‘‘Evolution Today’’ (G. G. E. Scudder and J. L. Reveal, Eds.), pp. 301–335, Carnegie-Mellon Univ., Pittsburgh, PA. Sibley, C. G., and Ahlquist, J. E. (1990). ‘‘Phylogeny and Classification of Birds,’’ Yale Univ. Press, New Haven, CT. Smith, A. B. (1994). Rooting molecular trees: Problems and strategies. Biol. J. Linn. Soc. 51: 279–292. Smith, M. F., Thomas, W. K., and Patton, J. L. (1992). Mitochondrial DNA-like sequence in the nuclear genome of an akodontine rodent. Mol. Biol. Evol. 9: 204–215. Sokal, R. R., and Rohlf, F. J. (1981). ‘‘Biometry,’’ Freeman, New York. Springer, M. S., Davidson, E. H., and Britten, R. J. (1992). Calculation of sequence divergence from the thermal stability of DNA heteroduplexes. J. Mol. Evol. 34: 379–382. Springer, M. S., and Kirsch, J. A. W. (1991). DNA hybridization, the compression effect, and the radiation of diprotodontian marsupials. Syst. Zool. 40: 131–151.

Springer, M. S., Kirsch, J. A. W., Aplin, K., and Flannery, T. (1990). DNA hybridization, cladistics, and the phylogeny of phalangerid marsupials. J. Mol. Evol. 30: 298–311. Springer, M. S., and Krajewski, C. (1989a). Additive distances, rate variation, and the perfect-fit theorem. Syst. Zool. 38: 371–375. Springer, M. S., and Krajewski, C. (1989b). DNA hybridization in animal taxonomy: A critique from first principles. Quart. Rev. Biol. 64: 291–318. Swofford, D. L. (1991). When are phylogeny estimates from molecular and morphological data incongruent? In ‘‘Phylogenetic Analysis of DNA Sequences’’ (M. M. Miyamoto and J. Cracraft, Eds.), pp. 295– 333, Oxford Univ. Press, New York. Swofford, D. L. (1993). ‘‘PAUP (Phylogenetic Analysis Using Parsimony), version 3.0,’’ Illinois Natural History Survey, Champaign, IL. Swofford, D. L., and Olsen, G. J. (1990). Phylogeny reconstruction. In ‘‘Molecular Systematics’’ (D. M. Hillis and C. Moritz, Eds.), pp. 411–501, Sinauer, Sunderland, MA. Swofford, D. L., Olsen, G. J., Waddell, P. J., and Hillis, D. M. (1996). Phylogenetic inference. In ‘‘Molecular Systematics’’ (D. M. Hillis, C. Moritz, and B. K. Mable, Eds.), pp. 407–514, Sinauer, Sunderland, MA. Vanden Berge, J. C. (1970). A comparative study of the appendicular musculature of the order Ciconiiformes. Am. Midland Nat. 84: 289–364. Verheyen, R. (1959). Contribution a` l’anatomie et a la syste´matique de base des Ciconiiformes (Parker 1868). Bull. Inst. Roy. Sci. Nat. Belgique 35(24). Werman, S. D., Springer, M. S., and Britten, R. J. (1990). Nucleic acids I: DNA–DNA hybridization. In ‘‘Molecular Systematics’’ (D. M. Hillis and C. Moritz, Eds.), pp. 204–249, Sinauer, Sunderland, MA. Wheeler, W. C. (1990). Combinatorial weights in phylogenetic analysis: A statistical parsimony procedure. Cladistics 6: 269–275. Williams, P. L., and Fitch, W. M. (1989). Finding the minimal change in a given tree. In ‘‘The Hierarchy of Life’’ (B. Fernholm, K. Bremer, and H. Jo¨rnvall, Eds.), pp. 453–470, Elsevier, Amsterdam. Wood, D. S. (1983). Phenetic relationships within the Ciconiidae (Aves). Ann. Carnegie Mus. 52: 79–112. Wood, D. S. (1984). Concordance between classifications of the Ciconiidae based on behavioral and morphological data. J. Ornithol. 125: 25–37. Yoder, A. D., Cartmill, M., Ruvolo, M., Smith, K., and Vilgalys, R. (1996). Ancient single origin for Malagasy primates. Proc. Natl. Acad. Sci. USA 93: 5122–5126.