2.2 DNA extraction and amplification

Genomic DNA was extracted from 2–4 mm³ of tissue sample using a Puregene® DNA isolation kit (Gentra Systems, Minneapolis, MN), according to the Animal tissue DNA isolation Protocol. Briefly, tissue was homogenized in a low molarity salt buffer and digested overnight with proteinase K. After the addition of a protein precipitation solution, the sample was centrifuged to remove cellular material from the supernatant. DNA was then precipitated at room temperature in 100% isopropanol with 1.5µl of 20 mg ml⁻¹ glycogen. We pelleted the DNA by centrifugation on high speed followed by washing in 70% ethanol. DNA was rehydrated in Puregene® DNA hydration solution and then stored at 4°C.

Aliquots of genomic DNA isolates were used as templates for PCR to amplify double-stranded DNA products from the seven gene fragments. Each PCR had a reaction volume of 25 µl and contained 1 µl of DNA stock regardless of stock concentration (diluted 1:10 in some cases), 2.5 µl 10X reaction buffer (100 mM Tris-HCl, 500 mM KCl, 15 mM MgCl₂, pH 8.3; Roche, Mannheim, Germany), 1.5 µl of 2 mM each deoxynucleotide triphosphates, 0.1 µl of 5 unit µl⁻¹ Taq polymerase (Roche, Mannheim, Germany), and 1.0 µl of each oligonucleotide primer, each at 10 µM concentration. Primers used for amplification and sequencing are given in Table 1. Some samples amplified weakly using standard PCR, and were reamplified with Ready-to-go™ PCR beads (Amersham Pharmacia Biotech Inc., Piscataway, NJ).

Each 16s, 12s and cytb PCR included an initial denaturation step at 94° C for 2 min, followed by 30 cycles of PCR, each cycle including denaturation at 94° C for 30 s, annealing at 49 or 60° C for 30 s, and extension at 72° C for 55 s. After a final extension step at 72° C for 2 min the PCR products were held at 4° C. Each nuclear gene PCR included an initial denaturation step at 95° C for 1 min, followed by 35 cycles of PCR, each cycle involving denaturation at 95° C for 30 s, annealing at 55° C for 1 min, and extension at 72° C for 90 s. The final extension step was at 72° C for 5 min. PCR reactions were performed using an MJ Research PTC-200 Peltier Thermal Cycler (MJ Research Inc., Watertown, Mass.).

PCR products were electrophoresed on 1% low melting agarose gels stained with ethidium bromide; bands were visualized under ultraviolet light then excised, melted at 70° C, and incubated with 1µl GELase™ Agarose Gel-Digesting Preparation (Epicentre, Madison, WI) for at least 3 hours at 45° C. PCR products were cycle-sequenced using an ABI PRISM® Big Dye™ Terminator Cycle Sequencing Ready Reaction Kit, version 3.0 or 3.1 (Applied Biosystems, Foster City, CA). Each reaction included 3.0 µl of buffer, 0.6 µl of primer, 1.0 µl of Big Dye™reagent, 0.7–2.0 µl of PCR product template DNA, and enough water to bring the total reaction volume to 10 µl. The cycling protocol was 25 cycles of denaturation at 96° C for 10 s, annealing at 50° C for 5 s, and extension at 60° C for 4 min. Sequences of both strands were generated on an ABI PRISM® 3100 or 3730 Genetic Analyzer (Applied Biosystems, Foster City, CA). All sequences have been deposited in GenBank with the following accession numbers: 12s EU601178–EU601227, 16s EU601228–EU601262, tmo4c4 EU601263–EU601306, Rag2 EU601307–EU601356, cytb EU601357–EU601405, otx1 EU601406–EU601455, bmp4 EU601456–EU601505, dlx2 EU601506–EU601554.

2.3 Sequence alignment

Sequences from all genes were trimmed to the size of the smallest fragment to minimize the amount of missing data in the data matrix. Sequences for tmo4c4, rag2, cytb, otx1, bmp4, and the exonic regions of dlx2 were aligned using Sequencher 4.2 (Gene Codes Corp). 12s and 16s sequences were aligned manually to a secondary structure for labrids (Westneat and Alfaro, 2005), using Se-Al v1.0a1 (Andrew Rambaut, University of Oxford). Gaps and ambiguously alignable regions were excluded from analysis. Se-Al was also used to align the intron sequences of dlx2; regions that could not be reliably aligned between taxa were also excluded. After excluding unalignable regions and trimming sequence ends, the character count for each gene was 12s (743), 16s (399), tmo4c4 (485), rag2 (715), cytb (477), otx1 (666), bmp4 (488) and dlx2 (514), for a total of 4487 nucleotide characters.

