Demographic history of the fragmented yellowthroated bulbul (Pycnonotus xantholaemus) population in the Deccan Peninsula, India

The yellow-throated bulbul (YTB) is an endemic passerine restricted to scrub forests along hill slopes with exposed rocky outcrops in the Deccan Peninsula, India. It is found in small, discontinuous populations and is vulnerable to extinction due to ongoing habitat loss and subsequent population decline. To assess the genetic connectivity and past demography, we sequenced 1050 nucleotide base pairs of the mitochondrial control region of 60 individuals that represent distinct populations in the geographic range of the species. We recovered 39 haplotypes defined by 81 variable sites. Haplotype diversity was high with low nucleotide diversity, suggesting rapid population growth from a founder population with a small effective population size. The negative values of Tajima’s D and Fu’s Fs and small positive value of Ramos-Onsins and Rozas’ R2 suggest deviation from neutrality and population expansion. The haplotype network and demographic ex pansion parameters further suggest historical population expansion. Mismatch analysis statistics and Bayesian skyline plots estimate population expansion during the late Pleistocene. Al though the species presently occurs in small, disconnected we found no structuring of the population. Dispersal events are the most likely explanation for the absence of genetic structuring in the YTB population. These results represent important data for the design of a conservation plan for this endemic and globally threatened species.


