Computational protocol: Genetic and morphological support for possible sympatric origin of fish from subterranean habitats

Similar protocols

Protocol publication

[…] DNA was extracted using the Chelex 100/200 method and the salt extraction method. The primers FishF1 and FishR1 in Ward, et al. were used for COI amplification. The PCR reactions were performed in 25-µl volumes containing 18.5 µl H2O, 2.5 µl 10 × buffer, 0.5 µl MgCl2 (50 mM), 0.5 µl of each primer (10 mM), 0.5 µl dNTP, 0.5 µl Taq DNA polymerase, and 2 µl DNA solution. Thermal cycles included 1 cycle of 94 °C for 5 min; 35 cycles of 94 °C for 1 min, 58 °C for 1 min, and 72 °C for 1 min; and a final extension at 72 °C for 5 min. PCR products were sequenced using a forward primer and a Prisma 3130 sequencer following the protocols provided by the manufacturer ( The sequences were aligned and edited visually using Bioedit. [...] Sequencing adapters were removed using cutadapt. Sequence quality of the first 10,000,000 reads was assessed using FastQC. Libraries were de-multiplexed using process_radtags in STACKS V.1.35. Reads were trimmed to 80 bp and shorter reads were discarded. We used the STACKS v1.35 analysis pipeline to score genotypes at 51,836 SNPs for our samples. Then these results were filtered using the scripts included in stacks_workflow ( The filtration parameters were: -l 0; -I 8; -p 70; -a 0.05; -A 0.05; -H 0.5; -f −0.5; -F 0.6; -s 10. A final set of 11,257 SNPs located on 7,048 reads were retained following the filtering steps.The sequences of the 7,048 loci were concatenated and, finally, a 563,840-bp sequence per individual was produced for phylogenetic analyses. When heterozygous, the different SNPs were named using IUPAC symbols. On average, 6% of locus information was missing per individual, and missing alleles were imputed by the most frequent allele in each species for each locus.To infer the historical demography based on Joint Allele Frequency Spectrum (JAFS) (see below), a dataset was prepared from the original VCF file. Several filtering steps – aimed at removing miscalled and low-quality SNPs, as well as false variation induced by merging paralogous loci – were performed using VCFtools v0.1.13. SNPs with more than 90% of missing genotypes in all individuals were removed, but a lower exclusion threshold (50%) was applied to the out-group to retain a maximum of orthologous loci. After filtering for Hardy-Weinberg disequilibrium for each population (p-value exclusion threshold of 0.01), the filtered datasets were merged. Finally, the most parsimonious ancestral allelic state was determined by keeping monomorphic loci in the out-group, but polymorphic in the complex G. typhlops-G. lorestanensis, aiming to infer the divergence between species, of the studied complex. These result in 5,890 oriented SNPs used to build the unfolded JAFS. [...] The admixture proportions among samples were inferred using the Bayesian clustering method implemented in the program STRUCTURE V.2.3.4. The structure was evaluated for K = 1–5 using admixture model with correlated allele frequencies. The MCMC chains were ran for 100,0000 generations. The support for different values of K was assessed from the likelihood distribution (lowest cross-validation error) and visual inspection of the co-ancestry values for each individual. In addition, two supplementary K-means clustering analyses, the Bayesian Information Criterion (BIC; Schwarz) and the Calinski–Harabasz pseudo-F-statistic, were performed on individuals using the GENODIVE v.2.0b25 program. For these K-mean clustering analyses, a simulated annealing method was used, where the optimal K value was determined via checking K values ranging from 1 to 5 for 5,000 permutations. A PCA implemented in the ade4 package was performed and the first three principal components were visualized using ggplot2 package available in R. [...] The mtDNA and nucDNA phylogenetic trees were reconstructed using Neighbour-Joining (NJ) and Maximum Parsimony (MP) methods using MEGA7 and Maximum Likelihood (ML) using RaxMLGIU 1.5b2. The best-fit model(s) of mtDNA sequence evolution were selected using the online ModelTest in the HIV sequence database ( As there may be incongruence between the gene trees and species trees due to different factors like incomplete lineage sorting (ILS),, species tree was reconstructed for the sequences of 5,843 loci, using SVDquartets + PAUP* (implemented in PAUP*4) software. SVDquartets inference of the species tree was performed using the multispecies coalescent tree model and QFM algorithm of quartet assembly. The branch supports in SVDqurtets method were calculated via implementing 1000 bootstrap repeats. Garra mondica was also included as an out-group. The K2P sequence distances for both data types were calculated using MEGA7. As the net molecular divergence (Da) has been reported to be a predictor of ongoing gene flow, this distance was also calculated using formula () (Camille Roux personal communication).2Da=(Δs−(piA+piB2))/nwhere ∆s is the average number of pairwise differences between sequences from species A and species B and Pi is the nucleotide diversity in each species, and n is the length of the concatenated sequence. The nucleotide differences and diversity indices used for the calculation of Da were calculated using MEGA7. […]

Pipeline specifications