2.4 Phylogenetic analysis

Phylogenetic analyses were conducted using PAUP* (4.0b10) or Mr. Bayes 3.1 (Ronquist and Huelsenbeck, 2003) on a quad-processor Apple Macintosh G5 computer or a distributed computer network (Illinois BioGrid). Heuristic parsimony searches to find the most parsimonious tree(s) were performed using TBR branch-swapping. One thousand random addition replicates were used to minimize the chance of finding only locally optimal trees (Maddison, 1991). All sites were equally weighted and gaps treated as missing characters. We used nonparametric bootstrapping (Felsenstein, 1985) to measure support of clades with 1000 pseudoreplicates and branch-swapping with two random addition replicates per pseudoreplicate.

For the maximum likelihood analysis the TRN+Γ+I substitution model (rmat = 1.0000 5.4800 1.0000 1.0000 7.5325) with invariable sites (pinvar = 0.5094) and among-site rate heterogeneity (shape=0.5298) was selected using hierarchical likelihood ratio tests implemented in Modeltest 3.7 (Posada and Crandall, 1998). A heuristic search with 100 random sequence additions was used to find the most likely tree.

We used MrBayes 3.1 to perform 4 independent runs of a Markov Chain Monte Carlo analysis with 6 chains running for 10 million generations on the Illinois BioGrid distributed computer network through DePaul University. The eight gene partitions were assigned separate parameters with their own model of sequence evolution, determined using ModelTest 3.7. We saved trees every 1000 generations for each run and used default (flat) priors for the transition/transversion rate ratio (Beta 1,1), branch length (Exp 10), alpha parameter of the gamma distribution for rate heterogeneity (Uniform 0.05–50), proportion of invariant sites (uniform 0–1), base frequencies (Dirichlet 1,1,1,1) and tree topology parameters (uniform over all possible topologies). Visual inspection of tree length, log-likelihood, and model parameters revealed that all parameters appeared to reach stationarity within 500,000 generations. To help ensure that stationarity was reached, we discarded the first 2 million generations (2000 sampled trees) from each run as burn-in and used the remaining 8 million generations (8000 sampled trees for each of the 4 runs) in all subsequent analysis. A majority rule consensus tree calculated from the post burn-in trees was constructed and used to determine the posterior probabilities of clades.

3.1 Molecular characterization and sequence divergence

The data set consisted of 4903 nucleotides of which 4487 were alignable, 1496 of the alignable bases were variable and 1121 were parsimony-informative. Base composition across the data set was fairly uniform (24.5% A, 28.9% C, 24.4% G, and 22.2% T). Within the dlx2 intron, we found that Calotomus carolinus and C. spinidens share a 5 bp insertion, and Scarus flavipectoralis and S. hypselopterus share a 14 bp deletion.

Character optimizations of gene fragments onto the total evidence Bayesian topology inferred from the concatenated sequences (Fig. 3–Fig. 5) illustrate the total amounts of nucleotide substitution along various branches of the tree. Sequence divergence among the parrotfishes and labrid outgroups was generally high (Fig. 3–Fig. 5), with the exception of comparisons within Chlorurus and Scarus. Peak 12s divergence (Fig. 3A) was 16.4%, occurring between Nicholsina usta and Pseudodax mollucanus, and within parrotfishes it was 13.1%, between Nicholsina and Scarus frenatus. Peak 16s divergence was lower, 10.3% and 8.6% between scarines and outgroups and within scarines, respectively. Cytochrome b divergence (Fig. 3B) was highest (25.2%) between Centrolabrus exoletus and Chlorurus capistratoides, and was maximally 21.6% within parrotfishes. Resolution of the phylogeny among Chlorurus and Scarus at the species level was aided most by variability in cytb data set within the genera (Fig. 3B). Maximum rag2 divergence (Fig. 4A) was between Lachnolaimus maximus and Leptoscarus vaigiensis (20.9%) and within scarines was 10.6% between Cryptomus and Scarus coelestinus. Peak sequence divergence of tmo4c4 was 17.1% between Epibulus and Hipposcarus, and within parrotfishes was maximal (11.5%) between Cryptomus and Hipposcarus.

