Demographic history of wood stork (Mycteria americana) Brazilian Pantanal colonies revealed by mitochondrial DNA analysis

We used mitochondrial DNA (mtDNA) sequences to investigate the demographic history of the wood stork (Mycteria americana) populations in the Brazilian Pantanal. Sequences of 390/460 bp fragment of the mtDNA control region were analyzed in 62 wood stork specimens from 8 colonies using neutrality tests, phylogeographic, and coalescent analyses. Population expansion was supported by the significantly negative values of Tajimas (D = -2.071) and Fus (Fs = -14.544) statistics and the unimodal pattern of mismatch distribution. Nested clade analyses indicated a historic range expansion event and recurrent gene flow that was restricted by isolation by distance as explanations for the haplotype distribution among the sampled colonies. High genetic diversity and the strictly unidirectional gene flow pattern emphasized the conservation importance of preserving the southern Pantanal colonies. Coalescence analyses suggested that northern and southern colonies diverged approximately 6,250 years before the present (YBP), and that their most recent common ancestor was approximately 18,900 YBP. Our results suggest that the contemporary wood stork Pantanal population originated from a more geographically limited founder population. Potential source populations may have occurred in the southern Pantanal or ancestry may reside in populations inhabiting the Brazilian central plateau or areas closer to the equatorial region.


