Phylogeography and diversification of Oriental weaverbirds (Ploceus spp.):A gradual increase of eurytopy

Abdul Razaq, Govann Forcna, Urban Olsson, Qan Tan, Robert Tzard,Nan Ln, Nla Pwnt, Aleem Ahmed Khan

aInstitute of Pure & Applied Biology, Zoology Division, Bahauddin Zakariya University, Multan, Pakistan

bUniversidad de Alcal(UAH),Global Change Ecology and Evolution Research Group(GloCEE),Departamento de Ciencias de La Vida,28805,Alcala de Henares,Spain

cCIBIO, Centro de Investigação em Biodiversidade e Recursos Genticos, InBIO Laboratrio Associado, Universidade do Porto, Campus de Vairão, 4485-661, Vairão,Portugal

dBIOPOLIS Program in Genomics, Biodiversity and Land Planning, CIBIO, Campus de Vairão, 4485-661, Vairão, Portugal

eDepartment of Biology and Environmental Science, University of Gothenburg, Box 463, SE-405 30, Gothenburg, Sweden

fGothenburg Global Biodiversity Centre, Box 461, SE-405 30, Gothenburg, Sweden

gDepartment of Biological Sciences, National University of Singapore, 16 Science Drive 4, Republic of Singapore

hWildlife Conservation Society, Global Conservation Program, Bronx, NY, USA

iNature and Wildlife Conservation Division, Ministry of Natural Resources and Environmental Conservation, Naypyidaw, Myanmar

Keywords:Avian phylogenetics Indian subcontinent Indochinese region Multilocus analyses Ploceidae

ABSTRACTWeaverbirds are a speciose group of colorful passerines inhabiting the Old World Tropics.Nevertheless, the Oriental weaverbirds(Ploceus spp.),widespread across southern Asia,are much less diverse and restricted to a few ecological niches compared to their African counterpart.To investigate their phylogeography, we retrieved 101 samples of Baya Weaver (P.philippinus), Streaked Weaver (P.manyar), Black-Throated Weaver (P.benghalensis)and Asian Golden Weaver(P.hypoxanthus)along with GenBank sequences of Finn's Weaver(P.megarhynchus).We reconstructed the first molecular phylogeny based on a dataset consisting of both mitochondrial and nuclear genes,dating the most recent common ancestor of Oriental Ploceus to～11 mya.Subsequent speciation appears to have been a combination of divergence within the Indian subcontinent and dispersal across a barrier situated between the Indian subcontinent and the Indochinese region, which provided habitats with a varying degree of isolations and ultimately promoted divergences in allopatry.Two descendants of the earliest nodes,P.megarhynchus and P.hypoxanthus,are both rare and local,often found near large river systems,which perhaps reflects niche conservatism and a lack of adaptive potential.The three smaller species are all widespread,common and less habitat specific.The most recent divergence, between western and eastern P.philippinus populations, is supported by both phylogenetic and morphological evidence, pointing toward limited gene flow between them.However, a zone of intergradation may exist in Myanmar and Brahmaputra flood plains, thus preventing a recommendation for species level recognition without further study.

1.Introduction

Weaverbirds (Aves: Ploceidae) are predominantly granivorous small and medium-sized passerines ranging from Sub-Saharan Africa, the southeasternmost part of the Arabian Peninsula and Indian Ocean islands,where—depending on species limits—110–119 species occur,to tropical Asia,which is home to just five species(Dickinson and Christidis,2014; del Hoyo and Collar, 2016; Gill et al., 2021).The latter are confined to the Oriental region, with the westernmost and easternmost boundaries marked by the Indus River in Pakistan and the Wallace line,respectively, and all belong to the genus Ploceus: Baya Weaver(P.philippinus), Streaked Weaver (P.manyar), Black-throated Weaver(P.benghalensis), Asian Golden Weaver (P.hypoxanthus), and Finn's Weaver(P.megarhynchus)(Fig.1).While Afrotropical ploceids inhabit a variety of habitats from sea level to forests and even alpine highlands across the Afrotropical region (Crook, 1963; del Hoyo et al., 2010; Sinclair et al.,2010),their Oriental counterparts are restricted to grasslands and swampy savannas (Ali and Ripley, 1983; Robson, 2000).Weaverbirds have long raised the interest of researchers from different areas due to their astonishing heterogeneity in terms of morphological,behavioral and socioecological features (e.g., Covas et al., 2002, 2004a, 2004b,2006; Brown et al., 2003; Covas and du Plessis, 2005; Spottiswoode,2009),one of the most peculiar being the capability of several species to weave their nests from leaf blades and foliage material(Crook,1963;Fry and Keith, 2004;del Hoyo et al.,2010).

The family Ploceidae was introduced by Sundevall (1836).The taxonomy of this group has been reviewed many times over the last century,with a seminal reappraisal conducted by Chapin (1917) on the basis of the size of the tenth and outermost primary preceding a comprehensive reassessment by Moreau (1960), which included additional morphological characters as well as eco-behavioral traits.According to most current literature (e.g., Dickinson and Christidis, 2014; del Hoyo and Collar,2016;Gill et al.,2021),the genus Ploceus contains 64–67 species,but in recent molecular studies based on mitochondrial and nuclear DNA markers, three Oriental species - P.benghalensis, P.manyar and P.philippinus-form a clade nested within(P€ackert et al.,2016)or sister to(De Silva et al., 2017, 2019) the one including the Afrotropical taxa Quelea and Foudia.De Silva et al.(2017)proposed that the genus Ploceus should thus be restricted to include only the five Oriental members,while their Afrotropical congeners as well as Anaplectes spp.should be transferred to the genus Malimbus.