Open in a separate window

Figure 3

Phylogram of rates of molecular evolution in the (A) 12s and (B) cytochrome b data partitions illustrated on the Bayesian consensus tree of all concatenated sequences. Numbers on selected branches indicate that the branchlength is N times the length shown.

Open in a separate window

Figure 4

Rates of molecular evolution in the (A) rag2 and (B) bmp4 data partitions illustrated on the Bayesian consensus tree of all concatenated sequences. Numbers on selected branches indicate that the branchlength is N times the length shown.

Open in a separate window

Figure 5

Rates of molecular evolution in the (A) dlx2 exon and (B) dlx2 intron data partitions illustrated on the Bayesian consensus tree of all concatenated sequences. Numbers on selected branches indicate that the branchlength is N times the length shown.

Regulatory genes showed similar, though slightly lower, levels of sequence divergence, with a maximal p-distance for bmp4 of 15.4% (Fig. 4B), between Lachnolaimus maximus and Leptoscarus vaigiensis, and a intra-scarine maximal divergence of 7.5%. Otx1 divergence was 10.4% between Tautoga onitis and Nicholsina usta, and within parrotfishes showed a peak between Nicholsina and Leptoscarus of 8.6%. Maximal overall dlx2 divergence (Fig. 5) was 13.7%, between Centrolabrus and Chlorurus japonicus, and was 7.1% between the scarine taxa Cryptomus and Scarus coelestinus. Most of this divergence is intronic however (Fig. 5), as the intron had a peak of 40% divergence (Epibulus to Calotomus spinidens) and a within parrotfish maximum of 23.4%, whereas the exonic regions were just 5.2% divergent overall (Tautoga-Sparisoma) and peaked at 2.5% within parrotfishes. Regulatory genes showed sequence divergences of zero or nearly zero between species pairs of Chlorurus and/or Scarus, and this was especially true of the dlx2 exon (Fig. 5A).

3.2 Phylogeny of the parrotfishes

A diverse DNA sequence data set including both mitochondrial and nuclear loci, intron and exon regions, and protein-coding and ribosomal regions, provided excellent resolution for the topology of the scarine phylogeny (Fig. 1). Bayesian, parsimony and likelihood methods produced nearly identical topologies of relationship among the parrotfishes. Our sample of species showed that the genera Scarus and Chlorurus were monophyletic, and together they form a well-supported clade. The predominantly seagrass-associated genera (Calotomus, Cryptotomus, Leptoscarus, Nicholsina, and Sparisoma) are shown to be distinct from a reef-associated clade (Bolbometapon, Cetoscarus, Chlorurus, Hipposcarus, and Scarus). Three species pairs in the genus Chlorurus received strong support: C. sordidus + C. japanensis, C. bleekeri + C. bowersi, and C. capistratoides + C. microrhinos. The sister species to the rest of the Chlorurus included in the tree was C. oedema. Four multi-species clades within Scarus can be identified with high Bayesian posterior probability and greater than 80% maximum parsimony bootstrap support (Fig. 1, Scarus clades 1–4). Scarus clade 1 includes four Western Pacific species, including the sister pair S. flavipectoralis + S. hypselopterus. Scarus clade 2 is composed of 7 of the sampled species, with the pair S. globiceps + S. oviceps closely related to S. rivulatus and S. quoyi, the pair S. chameleon + S. festivus sister to them, and, S. spinus as the sister to the rest of clade 2. Clade 3 is a strongly supported subclade of the four Caribbean Scarus taxa included in the study (S. coelestinus, S. guacamaia, S. iseri, and S. taeniopterus), with S. ghobban positioned as tentative sister to the Caribbean clade. Scarus clade 4 is a group of 8 species in the tree, arranged in two subclades (Fig. 1), including mostly IndoPacific taxa and a single member that ranges to the eastern Pacific Ocean (S. rubroviolaceus).

3.3 Phylogenetic utility of regulatory genes

