Geographical Structure of Genetic Variation in the Malagasy Scops-Owl Inferred from Mitochondrial Sequence Data

Abstract

Resumen

Introduction

Madagascar's peculiar and highly endemic fauna is the result of long isolation and can be explained by two different types of colonization events: vicariance when Madagascar was separated from the African landmass or over-ocean (“rafting”) dispersal from Africa and Eurasia. The origin, timing, and mode (vicariance versus rafting) of colonization by different animals are the foci of several research programs (Cibois et al. 2001, Vences et al. 2003, Yoder et al. 2003, Marks and Willard 2005, Fuchs et al. 2006). However, given the estimated age of the separation of Indo-Madagascar from Africa (ca. 165 mya; Coffin and Rabinowitz 1992, Storey et al. 1995), specifically with regard to the earliest known paleontological evidence of modern lineages of vertebrates, the majority of extant groups presumably arrived on Madagascar via over-ocean dispersal.

In contrast, phylogeographic studies associated with in situ differentiation on Madagascar, after colonization successfully took place, are uncommon and speciation patterns are still poorly understood for birds. Yamagishi et al. (2001) studied the species-level relationships among the Vangidae, an endemic radiation of Malagasy corvoids that probably dates back to the Oligocene (Fuchs et al. 2006). Marks and Willard (2005) recently addressed the geographic origin of the Madagascar Pygmy Kingfisher (Ispidina madagascariensis). They also suggest, using mitochondrial sequences from six individuals, that substantial genetic differentiation (1.4% divergence) may have occurred on Madagascar. Recently, a few studies have revealed the existence of cryptic vertebrate lineages on Madagascar (Yoder et al. 2005), corresponding, for example, to an east-west separation for mantelline frogs (Vences and Glaw 2002), apparent north-south differentiation among clades of mouse lemurs (Yoder et al. 2000), or altitudinal replacement for shrew tenrecs (Olson et al. 2004). However, avian genetic diversity below the species level has not been investigated to date for any Malagasy taxon using comprehensive geographic sampling.

The Malagasy Scops-Owl (Otus rutilus sensu lato) is a common endemic species across much of Madagascar. It inhabits a variety of ecosystems, from heavily degraded habitats to natural forests between 0 and 2050 m asl (Schulenberg 2003). As for many members of this genus, there is considerable variation in plumage coloration within and among populations (Rasmussen et al. 2000, Schulenberg 2003). Patterns of feather pigmentation in this taxon have been proposed to correlate with variation in vocalizations and geography (Benson et al. 1976, Rasmussen et al. 2000), suggesting that more than one species of scops-owl occurs on Madagascar. Rasmussen et al. (2000) resurrected the name madagascariensis for the form occurring in the west and proposed the following distinctions to separate the two endemic species: 1) rutilus sensu stricto (east Madagascar, hereafter rutilus) tends to have a deeper rufous and dark brown plumage, its song is a short flat hoot, and it is found in humid forest regions, which for the most part are in the east, and 2) madagascariensis has a more streaked gray or yellow-brown plumage, its song is a rapid series of two or more very short notes, and it occurs in the drier forested portions of the island, specifically the west. These differences in plumage and vocalization are also correlated with some subtle biometric characters: rutilus tends to have longer ear tufts, a shorter tail, and smaller claws than madagascariensis (Rasmussen et al. 2000). In some areas of the island, such as the extreme southeast, a zone of pronounced environmental shifts, morph replacement seems to occur across a distance of only 30 km.

However, the reliability of the characters proposed by Rasmussen et al. (2000) to distinguish these two taxa at the species level may be questionable. The type of discriminant plumage characters they proposed are known to show intraspecific variability, that are at least in part associated with abiotic factors. For example, Gloger's rule hypothesizes that populations in drier areas will be less heavily pigmented than populations inhabiting humid areas. Birds from drier areas generally have melanin granules primarily in the distal parts of their feathers to better protect feathers against abrasion (Barrowclough and Sibley 1980). The differences in plumage coloration they cited between rutilus and madagascariensis may thus only reflect adaptive phenotypic plasticity to local bioclimatic conditions.

