Contrasting Trends of Population Size Change for Two Eurasian Owlet Species—Athene brama and Glaucidium radiatum From South Asia Over the Late Quaternary

Climatic oscillations over the Quaternary have had a lasting impact on species’ distribution, evolutionary history, and genetic composition. Many species show dramatic population size changes coinciding with the last glacial period. However, the extent and direction of change vary across biogeographic regions, species-habitat associations, and species traits. Here we use genomic data to assess population size changes over the late Quaternary using the Pairwise Sequential Markovian Coalescent (PSMC) approach in two Eurasian Owlet species—the Spotted Owlet, Athene brama, and the Jungle Owlet, Glaucidium radiatum. While Spotted Owlets are typically associated with open habitats, Jungle Owlets are found in deciduous forests and scrublands. We find that the effective population size for the Spotted Owlet increased after the Interglacial period till the Last Glacial Maxima and subsequently declined toward the Mid-Holocene. On the other hand, effective population size estimates for the Jungle Owlet increased gradually throughout this period. These observations are in line with climatic niche model-based predictions for range size change for both species from a previous study and suggest that habitat associations at the local scale are important in determining responses to past climatic and vegetational changes. The Spotted Owlet result also aligns well with the expectation of open habitat expansion during the arid Glacial Maxima, whereas for the Jungle Owlet the contrasting expectation does not hold. Therefore, assessing the impacts of glacial history on population trajectories of multiple species with different habitat associations is necessary to understand the impacts of past climate on South Asian taxa.