The partitioned analysis, in which the regulatory genes and non-regulatory genes were analyzed separately, show that the regulatory genes expressed a similar phylogenetic signal to the more commonly used mitochondrial and nuclear characters, with both regulatory and non-regulatory partitions capturing the topology of the combined data tree (compare Fig. 1 and Fig. 2). Both partitions display the same pattern of distinct coral and seagrass clades, of the monophyletic association of Chlorurus and Scarus, of the monophyly of Sparisoma, Scarus, and Chlorurus, and of the main Scarus sub-clades seen in the combined analyses. Strongly supported nodes in the two partitioned analyses (Fig. 2) are overwhelmingly in agreement with one another and with the combined analyses. A number of the primary nodes along the backbone of the tree are common between the combined and non-regulatory partitions, but are not strongly supported by the regulatory partition. Only one well-supported node in the non-regulatory partition and three nodes in the regulatory partition (Fig. 2) are in conflict with strongly supported nodes in the combined tree (Fig. 1).

3.4 Rates of molecular evolution: synonymous and non-synonymous change

Phylogenetic maps of gene substitutions (Fig. 3–Fig. 5) illustrate that substitutions have accumulated most along the root branch of the scarine and cheiline clades, least along Scarus and Chlorurus tip branches, and that the different genes vary in the degree to which they contribute to branch lengths at particular locations in the phylogeny. For the protein coding genes, these substitution patterns were broken down into the rates of synonymous and non-synonymous changes within the phylogeny (Table 2). All six coding genes showed considerably higher rates of synonymous than non-synonymous change, resulting in dN/dS ratios well below 0.5, and usually below 0.1. A single exception is the rate ratio for tmo4c4 at the reef clade of parrotfishes, which is greater than 1.0 (Table 2). The branches leading to the cheiline and scarine subfamily nodes had the highest rates of synonymous substitution, and rates tended to decrease at finer levels in the phylogeny. Many of the shorter branches among species had no non-synonymous changes, preventing the computation of rate ratios.

The dS and dN phylogenetic maps for cytb (Fig. 6), rag2 (Fig. 7) and bmp4 (Fig. 8), illustrate the various synonymous and non-synonymous nucleotide substitution patterns in different genetic markers. The phylogenetic information content of cytb was high for the species relationships among Chlorurus and Scarus, but that information is almost entirely due to silent substitutions at the codon level (Fig. 6A), with low branch lengths for changes that result in amino acid changes (Fig. 6B). For rag2 and bmp4 (Fig. 7, Fig. 8), dS is also much greater than dN, although the dN trees shows some nonsynonymous change among Scarus species, but low variability within Chlorurus. The dS and dN character map for tmo4c4 was similar to that of bmp4, whereas those for otx1 and dlx2 show low dS and mostly zero dN branch lengths.

Open in a separate window

Figure 6

Synonyomus (A) and non-synonymous (B) substitution rate trees for cytochrome b, calculated using the Bayesian consensus topology and the Muse-Gaut (1994) codon substitution model implemented in HyPhy (Kosokovsky Pond et al., 2005).

Open in a separate window

Figure 7

Synonyomus (A) and non-synonymous (B) substitution rate trees for rag2, calculated using the Bayesian consensus topology and the Muse-Gaut (1994) codon substitution model.

Open in a separate window

Figure 8

Synonyomus (A) and non-synonymous (B) substitution rate trees for bmp4, calculated using the Bayesian consensus topology and the Muse-Gaut (1994) codon substitution model.

Overall patterns of molecular evolution displayed by the developmental regulatory genes differed from one another, with bmp4 usually having a rate similar to that of cytb, rag2, and tmo4c4 in most regions of the tree, whereas otx1 and exonic dlx2 had much lower rates of both synonymous and non-synonymous changes. Non-synonymous changes in otx1 and dlx2 were rare, found only in the branch leading to the scarines. However, the dN/dS rate ratio of dlx2 (0.33) is the highest of all six genes along the branch leading to the parrotfishes, significantly higher than the rate ratio of the rest of the phylogeny (Table 2). Four branches that were tested had accelerated dN/dS rate ratios (cheilines for rag2, scarines for dlx2, seagrass clade for bmp4, and reef clade for tmo4c4) and two had significantly lower rate ratios than the rest of the phylogeny (scarines for bmp4 and tmo4c4).

4.1 Phylogeny of the parrotfishes