Most of the song samples studied by Rasmussen et al. (2000) came from the extreme southwest and east. Only a few recordings were available from the central highlands (Schulenberg 2003), and one of these, a recording from Ihosy (south-central Madagascar; Fig. 1), is intermediate between the typical songs of rutilus and madagascariensis. Hence, it is not possible to exclude the possibility that vocal variations are clinal or associated with ecological factors. More precisely, six of the eight vocal parameters used by Rasmussen et al. (2000), for which statistically significant measurements were detected between rutilus and madagascariensis, involved frequencies. The humid forest form (rutilus) had lower frequencies for all parameters, suggesting that ecological factors related to habitat, rather than fixed interspecific differences, may be involved in geographical variation of scops-owl songs. Indeed, it has been shown for certain bird species that populations inhabiting dense forests use lower frequencies in their songs than populations of more open habitats (Waas 1988, Tubaro and Lijtmaer 2006).

Figure 1

Map of Madagascar indicating sampling localities for Malagasy Scops-Owls with an overlay of Cornet's (1974) bioclimatic zones.Otus rutilus (rufous morph) is found in humid and subhumid areas while O. madagascariensis (gray morph) is found in subarid and dry areas. Using the morphological and vocal characters outlined in Rasmussen et al. (2000), the Field Museum of Natural History (FMNH) specimen numbers in italics are referable to O. madagascariensis and those in standard script to O. rutilus. Numbers in parentheses indicate haplotype number.

Open in new tabDownload slide

Given that several of the characters used to recognize a second form of Otus on Madagascar may be ecologically linked, the use of independent data, such as nucleotide sequences, is needed to test if there are indeed two species on the island. To address this question, we sequenced two mitochondrial genes (ND2 and Cytochrome b) and investigated the phylogeographic structure of the Malagasy Scops-Owl. We recognize the potential for nonmonophyletic gene trees to be embedded within reproductively cohesive species due to incomplete lineage sorting (Pamilo and Nei 1988, Funk and Omland 2003) or recurrent gene flow (Slatkin and Maddison 1989). However, mtDNA haplotypes are expected to coalesce much more rapidly relative to nuclear markers due to their smaller effective population size (Moore 1995), and thus serve as relatively early indicators of reproductive isolation.

Methods

Laboratory Procedures

We obtained tissue samples from 34 Malagasy Scops-Owls throughout the distribution of this genus on the island (Fig. 1), including multiple representatives of the putative taxa Otus rutilus (Pucheran) (23 individuals) and O. madagascariensis (Grandidier) (11 individuals). Individuals were assigned to each species according to their geographical localities as proposed by Rasmussen et al. (2000).

Total DNA was extracted from frozen or alcohol-preserved tissues using a CetylTrimethylAmmonium Bromide-based protocol (Winnepenninckx et al. 1993) with an overnight Proteinase K (0.1 mg ml⁻¹) digestion. We sequenced portions of two mitochondrial genes (ND2 and Cytochrome b) that are routinely used in avian phylogeographical studies (Zink et al. 2002, Drovetski et al. 2004, Pavlova et al. 2005). ND2 and a portion of Cytochrome b were amplified with PCR using primer pairs L5219Met and H6313Trp (Sorenson et al. 1999) for ND2, and L14967 (5′-CATCCAACATCTCTGCTTGATGAAA-3′) with H15487 (5′-ATGAAGGGATGTTCTACTGGTTG-3′) and L15424 (5′-CGATTCTTCGCTTTACACTTCCTCC-3′) with H15916 (5′-ATGAAGGGATGTTCTACTGGTTG-3′) for Cytochrome b. Numbers refer to the 3′ position of the primers in the chicken (Gallus gallus) mitochondrial genome (Desjardins and Morais 1990). The amplification protocol was: 2 min at 94°C, followed by 36 cycles of 94°C for 40 sec, 54–56°C for 45 sec, 72°C for 40 sec, with a final extension at 72°C for 5 min. Three microliters of the amplification products were electrophoresed on a 1.5% agarose gel and revealed under UV light with ethidium bromide to check for correct fragment size and to control for the specificity of the amplifications.

