Sampling
A total of 19 feather samples of the Eurasian Collared Dove used in the present study were collected from the field at seven locations in Saitama Prefecture, Japan, and another 20 feathers were collected from individuals at five zoos in Saitama Prefecture, Japan (Fig. 1). One eggshell of the Eurasian Collared Dove was also obtained from the field in Saitama Prefecture. Wild samples were collected from nests or from the ground within a 3‒4 m radius from just under the nests, and only one sample was collected from each nest to avoid duplication of the same individual. Seven samples were collected in Kuki, six collected in Kawashima, two in Kasukabe and Misato, and one in Kounosu, Satte, and Saitama, respectively. From the zoo, one to nine feather samples were obtained from five zoos. A Ringneck Dove (Streptopelia risoria) sample was also collected from the zoo. All samples were collected between 2017 and 2019, and were sheltered from the sun immediately after collection and were frozen to avoid DNA degradation.
DNA extraction, PCR and sequencing
DNA was extracted from the feathers through the spin-column method (NucleoSpin Tissue Kit; MACHEREY–NAGEL, Düren, Germany). The feathers containing soft tissues were cut into small pieces using sterilized scissors, placed in 200 μL of DNA extraction buffer and incubated overnight at 56 °C, resulting to an eluted final volume of 50 μL. For the eggshell sample, the inside of the eggshell was wiped with a cotton swab soaked in DNA extraction buffer and incubated. A partial mtDNA fragment of COI and the Control Region (CR) was amplified through polymerase chain reaction (PCR) using the following primers specifically designed for this study; Sd-COI-F (5′-GTGACCCTAATCAATCGATG-3′) and Sd-COI-R2 (5′-TATGTAGCCGAAGGGTTC-3′) were used for COI, and SD-CR-F15665 (5′-CCCTGCATCTGTGTCCTATG-3′) and SD-CR-R3 (5′-CATTATTAATGGTTTGTCAGCGG-3′) were used for CR, respectively. The mtDNA control region has been widely applied to the population genetics of various animal species including birds because its mutation rate is approximately 2.8–5 times faster than other sequence segments (Avise et al. 1987; Nagata et al. 1998; Avise 2000).
Each amplification process was conducted using a reaction mixture of 20 μL containing 1 μL template DNA, 10 mM deoxynucleotides (2.5 mM each), 5 nM of each primer, 2 μL buffer, and 1 U of Ex Taq polymerase (TaKaRa, Shiga, Japan). The thermal cycling conditions involved an initial denaturation step at 94 °C for 5 min, 45 cycles of denaturation at 95 °C for 10 s, annealing at 57 °C for COI and 62 °C for CR for 30 s, elongation at 72 °C for 60 s, and final elongation at 72 °C for 2 min. The PCR products were confirmed by 2% agarose gel electrophoresis and purified using Nucleo Spin Gel and PCR Clean-up kit (MACHEREY–NAGEL). The purified product was eluted with 20 μL elution buffer. All the amplified fragments were sequenced with both forward and reverse primers. The PCR products were used as the template for the 10-μL cycle sequencing reactions using BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA). The sequences were generated using a four-capillary 3500 Genetic Analyzer (Applied Biosystems) and the specific primer.
Population genetic analyses
Multiple sequence alignment was then performed using CLUSTAL W (Thompson et al. 1994). The sequence waveform data were checked manually to avoid mistakes due to sequence reactions. Haplotype and nucleotide diversities were calculated based on the 1140-bp fragment of the CR using DnaSP version 6.12 (Rozas et al. 2017). The haplotype network was constructed using TCS 1.21 (Clement et al. 2000). The mtDNA CR sequences were deposited in the DNA Data Bank of Japan (DDBJ), European Molecular Biology Laboratory (EMBL), and GenBank under the accession numbers LC529745–LC529747, and LC530631, respectively.
Mismatch distribution analysis (Rogers and Harpending 1998) was performed with DnaSP 6 and Arlequin 3.5 (Excoffier and Lischer 2010) to test for signatures of recent population expansion. DnaSP 6 was used to develop figures and Arlequin 3.5 was used to calculate the values. Tajima’s D (Tajima 1983) and Fu’s FS tests (Fu 1997) of selective neutrality were performed using Arlequin 3.5. Negative Tajima’s D values suggest that a population has recently increased in size. The fit to a sudden expansion model was judged using the value of the sum of squared deviations (SSD). The lack of significant differences in the SSD values indicated that the population assumed an expansion model. The Harpending’s raggedness index (rg) was used to distinguish between expanded and stationary populations with the observed mismatched distributions (Harpending et al. 1993). The value of rg is typically low and inconsequential in expanding populations but high and significant in stationary populations (Harpending 1994). The phylogenetic relationships among the sequences were traced through the maximum likelihood (ML) method (Yang 1994) with MEGA 7 software (Kumar et al. 2016). The ML analysis was performed using the GTR + I + G model. The model was selected with the highest parameter of maximum likelihood fits using the Find Best-Fit Substitution Model of MEGA 7. A closely related S. chinensis and S. orientalis sample (DDBJ/EMBL/NCBI Accession number: KP636801 and NC031447) were used as the outgroup and the previously reported mtDNA of Chinese Eurasian Collared Dove was used for analysis in this study (KX372273). Their corresponding CR sequences were aligned for the phylogenetic analysis.