The higher-level pattern of parrotfish relationships resolved here agrees closely with that of Streelman et al. (2002) in regard to the split between seagrass and reef clades of scarines, and the ancestral splits of the genera Bolbometopon, Cetoscarus, and Hipposcarus at the base of the reef clade (Fig. 1). In addition, the present topology reinforces the morphological diagnosis of Bellwood (1994) regarding the genera Chlorurus and Scarus as being two distinct monophyletic lineages. The primary phylogenetic contribution of the present study is thus to advance our understanding of species relationships among the genera Chlorurus and Scarus. Most previous suggestions of relatedness among species have been based on color pattern and distribution information, so we summarize those characters that may be of interest in view of the current phylogenetic tree.

The earliest species to branch from the Chlorurus clade, C. oedema, has a relatively simple, dark color pattern and, like many Chlorurus, a somewhat restricted distribution, from Sri Lanka to the Philippines, and north to the Ryukyu Islands. The main Chlorurus clade sampled here, forming the sister-group to C. oedema, are more brightly and variably colored, and all have a cheekpatch of color in the terminal male coloration that C. oedema lacks. The species pair C. sordidus and C. japanensis share a disjunct color pattern on the caudal peduncle that abruptly changes from yellow (C. japanensis) or brownish (C. sordidus) to a light blue peduncle flecked with yellow. This species pair was suggested to be closely related by Choat and Randall (1986), and it also shares a pattern of the dorsal and anal fins in which the base and tips of the membranes are blue, with orange between, interrupted by a central blue band (also shared by C. capistratoides). C. sordidus is one of the most broadly distributed species of parrotfishes, ranging from the Red Sea south to Natal, South Africa and east to the Hawaiian islands, north to the Ryukyu Islands, and south to Perth, Australia, whereas C. japanensis is restricted to the region between the Rykyu islands and Australia.

The sister pair C. bleekeri + C. bowersi in our analysis share particularly distinctive cheek patches of green/white (C. bleekeri) to solid green/blue (C. bowersi), a double chin strap of blue, and nearly identical fin colorations. These two species have disjunct Western Pacific distributions, with C. bleekeri occurring in the Moluccas north to the Marshall Islands, and south to the Great Barrier Reef and Vanuatu, whereas C. bowersi occupies the region from Indonesia, to the Philippines, Ryukyu Islands, and Palau. Chlorurus troschelii is also likely a member of this clade (Choat and Randall, 1986). Finally, the sister-group C. capistratoides + C. microrhinos, is well supported by molecular data but few color pattern characters confirm this relationship in an obvious way. Cetoscarus microrhinos is thought to be part of a species complex of C. gibbus in the Red Sea, C. strongylocephalus in the Indian Ocean and C. microrhinos in the west-central Pacific (Myers, 1999), so inclusion of all members of this group in future work may clarify the evolutionary history of this subclade.

The first clade to diverge from the Scarus radiation is Scarus clade 1 (Fig. 1), represented in this study by four relatively widespread parrotfish species. Scarus psittacus is the sister to the rest of this Scarus clade in our tree, and is one of the most broadly distributed species among all parrotfishes, occurring from the Red Sea to South Africa and east to the Hawaiian, Marquesan, and Tuamoto islands. The other three species, S. schlegeli, S. flavipectoralis and S. hypselopterus are primarily Western Pacific with a distribution center around the Philippines, although S. schlegeli reaches the eastern Indian Ocean. The species pair S. flavipectoralis and S. hypselopterus share a TP color pattern in which the anterior half (slightly more in S. hypselopterus) of the body is orange and the caudal half is blue.

Scarus spinus is the sister species to the rest of Scarus clade 2 in our phylogeny, a relatively rare species distributed in the Western Pacific and eastern Indian Ocean. The species pair S. chameleon + S. festivus is next up the backbone of clade 2, two species with largely overlapping distributions from the Ryukyus to Lord Howe in the Pacific and west to the Indian Ocean, where S. festivus reaches to the African coast. S. spinus, S. chameleon and S. festivus share a characteristic deep head profile with a moderate or pronounced bulbous lump above the eye and a short snout that is particularly parrot-like. Choat and Randall (1986) identified the close relationship between S. chameleon and S. festivus based on meristics and color pattern. The four species remaining in Scarus clade 2 (S. quoyi, S. rivulatus, S. globiceps and S. oviceps) all have dark green-blue terminal males with complex facial lines and patterns, and the biogeography of this clade is intriguing because all four species have widespread distributions ranging from Indian Ocean to Western Pacific. The primary phases of S. globiceps and S. rivulatus are difficult to distinguish from one another. S. globiceps differs in that it is the only one with 3 scale rows and 5 or more scales in the ventral row. Previous workers (Choat and Randall, 1986) have remarked on the resemblance between S. oviceps and S. dimidiatus (here placed in Scarus clade 4) in both IP and TP color phases, but if our phylogeny is correct this would represent color pattern convergence.