INTRODUCTION
The yellow-throated bulbul Pycnonotus xantholaemus (YTB) is a habitat specialist, threatened passerine endemic to peninsular India. It is predominantly associated with scrub forests on hill slopes with exposed rocky outcrops and has patchy distribution across its geographical range (Ali & Ripley 1987, Subramanya et al. 2006. It is known from about 100 locations spread across peninsular India that include the Eastern Ghats mountain range, inland hillocks of the Deccan Plateau and eastern slopes of the southern West-ern Ghats mountain range. YTB habitat faces a severe threat from the mining industry. Quarrying of the rocky outcrops for granite and construction aggregate causes the destruction of associated scrub forest and subsequent habitat loss. India is a leading exporter of granite, and quarrying for granite is expected to expand in the region (Indian Bureau of Mines 2018). The YTB is classified as Vulnerable by the IUCN, due to its fragmented population and on going habitat loss in the region. Despite concern over the deficiency of data on abundance, population trends suggest a decline (BirdLife International 2016).
Molecular markers provide important measures to explore population genetic structure and geographic differentiation (Avise 1994). They have been employed to understand evolutionary processes and suggest appropriate conservation measures (Avise 1994). The non-coding mitochondrial control region (mtCR) evolves 3-to 5-fold faster than the remainder of the mtDNA, and it has been extensively employed to infer population genetic structure (Baker & Marshall 1997, Avise 2000. There have been few studies on Pycnonotus bulbuls utilizing variation in mtDNA sequences to understand the genetic structure in the population (olive-winged bulbul P. plumosus; Tang et al. 2016, and light-vented bulbul P. sinensis;Song et al. 2013).
Genetically distinct populations of species are likely to undergo extinction 3 orders of magnitude faster than the entire species (Hughes et al. 1997). Consequent loss of intraspecific genetic diversity (Manel & Holderegger 2013) makes populations susceptible to extinction. Given the threatened status of YTB, the persistence of small, isolated populations in a landscape that is fast transforming due to habitat destruction (Ramachandran et al. 2018) is a key factor for their survival. Therefore, in this study, we re port the first investigation of mtDNA sequence variation at the non-coding hypervariable control region in YTB to (1) describe the spatial patterns of genetic diversity between the populations, and (2) explore the signatures of any historical demographic changes.

Sampling
Mistnetting was carried out during March 2017 to September 2018 at 25 sites in 10 different locations covering the entire distribution range of YTB (Fig. 1). We employed 6 to 12 m long, 16 mm nylon mistnets and used call playback to capture YTB. Mistnetting was done within 3 h of sunrise. In all, 245 mistnetting hours were spent in 129 mistnetting sessions involving 150 field days. A total of 69 YTB were caught, of which 63 were unique individuals. Most of the recaptures were on the same day or the very next day at the same location. One YTB pair was recaptured at the exact location of original capture (13°21' 57" N, 77°11' 32" E) after 5 mo. Each captured bird was ringed, and its morphometric meas- Names of the states are in bold; the number shown near the sampling location denotes the total sample obtained from that particular region. Inset shows the location of mistnetting sites and the Deccan Peninsula in India. (b) Yellow-throated bulbul (YTB) in its natural habitat. (c) Typical YTB habitat: hillock with exposed granitic boulders and scrub vegetation. In the lower part of the hill, there is extensive habitat destruction, mining of soil for brick making and quarrying for construction aggregate urements viz. wing length, tail length, beak length (bill tip to notch in skull), beak (gape) width and beak depth were taken before it was released (Winker 1998). Aluminum leg rings, each engraved with a unique code, ensured that no individual was sampled twice. We either collected a blood sample from the brachial vein or plucked the outermost 2 tail feathers for DNA ana lysis. The bra chial vein was punctured using a fine sterile syringe, and the drop of blood obtained was adsorbed on Watmann filter paper strips and stored in sterile vials at room temperature. Similarly, the calamus portions of the plucked feathers were cut and stored in sterile vials at room temperature until DNA isolation. In all, blood samples from 6 individuals and feather samples from 57 individuals were used for the study.

DNA extraction, PCR amplification and sequencing
Total genomic DNA was extracted from all the samples using the standard method of phenolchloroform DNA isolation and precipitation (Sambrook & Russell 2001). Dithiothreitol (5 mM) was added to a lysis buffer for effective digestion of the calamus tissue. Primers used to amplify the control region in previous studies were not bulbul specific. Universal primers preferentially amplify the nuclear sequences of mitochondrial origin, or numts. These sequences evolve slowly and are similar to the ancestral mtDNA sequences (Sorenson & Fleischer 1996). To avoid cross amplification of numts, we designed bulbul-specific control region primers by aligning available mito chon drial genomes and control region sequences of Asian bulbuls and selecting the most conserved regions as primer binding sites. Newly generated primer sequences were analyzed using a sequence manipulation suite available online (Stothard 2000). We aimed at amplifying the partial control region of 1200 bp length and designed 2 pairs of overlapping primers: CRF1 (5'-GGT CTT GTA AAC CAA AGA TTG AAG-3'), CRR1 (5'-CTT CAA GTG CGT AGC AGG AGT T-3'), CRF2 (5'-GAA CCG AGC TAC TCA ACG TCA G-3') and CRR2 (5'-CCA ACT CCC AAA GCT GGC ATT TT-3'). Molecular sexing was not performed as the DNA yield was not sufficient in a few samples. Adult YTB occur in pairs, and we often caught the breeding pair in the same mistnet. Hence, we are confident that our samples are not sex biased.
The DNA amplification for the control region was performed using EmeraldAmp GT PCR Master Mix from Takara Bio. The PCR program included an initial denaturation at 95°C for 3 min, followed by 30 cycles of denaturation at 94°C for 30 s, annealing at 62°C for 30 s, extension at 72°C for 1 min and a final extension at 72°C for 5 min. PCR products were run on 1% agarose gel to confirm amplification of the desired amplicon: 746 bp for CRF1-R1 and 670 bp for CRF2-R2. PCR products were sequenced in the automated ABI 3730 DNA Sequencer (Applied Biosystems) using the same set of primers (all the samples were sequenced thrice for confirmation of variable sites). The sequences were analyzed and assembled in MEGA v.6.06 (Tamura et al. 2013) and Codon Code Aligner v.3.7.1 and were deposited in GenBank (accession numbers MT423843−MT423902).

Genetic diversity
Alignment and manual editing of the nucleotide sequences were performed in MEGA v.6.06 and Codon-Code Aligner v.3.7.1. Each variable site was visually confirmed. We grouped 63 samples in 10 populations based on geographic distance (> 50 km) and habitat continuity between the sampling sites ( Fig. 1). DNA poly morphism was quantified in DnaSP 6.12.03 (Librado & Rozas 2009), and the following diversity indices were estimated: number of segregating sites (s), number of haplotypes (h), haplotype diversity (hd), average number of pairwise nucleotide differences within populations (k) and nucleotide diversity (π). Among the 10 populations, Tirumala, Ponnuthu and Achampet had fewer than 3 samples; hence, population-level diversity indices, neutrality indices and demographic indices were not estimated for these 3 populations. Sequences from the 10 populations were pooled for the analysis of the entire dataset.

Historical demography
We employed 3 different methods -neutrality tests, mismatch distribution and Bayesian skyline plot ana lysis -to explore the historical demographic changes in the YTB populations. We chose these tests because of their statistical power in detecting significant changes in population size (Ramos-Onsins & Rozas 2002).
The deviation from neutrality of mtCR sequences was ascertained in DnaSP v.6.12.03 by estimating neutrality statistics -Fu's F s (Fu 1996), Tajima's D (Tajima 1996) and Ramos-Onsins and Rozas' R 2 (Ramos-Onsins & Rozas 2002). The neutrality test identifies whether a locus has evolved neutrally or under the influence of natural selection (Fu 1996). The significant negative values of D and F s and positive value of R 2 indicate an excess of low-frequency polymorphisms relative to expectations under the standard neutral model, implying population size expansion. Neutrality statistics (F s , D and R 2 ) were estimated separately for the populations and for the entire dataset using the empirical population sample size and the observed number of segregating sites. The statistical significance of the neutrality statistics was evaluated by comparing observed values with null distributions generated using 10 000 iterations in DnaSP v.6.12.03.
Mismatch distribution of the frequency of pairwise haplotype differences was implemented in Arlequin v.3.11 (Excoffier et al. 2007). Mismatch distribution ana lyzes the distribution of the number of site differences between pairs of sequences (Slatkin & Hudson 1991, Rogers & Harpending 1992. The expected values for a model of constant population size were plotted against the observed values. A unimodal pattern is indicative of exponential population growth in the recent past, while the multimodal pattern suggests long-term population stability via demographic equilibrium (Slatkin & Hudson 1991, Rogers & Harpending 1992. We validated the sudden demographic expansion model by estimating Harpending's raggedness index (HRI) and the sum of squared deviations (SSD). HRI quantifies the smoothness of the mismatch distributions and distinguishes between population expansion and stability (Harpending 1994). This index is generally expected to assume a large value in relatively stable populations. The SSD is an estimate of the fit between the observed and expected mismatch. We used demographic parameters (τ, θ 0 and θ 1 ) to estimate population expansion from an initial θ 0 to θ 1 in τ units of mutational time. The time in generations since expansion was calculated using the equation t = τ/2μ, where μ is mutation rate per generation for the sequence studied (Rogers & Harpending 1992). The time since expansion was estimated for YTB using the 1050 bp sequence of mtCR with a mutation rate of 1.6 to 2.8% per base pair per million years (Lerner et al. 2011). Since life history parameters were not available for this species, we used a generation time of 1 yr based on the estimates for the olivewinged bulbul Pycnonotus plumosus (Tang et al. 2016) and the Philippines bulbul Hypsipetes philippinus (Kennedy et al. 2000).
The Bayesian skyline plot (Drummond et al. 2005) method implemented in BEAST 2 (Bouckaert et al. 2014) was used to further explore the demographic history of YTB and to estimate the timing of the popu-lation expansion. We used 39 haplotype sequences to run the analysis using the coalescent Bayesian skyline tree prior. This coalescent-based approach estimates the posterior distribution of effective population size at intervals along a phylogeny, thereby allowing one to infer population fluctuations over time. Thirty-nine haplotype sequence alignments were imported in the program BEAUti2 (included in BEAST 2) in nexus format. The evolutionary model was set to Generalised Time Reversible, as estimated in jModeltest v.2.1.10 (Darriba et al. 2012) using Akaike's information criterion (Akaike 1973). We implemented a strict clock model, with an average substitution rate of 2.2% per base pair per million years (Lerner et al. 2011). The output file from BEAUti was run in BEAST 2. Two independent chains were run for 10 million generations each, sampling every 5000 generations, and the first 20% was discarded as burn-in. Convergence to stationarity was confirmed in TRACER v.1.6 (in cluded in BEAST 2) when an effective sample size > 200 was achieved for all the parameters. The Bayesian skyline plot was constructed in TRACER v.1.6 and exported to Inkscape v.0.91 for further annotation.

Population structure
The genetic structure of the haplotypes was plotted in NETWORK v.4.6 (www.fluxus-engineering.com). A haplotype network generated from intraspecific genetic data allows visualization of the strength of the relationship between haplotypes, and it also indicates the missing mutational connections (median vectors) (Bandelt et al. 1999). Aligned haplotype se quences (variable sites only) were imported in NETWORK v.4.6 in Roehl data format. The network parameters were generated using the median-joining algorithm. No pre-processing (star contraction) or post-processing (maximum parsimony calculation) was done. The haplotype diagram was drawn and annotated using NETWORK v.4.5.1.0 and Inkscape v.0.91.
Pairwise F ST values for all populations were calculated, and an analysis of molecular variance was performed in Arlequin v.3.11 to estimate the degree of genetic difference between disparate populations. Their significance was assessed by 1000 permutations. Geographic distances between populations were generated in ArcGIS 10.3 using the geographic coordinates of individual populations. We calculated the mean of the sample locations within a population to measure interpopulation distance. A Mantel test was conducted in R v.3.4.3 (R Core Team 2017) using the package ade4 (Dray & Dufour 2007) to test for isolation by distance, i.e. the association between matrices of pairwise geographic distance (in km) and genetic distance (F ST /(1 − F ST )) with 10 000 randomizations. A positive correlation between geographic distance and genetic distance would lend support to the hypothesis of isolation by distance. Isolation by distance occurs when the migration−drift equilibrium is reached over a long time after demographic changes (Crispo & Hendry 2005).

Sequence polymorphism
The analysis of the mtCR of YTB produced 60 sequences of 1100 to 1200 bp. Three individuals did not yield a good quality sequence and hence were not included in the analyses. After the alignment of 60 sequences in MEGA v.6.06, sequence extremes were trimmed off based on the quality of the readout. No gaps/ indels were found in the alignment (Table S1 in the Supplement at www. int-res. com/ articles/ suppl/ n043 p199_ supp. xlsx). A final set of 1050 bp of the mtCR from 60 YTB samples belonging to 10 different locations was used for the population genetics analysis. The total G+C content of the sequences was 0.474, there were 44 singleton variable sites (37 of which were parsimony informative), the total number of mutations was 85, and 81 polymorphic sites defined 39 haplotypes of varying frequencies and 12 median vectors. Overall, the YTB population showed high genetic diversity (hd > 0.5) and low π (< 0.005) ( Table 1).

Historical demography
We retrieved significant negative values for Tajima's D (−0.155, p < 0.05) and Fu's F s (−0. 60, p < 0.05) and a small positive value for Ramos-Onsins and Rozas' R 2 (0.10, p < 0.05) for the entire dataset consisting of 60 samples from 10 populations. The F s statistic was significantly negative (F s = −4.52, p < 0.05) only in the Ramanagara population. D, F s and R 2 values were not significant for the rest of the populations (Table 1). The negative values of D and F s and positive value of R 2 indicate population expansion.
The mismatch distribution graph was predominantly unimodal. There was a prominent peak at the 7 th pairwise difference and a weak peak at the 42 nd pairwise difference position (Fig. 2) Table 1. Genetic diversity for yellow-throated bulbul populations. Parameters for analyses at the intrapopulation level estimated under the sudden expansion model in Arlequin.
n: number of sequences; s: number of segregating (variable/polymorphic) sites; h: number of haplotypes; hd: haplotype/gene diversity (± SE); k: average number of nucleotide differences; π: nucleotide diversity (per site); τ: expansion parameter before (θ 0 ) and following (θ 1 ) expansion; SSD: sum of squared deviations; HRI: Harpending's raggedness index; na: not available is indicative of demographic stability. However, the values of the expansion parameters, θ 0 and θ 1 , were significantly different, which was indicative of population expansion ( Table 1). The HRI and SSD values were non-significant for all populations, except for the Tumkur population. These values lend support to population expansion rather than population stability (Table 1). We estimated the time since expansion for YTB based on the τ value of 5.8 for the entire population (Table 1). We estimate that the YTB population expanded in peninsular India during the Pleistocene, between 172 619 and 98 639 yr ago. The Bayesian skyline plot estimated the timing of the onset of population expansion of YTB at approximately 0.1 million years ago (mya) (Fig. 3). The estimated time of demographic expansion corroborates well with the estimate obtained using mismatch distribution statistics.

Population structure
The median-joining haplotype network consisted of 39 haplotypes and 12 median vectors. Among the 39 haplotypes, 29 were singleton (haplotype seen only once in a sample, i.e. an unshared haplotype) ( Table S1). The topo logy of the haplotype network showed that there is no locality-wise segregation of haplotypes, and there was no structure in the populations. Various haplotypes had mostly single nucleotide differences between them (Fig. 4). The abundance of singletons suggests a recent population bottleneck followed by population expansion (Taji ma 1989, Slatkin & Hudson 1991). Only 3 of 45 locality pairs had statistically significant F ST values (Table 2) Lines on bars that separate haplotypes represent the mutational steps between the haplotypes. Empty white circles represent haplotypes that were not sampled, also known as median vectors acting independently on the populations. Withinpopulation variance accounted for 98.82%, while among-population variance accounted for 1.16% of the total variation (Table 3). Almost all of the variation was observed within populations rather than among populations, which suggests that there was no critical historical disruption of gene flow. There was no significant relationship between genetic distance and linear geographic distance (Mantel's r = −0.064, p = 0.798, 9999 runs).

DISCUSSION
The YTB partial control region sequence had slightly higher GC content than the average Pycnonotus mitochondrial GC composition (0.4485− 0.4601, Ren et al. 2016a,b). It had high hd and low π. Similar results have been reported for other bulbul species: olive-winged bulbul Pycnonotus plumosus (Tang et al. 2016), light-vented bulbul P. sinensis (Song et al. 2013), collared finchbill Spizixos semitorques (Gao et al. 2011) and Philippines bulbul Hypsipetes philippinus (Silva-Iturriza et al. 2010). High hd and low π reflect that the population has undergone rapid expansion from a small ancestral population, with a small effective population size (N e ) (Rogers & Harpending 1992). As such, the time since population expansion is sufficient for the recovery of haplotype variation via mutation, resulting in a high hd value, but too short for the accumulation of a large difference in the sequences, hence the low π value (Grant & Bowen 1998, Avise 2000. The lack of an isolation by distance relationship suggests that mutation− drift equilibrium post-demographic expansion was not evident in the YTB populations. Further evidence of demographic expansion was provided by the neutrality tests. The overall negative values of Tajima's D and Fu's F s and positive value of Ramos-Onsins and Rozas' R 2 indicate an excess of rare mutations in the populations, implying recent population expansion. Alternatively, these values could have resulted from balancing selection on a nearby locus. Further work using neutral nuclear DNA markers might resolve this. The YTB mismatch distribution graph was weakly bimodal, which is indicative of population substructuring and mutation rate heterogeneity (Aris-Brosou & Excoffier 1996). The second peak in the mismatch distribution was contributed by Hap_1, the most diver gent of the 39 haplotypes; this came from 1 bird from the Tumkur population. Hence, the ob served bi modality of the mismatch distribution graph could be the result of a single highly divergent se quence from the Tumkur population, rather than due to demographical stability and/or population substructuring. This explanation is supported by the fact that the SSD and HRI values under the demographic expansion model for all the populations but Tumkur were non-significant, which indicates that the data are a relatively good fit for the demographic expansion model (Harpending 1994 : 199-207, 2020 with spells of aridification in the region as the climate oscillated between semi-arid and humid (Hora 1952 Our sampling was not equal across the geographic range of the species. About 70% of the samples were from 4 locations in the state of Karnataka, while the remaining 30% of the samples were from 6 locations in 3 states. The median vectors of the haplotype network represent unsampled genotypes or extinct ancestral sequences (i.e. representing missing intermediates) (Bandelt et al. 1999). Since the YTB population appears panmictic and our samples were not spatially uniform, we expect that several haplotypes in the populations were not captured by our dataset. The low sample sizes can yield false positive results (population structure observed when none exists). In our case, increasing the sample size will yield a resolved haplotype network with fewer or no median vectors. However, it will not change the inference drawn in this study that YTB exists as a panmictic population and underwent recent demographic expansion.

CONCLUSIONS
YTB occurs in small numbers in restricted habitats throughout its geographic range (Ali & Ripley 1987, Subramanya et al. 2006). It has been described as a rare and sedentary species (Subramanya et al. 2006). Such species are likely to show population structuring due to limited genetic exchanges between populations. They also face the threat of stochastic extinction due to small population size effects. Dispersal events can override these limitations and avert stochastic extinction. Our data show that YTB has high genetic diversity and there is no structure in the populations. This is most likely due to dispersal events. We hypothesize that dispersal events from natal areas are essential for long-term persistence of the isolated populations of YTB. Suboptimal habitats spread between the optimal habitats might serve as corridors for the dispersal of YTB. Ongoing mining activities to meet the growing demand for construc-tion aggregate are quintessential threats to YTB habitat and dispersal corridors. Because YTB habitats are not well represented in the protected area network in the region (Rahmani 2009) and dispersal events have shaped the YTB population, conservation efforts should focus on the landscape level. To leverage support from local communities and managers, we propose that YTB be considered a non-traditional flagship species (Entwistle 2000) for isolated inland hillocks and their associated scrub forests in peninsular India.