Secondary contact of recently diverged species may have several outcomes, ranging from rampant hybridization to reinforced reproductive isolation. In plants, selfing tolerance and disjunct reproductive phenology may lead to reproductive isolation at contact zones. However, they can evolve under both allopatric or parapatric frameworks and originate from adaptive and/or neutral forces. Inferring the historical demography of diverging taxa is thus a crucial step to identify those factors that may lead to putative reproductive isolation. We explored various competing hypotheses to account for the rapid divergence of a fir species complex (Abies flinckii - A. religiosa) distributed in ‘sky-islands’ across central Mexico (i.e., along the Trans-Mexican Volcanic Belt; TMVB). Despite co-occurring in two independent sympatric regions (west and center), these taxa rarely interbreed because of disjunct reproductive phenologies. We genotyped 1,147 SNPs, generated by GBS, across 23 populations, and compared multiple demographic scenarios based on the geological history of the TMVB. The best-fitting model revealed one of the most rapid and complete speciation cases for a conifer species-pair, dating back to ~1.2 Ma. Coupled with the lack of support for stepwise colonization, our coalescent inferences point to an early cessation of interspecific gene flow under parapatric speciation; ancestral gene flow during divergence was asymmetrical (mostly from western firs into A. religiosa) and exclusive to the most ancient (i.e., central) contact zone. Factors promoting rapid reproductive isolation should be explored in other slowly-evolving species complexes as they may account for the large tropical and subtropical diversity.

Please, see Materials and Methods.

GBS_Abies contains 11 subdirectories:

-ipyrad_Abies/: This directory has scrips/command lines and vcf files used to identify the optimal parameter values to perform a denovo SNP assembly using ipyrad 0.9.33 and VCFtools 0.1.16.

-out/: This directory contains a vcf file assembled with ipyrad after parameter optimization. This file is further filtered with scripts saved in bin/

-bin/: This directory holds most of the scripts used to perform population genomics computations, and genetic structure inferences.

-meta/: This directory has a metadata file and tables with the geographic/genetic information for each analized sample.

-out_TMVB/: This directory contains three vcf files (each with different level of missing data), which were used to asses the genetic diversity, structure and the past demography of fir populations from central Mexico. Outputs are also stored in this directory.

-out_religiosa/: This directory holds a vcf file with SNPs discovered only for A. religiosa. This files was used to perform a DAPC analysis.

-out_religiosa_NoMan/: This directory contains a vcf file with SNPs discovered only for A. religiosa, except the W genetic pool. We infered the genetic structure for these samples using ADMIXTURE v. 1.3.0.

-out_Megas/: It holds a filtered .vcf file with SNPs discovered from a pair of Megagametophyte samples (i.e., haploid tissue).

-fastsimcoal/: This directory contain scripts and input files (.tpl, .est, .obs) to perform demographic inferences using fastsimcoal2 v. 2.6.0.3-14.10.17

-DnaSP_results/: This directory holds input files to perform population genetic computations using DNA v. 6.12.03. Outputs are also saved.

-data_raw/: This directory was used to store non-demultiplexed samples (i.e., raw reads).