Bayesian discrete phylodynamics and phylogeography analyses

All sequences were aligned using MUSCLE implemented in MEGA-X (Kumar et al. 2018). Regression of root-to-tip distances of sequences was performed using TempEst software (v.1.5.3) (Rambaut et al. 2016). After removing sequences that affected the temporal signal, we used RDP4 (v.4.101) to confirm the absence of recombination signals in the dataset (Martin et al. 2015). The best-fit HKY + I + G nucleotide substitution model was determined based on the value of LnL in jModelTest (v.2.1.10) (Darriba et al. 2012). Time-scaled phylogenetic trees were constructed in BEAST (v.1.8.4) using tip dates (Suchard et al. 2018), a strict clock model, and a constant size coalescent model with a Markov chain Monte Carlo (MCMC) sample chain (10⁸ steps with sampling every 1,000 steps). The convergence of parameters was tested using Tracer (v.1.6) software (Rambaut et al. 2018), using an effective sample size (ESS) of greater than 200 as an acceptance criterion. A maximum clade credibility tree was constructed using Treeannotator (v.1.8.4) (Drummond and Rambaut 2007), with burn-in of the first 10% of samples, and visualized with FigTree (v.1.4.4) (http://tree.bio.ed.ac.uk/software/figtree). The discrete phylogeography analysis was performed in SpreaD3 (v.0.9.6) software (Bielejec et al. 2016). After converting to a keyhole markup language file, we calculated Bayes factors (BFs) in SpreaD3 using a Bayesian stochastic search variable selection (BSSVS) file to obtain statistically significant migration routes. The reliability of the analysis was verified by running the BSSVS independently 3 times, and posterior probability and BF cutoffs were used to define significance.