The third major Scarus clade recovered in our analysis is a mostly Caribbean clade of parrotfishes, but is anchored by an extremely widespread species, S. ghobban, that does not itself range into the Caribbean. S. ghobban extends from the Red Sea and South Africa, north to southern Japan, south to Perth, New South Wales, and across to the Eastern Pacific, ranging from the Gulf of California to Ecuador. Due to its widespread distribution, its occupation of non-reef and deeper water habitats (Randall, 1986), and its apparent sister-group position to the Caribbean lineage, the population biogeography of S. ghobban would be of considerable interest. Scarus taeniopterus and S. iseri are closely related and form the next branches up the topology of clade 3. These two species share a striped IP color pattern and similar TP coloration, as well as almost identical distributions in the Western Atlantic and Caribbean. The species pair S. coelestinus + S. guacamaia are large, often blotchily colored parrotfishes that also share Western Atlantic distributions that extend far down the coast of South America.

Scarus clade 4 is composed of eight Indo-West Pacific taxa that often have distinctive color patterns in both the initial and terminal phases. The clade formed by S. tricolor + (S. forsteni + S. rubroviolaceus) are a case in point, in which the initial phases of all three species are highlighted by a rosy orange/red on tail and anal fin or across the body and fins, contrasted by dark scales in black, blue or green colors at midbody. The terminal phase of S. forsteni and S. rubroviolaceus share a characteristic dark blue saddle anteriorly, with a paler green body posteriorly. Although to our knowledge S. rubroviolaceus has not been proposed to be related to S. forsteni or S. tricolor, Choat and Randall (1986) reviewed the taxonomic confusion between S. forsteni and S. tricolor, indicative of their close relationship. Scarus dimidiatus and S. scaber (not sampled here) are thought to be closely related (Choat and Randall, 1986), and we recovered the species pair S. frenatus + S. dimidiatus. This is biogeographically suggestive as S. frenatus occurs from the Red Sea to Western Australia and throughout the Western Pacific, whereas its apparent close relatives S. dimidiatus and S. scaber are split geographically, with S. scaber in the Indian Ocean and S. dimidiatus restricted to Indonesia, Samoa, north to the Ryukyu Islands, and south to the Great Barrier Reef. To round out our clade 4, Scarus niger is the sister to the species pair S. altipinnis + S. prasiognathos, three species with relatively deep bodies and similarly shaped caudal fins with a nearly lunate profile in large males. S. niger and S. prasiognathos are widely distributed Indo-Pacific taxa, whereas S. altipinnis is restricted to the Western Pacific. S. altipinnis and S. prasiognathos are both large species that have similar initial phase color patterns, and the TP males have a distinctive blue-green patch on the chin and throat. Choat and Randall (1986) suggested that these two species together with S. falcipinnis of the Western Indian Ocean form a species complex.

4.2 Evolutionary trends in mitochondrial, nuclear, and nuclear regulatory genes

The higher level splits among labrid subgroups and the distinctions between parrotfish genera were supported by each of the 8 genes used here, including regulatory genes. However, the strength of phylogenetic signal at the species level was less balanced, with cytochrome b, a gene often used in population genetic studies in fishes (Carmona et al., 2000; Feulner et al., 2007), providing half of the sites informative for parsimony within the 24 Scarus taxa. At the species level, the “standard” phylogenetic markers (12s, 16s, tmo4c4, and rag2) and regulatory gene fragments (bmp4, dlx2, and otx1), showed similar resolving power and similar phylogenetic structure (Fig. 2). However, regulatory partitions showed a wide range of variability in substitution rates, with bmp4 contributing much more information at the species level than the otx1 or dlx2 exon partitions (Fig 3–Fig 4). This result is similar to comparisons of bmp4, otx1, and other regulatory genes in African cichlid fishes (Terai et al., 2002) in which bmp4 showed the highest rates of both synonymous and non-synonymous change.

