Aim: The Neotropics constitute the most biodiverse region of the world, yet its patterns of diversification and speciation differ among Neotropical areas and are not equally well understood. Particularly, avian evolutionary processes are understudied in the open habitats of temperate South America, where the role of glacial cycles is not clear. We analyzed the evolutionary history of a Neotropical widespread bird species as a case study to evaluate its continental-scale patterns and processes of diversification, with a focus on Patagonia.

Location: Open habitats of the Neotropics.

Taxon: Vanellus chilensis (Aves, Charadriiformes).

Methods: We obtained reduced representation genomic and mitochondrial data from the four subspecies of V. chilensis to perform a phylogenetic/phylogeographic analysis and study the evolutionary history of the species. We complemented these analyses with the study of vocalizations, a reproductive signal in birds.

Results: The initial diversification event within V. chilensis, approximately 600,000 years ago, split a Patagonian lineage from one containing individuals from the rest of the Neotropics. We found considerable gene flow between these two lineages and a contact zone in northern Patagonia and showed that genomic admixture extends to northwestern Argentina. Shallower divergence was detected between the two non-Patagonian subspecies, which are separated by the Amazon River. Vocalizations were significantly different between the two main lineages and were intermediate in their temporal and frequency characteristics in the contact zone.

Main conclusions: Patagonian populations of V. chilensis are clearly differentiated from those of the rest of the Neotropics, possibly as a consequence of Pleistocene glaciations. A secondary contact zone in northern Patagonia with extensive gene flow among lineages appears to be the consequence of post-glacial, northward expansion of the Patagonian populations. Future analyses focused on the dynamics of the contact zone will allow us to establish whether the species continues to diverge or is homogenizing. 

We performed double-digest restriction site-associated DNA sequencing (ddRADseq; Peterson et al., 2012) following the protocol detailed in Thrasher et al. (2018). Genomic DNA from fresh tissue and blood samples was isolated using the DNeasy tissue extraction kit (Qiagen, CA, USA). For each individual, we isolated ~500 ng of genomic DNA at a standardized concentration of 20 ng/µl (quantified using a Qubit fluorometer; Life Technologies, NY, USA). DNA was digested with the restriction enzymes SbfI High Fidelity (8 base bp recognition site; 5’-CCTGCAGG-3’) and MspI (4 base bp recognition site; 5’-CCGG-3’) (New England BioLabs, MA, USA). The 5' and 3' ends of the digested genomic DNA were ligated to P1 and P2 adapters, respectively. Sequences of these adapters are available in Peterson et al. (2012). Samples with unique P1 barcodes were pooled into different indexing groups postdigestion/ligation and cleaned up using 1.59 Agencourt AMPure XP beads (Beckman Coulter, CA, USA). Samples were ordered in these groups according to their DNA concentration, independently of their geographic location or taxonomy. We size-selected fragments between 450 and 700 bp using Blue Pippin (Sage Science, MA, USA). To incorporate the full Illumina TruSeq primer sequences and unique indexing primers into each library, we performed low cycle number PCR with Phusion High-Fidelity DNA Polymerase (New England BioLabs), with the following thermocycling profile: 98°C for 30 sec followed by 11 cycles at 98°C for 5 sec, 60°C for 25 sec, and 72°C for 10 sec with a final extension at 72°C for 5 min. We visualized the product of this amplification on a 1% agarose gel and performed a final 0.73 AMPure to eliminate DNA fragments smaller than 200 bp. Finally, the index groups were combined at equimolar ratios and sequenced on an Illumina HiSeq 2500 lane at the Cornell University Biotechnology Resource Center, obtaining single-end 150-bp sequences.

We demultiplexed raw reads and applied standard quality filters in ipyrad 0.7.28 (Eaton & Overcast, 2020). We only exported loci that were present in at least 80% of the samples to generate the final nuclear DNA dataset used for analyses.