Fig.1.Distribution of Oriental weaverbirds and sampling localities (indicated with numbers; see Appendix Table S1 for further details).(A–E)Range of each one of the five species.Credits for the bird pictures (not to scale; all modified from the original version) used in this and following figures are reported in the Acknowledgments.(F) Map showing the most relevant geomorphological features (i.e., mountain ranges, river basins and plateaus) mentioned in the paper.

According to De Silva et al.(2017), the age of the most recent common ancestor (MRCA) of both the Oriental clade (Ploceus sensu stricto,sensu De Silva et al., 2019) and the African clade (Malimbus, sensu De Silva et al., 2019) are very similar, yet the latter has undergone an adaptive radiation and contains a number of species an order of magnitude higher than their Oriental counterparts, which have diverged in a much slower manner.Several studies have been carried out on African species,whereas their Oriental relatives are less well-known,and the few pertinent phylogenetic studies so far performed have not included nuclear markers, except for single species (De Silva et al., 2017, 2019).Moreover,there are no biogeographic studies based on genetic data.This clade is interesting in several ways: two of the species, P.megarhynchus and P.hypoxanthus, are rare and local, while three smaller species are widespread and common, with one of them, P.philippinus, occurring in sympatry with the other species throughout most of its range(Gill et al.,2021).Interestingly,this species shows marked morphological differentiation between western and eastern populations, with the Northeast India-Myanmar area appearing to be a zone of intergradation (Hartert,1902; Ticehurst, 1932).For all species, several subspecies have been proposed based on morphology, but their taxonomic significance has never been validated with genetic tools.Perhaps the morphologically clearest differences are found within P.philippinus (Appendix Fig.S1),where western birds are pure yellow on crown and breast, with yellow stripes on the mantle and a blacker face mask,whereas eastern birds are more orange tinged on the crown,with buffier breast and stripes on the mantle(del Hoyo et al., 2010).

All five Oriental Ploceus are associated with wetlands to some degree,often in connection with major river basins, the emergence and developmental history of which played a crucial role in the evolution and reshuffling of fauna from the Indian subcontinent (Karanth, 2003) and,arguably,also of this clade.When examining the distribution and ecology of Oriental Ploceus, it is clear that, unlike their African counterparts, its species have not been equally adapted to drier savanna or forests,which may partly explain why a local adaptive radiation did not occur in Asia.Disentangling the historical biogeography of this avian group may help improve our understanding of the drivers which shaped the biodiversity of the Indo-Malayan region.

In this study, we analyzed a combined dataset of mitochondrial and nuclear DNA markers to: i) elucidate the phylogeography of the five Oriental members of the genus Ploceus, (ii) infer their ancestral ranges and biogeographical evolution, and (iii) shed light on possible paleoclimatic events that may have been involved in driving Ploceus diversification across southern Asia.

2.Materials and methods

2.1.Taxon sampling and molecular markers

For reasons of practicality,in this work we refer to the genus Ploceus according to the traditional taxonomy(i.e., including both Oriental and Afrotropical taxa).We sampled seven P.benghalensis,four P.hypoxanthus,53 P.manyar, and 43 P.philippinus from different areas of South and Southeast Asia,while DNA sequences of P.megarhynchus were obtained from GenBank(accession numbers KJ455582 and KJ455132).Moreover,four Afrotropical Ploceus species and 10 other Ploceid genera were used as outgroups along with the House Sparrow(Passer domesticus) (Appendix Table S1).Genomic DNA from feathers, muscle and toe pads was extracted using the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany).Feather and toe pad samples were incubated in 30 μL 0.1%DTT with Proteinase K and lysis buffer to increase the DNA yield.

We partially sequenced two mitochondrial genes: nicotinamide adenine dinucleotide dehydrogenase subunit 2(ND2)and cytochrome-b(Cyt-b).We also sequenced a fragment of one nuclear gene, the glyceraldehyde-3-phosphodehydrogenase gene intron 11 (G3P).Two datasets were produced: (I) only the mitochondrial genes (the first 813 bp and 879 bp of the Cyt-b and ND2 alignments,respectively),which was used for demographic and genetic diversity analysis in subsets of the Oriental ploceids; (II) the concatenation of all the mitochondrial and nuclear loci(Cyt-b:1041 bp;ND2:749 bp;G3P intron 11:292 bp),which was used for investigating the radiation of Afrotropical and Oriental Ploceus.We also reconstructed single-locus phylogenies of all taxa and built a ND2 network including only Oriental ploceids(n¼42;first 727 bp of the alignment).

2.2.PCR amplification and sequencing

PCR amplifications were performed in 12.5 μL reaction volumes consisting of 10 μL VWR Red Taq DNA Polymerase Master Mix(Tris-HCl(NH4)2SO4, 2.2 mM MgCl2, 0.11% Tween 20, 0.22 mM of each dNTP,0.11 units/μL VWR Taq DNA Polymerase,inert red dye),0.5 μL of each primer(Appendix Table S2),1 μL deionized H2O and 0.5 μL of DNA on an Eppendorf autothermal cycler (AG Eppendorf, Germany).Cycling conditions for Cyt-b were 5 min at 95 C,40 cycles of 40 s at 95 C,1 min at 45 C, and 2 min at 72 C, and a final extension of 6 min at 72 C.For some degraded samples, we sequenced short overlapping fragments,using specific primer pairs and nested PCR protocols (Appendix Table S2).Cycling conditions for the ND2 were as for the Cyt-b but with an initial denaturation of 2 min and an annealing temperature of 50 C that we increased up to 58 C for 40 s in the case of G3P, for which we used a final extension of 1 min.PCR amplicons were then visualized on 1.0%agarose gels stained with GelRed and purified following a standard ExoSAP protocol (Fermentas Life Sciences, Waltham, MA, USA).Sequencing of purified PCR products was performed by GATC Biotech AG(Cologne,Germany).