INTRODUCTION
Quaternary climatic fluctuations, characterized by repeated cycles of cool and arid glacial periods punctuated by warm "interglacials, " are notable for their role in shaping the distribution and evolutionary history of species globally (Hewitt, 2000(Hewitt, , 2004. Northern latitudes experienced far greater extremes with glacial periods being marked by an expansion of ice sheets, retraction of temperate and boreal forests, and expansion of tundra vegetation (Hofreiter and Stewart, 2009;Beyer et al., 2020). Lowered sea levels exposed a greater land area and created land bridges. Though less pronounced, the tropics too experienced changes in temperature and precipitation, with glacial periods being drier and cooler, leading to contraction and fragmentation of tropical forests, expansion of grasslands, or dry woodlands (Qian and Ricklefs, 2001;Ray and Adams, 2001). Coincidental with these shifts in vegetation and available habitat, species ranges changed dynamically through the Quaternary, varying in their extent and geographical location (Dynesius and Jansson, 2000), impacting community composition and species distribution (Fuchs et al., 2019). Along with this, changes in population size and isolation-migration patterns altered the genetic composition of populations (Hewitt, 2000).
In the tropics, there is considerable empirical support for the impact of forest fragmentation and open habitat expansion during the glacial periods on the population history of many species (Schneider and Moritz, 1999;Cabanne et al., 2007;Wurster et al., 2010;Piñeiro et al., 2017;Song et al., 2020). However, data contradicting this general pattern (Batalha-Filho et al., 2012;Rocha and Kaefer, 2019) underscore the importance of understanding regional patterns across taxonomic groups. Among tropical regions, South Asia remains poorly represented in studies that include genetic data (Reddy, 2014), even fewer that examine historical climatic oscillation-induced species range dynamics. Global reconstructions of vegetation indicate that the region saw an expansion of open habitats including grasslands and dry scrubland during the last glacial period (Beyer et al., 2020). Palynological and isotope analysis suggest that monsoonal activity was also reduced during this period while being stronger during the interglacial (Prabhu et al., 2004). These changes are likely to have impacted forest and open-habitat associated species in an opposing manner, but to varying degrees depending on the strengths of association with the habitat.
Here we examine the impact of the late Quaternary climatic events on two phylogenetically close (from the same family-Strigidae) but ecologically differentiated avian species from the Indian subcontinent. The Spotted Owlet, Athene brama and the Jungle Owlet, Glaucidium radiatum, are small owls, widely distributed, and sympatric over large parts of their current distribution. Nevertheless, they have distinct dietary habits and breeding seasons (Mehta et al., 2018), and despite their range overlap, at smaller spatial scales they are segregated by their habitat preferences. Spotted Owlets avoid dense forests and are primarily associated with more open habitats such as cultivated areas, semi-desert areas, and open woodlands (Ali and Ripley, 1983;Holt et al., 2020b). Jungle Owlets are typically found in forested areas, including dense and secondary dry/moist deciduous forests, scrublands, and riparian areas (Ali and Ripley, 1983;Holt et al., 2020a). Therefore, it is likely that suitable habitats for the two species expanded and contracted in opposition to each other, with the glacial periods being favorable for Spotted Owlets and the interglacials for Jungle Owlets. Genomic data allow us to assess such predictions by looking for temporal changes in effective population size (N e ) coinciding with major climatic events.
The Pairwise Sequentially Markovian Coalescent (PSMC) model has frequently been used to infer changes in N e in the context of past climate (Nadachowska-Brzyska et al., 2015;Vijay et al., 2018;Lucena-Perez et al., 2020). This model relies on the rate of coalescence along a single diploid genome and uses this to infer population size over time (Li and Durbin, 2012). The PSMC framework has been particularly useful in cases where acquiring large genomic datasets is still logistically prohibitive (Mays et al., 2018). Based on PSMC curves, several species show dynamic changes in N e during the past glacial cycles, suggesting a strong impact of glacial cycles on population history (Nadachowska-Brzyska et al., 2015;Kozma et al., 2016;Mays et al., 2018;Vijay et al., 2018;Feng et al., 2019). Such analyses have often been used to complement and corroborate evidence from climatic niche models (Kozma et al., 2018;Mays et al., 2018). Moreover, longterm N e trends help understand contemporary genetic variation (Xue et al., 2015) and extinction risk (Hung et al., 2014), thus providing an evolutionary context for current population trends.

OBJECTIVES
Here we test the hypothesis that the demographic history of species will be strongly affected by the ecological differences, specifically habitat preferences. We expect to see an increase in N e from LIG-LGM followed by a decline in the LGM-MDH period for the Spotted Owlet. We expect the reverse trend for the Jungle Owlet, with contraction in N e during the glacial period and expansion toward the Holocene. However, a previous study using climatic niche models suggests that the extent of Jungle Owlet habitat was reduced (relative to the Spotted Owlet) during the LIG and gradually increased toward the MDH (Koparde et al., 2019). We use the PSMC-based reconstruction of past N e from whole-genome sequence data to investigate past population size dynamics of the Spotted Owlet and the Jungle Owlet over the late Quaternary.

Sample Collection
We obtained one sample per species from Central (Jungle Owlet) and Southern (Spotted Owlet) Western Ghats, with due permits from the Forest Department and extracted DNA using the DNeasy Blood and Tissue kit following the manufacturer's protocol with slight modifications. Whole DNA concentration was quantified using a Qubit 4.0 fluorometer and DNA integrity was visualized using 1% Agarose Gel electrophoresis and Tapestation. Samples were sequenced on an Illumina HiSeq 10X platform (150 bp paired-end), targeting 30x genome coverage.

Read Mapping and Variant Calling
Raw sequences were inspected for overall quality, adapter content, and the number of reads to ensure that it met the required quality standards using the program TrimGalore! 1 .
Reads were processed using the GATK short reads pipeline (McKenna et al., 2010) through which read-group information was added and Illumina adapters were marked. For both study species, no reference assembly was available. Therefore, the Burrowing Owl (Athene cunicularia) assembly was used as a reference genome (RAG1 nuclear gene p-distance-the proportion of nucleotide sites at which the compared sequences are different-with A. brama 1.07%; with G. radiatum 2.25%; Supplementary Table S1.). The unmasked assembly (athCun1) was downloaded from Ensembl (Yates et al., 2020). Reads were mapped against the reference using the Burrows-Wheeler Aligner (Li and Durbin, 2009;Li, 2013) with the "mem" option under default settings. PCR duplicates were marked using Picard Tools 2 . At this stage, the average depth and genome coverage for each sample were summarized using the "Depth Of Coverage" option from GATK as a quality check before proceeding further. The Burrowing Owl assembly has an accompanying sequence report which lists chromosomal information for each scaffold and identifies them as occurring on autosomes, the Z chromosome, or being unplaced. Therefore, reads mapping to the 51 autosomal scaffolds from 28 chromosomes were selected. The scaffold N50 of the reference genome is 42,147,404 bp (42 Mb) and the shortest scaffold used was 121,167 bp (0.12 Mb) long, making it suitable for use with PSMC (Patton et al., 2019). A bed file identifying autosomal scaffolds was created and used to extract reads aligning only to autosomes using the program samtools (Li, 2011) with the view option, along with -b -L flags. The aligned file was filtered by setting a base quality threshold (q 20). Variants were identified using the bcftools (Li, 2011) mpileup (with the -C50 flag) and call (with the -c flag) options and they were converted to a consensus file with depth filters (minimum 10, maximum 100), as suggested in the PSMC manual (Li and Durbin, 2012). The bcftools/samtools variant-calling pipeline was used, as suggested in the PSMC manual. Genome coverage was calculated using samtools depth for different depth filters to see how the filtering impacted the data.

Demographic History Through PSMC
The consensus file was converted to a suitable input format for the PSMC pipeline following instructions in the PSMC manual. PSMC analysis was repeated with different values of the "-t" and "-r" parameters, controlling the maximum time to the most recent common ancestor (TMRCA) and the number of free atomic time intervals, respectively. For each value of "-t" (starting with 15), "-p" was modified till each time span was supported by at least 20 events after 25 iterations. The final set of parameters was chosen to be "-N25 -t9 -r5 -p '26 * 2 + 4 + 7 + 1'." Bootstrapping with 100 replicates was performed with the identified set of parameters using a block size of 0.5 Mb. The output was written to a text format (using the -R option with the psmc_plot.pl script) for plotting in R. Estimates of the mutation rate and the generation time are required to convert the scaled population size (lambda) to effective population size (N e ) and the time (in generations) to units of years. As the mutation rate for the study species is unknown, it was assumed to be 4 × 10 −9 , 2 http://broadinstitute.github.io/picard/ based on an average estimate of mutation rates from chickens to passerines (Mueller et al., 2018). Nevertheless, additional plots with alternative values of mutation rate based on estimates from other species were also made for comparison. Generation time was assumed to be 2 years, calculated as the age of sexual maturity/reproduction (1 year; Pande et al., 2007) multiplied by 2, as suggested by Nadachowska-Brzyska et al. (2015). The scripts used have been uploaded to GitHub 3 .

RESULTS
The Jungle Owlet and Spotted Owlet samples yielded 52.8 and 44.5 Giga bases, respectively (from 150 bp paired-end reads). After alignment, the Jungle Owlet sample had an average depth of 30X and the Spotted Owlet 26.5X. Genome coverage, depth filters, and average depth are expected to impact the outcome of PSMC and for reliable results, 75% of genome coverage with a 10X depth filter is recommended (Nadachowska-Brzyska et al., 2016). The genome coverage for both samples decreased with quality filtering and removal of non-autosomal scaffolds (Supplementary Figure S1) but was within the recommended limit (Nadachowska-Brzyska et al., 2016).
The parameters chosen for PSMC resulted in a time-span covering the LIG, the LGM, and the MDH i.e., the period of interest. Altering the mutation rate and the generation time based on other studies on owls and other birds (Nadachowska-Brzyska et al., 2013;Fujito, 2020) had predictable effects on the PSMC curve, as seen in other studies (Nadachowska-Brzyska et al., 2015;Vijay et al., 2018), but did not alter the shape (Supplementary Figure S2).
The two species show contrasting trends of population size (N e ) change over the past 500,000 years (Figure 1). The N e of the Jungle Owlet rose gradually from approximately 14,500 during the LIG (120,000 YA) to about 78,000 during the LGM (22,000 YA) and eventually to about 160,000 at the MDH (6,000 YA). In contrast, the Spotted Owlet had the highest N e , about 293,000, around the LGM. It rose from about 33,000 around the LIG and again fell after the LGM to about 45,000 around the MDH. It must be kept in mind that PSMC results are less robust on recent timescales (Li and Durbin, 2012;Patton et al., 2019) and this affects the N e estimates around the Mid-Holocene. However, the population trends leading up from the LGM to the MDH still indicate divergent trajectories for the Jungle and Spotted Owlets. PSMC plots for the two owlets without scaling are shown in Supplementary Figure S3 and show the same contrasting pattern.

DISCUSSION
In this study, we elucidate population size trends for the Jungle Owlet and Spotted Owlet over the major climatic events of the past-the LIG-beginning of the Last Glacial Period, the LGM, and the MDH. We find that both species show expansion from the LIG to the LGM, although the Jungle Owlet has a much greater degree of expansion. While the Jungle Owlet continues to expand into the MDH, the Spotted Owlet shows a declining trend after the LGM. The direction of these population trends is broadly consistent with predictions from the species distribution models from a previous study (Koparde et al., 2019), suggesting that their distinct niches had a significant impact on their population history. We note that we have sampled one individual of each species and our data may not be representative of trends across the range of these species. Therefore, we restrict our inferences to the sampled Western Ghats populations of these widely distributed species.

Temporal Trends in N e Over the Late Quaternary for the Two Owlet Species
The LIG was a warm and wet period in India (Prabhu et al., 2004), with a strong southwest monsoon. Both samples show reduced N e during this period, a pattern observed in multiple species from East-Asia (Dong et al., 2017). The LGM was characterized by a cold and dry climate, increased aridification, and reduced precipitation (Prabhu et al., 2004). In the Indian peninsula, the rain forests were replaced by grassland vegetation (Sukumar et al., 1995). These shifts are likely to have been favorable for the Spotted Owlet, which shows a considerable increase in population size toward the LGM. The Jungle Owlet also shows an increasing trend, albeit to a much lower degree, possibly due to the species' preference for moist-dry forest and scrublands in potential refugia. Following this, the transition toward the Holocene climate resulted in a rise in temperature and precipitation. During this time, much of peninsular India shifted to wetter, forest habitats (Prabhu et al., 2004). The Spotted Owlet shows a decline in this period, consistent with the prediction of retraction of grasslands and open habitats, while the Jungle Owlet continues to show an upward trend.
One major concern with the PSMC model is the presence of population structure which has a confounding effect on the inference of population size changes (Mazet et al., 2015(Mazet et al., , 2016. For this reason, it is also possible that the population decline seen in the Spotted Owlet after the LGM may also be due to changes in population structure, and this needs to be tested with more data. Disentangling the effects of population structure and demography remains challenging, particularly with sequentially Markovian Coalescent methods (Mather et al., 2020). Another technical constraint is the absence of a reference genome assembly of the study species, which would have allowed better mapping of reads. Nevertheless, sampling more individuals and inclusion of more populations from across the range, while not logistically feasible for this study, may provide more conclusive results in the future.

The Influence of the Last Glacial Period on Species' History in the Tropics
The Last Glacial Period has been a time of major population size changes for many avian species (Nadachowska-Brzyska et al., 2015). However, the direction and extent of the impact may differ across species within the same region (Cabanne et al., 2016) and within species across different regions (Nadachowska-Brzyska et al., 2016). Moreover, the impacts of the Last Glacial Period were variable within the tropics. While African forests retracted into refugial fragments in favor of grassland expansion during the dry, cool LGM period, Amazonian forests may have been relatively stable, even though the forest composition itself was altered (Prentice et al., 2000;Lessa et al., 2003;Malhi and Phillips, 2004;Bush and de Oliveira, 2006;Piñeiro et al., 2017). A study on 11 East-Asian birds, on the other hand, shows a greater impact of the LIG, where several species show range contraction, as opposed to the glacial period, which is a period of expansion (Dong et al., 2017).
In the Indian subcontinent, changes in diversification rates and population history have been shown to coincide with paleoclimatic shifts (Robin et al., 2010;Agarwal and Ramakrishnan, 2017). A high-altitude forest specialist bird, Sholicola albiventris, shows post-glacial population expansion, aligning well with the idea of forest fragmentation during the LGM (Robin et al., 2010). On the other hand, the Indian peafowl Pavo cristatus-a generalist (occupying moist-dry deciduous forests, to open forests and woodlands, and cultivated areas and orchards), shows a stable population size throughout the Last Glacial Period, despite earlier bottlenecks (Jaiswal et al., 2018). The Spotted Owlet in our study, in line with observations for many other species, shows a cyclic expansion and contraction of population size around the LGM, whereas the Jungle Owlet does not. Such differences could result from differences in the extent of habitat change and the strength of the habitat association. Moreover, our Jungle Owlet sample is from the more widely distributed subspecies Glaucidium radiatum radiatum, collected from central India. The subspecies from Southern India-Glaucidium radiatum malabaricum, known to be associated with more moist forests (Ali and Ripley, 1983), may show a population history more consistent with other forest-associated species.

N e Trends in Understanding Long-Term Population History
The demographic reconstruction of population history of the two Owlet species appears to be broadly consistent with predictions from previously generated climatic niche models (Koparde et al., 2019) in terms of the direction of change. PSMC-based demographic inference coupled with Climatic niche models can be an effective means of understanding the population history of species (Kozma et al., 2018;Mays et al., 2018). Apart from indicating past population size trends, they are an important indicator of species vulnerability. Long-term declines have been observed in species that are now of conservation concern (Nadachowska- Brzyska et al., 2015) and are often considered to indicate reduced genetic variation and increased genetic load (Xue et al., 2015). While inferences from N e curves alone are valuable, incorporating niche models allows correlation with past climate and vegetation. Moreover, population trends may not be a simple reflection of past climatic and vegetation changes. Other factors such as the dispersal ability of a species or the rate of change of suitable habitat may interact with the prevailing climatic and vegetation conditions to modify the genetic consequences of range contraction (Arenas et al., 2012). Therefore, genomic data can be highly valuable when combined with climatic niche models to test predictions about the impact of past climatic changes on species history.
The two owlet species in our study show evidence for differential impacts of Quaternary climatic events, likely driven by their differing ecological niches, despite the phylogenetic similarity. Given our limited sampling, we highlight the need for a larger investigation of the impacts of past climate on South Asian taxa to understand the community and species-specific impacts of past climate change.

DATA AVAILABILITY STATEMENT
The raw illumina sequences for the two species have been submitted to NCBI and are available from the bioproject PRJNA665335.

ETHICS STATEMENT
Permit requests for this animal study were reviewed, approved, and issued by the State Forest Departments of Maharashtra and Kerala.

ACKNOWLEDGMENTS
We are grateful to the State Forest Departments of Maharashtra and Kerala for permits (permit letter numbers Desk-22(8)WL/CR-59(17-18)631/2019-20 and WL10-16346/2017, respectively); the Scientific Computing facility at IISER Tirupati and members of IT Department for HPC access and help with trouble-shooting; The Director, along with the Library, Administrative and Finance staff of SACON for logistic support; Naman Goyal, Ashwin Varudkar, and C. K. Vishnudas for help with sample collection. We are also grateful to the Ecology Evolution labs at IISER Tirupati for feedback on the manuscript.