Secondary contact rather than co-existence – Erebia butterflies in the Alps

Secondary contact zones are ideal systems to study the processes that govern the evolution of reproductive barriers, especially at advanced stages of the speciation process. An increase in reproductive isolation resulting from selection against maladaptive hybrids is thought to contribute to reproductive barrier buildup in secondary contact zones. While such processes have been invoked for many systems, it remains unclear to which extent they influence contact zone dynamics in nature. Here, we study a very narrow contact zone between the butterfly species Erebia cassioides and E. tyndarus in the Swiss Alps. We quantified phenotypic traits related to wing shape and reproduction as well as ecology in order to compare the degree of intra- and interspecific differentiation. Even though only very few first-generation hybrids occur, we find no strong indications for current reinforcing selection, suggesting that if reinforcement occurred in our system, it likely operated in the past. Additionally, we show that both species differ less in their ecological niche at the contact zone than elsewhere, which could explain why co-existence between these butterflies may currently not be possible.

Code for linear and linear mixed-effects models, as well as cline analyses can be found in 
Erebia_models_clines_and_plots.R
input file = phenotypics_FINAL.txt for models and clines
input file = Revision_FINAL.eigenvec, Revision_FINAL.eigenval, classifiers.txt for PCA based on genomic data
input file = Admixture_2.Q, Indvdata.txt for processing the results of the k=2 ADMIXTURE analysis.

Code for P matrix analysis (input file = phenotypics_FINAL.txt):
R script calculating theta and implementing a bootstrap approach (from Lucek et al. 2014 JEB) + bootstrap and calculation of the Mahalanobis distance between subsets:
Functions_for_Pmatrices_and_Mahalanobis_distances.R

R script for all P analyses for the allopatric comparisons:
P_matrices_and_Mahalanobis_distances_allopatric.R

R script for all P analyses along the zone of secondary contact:
P_matrices_and_Mahalanobis_distances_contact_zone.R

Output files from the P matrix analyses and Mahalanobis distance calculations:

Output files/wing_allo_pmax.txt                    => Output for P matrix of wing shape of allopatric populations
Output files/mahal_wings.txt                     => Output of Mahalanobis distances between allopatric populations for wing shape
Output files/wing_pmax.txt                        => Output for wing shape along the secondary contact zone
Output files/wings_pmax_eigenvalues.txt            => Eigenvalues of individual p matrices for wing shape along the secondary contact zone
Output files/wing_mahal.txt                     => Output of Mahalanobis distances for wing shape

Output files/orange_allo_pmax.txt                    => Output for P matrix of orange spot of allopatric populations
Output files/mahal_orange.txt                     => Output of Mahalanobis distances between allopatric populations for orange spot
Output files/orange_pmax.txt                        => Output for orange spot along the secondary contact zone
Output files/orange_pmax_eigenvalues.txt            => Eigenvalues of individual p matrices for orange spot along the secondary contact zone
Output files/orange_mahal.txt                     => Output of Mahalanobis distances for orange spot

Output files/genit_allo_pmax.txt                    => Output for P matrix of genital shape of allopatric populations
Output files/mahal_genit.txt                     => Output of Mahalanobis distances between allopatric populations for genital shape
Output files/genit_pmax.txt                    => Output for genital shape along the secondary contact zone
Output files/genit_pmax_eigenvalues.txt            => Eigenvalues of individual p matrices for genital shape along the secondary contact zone
Output files/genit_mahal.txt                     => Output of Mahalanobis distances across the contact zone for genital shape

Code for ecological variable selection and PCA = Erebia_vars_sel_sub.R, input files = phenotypics_FINAL.txt, output files = z_pca_all.txt, dat_niche_sel.txt
Code for niche overlap calculations = Erebia_niche_overlap_sub.R, input files = dat_niche_sel.txt