PCR products were purified using the QiaQuick PCR Purification Kit (Qiagen, Holden, Germany). Cycle-sequencing reactions were performed using the CEQ Dye Terminator Cycle Sequencing kit (Beckman Coulter, Inc., Fullerton, California) or the Big Dye (Applied Biosystems, Inc., Foster City, California) terminator chemistries kit employing the same primers that we used for PCR amplifications. For some samples, two internal primers for the ND2 gene were specifically designated (L5709: 5-ATRGCYRTCTCCTCARCAGC-3′; H5808: 5′-TTCAGCCYAGRTGGGAGATGG-3′) for the cycle-sequencing reactions. Both DNA strands were sequenced for all taxa on automated sequencers (CEQ2000 DNA Analysis System, Beckman Coulter, Inc., or ABI3100, Applied Biosystems, Inc.).

Sequences were edited, assembled, and aligned using Sequencher 4.1 (Gene Codes Corporation, Ann Arbor, Michigan). Every substitution was checked on the electropherograms to control for accuracy of base-calls. No insertions, deletions, or stop-codons were detected in the reading frames of ND2 and Cytochrome b, suggesting that our sequences were of mitochondrial origin and not nuclear pseudogenes (Sorenson and Quinn 1998). Newly generated sequences have been deposited in Genbank under the accession numbers EF198290–EF198323 (ND2) and EF198256–EF198289 (Cytochrome b).

Data Analyses

Many of the underlying assumptions of traditional phylogenetic tree-building methods (fully bifurcating trees, complete lineage sorting) are violated in intraspecific studies (Posada and Crandall 2001). Therefore, we employed alternative analytical methods to resolve spatial patterns among closely related lineages and used TCS version 1.21 (Clement et al. 2000) to reconstruct the minimum spanning network among haplotypes. We used Arlequin version 2.0 (Schneider et al. 2000) to calculate the number of haplotypes, haplotype diversity (H), nucleotide diversity (π), Tajima's D (Tajima 1989) and Fu's tests (Fu 1997) of selective neutrality, and the mismatch distributions within the O. rutilus and O. madagascariensis genetic datasets. The empirical mismatch distributions were tested against a model of sudden expansion using the generalized nonlinear least-squares approach (Schneider and Excoffier 1999). The overall genetic structure of populations was investigated with an analysis of molecular variance (AMOVA) using Tamura and Nei's pairwise distances (Excoffier et al. 1992), as calculated in Arlequin version 2.0 (Schneider et al. 2000). We used MrModeltest2 version 2.2 (Nylander 2004) and Akaike's information criterion to select the model that best fit our data. We performed a Mantel test to investigate if there was a correlation between genetic and geographic distances. We used program MDIV (Nielsen and Wakeley 2001, Nielsen 2002) to simultaneously estimate the posterior distributions of the parameters theta (θ, effective population size scaled to the mutation rate), M (migration rate between populations since divergence), and t (time of population divergence) between putative O. rutilus and O. madagascariensis populations. A series of simulations were conducted using different analytical parameters and prior combinations: number of generations between 10 and 50 million with a burn-in period corresponding to one-tenth of the sampled generations, M_max between 5 and 10, and t_max between 1 and 15. A finite sites model (Hasegawa-Kishino-Yane, or HKY, model) was used in all MDIV analyses. Values for θ, M, and t were plotted, and the mode of the posterior distribution was considered the best estimate. We calculated a 95% Bayesian confidence interval (CI) for θ and M. Values for t were converted into years before present (BP) using the formula:

$$

where L is the sequence length (1449 bp) and μ is the mutation rate per site per generation (Brito 2005, Steeves et al. 2005). Confidence intervals for t were derived from the 95% CI of θ. We used a mutation rate (μ) of 7 × 10⁻⁸ substitution per site per year, assuming a one-year generation time, which is appropriate for members of this genus. We obtained this mutation rate from our ongoing study of the origin and timing of colonization of Indian Ocean Otus spp. (JF et al., unpubl. data). The calibration point we used was the divergence between O. pauliani (endemic to Grande Comore) and O. moheliensis (endemic to Mohéli), and corresponded to the oldest age that would have allowed the colonization of Grande Comore by O. pauliani (0.5 million years BP; Nougier et al. 1986). Based on our unpublished results, we only included taxa that were close to O. rutilus and that putatively diverged from O. rutilus during the last two million years. This decision was made to have an accurate estimate of the mutation rate during the period of 0–2 million years BP (Ho et al. 2005, Ho and Larson 2006). We assumed a molecular clock, i.e., a constant rate of molecular evolution, and used the GTR + γ model in program PAML 3.14 (Yang 2004) to derive the mutation rate.