2.3.Molecular diversity and demographic analysis

Sequences were aligned using GeneStudio Professional 2.2.0.0(GeneStudio, Suwanee, GA, USA); however, some manual adjustments were necessary for G3P sequences.We translated the nucleotide sequences of mitochondrial genes to amino acids and aligned them using ClustalW in MEGA 7.0 (Kumar et al., 2016).Substitution models were estimated for each locus separately based on the Bayesian Information Criterion (BIC) as implemented in JModelTest 2.1.10 (Darriba et al.,2012).

A ND2 haplotype network was constructed using the median-joining(MJ)method(Bandelt et al.,1999)and setting the ε-tolerance parameter to 0 in POPART 1.7(Leigh and Bryant,2015).We estimated the number of variable and informative sites for each locus across the focal species(Appendix Table S3).Standard genetic diversity indices of Oriental Ploceus were calculated in DnaSP 6.10 (Rozas et al., 2017) using the two mitochondrial DNA loci separately(n¼64;666 bp and n¼34;879 bp for Cyt-b and ND2, respectively).Coverage by loci was uneven among samples,and only part of the included samples were sequenced for both mitochondrial loci.We divided the samples of western P.philippinus and P.manyar into two separate batches, one including all samples of respective taxon,and one with only those two samples that had data for two loci.Specifically, we computed the number of haplotypes (h),number of segregating sites(S),haplotype diversity(Hd),and nucleotide diversity(π)for the whole sample as well as for species and clades within species(P.philippinus;see Results).We also used DnaSP for various tests of demographic equilibrium:Tajima's D(Tajima,1989),Fu's FSneutrality test(Fu,1997), R2test(Ramos-Onsins and Rozas,2002), and mismatch distributions r (Rogers and Harpending, 1992).Coalescent simulations were run to test statistical significance.For populations showing evidence of deviation from demographic equilibrium,we included tau and calculated the time for the onset of population expansion(texp)based on a 2.1% per million years divergence rate for Cyt-b, using the mismatch calculator provided by Schenekar and Weiss(2011).Generation time was based on Bird et al.(2020).There is less consensus regarding the divergence rate of ND2,so for this locus texpwas not calculated.

Changes in effective population size over time based on single and multiple loci were inferred by Coalescent Bayesian skyline plot (CBSP)and Coalescent Extended Bayesian Skyline plot(CEBSP),respectively,as implemented in BEAST2 (Bouckaert et al., 2014).Specifically, five single-locus CBSP runs were conducted for all P.philippinus(Cyt-b,ND2),western P.philippinus(Cyt-b,ND2)and P.manyar(Cyt-b)individuals with 107states of Markov chain Monte Carlo(MCMC)and sampling every 105states.In view of the marked distinctiveness between the two clades,we additionally performed two separate multilocus CEBSP runs for comparative purposes,one for all P.philippinus and one for only western P.philippinus individuals, respectively, with 109states of MCMC and sampling every 107states.We could not process the eastern individuals alone due to their limited sample size.For all runs,we used GTR þ G þ I substitution model for ND2 and HYK þ G þ I for Cyt-b.Molecular clock rate was set to 0.021 substitution per site per million years with strict clock model.A burn-in of the first 10%of each run was discarded in both analyses.

The evaluation of independent runs in terms of convergence and mixing was carried out by observing and comparing traces of each statistics and parameter in Tracer 1.6 (Rambaut and Drummond, 2013).Effective sampling size (ESS) values > 200 were considered sufficient indicators of parameter mixing.

2.4.Phylogenetic analyses and divergence time estimates

Molecular phylogenies based on Maximum Parsimony (MP) and Maximum Likelihood(ML)phylogenetic reconstructions were made with PAUP 4.0b10a(Swofford,2002).The models that best fit the multilocus analysis(MLA)of Cyt-b þ ND2 þ G3P are shown in Appendix Table S3.Trees were estimated by analyzing all three loci separately and all genes concatenated(dataset II),and by Bayesian inference with BEAST 1.10.4(Drummond et al., 2012).The xml files were created in BEAUti 1.10.4.Analyses were run under a birth-death incomplete sampling tree prior and a relaxed lognormal clock with two different approaches to dating the phylogeny.One used a normal prior (ucld.mean) for Cyt-b with a mean of 0.0105 and a standard deviation of 0.000567,corresponding to a substitution rate of 2.1% per million years or 0.0105 substitutions/site/lineage/million years(Weir and Schluter,2008);for ND2 and G3P we used a relaxed lognormal clock model and a CTMC rate reference prior (ucld.mean).The other approach fixed the age of the MRCA of the least inclusive clade containing both Passer and Ploceus to 20.7 mya (million years ago) based on Olsson and Alstr€om (2020),applying a normal prior and a standard deviation of 0.05, which produced a 95%HPD of 20.2–21.2 mya;for Cyt-b,ND2 and G3P we used a relaxed lognormal clock model and a CTMC rate reference prior (ucld.-mean).Other priors were applied with default values.Four Markov Chain Monte Carlo (MCMC) chains were run for 5 109generations, with sampling every 104generations.The logfile was analyzed in Tracer to assess whether effective estimates of the posterior distribution of the parameters had been acquired.The first 25% of the generations was discarded as burn-in.Maximum clade credibility trees with median node heights were built using TreeAnnotator 1.8.4(Rambaut and Drummond,2016), and visualized with FigTree 1.4.3(Rambaut,2016).

2.5.Ancestral range reconstruction