Introduction
Combinations of historical and contemporary processes determine the genetic structure and geographical distribution of genetic diversity among populations (Templeton et al., 1995). Genetic markers can be used to track population demographic responses to climatic shifts and other historical processes (Lessa et al., 2003;Templeton et al., 1995) while phylogeographic studies allow inferences to be made about the history of population divergence based on associations between the geographical distribution of mitochondrial DNA haplotypes and their genealogical relationships (Avise, 2000).
The wood stork (Mycteria americana, Ciconiidae, Ciconiiformes) is a colonial nesting wading bird adapted to the fluctuating water levels found in wetlands. Its breeding range extends from the southeastern United States (US) to northern Argentina, its distribution being is almost contiguous where suitable habitat exists (American Ornithologists' Union, 1998). Movement patterns of wood stork have been linked to the abundance and availability of food (Bryan Jr., 2001) and in the freshwater wetlands of southern Florida and the Georgia and South Carolina coastal marshes these birds have been observed throughout the year. Like other species in the family Ciconiidae that are threatened by habitat loss, the wood stork US population is endangered because of its sensitivity to environmental disturbances such as anthropogenic water level manipulations (Coulter et al., 1999;Wetlands International, 2002). In the US, wood stork has been listed as an endangered species after declining from an estimated population of 60,000 in the 1930s to about 5,000 in 1978 (USFWS, 1984;Ogden et al., 1987). The species also occurs in the Brazilian Pantanal, a region characterized by wet/dry cycles which result in large fluctuations in the availability of prey that limits the presence of wood stork during the dry season.
Waterbird colonies are established annually at the same sites and regional genetic structuring may be determined by evaluating the distribution of genetic diversity among colonies (Friesen, 1997). Previous studies investigated the spatial patterns of variation and genetic differentiation among wood stork colonies using nuclear DNA markers (Stangel et al., 1990;Van Den Bussche et al., 1999;Del Lama et al., 2002;Rocha et al., 2004). Allozyme (Stangel et al., 1990), microsatellite, and DNA fingerprint-ing analyses (Van Den Bussche et al., 1999) have indicated low levels of genetic differentiation among wood stork colonies in the US, and the authors speculated that these patterns were the result of recent colonization and high gene flow. The genetic structure of wood stork colonies in the Brazilian Pantanal has been assessed using the same allozyme and microsatellite loci described in the US population. No differentiation was detected and the survey indicated high gene flow among members of neighboring colonies (Del Lama et al., 2002;Rocha et al., 2004).
Mitochondrial DNA (mtDNA) is an appropriate genetic marker for evaluating phylogeographic patterns and the recent evolutionary history of species (Avise, 2000). The control region is a hypervariable, non-coding portion of the mitochondrial genome that is often used in studies of genetic population structuring, phylogenies of closely related taxa and phylogeography (Avise, 2000;Ruokonen and Kvist, 2002).
In the study described in this paper we hoped to shed light on the natural history and evolution of wood stork colonies in the Brazilian Pantanal using mtDNA control region sequences, neutrality tests, phylogeographic analyses and coalescent theory to investigate the historical processes affecting the genetic structure of this population. Additionally, the results presented can be used as a basis for future phylogeographic studies of other wood stork populations.

Study area and sampling
The 150,000 km 2 Pantanal plain region of centralwestern Brazil is one of the world's largest freshwater wetlands which was formed approximately 2.5 million years ago between the upper Pliocene and lower Pleistocene by the westward migration and compression of the Brazilian shield by the uplift of the Central Andes (Ussami et al., 1999). Topographically, the Pantanal is from 80 to 200 m above sea level with a low inclination gradient of about 1-2 cm km -1 from north to south and 6 to 8 cm km -1 from east to west. The climate is characterized by a wet season during October to March (spring through summer in the southern hemisphere) and a dry season from April to September (autumn through winter) and it is an important avian migratory stopover and wintering ground and the most significant South American waterbird breeding area during the dry season (Alho et al., 1988).
For our study we sampled eight historically recognized wood stork colonies from the Brazilian Pantanal over three breeding seasons (spring through summer) during 1997, 1999 and 2000 (Table 1, Figure 1). Approximately 0.5 ml of blood was collected from wood stork nestlings in each of 62 different nests, using EDTA as anticoagulant. Only non-related nestlings were included in the population genetic analyses.

Genetic analysis
Total genomic DNA was extracted from the red cells of wood stork nestlings sampled using the method of Sambrook et al. (1989) and the concentration of DNA estimated by electrophoresis in 1% (w/v) agarose gel.
All the wood stork control region, part of the NADH dehydrogenase subunit 6 (ND6) gene, the complete tRNA-Glu and tRNA-Phe gene sequences, and the partial sequence of the 12S ribosomal RNA gene were amplified and sequenced using two sets of primers, the H1251 and L16525 primers (Sorenson et al., 1999) and two primers developed in this study, LCBS-1 (5'-TTTAGTGAATGC TTGTCGGA-3') and HCBS-1 (5'-TCCGACAAGCATT CACTAAAT-3'). Based on the wood stork sequences a new set of primers, the previously tested LCBS-1 primer and the new primer HMYC (5'-GGTCATTGGGGGAA GAAGCTC-3'), were used to amplify a fragment of 390 or 460 base pairs (bp) of the third domain mitochondrial control region. 242 Lopes et al.  Polymerase chain reaction (PCR) amplification was performed in a final reaction volume of 20 µL containing approximately 320 ng of genomic DNA, 0.5 mM of each primer, 0.03% (v/v) Tween, 0.01 M Tris-HCl (pH 8.3), 1.5 mM MgCl 2 , 0.03% (v/v) DMSO, 0.15 mM dNTP and 1 unit of Taq polymerase (Biotools). The touchdown program was: 5 min at 94°C; five cycles of 20 s at 94°C, 20 s at 60°C (annealing, with a decrease of one degree per cycle) and 20 s at 72°C; 15 cycles of 20 s at 94°C, 20 s at 55°C and 30 s at 72°C; 20 cycles of 20 s at 94°C, 20 s at 52°C, and 30 s at 72°C; and 10 min at 72°C. One microlitre of the PCR product was eletrophoresed on a 1% ethidium bromide stained agarose gel to assess the PCR product quality. The PCR products with more than 60 ng µL -1 were purified by incubating 5 µL of PCR product with 1 unit of shrimp alkaline phosphatase (SAP; Pharmacia) and 10 units of exonuclease I (Exo I; Pharmacia) at 37°C for one hour followed by 10 min at 80°C (Werle et al., 1994).
Sequencing reactions were carried out using the Big-Dye Terminator Cycle Sequencing kit (Perkin Elmer) as per the manufacturer's instructions, along with approximately 400 ng of purified DNA template and 3 pmol of LCBS-1 primer in an ABI 377 automated sequencer. For accuracy, the light strand was sequenced twice per sample because the primer used to sequence the heavy strand did not produce a satisfactory sequence.

Statistical analyses
Sequences were aligned by MULTALIN (Corpet, 1998) in the DNA 5-0 option and visually verified. Haplotypes had four or five 70 bp repeats, which required gaps to align sequences, these gaps being treated as missing data. Positions with missing data were excluded from the variability analysis but not from the nested clade analyses.
Haplotype diversity (h) and nucleotide diversity (π) (Nei, 1987) were calculated with the ARLEQUIN v. 2.0 software package (Schneider et al., 2000) and used to describe the amount of genetic variation in a population. Theta (θ) values based on the observed number of segregating sites were calculated using the DnaSP program v. 4 (Rozas et al., 2003). Deviations from selective neutrality were assessed using Tajima's D (Tajima, 1989a) and Fu's Fs statistic (Fu, 1997) in DnaSP v. 4. Values of Tajima's D and Fu's Fs have also been used as indicators of a recent change in population sizes because they are sensitive to departures from demographic equilibrium (Tajima, 1989a, b;Aris-Brosou and Excoffier, 1996;Fu, 1997). Mismatch distributions were calculated for northern (MI, TU, PF, FI, and BG), southern (BB, RV, and FR) and all Pantanal colonies as an additional test for demographic expansion (Rogers and Harpending, 1992). The fit between the observed and expected distributions was evaluated by the sum of square deviations (SSD; Schneider and Excoffier, 1999). Total control region sequences were analyzed using the BLAST program (McGinnis and Madden, 2004).
Population substructure was investigated using analysis of molecular variance (AMOVA; Excoffier et al., 1992). AMOVA estimates a conventional F-statistic based on analysis of variance of gene frequencies. However, it takes the number of mutations between haplotypes into consideration and performs a hierarchical analysis of variance to partition variation into components due to intra-and inter-individual differences and inter-population differences (Schneider et al., 2000).
Using phylogenetic information, the nested clade analysis (Templeton et al., 1995) tests the phylogeographic association of haplotypes while discriminating between the effects of recurrent but restricted gene flow (particularly gene flow restricted by isolation by distance; Wright, 1943) vs. non-recurrent historical events (e.g., range expansion, colonization, or past fragmentation) (Templeton et al., 1992;Templeton et al., 1995;Templeton, 1998). Aligned DNA sequences were inputted into the TCS program v. 1.13 (Clement et al., 2000) to generate a nested haplotype network using the procedures described in Templeton et al. (1987). The GEODIS program  was used to test for significant association between the resulting nested clade design and the geographical distances between sample areas. When the null hypothesis of`no geographical association' was rejected, an inference key (Templeton, 2004) was used to interpret the outcome of the geographical association analysis.
To estimate gene flow between colonies located in the north and south Pantanal we used the MIGRATE program v. 0.9.10 (Beerli, 2000) that calculates N ef m (N ef = female effective population size, m = migration rate) maximumlikelihood estimates based on the coalescent theory. Estimation of migration rates based on coalescent theory allows detection of asymmetries in migration rates and differences in population sizes (Beerli, 1998;Beerli and Felsenstein, 1999). The GENETREE program (Griffiths and Tavaré, 1994) produces likelihood estimates for DNA sequence data under the infinite sites mutation model and this was used to estimate the likelihood of migration matrix parameters between northern and southern Pantanal populations.
The MDIV program (Nielsen and Wakeley, 2001) was used to estimate the scaled coalescent time (T), using T = t/N ef , where t is the time in years between northern and southern Pantanal colonies' divergence and the most recent common ancestor to the Pantanal population. The parameter N ef was estimated from N ef = θ/2µ, where µ is the mutation rate per sequence per generation. The program was run under the HKY finite sites model assuming the possibility of multiple hits, differences in nucleotide frequencies and the presence of transition/transversion bias (Nielsen and Wakeley, 2001). The divergence time was calculated using 389 bp (number of base pairs not considering the gap region), the generation time of 4 years (Coulter et al., 1999), the theta value from MDIV (4.5) and a mutation rate of 10% per million years. The conventional avian mtDNA clock suggests 2% sequence evolution between a pair of lineages per million years, but for control region sequences a ten-times faster rate has been employed for calibration (Quinn, 1992;Baker and Marshall, 1997). The estimation of the most recent common ancestor in years was obtained by multiplying the most recent common ancestor coalescent units by N ef (Harris and Hey, 1999).

Results
Three amplification patterns were observed, including a 390 bp pattern, a 460 bp homoplasmic pattern and a 390/460 bp heteroplasmic pattern. Heteroplasmic specimens, characterized by the presence of different mtDNA copies in their genome, were excluded from the population analyses and three family samples collected from siblings in the same nest were used to investigate the heritage of these patterns. The observed heteroplasmic pattern was explained by an insertion/deletion (indel) of a 70 bp repeat that results in simultaneous amplification of 390 and 460 bp fragments in the same individual. The fragment amplification of three heteroplasmic individuals and their siblings confirmed the maternal transmission of the heteroplasmic pattern (data not shown).
When BLAST analyzed, total control region sequences produced significant alignments solely with other avian mtDNA control region sequences but not with sequences from other taxa. Twenty-two variable nucleotide positions were observed in the 62 wood stork specimens investigated. Nineteen haplotypes were identified, excluding variations detected inside the gap region (Table 2) which contained three more mutation points: a G insertion at nucleotide position 87 and two substitutions at positions 116 and 144. Among 19 haplotypes, 15 were observed in only one colony (private or endemic haplotypes), three were shared by two or three colonies, and haplotype 1 was observed in 38 specimens distributed among all colonies. The majority of the private haplotypes corresponded to tips (i.e. haplotypes that have only a single mutational connection to the other haplotypes within a cladogram) in the nested clade analysis cladograms.
Haplotype diversity averaged 0.624 (±0.073) and ranged from 1.000 (RV) to 0.464 (TU, BG, and FR) ( Table  3). The estimated haplotype diversity in northern colonies was 0.568 (±0.095, 38 sequences, 11 haplotypes) and 0.757 (±0.094, 24 sequences 12 haplotypes) in southern colonies. Nucleotide diversity for each colony ranged from 0.010 (MI) to 0.002 (TU, BG and FR) and the average nucleotide diversity value for the total sample was 0.075. The estimated average theta based on the number of segregating sites was 0.012 and the value for each colony ranged from 0.010 (RV) to 0.003 (TU, PF, BG, and FR).
The neutrality tests applied to the wood stork Pantanal population showed significantly large negative 244 Lopes et al.  values, indicating population expansion (Kvist el al., 1998;Moum and Árnason, 2001;Lessa et al., 2003;Uimaniemi et al., 2003;Bowie et al., 2004;Pearce et al., 2004;Gottelli et al., 2004; and others). Tajima's D did not differ significantly from zero for any single colony, but it was significantly negative when all colonies were combined (D = -2.071, p < 0.05) and when only southern colonies were considered (D = -2.098, p < 0.05). Significant negative values of Fu's Fs statistics were detected in the RV colony (Fs = -5.899, p < 0.01), northern (Fs = -4.608, p < 0.01) and southern (Fs = -8.616, p = 0.00) regions and when all Pantanal colonies were considered (Fs = -14.544, p = 0.00) ( Table 3). The southern, northern, and all Pantanal colonies mismatch distribution analyses showed a unimodal lefthand peak, which supports a population expansion model (data not shown). According to Rogers and Harpending (1992), long-term stable or slowly declining populations usually generate a multimodal mismatch distribution while a unimodal peak on the left side of the graph indicates a very recent expanding population. The sums of squared deviations of mismatch distribution under the population expansion model were significant for the southern colonies (SSD = 0.117, p = 0.04) and marginally significant for the northern colonies (SSD = 0.126, p = 0.05), but no differences resulted from an evaluation of all Pantanal colonies (SSD = 0.008, p = 0.75).
The AMOVA showed no significant genetic differentiation among colonies (global Fst = 0.04, p > 0.05). Because the polymorphism in the gap region was included in the nested clade analysis, the haplotype 15 was subdivided into haplotype 15 and 15' resulting in 20 haplotypes. Since it was possible to insert haplotype 19 at two different points of the network two possible cladograms were generated ( Figure 2). Results of the nested contingency analysis are showing significant geographical association in the total cladogram for network A (χ 2 = 13.148, p < 0.05). In both cladograms, nested clade analysis indicated seven instances in which the null hypothesis of high gene flow levels (panmixia) was rejected. Interpretations of significant results for both cladograms are given in Table 4, following the inference key (Templeton, 2004). The two cladograms generated by the different position of haplotype 19 produced basically the same nested clade analysis results, so only one of the analyses is shown.
To elucidate the historical processes that shaped the evolution of the wood stork Pantanal population we used nested clade analysis. In this type of analysis there are three reasons why the null hypothesis of no association between haplotypes and geographical distribution may not be rejected, i.e. inadequate sample sizes, panmixia or insufficient genetic variation (Templeton et al., 1995). Our analysis detected continuous range expansion in clade 1-2 composed of the most commonly sampled haplotype (haplotype 1), which was present in all eight colonies, with five tip haplotypes that were present in four colonies (BG, BB, RV, and FR). The results suggest restricted gene flow with isolation by distance (IBD) affecting clade 2-3 but our sampling design did not allow discrimination between fragmentation and isolation by distance in clade 1-7. The only difference between the two cladograms is that shown by the results for clade 2-4, for which cladogram B suggests restricted gene flow with isolation by distance. However, we used the conservative results in cladogram A, which failed to distinguish between fragmentation and isolation by distance. To summarize, nested clade analysis suggested that geographical association of haplotypes in the wood stork Pantanal wood stork demography 245  (Figure 3) were the product of historical range expansion combined with gene flow restricted by isolation by distance and/or fragmentation. Inadequate geographical sampling and inconclusive outcome (due to absence of tip/interior status) results are included in Table 4. Migration estimates between colonies located in northern and southern Pantanal were N ef m = 81.6 from south to north and N ef m = 0.00 from north to south according to the maximum-likelihood method, using the MIGRATE program. These values were similar to those es-timated using the GENETREE program (i.e. N ef m = 33.0 from south to north and N ef m = 0.01 from north to south). The value of the scaled coalescent time for southern and northern Pantanal wood stork colonies was 0.432 and a divergence time estimate from coalescence analysis was approximately 6,250 years before present (YBP). The most recent common ancestor estimated was 1.307 coalescent units for the total sample and the estimated time in years based on the most recent common ancestor coalescent units was 18,900 YBP.

Discussion
Intra-individual sequence length variation in the mtDNA control region has been described in several bird families, including Laridae (Berg et al., 1995;Kidd and Friesen, 1998;Crochet and Desmarais, 2000), Scolopacidae (Berg et al., 1995) and Paridae . Although mtDNA is predominantly transmitted through maternal lines to offspring in most species, including birds, 246 Lopes et al.   Geographical sampling inadequate to discriminate between fragmentation and isolation by distance heteroplasmy can be generated by relatively common somatic mutations in the oocytes. This can occur within an individual or by paternal leakage of mitochondria during egg fertilization . The mitochondrial control region is a highly variable genetic marker, frequently showing insertion/deletion (indel) mutations that enable detection of heteroplasmy. Our results from wood stork siblings suggest that the heteroplasmic pattern was matrilineally transmitted. To our knowledge, this is the first description of a heteroplasmic mtDNA pattern in the family Ciconiidae.
This study also revealed a pattern of gene flow restricted by distance, which has not been detected with nuclear data. Previous studies of the genetic structure of the wood stork Pantanal breeding colonies using nuclear genetic markers showed no genetic differentiation and a high level of gene flow among breeding colonies (Del Lama et al., 2002;Rocha et al., 2004) which led investigators to conclude that all the wood stork Pantanal colonies belong to a single panmictic population and that distance did not represent a significant barrier to gene flow among colonies in the same region. Our AMOVA showed no differentiation among wood stork colonies in the Pantanal. While our nested clade analysis results also showed recurrent gene flow among wood stork Pantanal colonies, they also indicate distance-associated restrictions to gene flow among colonies located in the Pantanal area (clade 2-3). Our results showed a single widely distributed haplotype accompanied by other lower-frequency localized haplotypes, which is a pattern indicative of either gene flow among colonies or recent demographic expansion. Expansion in the Pantanal wood stork population is further evidenced through significant negative Tajima's D and Fu's Fs values for the total cladogram, the unimodal mismatch distribution pattern and the nested clade analysis (clade 1-2). Our mtDNA analyses suggest that the genetic homogeneity previously detected may be indicative of both a recent population expansion in the Pantanal area and gene flow among neighbor colonies.
The nested clade analysis did not discriminate between population fragmentation and isolation by distance in two clades (1-7 and 2-4), including haplotypes from colonies on the extremities of the sampling area ( Figure 3). The classification of the wood stork Pantanal colonies into northern and southern colonies was based on their geographical location, but there were colonies between these areas that were not sampled. According to Templeton (1998), gaps in sampling prevent the distinction between gene flow restricted by isolation by distance (a recurrent event) and fragmentation (an historical event). Based on this information, we expect that sampling geographically intermediate colonies would confirm gene flow restricted by isolation by distance. These results should help direct future sampling in the Pantanal, as well as in other wood stork breeding areas.
Results showing a strictly unidirectional migration from southern to northern Pantanal may express an historical expansion that occurred during colonization of the area by wood stork. The elevated wood stork haplotype diversity in the southern Pantanal colonies and the apparent south-to-north gene flow suggest that founders from the southern Pantanal region expanded their range northward. The topography of the Pantanal enables southern regions to remain flooded longer during the wet season, so if rainfall in the Pantanal was lower in the past, wood stork may have colonized the southern lowlands first because they were flooded while northern areas remained dry.
Our coalescence data suggest a recent divergence time of 6,250 YBP between southern and northern colonies and an older value of 18,900 YBP for the most recent common ancestor. Differences between these estimates are to be expected since ancestor divergence is initiated before the population splits (Edwards and Beerli, 2000;Arbogast et al., 2002). The most recent common ancestor is an estimate of the time at which all the haplotypes in the current population coalesce and indicates the degree to which the sample sequences are related. The most recent common ancestor estimates (18,900 YBP) suggests that the colonization of the Pantanal occurred at the end of the Late Pleistocene. However, the ancestral origin of the current wood stork Pantanal population remains unclear because the climatic conditions in the Pantanal at that time are unknown. Temperature and water level may have influenced distribution and breeding activity of wood stork in South America during the Late Pleistocene because these birds are tactile feeders that depend on open shallow wetlands for foraging. Climatic instability has been shown to influence the behavior of wood stork with, nest abandonment having been observed after heavy rainfall and during cold weather in US colonies (Kahl, 1964;Coulter and Bryan Jr., 1995;Mitchel, 1999). Reproduction has also failed when the Pantanal received more rainfall than usual during an El Niño event in 1998 (personal observation by I.F. Lopes). Furthermore, anthropogenic hydropattern changes in the Florida Everglades altered the location and timing of wood stork breeding colonies (Ogden, 1991;Mitchel, 1999;Coulter et al., 1999).
Despite the well-documented patterns in the northern hemisphere, there is debate regarding the changes in climate and vegetation that occurred in tropical areas during the Pleistocene (two million to twelve thousand YBP). Authors generally agree that during the Pleistocene the tropics were characterized by average temperatures which were 4 to 12°C lower than at present but questions remain about rainfall (Ab'Saber, 1977;Prance, 1982;Clapperton, 1993;Stute et al., 1995;Colinvaux et al., 2000;Assine and Soares, 2004). Some authors asserting that the climate was drier during the Late Pleistocene (Damuth and Fairbridge, 1970;Ab'Saber, 1977;Clapperton, 1993) while others maintain that it was wetter (de Oliveira, 1992;Ledru et al., 1996;Ferraz-Vicentini and Salgado-Laboriau, 1996;Colinvaux et al., 2000). The climatic conditions in the Pantanal during the Late Pleistocene have also been debated, with some authors proposing that at this time the Pantanal consisted of sand dunes and that the climate was arid (Ab'Saber, 1977;Clapperton, 1993), while others suggest a semi-arid climate with occasional heavy rains (Boggiani and Coimbra, 1995) or conclude that the data from the Pantanal are too debatable to allow any clear inference (Colinvaux et al., 2000). However, studies in areas neighboring the Pantanal indicate that during the Late Pleistocene the climate was moister than previously suggested (De Oliveira, 1992;Ledru et al., 1996;Ferraz-Vicentini and Salgado-Laboriau, 1996;Assine and Soares, 2004). The formation of the Pantanal wetland has been attributed to a change to a more humid climate with the individualization of isolated lakes occurring at the Pleistocene/Holocene transition between 10,200 and 5,190 YBP (Bezerra, 1999;Assine and Soares, 2004). Other studies assign the establishment of the present vegetation of the Brazilian southeast and central-west regions to about the same period, i.e. about 10,000 to 6,000 YBP (De Oliveira, 1992;Barbieri et al., 2000;Behling, 2002;Assine and Soares, 2004).
Based on the evidence suggesting wetter climatic conditions we hypothesize that the wood stork Pantanal breeding population was not established before the end of the Pleistocene and that when colonization occurred it is more probable that it was in the southern area which is topographically more likely to have provided wetland habitats. This hypothesis is supported by the southern ancestry and south to north migration pattern indicated by our data. The time estimated for the most recent common ancestor supports wood stork colonization of the Pantanal during the peak of the last northern hemisphere glaciation event about 18,000 YBP (Pielou, 1991). However, the most recent common ancestor may not necessarily coincide with the time that these wood stork populations arrived in the Pantanal. Alternatively, the ancient wood stork haplotype may have originated elsewhere and colonization of the Pantanal may have occurred during the Holocene following changes in the Pantanal landscape, increased humidity and a warmer climate between 10,200 and 5,190 YBP (Bezerra, 1999;Assine and Soares, 2004). If this is the case, the ancient haplotype may have originated in the Brazilian central plateau (Barbieri et al., 2000;Ferraz-Vicentini and Salgado-Labouriau, 1996) or areas closer to the equator that remained climatically stable (Colinvaux et al., 2000). Expansion of this analysis using samples from regions neighboring the Pantanal and from other regions in South America will be helpful in elucidating the recent colonization and expansion pattern suggested here.
In conclusion, this first approach using mtDNA to analyze the genetic structure of the wood stork Pantanal population supports the existence of contemporary gene flow among colonies and also shows demographic patterns not previously revealed by nuclear markers. Our results suggest that the current wood stork Pantanal population may have originated from a more geographically limited founder population that expanded recently. Distance as a limiting factor to gene flow in wood stork is also a pattern not previously detected with nuclear markers. We also found that wood stork in southern Pantanal colonies show higher genetic variability than those in the northern colonies, which may be an important aspect to consider when establishing conservation guidelines. The most recent common ancestor estimate is compatible with the establishment of wood stork breeding colonies in the Pantanal after the last glaciation maximum, and the same estimates based on nuclear genes should be able to confirm this. Mitochondrial DNA haplotypes, nested clade analysis and coalescence analyses are useful tools for revealing demographic historical processes affecting Pantanal colonies and can be applied to future studies of wood stork breeding populations in other regions. Additionally, our methodology may be an important tool for investigating the comparative phylogeographic pattern of co-distributed waterbird species and testing predictions of the climatic influence on Pantanal biota.