Results

We analyzed 1449 bp of mtDNA (1036 bp from ND2 and 413 bp from Cytochrome b) obtained from 34 individuals, corresponding to positions 5246 to 6281 (ND2) and 15022 to 15434 (Cytochrome b) of the Gallus gallus mitochondrial genome sequence (Desjardins and Morais 1990; GenBank accession number X52392). Among the 1449 sites, 22 were variable (20 transitions and two transversions) and 10 were parsimony informative. Most of the variation was at third codon positions (15 transitions and one transversion), followed by first codon positions (three transitions and one transversion), and second codon positions (two transitions). These substitutions yielded five nonsynonymous amino acid changes (AA160: Ala-Thr, AA173: Ile-Thr, AA269: Val-Ile for ND2; AA356: Ala-Thr and AA357: His-Asn for Cytochrome b). Nineteen haplotypes were detected among the 34 sequenced individuals; four haplotypes were shared by more than one individual. Maximum sequence divergence was 0.6% (mean  =  0.24%) and was observed between haplotypes H1 and H11. The most common haplotype (H8) was shared by 10 individuals originating from eastern and western localities and different bioclimatic zones (Fig. 1). Conversely, 12 haplotypes were exclusive to individuals classified as O. rutilus and five to those assigned to O. madgascariensis. The 95% minimum spanning tree among all the haplotypes obtained with TCS version 1.21 is depicted in Figure 2. The minimum spanning tree indicated a low level of haplotype divergence with several haplotypes differing by only one substitution from the most common haplotype. Number of haplotypes, haplotype diversity (H), nucleotide diversity (π), Tajima's D, and Fu's statistics for each population are presented in Table 1. MrModeltest2 version 2.2 (Nylander 2004) selected the General Time Reversible (GTR) model incorporating a proportion of invariant sites (I) as the best-fit model to our dataset. Analysis of molecular variance was thus performed with the Tamura and Nei model, as it is the closest model to the GTR model implemented in Arlequin. Results of the AMOVA indicated significant partitioning of the genetic variability between O. rutilus and O. madagascariensis (F_1,32  =  57.2, P  =  0.006), with 17% of the variance explained by differences between the two supposed species. The mismatch distributions for O. madagascariensis and O. rutilus appeared bimodal (Fig. 3), although a model of sudden population expansion could not be statistically rejected (O. madagascariensis: test of goodness-of-fit, sum of squared deviation  =  0.04, P  =  0.30, Harpending's raggedness index  =  0.08, P  =  0.55; O. rutilus: test of goodness-of-fit, sum of squared deviation  =  0.02, P  =  0.85, Harpending's raggedness index  =  0.04, P  =  1.0). The Mantel test was not statistically significant (r  =  0.04, P  =  0.40), indicating that there was no correlation between genetic and geographic distances.

Figure 2

Unrooted haplotype network of Malagasy Scops-Owl obtained with program TCS version 1.21 (Clement et al. 2000). Area of circles is proportional to number of individuals with that haplotype. Colors refer to Otus rutilus (white) and O. madagascariensis (gray). The smallest circles represent unsampled or extinct haplotypes.

Open in new tabDownload slide

Table 1

Open in new tab

Number of haplotypes, haplotype diversity (H), nucleotide diversity (π), and Tajima's D and Fu's statistics of selective neutrality obtained for each taxon of Malagasy Scops-Owl (Otus rutilus and Otus madagascariensis) as well as for all samples combined. Otus rutilus refers to populations found in humid and subhumid forests of eastern Madagascar and Otus madagascariensis refers to populations inhabiting arid areas and dry forests from western Madagascar (see Fig. 1 for exact sampling localities). Fu's D and F and Tajima's D-values reveal that the mtDNA sequences conform to neutral expectations.

Table 1

Open in new tab

