Population structure and phylogenetic analysis based on genomic DNA
We demultiplexed raw reads and applied standard quality filters in Ipyrad 0.7.28 (Eaton and Overcast 2020). We discarded sequences when a single base had a Phred quality score below 10 or more than 5% of bases had a Phred quality score below 20. Additional filtering was applied to only retain reads that did not have enzyme cleavage sites, adaptor sequences or index sequences and included up to 3 ambiguous sites. Finally, we only kept loci that were present in at least 80% of the samples. The ddRAD data is available in Dryad (doi: xxx).
We performed a principal component analysis (PCA), using the gdsfmt and SNPRelate packages (Zheng et al. 2012) in R 3.5.1 (R Core Team, 2018), to visualize possible clusters in the data. We used all the SNPs of the dataset (including multiple SNPs per locus), but because the PCA is sensitive to missing data, SNPs missing from at least one individual were removed. We also assigned individuals to genetic clusters (K) using Structure 2.3.4 (Pritchard et al. 2000). For this analysis we used a single random SNP from each RAD locus. We implemented the admixture ancestry model with correlated allele frequencies and an allele frequency prior of λ = 1. We conducted 10 runs for each value of K = 1-4, and each run consisted in 500,000 generations following a burn-in of 100,000. The most likely value of K was determined following the ΔK method described by Evanno et al. (2005) and implemented in Structure Harvester 0.6.94 (Earl and vonHoldt 2012). We averaged results across the 10 runs using the greedy algorithm in the program CLUMPP 1.1.2 (Jakobsson and Rosenberg 2007) and visualized results using the conStruct package (Bradburd et al. 2018) in R.
We used RAxML version 8.2.4 (Stamatakis 2014) to infer a maximum likelihood phylogenetic tree. We implemented the ASC_GTRGAMMA model and the Lewis correction for ascertainment bias and we conducted 200 bootstrap replicates to assess node support. We included one specimen ofT. caerulescens as outgroup (we obtained ddRADseq data for this specimen from a previous study; Kopuchian et al. 2020; see Table S1).