The evolution of geographical ranges of Oriental ploceids were reconstructed using statistical dispersal-vicariance analysis (S-DIVA) as implemented in the package RASP (Reconstruct Ancestral State in Phylogenies: Yu et al., 2015) with default parameters.The distributional range of the taxa included in this study was subdivided into four bioregions(area code:A–D):(A)Indian subcontinent(Indian sub-region and Ceylonian sub-region); (B) Indochinese sub-region; (C) Indo-Malayan sub-region; (D) Afrotropical region.To account for uncertainty in the tree topology,all trees from the BEAST MLA analysis were used as input.The optimal S-DIVA reconstruction was summarized on the maximum clade credibility tree of the concatenated dataset.

3.Results

A total of 1790 bp were obtained for the joint mitochondrial DNA genes(1041 bp and 749 bp for the Cyt-b and ND2,respectively,without evidence—i.e., no stop codons or indels—of being pseudogenes), while the amplified G3P intron was 292 bp long.After alignment, the concatenated sequences were 2082 bp long from the five sampled species(Appendix Table S3).The sequences were deposited in GenBank (Appendix Table S1).

3.1.Genetic diversity of oriental weaverbirds

The DnaSP analysis identified a total of 32 and 34 haplotypes among the Cyt-b and ND2 sequences, respectively, based on 141 (32 transition and 5 transversion pairs) and 170 (46 transition and 12 transversion pairs)parsimony informative sites,respectively.For western populations of P.philippinus, both Cyt-b and ND2 were sequenced in 16 individuals.For P.manyar, both Cyt-b and ND2 were sequenced in 6 individuals(Appendix Table S3).Standard diversity indices computed for the different datasets and the individual loci are represented in Table 1.The MJ network built on the basis of 33 ND2 haplotypes showed them to be joined in a complicated pattern of well-separated clusters not mirroring the inter- and intraspecific relationships inferred by the phylogenetic analyses.Two of these clusters corresponded to the two P.philippinus clades (H6 and H16-H33 plus H11) and differed by 33 substitutions,which is comparable to the species level divergence between eastern P.philippinus (H6) and P.manyar (H12, H13 and H15).A remarkable diversity was found also within clade A (mostly including Pakistani representatives of P.philippinus,the only exception being the sample from Myanmar) made up of as many as 19 haplotypes differing from the locally most abundant haplotype(H20)by up to 5(H21 and H27)and 7(H23) substitutions.However, the highest intraspecific sequence divergence was found in P.benghalensis, with haplotypes H9 and H10(belonging to individuals from Pakistan and India,respectively)differing by 15 substitutions(Fig.2).

3.2.Demographic analyses

Significant Fu's FS,R2and Tajima's D values suggest recent population bottleneck or selective sweeps to have occurred in the western P.philippinus clade (Table 1).Such demographic dynamics is also observed in P.manyar (Fu's FS,R2and Tajima's D; both loci), but not in other studies species.Ploceus hypoxanthus exhibited significant R2value but only for Cyt-b.Mismatch distribution raggedness index r generally did not significantly deviate from the null hypothesis (which is population expansion for r) except for ND2-based results for western P.philippinus and P.hypoxanthus (but note that this estimator is considered the weakest, and in case of conflicting results should be awarded least confidence:Ramos-Onsins and Rozas,2002).

Coalescent simulations on P.philippinus suggests a population decline at 200–300 kya (thousand years ago) followed by a period of recovery,likely associated with the expansion of the western population, which started around 60–180 kya(Fig.3A,Appendix Figs.S2A and S2B).This recovery might have been mainly contributed by the onset of population expansion of western P.philippinus starting 124–145 kya (Fig.3B, Appendix Figs.S2C and S2D).In contrast, P.manyar exhibits an overall demographic stability over most time, with a slight recent increase starting approximately 40 kya(Appendix Fig.S2E).

3.3.Phylogenetic analyses

The MLA tree topology was consistent across all the methods tested(Bayesian inference,MP and ML)and well supported(Appendix Fig.S3).The five target species were all reciprocally monophyletic, with P.megarhynchus being sister to the rest of the Oriental ploceid clade and P.hypoxanthus to the three smaller species.Among these,P.benghalensis was sister to a clade containing P.manyar and P.philippinus, with the latter further divided into two strongly supported clades.Clade A included individuals from Pakistan and Myanmar, whereas clade B consisted of individuals from Thailand and Singapore.Overall, the phylogenetic placement of taxa within this species complex met expectations based on morphological data, with the exception of our only P.philippinus sample from Myanmar, which clustered in clade A withwestern representatives while displaying the typical phenotype of the eastern subspecies(see Discussion).

Table 1Descriptive statistics on genetic variation in Oriental Ploceus, with significance determined using coalescent simulations; all calculated in DNASP 5.10.(Librado and Rozas, 2009).Coverage by loci was uneven among samples.In the table, all clades except Clade 1b and Clade 3b include all sequenced samples belonging to any particular geographical population.For western populations of Ploceus philippinus both Cyt-b and ND2 were sequenced for 16 individuals.The results for these are given under Clade 1b.For Ploceus manyar, both Cyt-b and ND2 were sequenced for 6 individuals.The results for these are given under Clade 3b.We included tau for populations indicated to deviate from demographic equilibrium, and calculated texpbased on Cyt-b data.

Fig.2.ND2 median-joining network of Oriental ploceids.(A) Overall structure.(B–F) Close-up on haplotype by species.The correspondence between haplotype number and sample ID is reported in Appendix Table S1.The number of mutational steps between haplotypes (i.e., pies) is reported along the line connecting them(when >1).