The utility of the regulatory genes for phylogenetic questions differs among hierarchical levels, precisely the pattern of information content that phylogeneticists may be looking for when addressing large-scale and finer scale questions in their analyses. Bmp4 showed dS and dN levels at many nodes of the tree similar to that of tmo4c4 (Table 2), a gene fragment which has been popular for family-level phylogenetic studies of fishes (e.g., Near et al., 2004, Sparks and Smith, 2004; Westneat and Alfaro, 2005) since its characterization by Streelman and Karl (1997). We conclude that bmp4 sequence represents a useful additional nuclear gene that will be informative in studies of family-level systematics in fishes. In contrast, otx1 accumulates mutations much more slowly than bmp4 and so may be more valuable in studies with a greater time depth, for example among families within the Perciformes, a problem gaining increased attention with a wide array of genetic and genomic data (Miya et al., 2003; Smith and Craig, 2007). Similarly, the two dlx2 exon fragments that we sampled evolve slowly, no doubt in part because almost half of our nucleotides are from sequence of the homeodomain. The intervening intron sequence, however, displays considerable variation among parrotfishes (Fig. 5B), but alignment with the other labrid outgroups was problematic. This segment is perhaps best suited to phylogenetic investigation among closely related species or perhaps even populations, but it remains to been seen if other taxa will display such significant divergence in this region.

Phylogenetic analysis of substitution patterns in protein coding genes can reveal regions of both the genes and the phylogenies in which an increase in nonsynonymous change serves as a signal for selection (Guindon et al., 2004; Nei, 2005; Seo et al., 2004). All of the genes examined here showed higher rates of synonymous than non-synonymous change, resulting in dN/dS ratios well below 0.5, and usually below 0.1 (Table 2), suggesting that a process of purifying or stabilizing selection has acted on these molecules. A single exception is the rate ratio for tmo4c4 at the reef clade of parrotfishes, which is greater than 1.0 (Table 2), indicating that there was a burst of nonsynonymous change along this branch that may be indicative of positive selection for protein evolution in this gene. Despite an overall low substitution rate for regulatory genes, it is suggestive that dS and dN were highly variable on the phylogeny, and that the branches leading to the cheiline and scarine subfamily nodes had the highest rates of nucleotide substitution, including both dS and dN. For example, the bmp4 dS and dN trees (Fig. 8) illustrate this variability, with long branch lengths to the cheiline and scarine clades. Also, it is noteworthy that the dN/dS rate ratio of dlx2 leading to the parrotfishes (0.33) is the highest of all six genes along this part of the phylogeny. An elevated level of variation in regulatory genes has often been shown to be associated with morphological diversity in plants (Aagard et al., 2006; Barrier et al., 2001). Recent exploration of skull function of labrid fishes in the context of a phylogeny (Westneat et al., 2005) has shown that the cheiline and scarine labrids are a hot-spot of structural and functional evolution in these fishes, and dlx2 is known to regulate formation of pharyngeal elements of the feeding apparatus during development (Qiu et al., 1997; Stock et al., 1996). The discovery of elevated dN/dS ratios at major nodes of the parrotfishes for dlx2 and bmp4 (Table 2) suggests that these phylogenetic branches may be an intriguing place to search for changes in developmental patterns in the skull and pharyngeal arches. If the present study can be broadened to include multiple clades of fishes with differential levels of morphological and mechanical divergence, we may discover whether evolutionary changes in functional morphology are correlated with higher rates of change in the developmentally important signaling molecules produced by regulatory genes.

Field work for specimen collection was aided by the Lizard Island Research Station, Bureau of Fisheries and Aquatic Research of the Philippines, Mark McGrouther, Jeff Leis, Randy Mooi, Tom Trinski, Jeff Williams, Kent Carpenter, Luz Regis, Jeff Janovetz, Brad Wright, Aaron Rice, Jim Cooper, and others on expeditions to the northern Great Barrier Reef, Vanuatu, Solomon Islands, and the Philippines. DNA sequencing was performed in the Pritzker Laboratory of Molecular Evolution at the Field Museum of Natural History. Thanks to Rob Vogelbacher and David Angulo for access to the Illinois BioGrid through DePaul University. This research was primarily funded by NSF 0235307 to M. Westneat and an NSF/Alfred P. Sloan Foundation Postdoctoral Fellowship (NSF DBI 9803946) grant to J. T Streelman, with support during analysis and manuscript preparation by NSF DEB 0445453 to M. Alfaro and by John D. and Catherine T. MacArthur Foundation support of the Encyclopedia of Life.