Number of haplotypes, haplotype diversity (H), nucleotide diversity (π), and Tajima's D and Fu's statistics of selective neutrality obtained for each taxon of Malagasy Scops-Owl (Otus rutilus and Otus madagascariensis) as well as for all samples combined. Otus rutilus refers to populations found in humid and subhumid forests of eastern Madagascar and Otus madagascariensis refers to populations inhabiting arid areas and dry forests from western Madagascar (see Fig. 1 for exact sampling localities). Fu's D and F and Tajima's D-values reveal that the mtDNA sequences conform to neutral expectations.

Figure 3

Mismatch distribution for the Otus rutilus and O. madagascariensis populations obtained with program Arlequin version 2.0 (Schneider et al. 2000). Gray bars represent observed data, lines represent simulated data used to test goodness-of-fit to the sudden population expansion model. A model of sudden population expansion could not be rejected for either taxon (O. rutilus: test of goodness-of-fit, sum of squared deviation  =  0.02, P  =  0.85; O. madagascariensis: test of goodness-of-fit, sum of squared deviation  =  0.04, P  =  0.30).

Open in new tabDownload slide

Estimates of parameters associated with the coalescent-based analyses (θ, M, and t) were 4.96 (95% CI  =  2.99–9.55), 1.18 (95% CI  =  0.46–9.55), and 0.33, respectively (Fig. 4). These estimates were not dependent on the assigned prior values for M and t. Our data suggest that populations associated with the names O. rutilus and O. madagascariensis exchange on average one individual (female) per generation and that the two populations started to diverge very recently, ca. 8070 years BP (95% CI  =  4860–15 540). The “time to most recent common ancestor” (TMRCA) estimate was 1.54, yielding a date of 36 300 years BP (95% CI  =  22 700–72 520).

Figure 4

Posterior distributions of θ (theta, population size scaled to mutation rate), M (migration rate between the two taxa after divergence), and t (population divergence time) for Malagasy Scops-Owls based on analysis of 1449 bp of mtDNA using program MDIV (Nielsen 2002). Two different simulations were conducted (10 and 50 million generations) and the result from the longer run is presented here. Priors were set as T_max  =  15 and M_max  =  10 and a Hasegawa-Kishino-Yane (HKY) model was used in all analyses. The posterior distribution of the migration rate indicates the occurrence of gene flow between populations assigned to Otus rutilus and O. madagascariensis.

Open in new tabDownload slide

Discussion

To test the two-species hypothesis and assess the phylogeographic structure of Otus on Madagascar, we sequenced two mitochondrial genes for 34 individuals, consisting of samples obtained across much of the distribution of these two putative species. Populations assigned to O. rutilus and O. madagascariensis, based on the designations of Rasmussen et al. (2000), were not supported as monophyletic units. The mean and maximum genetic divergences between haplotypes were low (0.24% and 0.60%, respectively), with the two most divergent haplotypes (H1 and H11) differing by only nine substitutions. This level of divergence is ten-fold lower than that observed between Otus on Madagascar and their closest relatives on western Indian Ocean islands (~5%–6%; JF et al., unpubl. data). The maximum amount of genetic divergence observed between haplotypes of Malagasy Otus is about half the level of that observed between several raptor sister-species (around 1.5%; Wink and Sauer-Gürth 2000). Yet, amount of genetic divergence generally is not a reliable indicator of species status of taxa; several studies have shown that only a few mutations on genes involved in premating or postmating isolation can have dramatic consequences for gene flow (Orr 1992, Turner et al. 2005). Despite the low level of overall genetic differentiation, the analysis of molecular variance of Malagasy Otus revealed that the two populations exhibited significant differences in haplotype frequencies. However, only a small part of the total genetic variance was accounted for by genetic divergence between eastern and western populations. Reciprocal monophyly of mitochondrial haplotypes has often been considered a major requirement to delimit species or subspecies boundaries (Zink 2004). Yet, several factors, such as time elapsed since population divergence, initial population size, and the stochastic process of haplotype fixation, can lead to situations in which two biological species do not show a pattern of gene tree monophyly. Indeed, coalescence theory predicts patterns of polyphyly at an early stage of speciation (Avise 1994, Funk and Omland 2003). The situation encountered here is similar to that of intermediate polyphyly and neotypy as defined by Omland et al. (2006). In the case of neotypy, only common ancestral haplotypes are shared and new haplotypes, created by mutations, are not common among populations. In the case presented here, two haplotypes (H5 and H8) were shared between the two supposed species. Although the most common haplotype (H8) occupies an internal position in the network, suggesting that it is a retained ancestral haplotype, one cannot exclude that the sharing of H5 between the two populations is the result of recent or ongoing gene flow. This latter hypothesis receives further support from our coalescent analyses, which suggest that the two populations exchange on average one individual (female) per generation. Interestingly, one of the two individuals with the H5 haplotype (FMNH 345654) originated from the humid forest of Marosohy, just east of the ecotone between the two habitats and a contact zone between the two putative taxa.