Tree topologies concerning the Afrotropical taxa were entirely consistent with results obtained by De Silva et al.(2017).The single-gene phylogenies based on ND2(Appendix Fig.S4),Cyt-b(Appendix Fig.S5)and G3P (Appendix Fig.S6), were not entirely consistent with the multigene phylogeny(Appendix Fig.S3), mainly concerning the position of Passer, but there were no strongly supported conflicts.Other minor discrepancies between the two mitochondrial and the nuclear marker should be interpreted on the basis of the relatively limited information content and low support of the G3P topology.

Fig.3.Historical dynamics of effective population size by Coalescent Extended Bayesian skyline plot on Cyt-b and ND2 for (A) P.philippinus and (B) western P.philippinus.CPD: Central Posterior Density.

3.4.Molecular dating

The median MRCA age of the clade including Foudia, Quelea and Oriental Ploceus was estimated to around 11.9 mya (Fig.4A, molecular divergence rate,hereafter Rate)or 11.3 mya(Fig.4B,node age calibration,hereafter Node),depending on calibration method(see Fig.4A and B for 95% HPD).The median age of the first divergence in the clade is estimated to have occurred around 9.49 mya(Fig.4A,Rate)or 9.0 mya(Fig.4B,Node).At this time lineages giving rise the P.megarhynchus on the one hand and all other Oriental Ploceus on the other hand diverged.The median age of the split between P.hypoxanthus and its sister taxa was estimated at about 8.25 mya(Fig.4A,Rate)or 7.82 mya(Fig.4B,Node).The median age of the MRCA of the three smaller species was inferred to around 5.97 mya (Fig.4A, Rate) or 5.68 mya (Fig.4B, Node), and the median age of the MRCA of P.manyar and P.philippinus at 3.72 mya(Fig.4A, Rate) or 3.56 mya (Fig.4B, Node).The median age of the divergence between western and eastern P.philippinus populations was estimated at 2.4 mya(Fig.4A,Rate)or 2.3 mya(Fig.4B,Node).Based on the estimates made by De Silva et al.(2017), the Oriental Ploceus split from their African sister clade 5.50–7.57 mya and their MRCA lived approximately 4.35–5.95 mya.

3.5.Ancestral area reconstruction

The ancestral area reconstruction in S-DIVA based on the four predefined bioregions (Fig.5A and B) suggested the Indian subcontinent to be the area where the MRCA of all extant Oriental ploceids lived.Its initial split from the African sister lineage was indicated to be attributable to vicariance(Fig.5C).The Indian subcontinent was pointed out as the region witnessing the major splits leading to the current species.The exceptions are nodes 2 and 5(Fig.5B),the inferred divergence between the ancestor of P.hypoxanthus and its sister lineage giving rise to all other Oriental weaverbirds, for which the range of the MRCA is inferred to have included both the Indian subcontinent and the Indochinese subregion, and the MRCA of the P.philippinus lineages, respectively.The divergence at node 2 (Fig.5B) was indicated to be a vicariance rather than a dispersal event(Fig.5C).The range of the MRCA of the currently widespread P.philippinus is ambiguous,but was most likely to the east of the Indian subcontinent.The species level splits at nodes 1, 3 and 4(Fig.5) pointed toward a MRCA most likely living in the Indian subcontinent.The split between P.manyar and P.philippinus(node 4,Fig.5C)is the divergence event most strongly inferred to be a result of dispersal.

4.Discussion

4.1.Systematic placement, origin and diversification of Oriental weaverbirds

Based on a multilocus approach and a more comprehensive sampling of Oriental Ploceus than previous studies, our phylogenetic results were nonetheless still consistent with their findings, namely supporting both the monophyly of this group and its sister relationship to the Afrotropical taxa Quelea and Foudia(Warren et al.,2012;P€ackert et al.,2016;De Silva et al., 2017, 2019) as well as the polyphyly of the genus Ploceus as currently circumscribed (Dickinson and Christidis, 2014; del Hoyo and Collar, 2016; Gill et al., 2021).We recovered six well-supported reciprocally monophyletic groups of Oriental Ploceus, corresponding to the currently recognized species in all phylogenetic reconstruction methods tested.Our results are congruent with those obtained in previous studies in terms of systematic relationships, but we also detected a deep intraspecific divergence between eastern and western P.philippinus.However,even though all populations of this species from Assam (northeastern India)and eastwards—assigned to P.p.infortunatus(Hartert,1902),P.p.burmanicus (Ticehurst, 1932), and P.p.angelorum (Deignan, 1956)—differ morphologically from the western ones in a similar way,our only sample from Myanmar belonged genetically to the western mitochondrial DNA clade.The nuclear marker did not resolve the topology at this shallow level, presumably due to the limited amount of information yielded in the short fragment (～200 bp) analyzed.As a matter of fact,most of the phylogenetic signal retrieved in our study is conveyed by mitochondrial DNA only.

The major lineages leading to Malimbus þ Anaplectes þ African Ploceus, Euplectes, and Oriental Ploceus þ Quelea þ Foudia split around 14.8–11.3 mya (cf.Olsson and Alstr€om, 2020).This timing covers the early stages of a shift in many Afrotropical ecosystems from forest to grassland during mid-Miocene,at a time when the C4 grasses are deemed to have first appeared (Jacobs et al., 1999).While divergence of these taxa occurred in the same temporal window which is typically associated with similar events in other Afrotropical avian groups (e.g., Crowe,1978),internal radiation in these lineages started later,presumably when the availability of open habitat niches had increased.

Interestingly, the RASP event density plot suggests that the divergence of Oriental Ploceus MRCA from its Afrotropical ancestor was a vicariance event, implying a once continuous range spanning parts of both the Afrotropical and the Indian subcontinent regions.We estimated that the range of the Oriental Ploceus MRCA(node 1,Fig.5),the starting point of diversification in Asia, was located in the Indian subcontinent where the ancestor of the clade may have lived along one or more of the major river systems,an ecoregion which today harbors all but one of the extant Oriental Ploceus species.The scenario emerging from our ancestral range reconstruction indicates that much of the subsequent Oriental Ploceus diversification involved dispersals and vicariance between the Indian subcontinent and the Indochinese region, centered right at the border between these areas where adverse geographical conditions acting as a barrier to dispersal are inferred to have been in place (cf.Barley et al.,2014;Chingangbam et al.,2015).

Fig.4.Molecular dating of ploceid multilocus phylogeny.(A)Dating based on a Cyt-b substitution rate of 2.1%per million years(Weir and Schluter,2008).(B)Dating based on the MRCA age of the least inclusive clade containing both Passer and Ploceus(Olsson and Alstr€om,2020).The color of branches refer to the distribution map of Oriental ploceids shown in Fig.1.2-Digit ISO codes are used for countries: PK: Pakistan; IN: India; MM: Myanmar; TH: Thailand; SG: Singapore.

The first major split internal to the Oriental clade appears to have been a result of an eastward range expansion,including dispersal through the Arakan Mountain Range into the Indochinese region,that culminated in an allopatric divergence(around 9–9.49 mya)between a presumably western population,the ancestor of P.megarhynchus,and the MRCA of all other Oriental Ploceus,which probably ended up occupying a range that included parts of both the Indian subcontinent (e.g., the Brahmaputra basin) and the Indochinese region (e.g., the Ayeyarwadi basin).At this time, the Ganges and Brahmaputra river basins were seemingly disconnected (Govin et al., 2018), and may have served as separate regions favoring allopatric divergence.This first split in the Oriental clade(node 1, Fig.5) was estimated to have occurred at roughly the same time as when the ancestor of Foudia diverged from Quelea.However, we estimated that the massive radiation of continental Afrotropical Ploceus(Malimbus, sensu de Silva et al., 2019) started 7.75–8.17 mya (Fig.3),perhaps a million years or more later than the initial branching event in the Oriental clade(nodes 1–3,Fig.5).

4.2.Paleoclimatic and paleogeological speciation drivers

During the time of the first inferred stages of the ancestral Ploceus becoming established in the Indian subcontinent, grasslands were expanding in both this region(Molnar et al.,1993;Copeland,1997)and eastern Asia (Molnar, 2005), spurred by a significant climate shift that led to the onset of present-day monsoon cycles around 8 to 9 mya.Interestingly,recent studies indicate that a significant monsoonal regime was already established in Southeast Asia as early as 40 mya(Licht et al.,2014a, 2014b), but both pollen and megafauna fossil records indicate that humid climate persisted during the Miocene before drying during the Late Pliocene (Kershaw et al., 2006; Morley, 2012).The establishment and spread of present-day local savanna biotas benefited from recurrent connections of Indochina to Sundaland during the Pliocene and Pleistocene,with savannas stretching from Myanmar to Java(Louys and Meijaard,2010; Suraprasit et al.,2014).Conversely,the Indian subcontinent had already experienced major paleoclimatic changes that led to the establishment of open habitats.Indeed, during late Miocene and Pliocene the northward shift of the Indian tectonic plate with the orogenesis of the Himalayas and the Western Ghats marked the onset of an arid climatic phase, causing the conversion of tropical evergreen forests into deciduous woodland and the appearance of grasslands over most of the Deccan Plateau (Meher-Homji, 1983).Evidence for this change in vegetation include that C4 plants were recorded isotopically as a minor dietary component for mammals in the Siwalik sequence of Pakistan at 9.4 mya, but increased extensively as a foraging resource about 5 mya (Morgan et al., 1994; Cerling et al., 1997).Furthermore,analyses of paleosol carbonates from Pakistan and Nepal (Quade et al.,1989,1995)and organic matter from the Bengal Fan(France-Lanord and Derry, 1994) indicate a shift from C3- to C4-dominated vegetation between 8.1 and 6.5 mya.This period also witnessed a significant faunal reshuffling in Pakistan,which is indicative of an increase in the extension of open habitats and ecological niches offering novel feeding opportunities to a plethora of avian taxa as well as boosting their potential for diversification and dispersal(Jacobs et al.,1999).

4.3.River basins and mountain ranges as drivers of diversification

The ancestor of P.hypoxanthus was estimated to have diverged from the incipient ancestor of P.benghalensis,P.manyar and P.philippinus in an inferred vicariance event between an Indian and an Indochinese subpopulation around 7.82–8.25 mya(node 2,Fig.5),probably separated by the Arakan Mountain Range.This eastern population may have initially become established in the Ayeyarwadi basin, and today occurs patchily from Myanmar to Java.The other population involved in this vicariance event was indicated to as having most likely remained in the Indian subcontinent.More specifically,it is plausible that at least a part of it was residing in the Brahmaputra basin,at this time still disconnected from the Ganges basin (Govin et al., 2018), being therefore separated to the east from the range of P.megarhynchus ancestor.The starting point of the next stage in Ploceus divergence may thus have been three allopatric populations,one in the Ganges basin,and one more recently divided between the Brahmaputra and the Ayeyarwadi basins.These events would have occurred in conjunction with an additional decrease of average temperatures leading to a more dramatic shrinking of C3 tropical forest patches as they were replaced by C4 grassland and savannas (especially in sub-Himalayas: Valdiya, 1998), which, starting from ca.7.4–7 mya onwards, increased to represent more than 90% of the biomass (Quade et al.,1989; Morgan et al.,1994).