Based on evidence from subfossil deposits, the modern Malagasy dry forest biomes, particularly those in the southern quadrant of the island, have undergone significant climatic change during the Quaternary period. Evidence from subfossil pollen deposits demonstrates that approximately 3000 to 2000 years ago the region underwent considerable desiccation (Burney 1993), which included, for example, substantial shifts in vegetational formations and the drying up of freshwater lakes. These vicissitudes coincided with notable changes in the regional vertebrate fauna and extinction events (Goodman and Rakotondravony 1996, Goodman and Rakotozafy 1997, Godfrey et al. 1999, Goodman et al. 2006). The timing of these shifts in habitat types roughly coincides with the period of divergence in Malagasy Otus, as estimated by our coalescent-based analyses (8070 years BP, 95% CI  =  4860–15 540 BP). Consequently, a parapatric ecological speciation model may best explain the observed pattern for Malagasy Otus. Indeed, several empirical (Schneider et al. 1999, Via 2001, Ogden and Thorpe 2002, Spector 2002) and theoretical studies (Doebeli and Dieckmann 2003) have suggested that speciation in parapatry can occur with moderate levels of gene flow, if selection is sufficiently strong and if the traits under selection are linked to reproductive isolation or survival. Differences in song characteristics between the eastern and western populations of the Malagasy Scops-Owl may act as a premating barrier. However, it is not known to what extent an individual Otus can recognize a song as conspecific or not, and whether these differences efficiently prevent gene flow among populations. Moreover, further studies are needed to clarify the impact of habitat on selection for plumage coloration.

The overall genetic difference between O. rutilus and O. madagascariensis is very weak, and of a degree typically observed at the subspecies level for other bird species. More importantly, our coalescent analyses suggest that gene flow occurs between the two “taxa.” Therefore, we conclude that the two-species hypothesis is not supported by mtDNA data. Nevertheless, our study provides molecular evidence for differences in haplotype frequencies between western and eastern population of Otus, resulting in significant geographical structure of genetic variation on Madagascar. Further bioacoustic studies, overlaid on ecological parameters, as well as additional genetic studies, including nuclear markers and larger sample sizes, are clearly needed to fully evaluate the level of reproductive isolation in populations of Otus on Madagascar and to ascertain their systematic status. This research might focus on the extreme southeastern contact zone, where morph replacement appears to be abrupt (Rasmussen et al. 2000).

Acknowledgments

Laboratory work at Muséum National d'Histoire Naturelle was supported by the ‘Service Commun de Systématique Moléculaire,’ IFR CNRS 101, MNHN, and by the Plan Pluriformation ‘Etat et structure phylogénétique de la biodiversité actuelle et fossile.’ We thank Annie Tillier, Céline Bonillo, and Josie Lambourdière for their assistance in the laboratory and Philippe Gaubert and Sarah Samadi for discussions on haplotype networks. A number of the samples from Madagascar were collected by David E. Willard and Thomas Gnoske. Lucienne Wilmé kindly prepared Figure 1. On Madagascar, we gratefully acknowledge the assistance of Association pour le Gestion des Aires Protégées (ANGAP) and the Direction des Eaux et Forêts for issuing permits to conduct this research, and Daniel Rakotondravony and Olga Ramilijaona, Département de Biologie Animale at the Université d'Antananarivo, for help with certain administrative details.

Literature Cited

Author notes