Of the two lineages that evolved from node 3 (Fig.5), namely from the ancestral population presumably inhabiting the Brahmaputra basin around 5.68–5.97 mya, one became the ancestor of P.benghalensis, a species that natively occurs only in the Indian subcontinent at present,while the other would in time give rise to P.manyar and P.philippinus.One possible scenario for this split could be that the ancestor of P.benghalensis expanded its range westwards into the Ganges basin, as contact between the latter and Brahmaputra basin was established sometime between 7 and 5 mya (Govin et al., 2018).Such westward expansion may have reached the range of P.megarhynchus, possibly establishing the first case of sympatry among the Oriental Ploceus.In this scenario,gene flow between the ancestor of P.benghalensis,now living in the Ganges basin,and the population in the Brahmaputra basin may have eventually ceased, giving rise to the ancestor of P.manyar and P.philippinus therein.In an alternative scenario, P.benghalensis ancestor may have inhabited both the Brahmaputra and Ganges basin,while part of the population expanded eastwards across the Arakan barrier into the Indochinese sub-region,where it subsequently diverged.

Diversification among the two lineages descending from node 4(Fig.5) to become present day P.manyar and P.philippinus seems more complicated, and may have involved all Asian regions.It was reconstructed as being characterized by initial dispersal, evidently creating large ranges and a very complicated phylogeographic pattern.The location of the ancestral range of P.manyar and P.philippinus MRCA(node 4, Fig.5) appears similar to the range reconstructed at node 3(Fig.5), possibly involving the Brahmaputra basin, but with a slightly higher probability of a more easterly distribution.Around the time when P.manyar and P.philippinus lineages diverged, the ancestors of P.megarhynchus and P.benghalensis were established in the Indian subcontinent.Ploceus manyar and P.philippinus lineages emanating from node 4(Fig.5) both appear to have expanded their ranges eastwards to some extent.It is possible that P.manyar ancestor evolved in a range including both Brahmaputra and Ayeyarwadi basin,whereas the range of P.philippinus ancestor might have stretched further east, away from the Indo-Myanmar barrier.Denser and geographically more complete sampling would be required to understand the divergence patterns in this clade.

Noteworthy,this process coincides with a period when Central China experienced intensified precipitations, approximately 4.5–2.7 mya, also coinciding with a period of decreasing temperatures (Nie et al., 2014).The split between P.philippinus and P.manyar around 3.56–3.72 mya follows on the second great upheaval of the Tibetan Plateau in the middle Pliocene and the appearance of very high jagged Himalayan Ranges(Roberts, 1992).This resulted in the emergence of large river systems periodically flooding the downstream lowlands with fluvial incision favoring hillslope denudation across their entire basin (Lave and Avouac,2001),thus offering a suite of new ecological niches facilitating the early stages of secondary sympatry,even though all evidence indicate that the actual diversification process was driven by allopatric isolation and drift rather than segregation in sympatry.Recurrent climatic oscillation cycles of dry and wet periods led to expansion and contraction of humid grasslands and open swamps which progressively retreated to disjunct patches in India(Karanth,2003),possibly creating an opportunity for the most recently evolved species to adapt to new niches,less dependent on permanent wetlands.Both the ancestors of P.manyar and P.philippinus,respectively, appear to have been less demanding in their habitat requirements, which probably facilitated their successful and relatively recent colonization of areas already inhabited by congeners as well as novel areas.These species are now found across large parts of both the Indian subcontinent and the Indochinese region, inhabiting a variety of wetland areas away from the major rivers.

Finally, a deep divergence within P.philippinus emerged approximately around 2.3–2.4 mya.At this time Pleistocene climatic oscillations,bringing about humid and warmer conditions in lowlands and cooling of the Himalayan region, led to an increase in forest cover and a fragmentation of the open habitats inhabited by weaverbirds.

Most speciation events seem to be the result of dispersal followed by isolation in allopatry,but gradual changes in niche availability resulting from paleoclimatic events likely contributed to conditions allowing for secondary sympatry.It is plausible that at least during the early stages of the diversification process the colonization of the Indian subcontinent followed the major trajectories along the Indus,Ganges and Brahmaputra basins, as wells as the Ayeyarwady basin, similar to other faunal radiations(e.g.,Forcina et al.,2018),offering a suite of favorable habitats such as forest belts and large marshlands.Also, it is possible that the intervening areas between these basins in the Indian subcontinent may have served as barriers to dispersal, thus promoting vicariance within the subcontinent, as the evolution of P.benghalensis may indicate.Our ancestral area reconstruction points to multiple range expansions including dispersal events between the Indian subcontinent and the Indochinese region across the Arakan Mountain Range, an area that for long periods seem to have served as a barrier to gene flow.Once this barrier was cleared,more indiscriminate expansions into both the Indian subcontinent and the Indochinese region seem to have been straightforward,particularly in the more recent cases.It is worth mentioning that Myanmar—the country of Southeast Asia bordering the Indian subcontinent, but separated by the Arakan Mountain Range—represents the easternmost distribution limit for a suite of avian species ranging across southern Asia (Rasmussen and Anderton, 2012).Therefore, these mountains appear to constitute either a real barrier to dispersal, or a transition to a region of unfavorable habitat for these species (which is possibly the most likely hypothesis in view of the obstacle the range represents to the southwestern monsoon rains that would otherwise reach Central Myanmar).Moreover, Miocene cooling events initiated a gradual aridification of central Asia(Favre et al.,2015),which probably limited a northward expansion of Oriental ploceids along with the barrier of the Himalayas,and thus restricted the total area available to this group of birds as occurred with avian clades (e.g., francolins: Forcina et al.,2012).

4.4.Eurytopy and niche filling

There is no evidence of sympatric altitudinal or habitat-related speciation in the Himalayas(Price et al.,2014),whereas all sister taxa of young age are widely allopatric.Even though many species have closer extralimital relatives, the most closely related species that occur in sympatry in the Himalayas are estimated to be on average 7.5 mya old(Price et al.,2014).This is older than the 2–3 mya often regarded as the time required for speciation in birds, suggesting that niche filling is limiting immigration of new populations, thus slowing down the local speciation process.In contrast, according to our data, the average time between significant cladogenesis events among Oriental Ploceus is 1.7 million years.It is unclear when the range of P.megarhynchus ancestor was first invaded by another Ploceus species,but P.benghalensis ancestor seems to have been expanding in the Indian subcontinent within 4.5 million years of its ancestral lineage diverging from its MRCA with P.megarhynchus.The split between the now sympatric P.manyar and P.philippinus happened around 3.5 mya.Both these estimates are younger ages than for the Himalayan examples,suggesting that there is no similar slowdown in this clade,in turn suggesting that potential niches for these birds may not be filled yet.

Ploceus philippinus hosts the highest number of morphological subspecies (five), with a distinct morphological differentiation between eastern and western populations.Additional genetic data from geographically isolated populations not included in this study—for instance Indonesia—may reveal even more recent divergence among this and other species.The shortest time span between major splits and subsequent secondary sympatry in the Oriental Ploceus clade(3.5 mya),suggests that the approximately 2.3–2.4 mya that separate the evolutionary trajectories of eastern and western P.philippinus lineages may not be enough divergence time in allopatry to successfully maintain their identity in case they came into sympatry with one another.This may actually be under current testing in the Brahmaputra basin and Myanmar,where birds morphologically belong to eastern populations,but our only sample was carrying a western haplotype.Further investigations of eastern and western P.philippinus lineages are needed to clarify the degree of gene flow, cryptic introgression and backcrossing in regions where molecular and morphological evidence is contrasting.

There is no doubt that this species complex deserves major interest in terms of microevolutionary studies at lower taxonomic resolution.It offers an example of ongoing speciation in a clade of birds that have closely related ecological and morphological counterparts in Africa, but which have evolved in a much different fashion.The latter experienced an adaptive radiation in grasslands and savanna, and also at one point colonized forests and woodlands, where a parallel adaptive radiation occurred.None of this happened in Asia based on the ecology and distribution of local species (cf.Ali and Ripley, 1983; Robson, 2000), and future studies may investigate possible genetic, ecological or morphological predisposition favoring or hampering the ability to adapt to new circumstances in these groups.

The evolution of P.philippinus may be an interesting case study, not only due to its intermediate status in the speciation process, but also regarding its adaptive genetic capacity.It is the only Oriental Ploceus to inhabit heavily urbanized areas in Southeast Asia such as Singapore,pointing to the capability of adapting to new habitats.

In light of these considerations,future studies on morphology,niche requirements and molecular ecology of Oriental Ploceus are needed to decipher the evolutionary history and provide evidence and data beneficial to the long-term conservation of this neglected avian group.

Authors’contribution

AAK, UO: conceptualization, methodology, data curation, formal analysis, funding acquisition, writing–review & editing.AR: data curation, formal analysis, investigation, writing–original draft.GF: data curation, formal analysis, visualization, writing–review & editing.QT:formal analysis,visualization,writing–review&editing.RT,NL,NP:data acquisition,writing–review&editing.All authors read and approved the final manuscript.

Funding

This work was additionally supported by fellowships of the Ministry of Universities of the Spanish Government (María Zambrano/Next Generation EU)and the Portuguese Foundation for Science and Technology(FCT,PTDC/BAA-AGR/28866/2017)to GF.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

Authors are grateful to Sharon Birks (Burke Museum, University of Washington), Joel Cracraft, Paul Sweet and Thomas J.Trombone(American Museum of Natural History), Jon Fjeldså and Jan Bolding Kristensen (Zoological Museum, University of Copenhagen), Staffan Andersson (University of Gothenburg), Frank Rheindt and colleagues(National University of Singapore), Per Alstr€om (Uppsala University),Shahbaz Muhammad, Ghulam Qammer-ud-din Niazi and Muhammad Shakeel(Bahauddin Zakariya University Multan,Pakistan)for providing samples.UO was supported by Mark&Mo Constantine.The authors are also thankful to a list of nature photographers for granting permission to use their photos: Shantanu Kuveskar and Trisha Shears (Ploceus philippinus); Dick Daniels (Sporopipes frontalis); Nik Borrow (Malimbus malimbicus); Derek Keats (Anaplectes rubriceps and Vidua chalybeata); Jaques Erard (Ploceus nigerrimus); Jean-Marie Gradot (Quelea and Foudia madagascarensis); Dibyendu Ash (Ploceus megharhyncus); Anthony Villaume,IBC1034370 (Quelea cardinalis); Valentina Storti (Plocepasser mahali);Michael Pazzani(Euplectes afer);Mike Barth(Ploceus vitellinus);Shantanu Kuveskar (Ploceus philippinus); Charles J.Sharp (Ploceus cucullatus);Savithri Singh(Ploceus benghalensis);Imran Shah(Ploceus manyar);Apisit Wilaijit (Ploceus hypoxanthus); Juan Emilio (Estrilda astrild).The picture of Passer domesticus was freely available at https://www.pngfuel.com/free-png/rgjrb.Further details of author affiliations and links are available in Appendix Table S4.The authors are indebted also to Prof.Muhammad Arshad (University of Sargodha) for his effort, advice and helpful discussion aimed at improving this work.

Appendix A.Supplementary data

Supplementary data to this article can be found online at https://do i.org/10.1016/j.avrs.2023.100120.