RESEARCH ARTICLE Open AccessPhylogeny with introgression inHabronattus jumping spiders (Araneae:Salticidae)Geneviève Leduc-Robert1 and Wayne P. Maddison1,2*AbstractBackground: Habronattus is a diverse clade of jumping spiders with complex courtship displays and repeatedevolution of Y chromosomes. A well-resolved species phylogeny would provide an important framework to studythese traits, but has not yet been achieved, in part because the few genes available in past studies gave conflictingsignals. Such discordant gene trees could be the result of incomplete lineage sorting (ILS) in recently diverged partsof the phylogeny, but there are indications that introgression could be a source of conflict.Results: To infer Habronattus phylogeny and investigate the cause of gene tree discordance, we assembledtranscriptomes for 34 Habronattus species and 2 outgroups. The concatenated 2.41 Mb of nuclear data (1877 loci)resolved phylogeny by Maximum Likelihood (ML) with high bootstrap support (95-100%) at most nodes, with someuncertainty surrounding the relationships of H. icenoglei, H. cambridgei, H. oregonensis, and Pellenes canadensis. Speciestree analyses by ASTRAL and SVDQuartets gave almost completely congruent results. Several nodes in the MLphylogeny from 12.33 kb of mitochondrial data are incongruent with the nuclear phylogeny and indicate possiblemitochondrial introgression: the internal relationships of the americanus and the coecatus groups, the relationshipbetween the altanus, decorus, banksi, and americanus group, and between H. clypeatus and the coecatus group. Todetermine the relative contributions of ILS and introgression, we analyzed gene tree discordance for nuclear loci longerthan 1 kb using Bayesian Concordance Analysis (BCA) for the americanus group (679 loci) and the VCCR clade (viridipes/clypeatus/coecatus/roberti groups) (517 loci) and found signals of introgression in both. Finally, we tested specifically forintrogression in the concatenated nuclear matrix with Patterson’s D statistics and DFOIL. We found nuclear introgressionresulting in substantial admixture between americanus group species, between H. roberti and the clypeatus group, andbetween the clypeatus and coecatus groups.Conclusions: Our results indicate that the phylogenetic history of Habronattus is predominantly a diverging tree, butthat hybridization may have been common between phylogenetically distant species, especially in subgroups withcomplex courtship displays.Keywords: Phylogeny, Introgression, Hybridization, Transcriptome, Jumping spiders, Salticidae, Salticinae, Plexippini,Harmochirina, Habronattus* Correspondence: wayne.maddison@ubc.ca1Department of Zoology, University of British Columbia, Vancouver, BC V6T1Z4, Canada2Department of Botany and Beaty Biodiversity Museum, University of BritishColumbia, Vancouver, BC V6T 1Z4, Canada© The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, andreproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link tothe Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.Leduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 https://doi.org/10.1186/s12862-018-1137-xBackgroundThe visually-acute jumping spiders include Habronattus,a clade of about 100 described species whose colourfuland diverse courtship displays are among the most com-plex found in arthropods [1–6]. The group is emerging asa model to study the role of sexual selection in divergence[4, 7–11], the evolution of sex chromosomes [12, 13], andarthropod visual systems [14–16].Although these studies of characters and diversifica-tion in Habronattus have been guided by our knowledgeof the group’s phylogeny, their ability to clearly traceevolutionary processes in this densely diverse group hasbeen limited by poor phylogenetic resolution. Severalsubclades were recognized by distinctive structures andbehaviours by Griswold [3] and confirmed by moleculardata from two gene regions (nuclear Ef1-α and mito-chondrial 16S-ND1, [5]). These include a large clade of42 species whose males have fringed first legs and modi-fied third legs (here called the VCCR clade, subdividedinto the viridipes, coecatus, clypeatus, and robertigroups), the americanus group, the dorotheae group,and several groups of robust-bodied, often shrub-dwelling species (agilis, amicus, tranquillus groups, col-lectively named here the AAT clade). Other well-supported groups, such as the decorus group, emergedonly with molecular data [5]. However, relationshipswithin and among these species groups are, for the mostpart, little resolved (see, e.g., the conservative tree of Fig-ure 4 in [12]).Difficulties in resolving Habronattus phylogeny maystem from the recency of its diversification — possiblyless than 5 million years (Figure 8 in [17]) — yielding in-sufficient sequence divergence or high rates of incom-plete lineage sorting (ILS, [18]). Maddison and Hedin’s[5] gene trees show several cases of distinct morphospe-cies intermingling when represented by multiple speci-mens (e.g., H. pyrrithrix, H. virgulatus, and othercoecatus group members). In addition, the recency oftheir divergence may leave Habronattus species suscep-tible to hybridization [19, 20], leading to introgressionand thus discordant phylogenetic signals.Indeed, signals of introgression have been noted inmitochondrial data between sympatric and closely re-lated species [21] and even between distant specieswhose courtship ornaments and genitalia differ markedly[5]. Hints are also seen in sexually selected ornaments inmales from divergent H. pugillis populations, whose pat-terns of convergence could be explained by introgression[4]. Hybrid zones are known among several pairs ofclosely-related species ([36], unpublished observations).The possibility of introgression is supported by behav-ioural studies, which have found that females from dif-ferent populations of H. pugillis have preferences forforeign males with divergent courtship displays, apossible result of antagonistic coevolution between thesexes [7]. If this has happened throughout Habronattus,then we may be faced with an unfortunate irony: thevery processes of sexual selection that make this groupso tempting to study may at the same time obscure thephylogenetic history we depend on for comparativeanalyses.To whatever extent a divergent phylogeny exists in thegroup, our goal here is to use genomic data to resolve it.We also seek to determine whether the previously-inferred mitochondrial introgression in Habronattusstands alone as in other taxa (e.g., [23]) or is accompan-ied by nuclear introgression. Such introgression coulddo more than confuse phylogeneticists; it could have in-troduced new genetic variation at a rate faster than pos-sible by mutation alone, leading to the sharing ofadaptive loci across lineages and facilitating diversifica-tion [24–26]. Although introgression of courtship traitsinto established systems of mate choice with such elab-orate signals might seem difficult, sensory bias, Fisherianrunaway, and antagonistic coevolution models could allpromote this dynamic [7, 9, 27, 28]. Determining howmuch genetic discordance in Habronattus can be attrib-uted to introgression may provide crucial insights intowhether hybridization may have had any substantial cre-ative role in the group’s diversification.The importance of introgression in animal evolution isuncertain [29], in part because distinguishing discordantsignals resulting from ILS and introgression is difficultwithout phylogenomic datasets [30–33]. We collectedtranscriptome data for 34 Habronattus species and twooutgroups, the first genomic dataset assembled for salti-cid spiders. With a well-resolved phylogenetic tree, wewere then able to use a phylogenetic approach tocharacterize nuclear and mitochondrial introgression inthe group. We focused on the americanus group and theVCCR clade, two of the most diverse clades within Hab-ronattus. To investigate nuclear introgression, we firstconducted a Bayesian Concordance Analysis [33, 34] toinvestigate discordance (caused by either ILS and/orintrogression) in gene trees, and we then applied Patter-son’s D statistic [32] and DFOIL [22] tests to explicitlytest allele patterns for introgression. We were able to re-solve most nodes of the phylogeny with high supportand identified several instances of hybridization in thegroup.MethodsTaxon samplingWe sampled a total of 36 species, including representa-tives from most major clades within Habronattus, and 2outgroups (Table 1). We sampled more extensively inspecies-rich groups (i.e., viridipes/clypeatus/coecatusgroup, and the americanus group), and prioritized theLeduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 2 of 23Table 1 Specimens of Habronattus and outgroups sequencedSpecies Voucher Locality, with latitude and longitudeOutgroupsEvarcha proszynskii Marusik & Logunov, 1998 GLR135 ♂ Canada: BC: Mission 49.166, − 122.409Pellenes canadensis Maddison 017 GLR106 ♂ Canada: BC: Mt. Baldy Road 49.1135, − 119.2103AAT CladeHabronattus conjunctus (Banks, 1898) GLR234 ♂ U.S.A.: AZ: Madera Canyon 31.7417, −110.8847H. hirsutus (Peckham & Peckham, 1888) GLR080 ♂ Canada: BC: Mt. Kobau Road 49.095, −119.610H. signatus (Banks, 1900) – – U.S.A.: CA: Ocotillo 32.7421, −115.9949H. ustulatus (Griswold, 1979) – – U.S.A.: CA: Boulder Oaks 32.7302, −116.4607americanus groupH. aestus Maddison 2017 GLR287 ♀ México: Sonora: Puerto Peñasco 31.418, −113.626H. americanus (Keyserling, 1885) GLR014 ♂ Canada: BC: Iona Beach 49.221, −123.214H. ophrys Griswold, 1987 GLR023 ♂ Canada: BC: Iona Beach 49.221, −123.214H. ophrys GLR015 ♀ Canada: BC: Iona Beach 49.221, −123.214H. sansoni (Emerton, 1915) GLR066 ♂ Canada: BC: Kelowna 49.954, −119.398H. tarsalis (Banks, 1904) GLR297 ♂ U.S.A.: AZ: Yuma 32.731, − 114.612DTB cladeH. altanus (Gertsch, 1934) GLR180 ♂ Canada: AB: Smoky Lake 54.112, −112.198H. chamela Maddison 2017 GLR352 ♂ México: Jalisco: Chamela 19.5038, −105.0334H. decorus (Blackwall, 1846) GLR132 ♂ Canada: BC: Mission 49.166, −122.409H. zapotecanus Griswold, 1987 GLR339 ♂ México: Jalisco: Chamela 19.5316, −105.0707roberti groupH. roberti Maddison 2017 JAL14-9281 ♂ México: Jalisco: Chamela 19.496, −105.042viridipes groupH. calcaratus maddisoni Griswold, 1987 GLR321 ♂ Canada: ON: Haileybury 47.45, −79.708H. jucundus (Peckham & Peckham, 1909) GLR320 ♂ U.S.A.: OR: Bolan Lake, 42.024, −123.461clypeatus groupH. aztecanus (Banks, 1898) GLR347 ♂ México: Jalisco: Puerto Vallarta 20.670, −105.274H. clypeatus (Banks, 1895) GLR227 ♂ U.S.A.: AZ: Mt. Hopkins Road 31.686, −110.975H. gilaensis Maddison & Maddison 2016 AS56 ♀ U.S.A.: New Mexico: Silver Citycoecatus groupH. borealis (Banks, 1895) GLR040 ♂ Canada: ON: Burlington 43.33, −79.8H. captiosus (Gertsch, 1934) GLR356 ♀ Canada: AB: Guy 55.4505, −117.1440H. empyrus Maddison 2017 GLR282 ♂ México: Sonora: Puerto Peñasco 31.293, −113.452H. festus (Peckham & Peckham, 1901) GLR094 ♂ Canada: BC: Hayne’s Lease 49.0813, −119.5181H. festus GLR088 ♀ Canada: BC: Hayne’s Lease 49.0813, −119.5181H. mexicanus (Peckham & Peckham, 1896) GLR353 ♂ México: Jalisco: El Tuito 20.341, −105.350H. pyrrithrix (Chamberlin, 1924) GLR304 ♂ U.S.A.: AZ: Yuma 32.731, − 114.612H. virgulatus Griswold, 1987 GLR205 ♂ U.S.A.: AZ: Mt. Hopkins Road 31.689, −110.975Other HabronattusH. cambridgei Bryant, 1948 GLR351 ♂ México: Jalisco: Puerto Vallarta 20.670, −105.274H. geronimoi Griswold, 1987 GLR267 ♂ U.S.A.: AZ: Miller Canyon 31.416, −110.276H. hallani (Richman, 1973) GLR209 ♂ U.S.A.: AZ: Arivaca 31.668, −111.245H. icenoglei (Griswold, 1979) GLR283 ♂ México: Sonora: Puerto Peñasco 31.273, −113.361H. luminosus Maddison 2017 GLR218 ♀ U.S.A.: AZ: Mt. Hopkins Road 31.6759, −110.9289H. oregonensis (Peckham & Peckham, 1888) GLR149 ♂ Canada: BC: Squamish 49.8465, −123.1452H. paratus (Peckham & Peckham, 1896) GLR363 ♂ Panama: Isla Colon 9.40376, −79.8635H. pugillis Griswold, 1987 GLR236 ♂ U.S.A.: AZ: Mt. Hopkins Road 31.689, −110.975Leduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 3 of 23coecatus group because of possible introgression [5].This is the first phylogenetic analysis to include the spe-cies H. aestus, H. empyrus, and H. luminosus [35]. Someof our sampled species were studied by Maddison andHedin [5] under different names (H. chamela, H. gilaen-sis, H. roberti; see [35, 36]). The outgroup Pellenes cana-densis is closely related to Habronattus within thesubtribe Harmochirina, while Evarcha proszynskii ismore distantly related, in the neighboring subtribe Plex-ippina [17, 37].Specimens were collected from 2012 to 2014 from thelocations listed in Table 1, following institutional and gov-ernmental regulations. Permits for specimens fromMexico were granted through the collaboration of Dr. TilaMaría Perez by the Secretaría de Medio Ambiente yRecursos Naturales (Semarnat), Mexico. Adult male speci-mens were preferred because they are easier to identifymorphologically to species. We resorted to adult femalesfor H. aestus, H. luminosus, H. captiosus, and H. gilaensisbecause males were not available, but in each casethere are no sympatric closely related species likely tobe confused with them, and males have beencollected from the same locations. Both a male and afemale specimen were included for H. ophrys and H.festus in an effort to assemble a more complete referencetranscriptome.DNA of H. paratus was preserved in 95% EtOH. Allother specimens were killed by submersion in RNAlaterfor RNA preservation. To maximize tissue exposure, thecephalothorax and abdomen were opened immediatelyupon submersion. All specimens were stored at − 20 °C.Legs and the male palps or the female epigynum werepreserved separately as vouchers (stored at the BeatyBiodiversity Museum at the University of BritishColumbia, Table 1).Molecular extractions and sequencingTotal RNA was extracted from whole specimens using acombination of TRIzol extraction (Life Technologies)and RNeasy Mini Kit (Qiagen) for RNA purification andDNAse digestion. DNA was extracted from the legs andabdomen of H. paratus using a QIAamp DNA Investiga-tor Kit (Qiagen), assessed for integrity on a 21,000 Bioa-nalyzer. Libraries were constructed with BIO-ONEXTflex Library Prep Kits (Bioo Scientific, Inc.) withinsert sizes averaging 220 bp for RNA and 300 bp forDNA, and sequenced as 100 bp paired-end reads on anIllumina HiSeq 2000 (Illumina, Inc.) at the BiodiversityResearch Centre Next Generation Sequencing Facilities(University of British Columbia). For economic feasibil-ity we chose a strategy of building a de novo assembledtranscriptome for one of two species that were deeplysequenced (H. ophrys and H. festus), then assembling theother species by reference to it, allowing them to be lessdeeply sequenced. For the latter species, we aimed for atleast 20 million paired reads per species before trim-ming, and achieved 5-30 million reads after trimming(see Additional file 1: Table S1 for sequencing summary).Marshal Hedin supplied sequence reads of transcrip-tomes of H. signatus and H. ustulatus, prepared and se-quenced as 50 bp paired end sequences, and MeganPorter supplied sequence reads of H. gilaensis, preparedand sequenced as 150 bp paired end sequences. Se-quence reads are deposited in the Sequence Read Arch-ive (SRA submission SUB3319693 [38]); accessionnumbers are in Additional file 1: Table S1.Sequence read filtering and trimmingAny sequence read with an average Phred score underQ = 30 was discarded. All remaining reads were qualitychecked with FASTQC V0.10.1 [39]. Terminal nucleo-tides were trimmed using fastq-mcf from ea-utils [40] ifthey had a score below Q = 30 or if they were sequen-cing adaptors. Reads that were 95% or more homopoly-mer were discarded and any suspected contaminants,detected from the GC content curve of FASTQC, weretrimmed using PRINSEQ-lite [41]. Any read shorterthan 33 bp after trimming was discarded.Reference transcriptomeTranscriptomes were assembled de novo for H. festus andH. ophrys in Trinity RNA-seq_v20140717 [42] with thecommand: “Trinity.pl –seqType fq –left leftreads.fq –rightrightreads.fq –JM 110G –CPU 12 –inchworm_cpu12 –bflyCPU 12 –min_contig_length 200 –-kmer_cov 2”.Both assemblies were similar in size and quality, so weselected H. ophrys as the reference transcriptome be-cause of predicted ease of obtaining this species for fu-ture studies.We filtered the approximately 100,000 transcriptsfrom H. ophrys prior to using its transcriptome as a ref-erence for subsequent assemblies. We determined tran-script abundance with RSEM v1.2.19 [43] and kept onlythe most abundant transcript variant per gene. Anyremaining redundant transcripts identified with CD-HIT-EST [44] with a similarity threshold of 95% were re-moved (176 transcripts removed). To decrease the likeli-hood of paralogous genes assembling on a referencetranscript during reference-based assemblies, we con-ducted an all-versus-all BLAST with all remaining H.ophrys contigs and removed any contig with a contigother than itself as a significant hit (blastx, evalue = 10−3; 34 transcripts removed). To set codon positions, thereference transcriptome was scanned for open readingframes using TransDecoder_r20131110 [42] and the lon-gest open reading frame of a transcript was chosen as itsprotein coding region. If multiple non-overlappingLeduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 4 of 23coding regions were found on a transcript, those tran-scripts were split between coding regions.We conducted a BLAST search of the entire H. ophrystranscriptome (evalue = 10− 3, min. HSP length = 33,max_target_seqs = 20) to the H. oregonensis mitochon-drial genome [45] and any significant hits were removedfrom the final reference nuclear transcriptome. The an-notated H. oregonensis mitochondrial genome [45], wasused instead as the reference for all mitochondrialassemblies.Reference-based assembly of transcriptomesFor every species (including reference species H.ophrys),sequencing reads were first mapped to the H. oregonen-sis mitochondrial genome and the remaining unmappedreads were mapped to the reference transcriptome usingCLC Genomics Workbench (CLC Bio), chosen becauseof its extensive facilities for making and visualizingreference-based assemblies. Reads of H. ophrys were alsoremapped to its own reference so as to obtain sequen-cing depth information and follow a trimming protocolcomparable to the other species. Assembly parametersused were: mismatch cost = 2, insertion cost = 3, deletioncost = 3, length fraction = 0.5, similarity fraction = 0.8.When consensus sequences were extracted from theread mappings, polymorphisms were retained as ambi-guity codes in the consensus sequence if the variant rep-resented 30% or more of mapped alleles, except for H.festus which used a combination of read count and sumof quality scores to resolve the base conflict. Nuclear se-quences with average sequencing depth less than 5×were discarded, and contigs were split into fragments atany region where sequencing depth was less than 5×.Following these steps, only contigs longer than 200 bpwere retained. Total raw sequencing reads, trimmed se-quencing reads, the number of reads assembled, totaltranscripts, and sequencing depth are summarized forevery species in Additional file 1: Table S1. Only mito-chondrial data were used for H. paratus, as most of itsnuclear sequences were incomplete or with inadequatesequencing depth.Alignment and filtering of lociWe converted species-based FASTA files into locus-based FASTA files and trimmed sequences of anyremaining poly-A or poly-T tails using custom pythonscripts. If a locus was fragmented into multiple se-quences for a species, and the lengths of those fragmentswhen added together met our length cut-off (200 bp),then those sequence fragments were scaffolded usingtheir positions aligned against the reference transcript,filling in with “N”s to represent the missing data.At this stage, there were 55,316 loci, most present onlyin the two species with a much higher number of reads,H. ophrys and H. festus. Of those, 7132 were present inat least 15 species, and 4188 were present in at least 25species. Each locus was aligned using MAFFT v7.058b[46] using L-INS-I and parameters –localpair –maxite-rate 100. Alignments of a few loci were manually cor-rected to resolve obvious reading frame misalignments.Nucleotides trailing from either end of an alignment(present in 30% or less of species) were trimmed and asequence was discarded if it was 30% or less than theaverage sequence length for that alignment. All align-ment, trimming and partitioning steps were completedwith the aid of Mesquite 3.02 [47], with a few steps in-volving the package Gataga (Maddison and Maddison,unpubl.).The aligned matrix for each locus was partitioned bycodon position or non-coding based on Transdecoderresults. To better understand the cause of gaps and am-biguous nucleotides in these alignments, we aligned asubset of loci with annotated sequences with knownprotein-coding regions matched using a BLAST searchto the SWISS-PROT database, to Ef1-α [5], and to theannotated mitochondrial genome [45]. All sites withgaps were caused by an insertion of a nucleotide in oneor a few sequences, implying highly unlikely frame-shiftsin conserved genes. These insertions occurred only inthe reference-based assemblies and never occurred inthe (higher quality) reference sequence. Therefore, weinferred that these insertions of gaps and nucleotideswere assembly errors rather than true insertions. As aresult, columns with gaps were excluded whenevercoding-region gaps were not present in the reference se-quence (to avoid any frameshift mutations). In noncod-ing regions there were also occasionally insertions insome species — in the vast majority of columns, in onlyone or two species, again suggesting assembly errors.However, a more relaxed exclusion criterion was usedfor non-coding regions: a column was excluded if the in-sertion was present in less than 50% of species.Some loci showed high levels of ambiguous sites,which might not represent true heterozygosity, but ra-ther multiple transcript variants or paralogs assembledon the same reference transcript. To be conservative,loci were excluded from analyses if ambiguous sites con-stituted more than 3% of the alignment. 236 loci werethereby removed, resulting in 1877 loci in our highestquality subset (present in at least 25 species, at least200 bp long, with a coding region identified by Transde-coder, and with less than 3% ambiguity).Because reference-based assemblies generate se-quences aligned to the same reference transcripts, se-quences assembled on the same reference are treated asorthologous. Filtering for paralogs in the H. ophrys refer-ence transcriptome, reducing reference transcriptomeredundancy, and removing sequences with high levels ofLeduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 5 of 23ambiguous sites are expected to have reduced the possi-bility of paralogs. To assess whether there were obviousparalogs, gene trees were constructed for all genes ofalignment lengths over 1 kb with a single search repli-cate in RAxML 7.7.9 [48]. None of the trees producedhad unusually long branches or phylogenetic structuresobviously indicative of paralogy. While this reassures usthat there were not deep paralogs, some more recentparalogs may have passed through our filters. Our ex-pectation is, however, that the effect of paralogs averagedover many loci would be to add noise rather than sys-tematic biases, as noted in the discussion about intro-gression and artifacts.Nuclear subsetsWe separated the nuclear loci into four groups. The firstgroup, the Primary Subset, included the 1877 highestquality loci as described above (2.41 Mb alignment; aver-age of 2,036,173 base pairs per species). These were usedfor our primary nuclear phylogenetic analyses. Amongthe remnant “low quality” loci, those matching the cri-teria of high quality except for being present in only 15to 24 species were treated as the second group, theMissing Species (MS) subset (1019 loci; average of548,107 bp per species). The other remnants without anidentifiable coding region but present in 25 or more spe-cies were the third group, the Noncoding Loci (NL) sub-set (236 loci; average of 92,567 bp per species). Thefourth group consisted of loci not meeting any of thesecriteria; they were discarded.Mitochondrial dataTo compose a complete mitochondrial alignment, se-quences for 16S RNA (1022 bp), 12S RNA (691 bp),ND1 (921 bp), ND2 (959 bp), ND3 (342 bp), ND4(1289 bp), ND4L (268 bp), ND5 (1638 bp), ND6(429 bp), ATP6 (666 bp), ATP8 (158 bp), Cytochrome B(1111 bp), COX1 (1542 bp), COX2 (666 bp), COX3(786 bp) were aligned, concatenated, and assigned codonpositions based on annotations from the H. oregonensismitochondrial genome [45]. This yielded theconcatenated mitochondrial matrix, with a total align-ment length of 12.33 kb.Maddison and Hedin’s [5] two-gene Habronattus phyl-ogeny includes many species absent in our transcrip-tome data. In order to get a denser perspective onmitochondrial introgression within the VCCR clade, wetook sequences from the 16SND1 region from our mito-chondrial data (1047 bp) for our 35 transcriptome spe-cies and DNA specimen H. paratus, and added to themdata for the same region obtained by Maddison and He-din [5], Masta and Boore [45], and Masta and Maddison[11]. Combined these yielded a matrix with 196 se-quences of 16SND1 from across the genus, including 37VCCR species (in comparison to 13 species in our tran-scriptome data).Phylogenetic analysesHabronattus species phylogeny from nuclear genes wasinferred by maximum likelihood from concatenatedalignments, as well as by coalescent-based methods (AS-TRAL, SVDQuartets).Maximum likelihood on concatenated matricesAll maximum likelihood phylogenetic searches on theconcatenated matrices were run in RAxML 7.7.9 [48]with 20 search replicates for the ML tree, and one searchreplicate for each of the 1000 bootstrap replicates. Weused PartitionFinder (v.1.1.1 [49]) to assess the substitu-tion models, using a greedy algorithm search and AICmodel selection, and considering both locus and codonpositions as possible partitioning criteria.In addition to the primary concatenated matrix of1877 genes, 12 other matrices were analyzed represent-ing different subsets of nuclear loci, to explore theconsistency of support. Eight of these were equal subsetsof our primary matrix (nuclear subsets 1-8; 302,200 bpeach). Locus order was randomized prior to subset div-ision and concatenation to ensure that each subset rep-resents a random sample of loci. In addition, theremnant (low quality) loci were used to make two otherdisjoint matrices (MS and NL, as defined above).To assess the strength of the mitochondrial phylogen-etic signal with less data, we divided mitochondrialrRNA (1.72 kb) from protein-coding sites and then sepa-rated the protein-coding alignment into 4 even subsets(mtDNA subsets 1-4; 2.53 kb each).In the analysis for the expanded 16SND1 matrix weconstrained the inference to enforce any node in theconcatenated mitochondrial tree that had at least 90%bootstrap support. We did this so as to take advantageof the strong resolution available for the 36 transcrip-tome taxa, and to determine how the additional se-quences fell on that skeletal constraint tree. Althoughwe performed the analysis across the whole genus, ourresults focus on the VCCR clade.Coalescent methods for the species treeUsing the primary subset of nuclear genes (1877 loci),two methods based on multi-species coalescent modelswere used to infer the species phylogeny, ASTRAL [50]and SVDQuartets [51]. ASTRAL version 4.7.12 [52] wasapplied to the 1020 loci from the primary nuclear subsetof alignment length at least 1000 bp. For each locus, asingle gene tree was reconstructed by a simple max-imum likelihood search by RAxML (model GTRGAM-MAI, unpartitioned), and the set of 1020 gene trees wasanalyzed by ASTRAL using default settings. For theLeduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 6 of 23SVDQuartets analysis, PAUP* version 4.0a149 (2016;[53]) was applied to the primary concatenated nuclearmatrix (2.41 Mb alignment), default search settings, with1000 bootstrap replicates.Analysis of introgressionTo assess signals of introgression in the nuclear data, weexamined nuclear signals discordant with the speciestree using Bayesian Concordance Analysis (BCA), andtested for introgression using D-statistics.Bayesian concordance analysisWe used BCA to assess the extent of nuclear discord-ance and the possibility of introgression [33, 34]. From aBayesian sample of gene trees for each locus, BCA de-rives a primary concordance tree representing the dom-inant phylogenetic signal (an estimate of the speciestree). and secondary concordance factors (CF) indicatingsubstantial support from some genes but discordant withthe dominant species tree [33, 34, 54]. While BCA doesnot test explicitly for introgression, if two conflictingsecondary CFs are unequal and their 95% credibilityintervals do not overlap (and therefore are significantlydifferent), introgression may be inferred as a possiblecause of discordance because under ILS the CFs forconflicting secondary clades would be expected to beequal [54].We conducted all Bayesian gene tree searches withMrBayes 3.2.2 [55] with 4 chains (3 cold, 1 hot), 2 runs,and 25% burn-in. To determine the number of MCMCgenerations required for convergence, 20 genes were rununtil convergence (standard deviation of split frequen-cies < 0.01; [56]) as a test. The convergence time rangedfrom 570,000 to 13,087,000 and averaged 4,025,500 gen-erations. Based on this, we set the number MCMC gen-erations conservatively at 20,000,000 generations pergene. Codon positions were used as partitions.BCA analyses were conducted using BUCKy 1.4.3 [57]with 2 runs and 4 chains per analysis. Due to computa-tional limitations, we analyzed only the americanusgroup and the VCCR clade, and included only genes lon-ger than 1 kb that were present in all species being ana-lyzed. There were 679 genes included for theamericanus group analysis (with H. signatus as out-group) and 517 genes included for the VCCR clade ana-lysis (with H. ophrys as outgroup).For the adjustable prior α, the a priori level of genetree incongruence, we tried values of α = 0.1, 1, 2, 5 and10. We found no substantial difference in results usingdifferent α for the americanus group, so kept α set to 1.Analyses involving the VCCR clade had difficulty con-verging at higher α, so α was set to 0.1 for that analysis,although there was very little difference between CFs de-pending on α for this clade. The americanus groupanalysis ran for 10,000,000 generations and the VCCRclade analysis ran for 30,000,000 generations due to thelonger time required for convergence.Patterson’s D statistic and DFOILTo distinguish incomplete lineage sorting (ILS) andintrogression patterns in SNP data, we conducted Patter-son’s D statistic tests [58] and the related test DFOIL [22].These tests compare patterns of shared SNPs across setsof 4 and 5 taxa, respectively. Under ILS in a 4-taxa bin-ary tree, it is expected that the number of alleles sharedby non-sister taxa (i.e. discordant with the species tree)would be equal for each possible non-sister pairing (i.e.,patterns ABBA and BABA) [32]. Introgression betweentwo species is inferred when they have significantly moreshared alleles than alternative discordant pairs (i.e., moreABBA than BABA or vice versa). The same principlecan be applied to a 5-taxa tree with structure (((species1, species 2),(species 3, species 4)), outgroup), with someadded complexity [22]. Whether Patterson’s D statisticor DFOIL was used depended on the structure of the spe-cies tree of the taxa being tested.We tested sets of species based on indications of intro-gression in previous studies [5], in the mitochondrialphylogeny, or in the BCA. There were 3 principal hy-potheses of introgression tested: (1) among species ofthe americanus group, (2) between H. roberti and theother groups of the VCCR clade and (3) between the cly-peatus group and coecatus group. In addition, the possi-bility of very distant introgression between the DTBclade (decorus/texanus/banksi groups), VCCR clade andthe americanus group was also explored.For the americanus group tests, we included all 5americanus group species. For the other tests, we se-lected a few species as representatives of each clade. Forthe H. roberti/VCCR D statistics analysis, we used thespecies with the most sequencing data per clade for thefour clades involved (H. roberti, and one representativeof each of the viridipes, clypeatus, coecatus groups). Forthe clypeatus group - coecatus group analysis, we did 6combinations of species, including 2 coecatus group and2 clypeatus group species per DFOIL test, and using H.ophrys as the fifth species (outgroup). H. pyrrithrix (coe-catus group) was included in each combination becauseits phylogenetic position is consistent in both theconcatenated nuclear and mitochondrial phylogenies; H.clypeatus (clypeatus group) was included because ofmitochondrial introgression detected in the mitochon-drial phylogeny and in other members of its speciesgroup by Maddison and Hedin [5]. By using more or lessdivergent sister taxa in different tests for comparison, wecan better approximate where in the coecatus and cly-peatus group introgression occurred. Even still, D-statistics and DFOIL cannot rule out the possibility thatLeduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 7 of 23detected introgression was with related lineages that areextinct or were otherwise not sampled (i.e., ghost line-ages) [59].DFOIL was run in mode dfoilalt to reduce noise fromsynapomorphic sites. We used a custom R script tocount allele patterns. All sites that included ambiguousnucleotides, gaps, or missing data were excluded fromthe analysis. DFOIL estimated divergence T-values wereverified against the assumption that T12 < T34 < T1234.Because each DFOIL and D-test is a chi-square binomialtest, we adjusted significance for multiple comparisons(62 in total, including all DFOIL and D tests) with a Bon-ferroni correction to a p-value lower than 0.0008 to indi-cate 95% significance.ResultsTranscriptome assemblies and data filteringThe unfiltered H. ophrys reference transcriptome in-cluded 117,859 transcripts (total 53,927,457 bases as-sembled), with an N50 (analogous to median contiglength; [60]) of 516, and an average sequencing depth of103×. Following filtering for redundancy, selection of asingle variant per gene, removal of possible paralogs,and the separation of connected transcripts, there were92,343 transcripts left.Additional file 1: Table S1 gives a summary of tran-scriptome assemblies, excluding unused sequences withlow (< 5×) sequencing depth. After reads were re-mapped, 51,143 H.ophrys transcripts had sufficient (5×)sequencing depth (average sequencing depth 111×). Forall other species, reference-based transcriptome assem-blies mapped on average 77% of trimmed reads to eitherthe nuclear reference transcriptome (average nuclear se-quencing depth = 67×) or mitochondrial reference genome(average mitochondrial sequencing depth = 13,640×).There was an average of 10,164 transcripts assembled perspecies, although numbers ranged widely (depending onthe number of reads) from 3746 for H. roberti to 28,846for H. festus. Aligned matrices for each of the partitions(Primary, Missing Species, Noncoding Loci, and mito-chondrial) are available in Additional file 2.Substitution model selectionFor the primary concatenated nuclear matrix, Partition-Finder was unable to analyze the approximately 7500partitions (codon positions for each of 1877 genes) be-cause of computational limits. Thus, we applied it to themitochondrial genes and on a sample of 20 nucleargenes, using it to assess models rather than choose parti-tions. GTR +G + I or GTR + G was chosen as the opti-mal substitution model using AIC for all mitochondrialpartitions. For the 20 nuclear genes tested, 33 partitionshad a GTR model selected, 22 had TVM, 10 had K81uf,5 had TIM, 4 had HKY, and 3 had TrN (these numberinclude model variations like +G or +I). We were unableto set a different model for each partition inconcatenated matrices due to computational limitations.Instead, we set GTR+G+ I as the substitution model in allmaximum likelihood (ML) and Bayesian analyses (nst = 6rates = invgamma) because it was the most commonlychosen and most widely applicable mode. We used 4 parti-tions based on codon position (position 1, position 2, pos-ition 3, and noncoding) for all figured analyses. To testwhether partitioning played a major role, we also ranan unpartitioned likelihood analysis on the 1877 locusconcatenated nuclear matrix, as well as one parti-tioned by locus.Nuclear phylogenyThe nuclear phylogenetic results are summarized inFigs. 1 and 2a; trees are available in NEXUS format inAdditional file 3. Most species relationships are resolvedwith strong support, and concordantly among theconcatenated nuclear and species tree analyses. The in-dependent “low quality” remnant matrices also supportmany of the clades. As suggested by Maddison andHedin’s [5] results from a few genes, Habronattus is di-vided into the AAT clade (amicus, agilis, and tranquillusgroups) and a large clade of the remaining species,within which H. geronimoi is sister to the rest. Theamericanus group and the VCCR clade are also resolved,and within the latter the viridipes, clypeatus and coeca-tus groups are each monophyletic.The phylogeny has a few notable regions of uncer-tainty, which are also points of disagreement among theanalyses. Although the previous morphological data (dis-cussed below) provides good support for the monophylyof Habronattus, only one of the 1/8th nuclear partitionssupports it here. The concatenated nuclear ML tree, theconcatenated bootstrap consensus, three of the 1/8thpartitions, and the ASTRAL tree all place Pellenes cana-densis as sister to the AAT clade; the remnant matricesplace P. canadensis as sister to the major Habronattusclade that excludes the AAT clade, while the SVDQuar-tets tree weakly places P. canadensis as sister to H. con-junctus. The clade including the VCCR clade plus H.hallani, H. pugillis and H. luminosus (“VCCR+” clade) iswell supported, as is the larger clade that adds H. cam-bridgei, H. icenoglei, and H. oregonensis. However, the re-lationships among the latter three are unstable: theconcatenated nuclear ML tree chooses (cambridgei,((icenoglei, oregonensis), VCCR+)), the concatenatedbootstrap consensus (i, ((c, o),V)), ASTRAL (c,(i,(o,V))),and SVDQuartets ((c,i),(o,V)), though the last with lowbootstrap support. A few discordant placements showup in some analyses (e.g., the SVDQuartets analysisstrongly places H. tarsalis sister to H. americanus andH. sansoni, and weakly places H. roberti as sister to theLeduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 8 of 23viridipes group). The unpartitioned ML analysis of theconcatenated nuclear matrix yielded the same topologyas in Fig. 1a, but partitioning by locus yielded the closealternative resolution (icenoglei, ((oregonensis, cambrid-gei), VCCR+). The primary concordance trees from theBCA analysis in the americanus group and the VCCRclade match those portions of the concatenated nuclearML tree except for the placement of H. virgulatus.Mitochondrial phylogenyThe mitochondrial transcriptome tree (Fig. 2b; Add-itional file 3) is broadly concordant with the nucleartree, concurring on the VCCR clade, the next largerclade adding H. hallani, H. pugillis and H. luminosus,and the next larger clade adding H. cambridgei, H.oregonensis, and H. icenoglei. The americanus groupis monophyletic, as is the AAT clade. Bootstrap sup-port values are generally high, and some key resultsare consistent across rRNA and four protein-codingmitochondrial subsets (Fig. 2b). The strongly sup-ported results, however, include several notable differ-ences with the nuclear phylogeny (Fig. 2a), includingthe placement of the H. clypeatus specimen and therelationship of the DTB clade and the americanusgroup. These will be discussed in the context ofintrogression.Pellenes canadensisEvarcha proszynskiiH. conjunctusH. hirsutusH. signatusH. ustulatusH. geronimoiH. aestusH. tarsalisH. ophrysH. americanusH. sansoniH. altanusH. decorusH. chamelaH. zapotecanusH. cambridgeiH. icenogleiH. oregonensisH. hallaniH. luminosusH. pugillisH. robertiH. c. maddisoniH. jucundusH. aztecanusH. gilaensisH. clypeatusH. mexicanusH. pyrrithrixH. empyrusH. virgulatusH. festusH. borealisH. captiosusaASTRALSVDQuartetsConcat. 1/8thsMS remnantNL remnantConcat. bootstrap100bH. aestusH. tarsalisH. ophrysH. americanusH. sansonidH. aztecanusH. gilaensisH. clypeatusH. mexicanusH. pyrrthrixH. empyrusH. virgulatusH. festusH. borealisH. captiosusH. robertiH. c. maddisoniH. jucundusH. gilaensisH. clypeatusH. mexicanusH. pyrrithrixH. empyrusH. virgulatusH. festusH. borealisH. captiosus100100VCCR clade97100100 10010010010010010010010010010072100100100100100100981001007998100100100100AAT cladeamericanus groupviridipesgroupclypeatusgroupcoecatus groupDTB cladeH. aztecanuscFig. 1 Species phylogeny from nuclear data in Habronattus (for branch lengths, see Fig. 2a), and main conclusions of nuclear introgression. aMaximum likelihood tree from 1877 concatenated nuclear loci (2.41mb alignment; RAxML); largely concordant with the ASTRAL tree (stars).Colours mark species groups as labeled. Legend for decorations at lower left. Bootstrap percentages for ML analysis of concatenated nuclearmatrix. Spots and vertical bars show presence of clade in the ML tree for various partitions of the data (black = clade present). Spots showpresence of clade in tree for each of the 1/8th portions of concatenated nuclear loci. Vertical bars show presence of clade in the trees from theremnant Missing Species (MS) and the Noncoding Loci (NL) matrices. Stars show clades present in ASTRAL analysis of 1020 ML gene trees ofalignment length at least 1000 bp. Diamonds show clades with (black) > 94% or (grey) 75-94% bootstrap support in SVDQuartets analysis.Discordant results shown by grey arrows with circle (concatenated ML bootstrap consensus) or star (ASTRAL). For instance, the concatenatedbootstrap consensus places H. cambridgei as sister to oregonensis, and H. icenoglei in a more basal position with respect to the VCCR clade. Fromb to d: Conclusions of introgression from D-statistics and BCA for b the americanus group, c H. roberti and the VCCR clade, and d the clypeatusand coecatus groups. See text for signals of other introgression eventsLeduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 9 of 23The 16SND1 phylogeny (Fig. 3a; Additional file 3) ofthe VCCR clade, based on a combination of Sanger se-quencing data and transcriptome data, generally resolvesthe groups but not the species. Of the eight species ofthe coecatus group represented by more than one speci-men, only two appear as monophyletic on the tree (H.ammophilus, H. festus). The clypeatus group is notmonophyletic, with three specimens (H. cf. arcalorus“CHIH” [HA292]; H. velivolus [HA659]; H. clypeatus[GLR227, transcriptome]) appearing within the coecatusgroup.Bayesian concordance analysesKey findings from the BCA are summarized in Fig. 4afor the americanus group and Fig. 5a for the VCCRclade. All additional concordance factors > 0.05 are listedin Additional file 4.The BUCKy analysis of 679 americanus group loci(Fig. 4a; Additional file 4) converged with an average SDof mean sample-wide CF of 3.24 × 10− 5 (all 105 topolo-gies represented among the 15 M trees sampled). Theanalysis supports the same americanus group phylogenyas the concatenated nuclear phylogeny seen in Fig. 1a.H. ophrys is linked to H. americanus and H. sansoni viaa CF (CF = 0.305, CI = 0.272 - 0.339) that is more thantwice the CF shared for the equivalent pairing H. tarsalisand H. sansoni/H. americanus (CF = 0.136, CI = 0.112 -0.162). The BCA also found significant asymmetric sup-port linking H. aestus with H. tarsalis (CF = 0.136, CI =0.112 - 0.162) compared to the conflicting clade H. aes-tus and H. ophrys (CF = 0.39, CI = 0.023 - 0.055).The BUCKy analysis with 517 loci for the VCCR clade(Fig. 5a; Additional file 5) converged with an average SDof mean sample-wide CF of 0.003 (4,795,750 topologiesand 8177 distinct splits represented among the 15 Mtrees sampled). The primary concordance tree (Fig. 3b)from the VCCR-clade analysis is in agreement with theconcatenated nuclear phylogeny except for the positiona bFig. 2 Comparison of nuclear and mitochondrial phylogenetic results. a Maximum likelihood nuclear tree from the concatenated 2.41mbalignment, as in Fig. 1, but with branch lengths (RAxML). Named groups shown with same colours as in Fig. 1. Numbers show bootstrappercentages; grey branches with < 95% bootstrap support. b Maximum likelihood mitochondrial tree from the concatenated 12.33 kb alignment,with bootstrap percentages (RAxML). Bars show presence of clade in the ML tree for each of five subdivisions of the concatenated matrix (rRNA,followed by 4 portions of coding regions; black = clade present)Leduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 10 of 23of H. virgulatus, though this placement is accom-panied by a very low concordance value (CF = 0.1-0.14). Two substantial conflicting (secondary) CFsplace H. roberti with the clypeatus/coecatus clade(CF = 0.196, CI = 0.174 – 0.222) and with the clypea-tus group (CF = 0.188, CI = 0.17 – 0.205). Both ofthese CF credibility intervals are significantly higherthan those for alternative conflicting clades: H.roberti with the viridipes group (CF = 0.052, CI =0.039 – 0.06) and H. roberti with the coecatus group(CF = 0.041, CI = 0.027 - 0.058). There are 35 verysmall but significant (not overlapping 0) CFs aver-aging 0.01 linking H. clypeatus with particular coe-catus group species (see Additional file 5).Fig. 3 Mitochondrial and nuclear results from the VCCR clade. a VCCR portion of maximum likelihood tree for the 16SND1 mitochondrial region(1047 bp alignment; RAxML; branch lengths proportional to change). 160 Sanger sequenced specimens and all transcriptome specimens (markedby “Transcriptome”) are included. Tree constrained to the phylogenetic structure of the concatenated mitochondrial transcriptome phylogeny(nodes with > 90% bootstrap support only). Symbols show species polymorphic for coecatus-group and clypeatus-group type mitochondria: ● = H.clypeatus; ▼ = H. velivolus; ■ = H. cf. arcalorus. b BUCKy Primary Concordance tree for the VCCR clade, based on 517 genes longer than 1 kb.Node values are the CF credibility intervals for each clade. Other CFs and their credibility intervals are indicated in Fig. 5a and inSupplemental ResultsLeduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 11 of 23DFOIL and D-statistics tests of introgressionAsymmetries of gene-lineage sharing in the americanusgroup and the VCCR clade are seen also in the DFOILand D statistics (ABBA-BABA) results. Allele counts forthese tests are listed in Additional file 6: Table S2 andAdditional file 7: Table S4. Results (D values and pvalues) are reported in Additional file 7: Table S4 andAdditional file 8: Table S3.DFOIL tests support introgression between H. ophrysand the common ancestor of H. sansoni and H. ameri-canus (Fig. 4b, DFO = 0.232 p < 10− 11, DIL = 0.230 p = 10− 11, DFI = − 0.036 p = 0.549, DOL = − 0.043 p = 0.05).ABBA-BABA tests detect an introgression signal be-tween H. aestus and H. tarsalis, whether the fourth spe-cies is H. ophrys (Fig. 4d, D= 0.160, p < 10− 6), H.americanus (D = 0.192, p < 10− 10), or H. sansoni (D =0.165, p < 10− 8).Introgression is detected using ABBA-BABA tests be-tween H. roberti and H. gilaensis (clypeatus group) whenthe third species used for comparison is H. jucundus(viridipes group) (Fig. 5b, D= − 0.405, p < 10− 12; inter-pretation in Fig. 5c) or H. festus (coecatus group) (Fig. 5dD = 0.377, p < 10− 12; interpretation in Fig. 5e). Introgres-sion is also detected between H. festus (coecatus group)abd e fcgFig. 4 Signals of nuclear introgression in H. americanus species group (see Fig. 1b for summary interpretation). a BCA analysis using BUCKy,showing concordance factors, using Bayesian sample of trees from 679 loci. Higher CF for the discordant grouping ophrys + americanus + sansonithan tarsalis + americanus + sansoni suggests introgression. b Total biallelic pattern counts for all DFOIL tests for introgression between americanusgroup species; * indicates significant difference at 0.05 level after Bonferroni correction. c Interpretation of DFOIL: non-neutral DFO and DIL values inthe D-statistic signature indicate introgression between the ancestor of H. sansoni/H. americanus and H. ophrys. d, f ABBA vs. BABA allele patternscounts for D statistic tests for introgression in americanus group species; * indicates significant difference at 0.05 level after Bonferroni correction.e, g Interpretation: Introgression is detected between H. aestus and H. tarsalis (e); none between H. aestus and H. americanus/H. sansoni (g).Species not participating in the particular test are greyedLeduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 12 of 23and H. roberti when H. gilaensis (clypeatus) group is excluded(Fig. 5f, D=− 0.177, p < 10−12; interpretation in Fig. 5g).In DFOIL analyses of the clypeatus and coecatus groups(Fig. 6, Additional file 8: Table S3), introgression is de-tected between H. clypeatus and H. pyrrithrix when spe-cies 1 is the closely related H. gilaensis, and species 4 iseither H. mexicanus (Fig. 6g, h) or H. borealis (Fig. 6i, j).The neutral rather than negative DFO and DFI values ofthese DFOIL signatures indicate introgression with uncer-tain or reciprocal direction (signatures are DFO = 0,DIL = −, DFI = 0, DOL = −).Introgression is also detected between the commonancestor of H. aztecanus and H. clypeatus and H. pyr-rithrix when the fourth species included is the distantcoecatus group member H. mexicanus (Fig. 6e, f ) andthe less distant H. borealis (Fig. 6k. l) because the sig-nature shifts from DFO = +, DIL = -DFI = 0, DOL = 0to DFO = +, DIL = −, DFI = 0, DOL = 0. With a neutralrather than positive DFO value, the signature is nolonger indicative of any single introgression event.However, it does hint at introgression between H.clypeatus and H. pyrrithrix. Support for introgressionbetween H. pyrrithrix and the clypeatus group disap-pears when H. empyrus is the fourth species (Figs. 6a,c), ruling out H. pyrrithrix-specific introgression.DFOIL and D-statistics results are presented in Fig. 7 andAdditional file 6: Table S2, Additional file 7: Table S4,Additional file 8: Table S3 for the study of deeperintrogression among the americanus group, DTBclade, and others. Several signals of introgression weredetected among parts of Habronattus phylogeny thatare now highly divergent.DiscussionUsing transcriptome sequence data from 1877 nuclearloci, we were able to reconstruct a phylogeny for 34species of Habronattus with high confidence at mostnodes. Far better resolved than previous results basedon just two loci [5], this robust phylogenetic frame-work will enable more refined interpretations of char-acter evolution, sexual selection, and hybridization.Differences in the mitochondrial and nuclearABBA-BABA testsH. ophrysH. robertiH. jucundus(viridipes group)H. gilaensis(clypeatus group)H. festus(coecatus group)cABBA BABAH. jucundus + H. robertiH. gilaensis+ H. robertiSNP Count010002000b*H. ophrysH. robertiH. jucundus(viridipes group)H. gilaensis(clypeatus group)H. festus(coecatus group)eSNP Count010002000ABBA BABAH. gilaensis+ H. robertiH. festus+ H. robertid*H. ophrysH. robertiH. jucundus(viridipes group)H. gilaensis(clypeatus group)H. festus(coecatus group)gSNP Count010002000ABBA BABAH. jucundus + H. robertiH. festus+ H. robertif*Bayesian Concordance Analysisa0.00.10.20.30.40.5Concordance Factor (CF)Discordant patternsConcordantH. ophrysH. robertiviridipes groupclypeatus groupcoecatus group0.267 0.23 0.052 0.196 0.188 0.041 CFFig. 5 Signals of introgression between H. roberti and other VCCR clade members (see Fig. 1c for summary interpretation). a BCA analysis usingBUCKy using Bayesian sample of trees from 517 loci, showing concordance factors. Higher CF for discordant grouping roberti + clypeatus groupthan alternatives roberti + viridipes group or roberti + coecatus group suggests introgression. b, d, f ABBA vs. BABA allele patterns counts forPartitioned D statistic tests for introgression between H. roberti and the VCCR clade. * indicates significant difference at 0.05 level after Bonferronicorrection. c, e, g Interpretation: Introgression is detected between H. roberti and either the clypeatus group or the coecatus group, but the latteronly when the clypeatus group is absent. This ghost lineage effect suggests direction introgression from the clypeatus group into H. roberti.Species not participating in the particular test are greyedLeduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 13 of 23phylogenies indicate several possible instances of mito-chondrial introgression, while nuclear introgression inseveral regions of the phylogeny is suggested by BayesianConcordance Analysis and D-statistics and DFOIL.Habronattus phylogenyThe consistency of clade support from different partitionsand methods suggests that the genetic history ofHabronattus is predominantly divergent, despite clearsigns of introgression. Except for the placement of Pellenescanadensis, the relationships of H. cambridgei, H. oregonen-sis and H. icenoglei, and the placement of H. virgulatus, thephylogenetic tree shown in Fig. 1a is solidly supported asthe dominant genetic history of Habronattus species.We are doubtful of the placement of P. canadensiswithin Habronattus because of morphological support(DFO) (DIL) (DFI) (DOL)+ H. clypeatus+ H. aztecanus+ H. clypeatus+ H. aztecanus+ H. empyrus+ H. pyrrithrix+ H. empyrus+ H. pyrrithrix2006001000SNP CountH. empyrus H. pyrrithrix H. clypeatus H. aztecanus2006001000SNP Count+ H. empyrus+ H. gilaensis+ H. clypeatus+ H. gilaensis+ H. clypeatus+ H. empyrus+ H. pyrrithrix+ H. pyrrithrixH. empyrus H. pyrrithrixH. clypeatusH. gilaensisH. ophrysH. aztecanusH. gilaensisH. clypeatusH. mexicanusH. pyrrithrixH. empyrusH. borealisH. ophrysH. aztecanusH. gilaensisH. clypeatusH. mexicanusH. pyrrithrixH. empyrusH. borealisH. ophrysH. aztecanusH. gilaensisH. clypeatusH. mexicanusH. pyrrithrixH. empyrusH. borealisH. ophrysH. aztecanusH. gilaensisH. clypeatusH. mexicanusH. pyrrithrixH. empyrusH. borealisH. ophrysH. aztecanusH. gilaensisH. clypeatusH. mexicanusH. pyrrithrixH. empyrusH. borealisH. ophrysH. aztecanusH. gilaensisH. clypeatusH. mexicanusH. pyrrithrixH. empyrusH. borealis+ H. mexicanus+ H. pyrrithrix+ H. mexicanus+ H. pyrrithrix+ H. aztecanus+ H. clypeatus+ H. aztecanus+ H. clypeatus2006001000SNP CountH. aztecanus H. clypeatus H. mexicanus H. pyrrithrix* *+ H. mexicanus+ H. pyrrithrix+ H. mexicanus+ H. pyrrithrix+ H. gilaensis+ H. clypeatus+ H. gilaensis+ H. clypeatus2006001000SNP CountH. gilaensis H. clypeatus H. mexicanus H. pyrrithrix* *+ H. borealis+ H. pyrrithrix+ H. borealis+ H. pyrrithrix+ H. gilaensis+ H. clypeatus+ H. gilaensis+ H. clypeatus2006001000SNP Count*H. gilaensis H. clypeatus H. borealis H. pyrrithrix* *+ H. borealis+ H. pyrrithrix+ H. borealis+ H. pyrrithrix+ H. aztecanus+ H. clypeatus+ H. aztecanus+ H. clypeatus2006001000SNP CountH. aztecanus H. clypeatus H. borealis H. pyrrithrix*a bdfhjlcegik *Fig. 6 Signals of introgression between the coecatus and clypeatus groups from DFOIL tests (see Fig. 1d for summary interpretation). On the left (a, c, e, g,i, k): total biallelic pattern counts for all DFOIL tests; * indicates significant difference at the 0.05 level after Bonferroni correction. On the right (b, d, f, h, j, l)are phylogenetic interpretations of the tests: Introgression is detected between H. pyrrithrix and H. clypeatus (h, j) or H. clypeatus/H. aztecanus (f) in somebut not all testsLeduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 14 of 23for the monophyly of Habronattus. Against its little-ornamented relatives Pellenes and Havaika, Habronattusappears well delimited by shorter first legs and a longand thin terminal apophysis (part of the male genitalia)that has a distinctive elbow on it, though the elbow issecondarily lost in the coecatus group [3, 5]. Despite theapparently clear synapomorphy of the elbowed terminalapophysis, the genus does not hold together as mono-phyletic in many of our analyses, with Pellenes canaden-sis falling inside, near H. conjunctus or the AAT clade.We are therefore faced with three possibilities: that theelbowed terminal apophysis arose twice, that the non-elbowed terminal apophysis of Pellenes canadensis rep-resents a reversal to the ancestral state, or that Pellenescanadensis and other members of its subgenus Pellenat-tus form a close sister group to Habronattus, withextremely short branches separating their early diver-gence, leading to difficulties in resolving the deeper rela-tionships, especially given that the outgroup Evarcha isphylogenetically distant, in a separate subtribe (the Plex-ippina [37]). The choice may be resolved by having bet-ter sampling of outgroups among the Harmochirina.Because Pellenes species are also distributed in Asia [61]and Europe [62], a broader global sample of Pellenesspecimens and other closely related groups (e.g., Harmo-chirus) should be included to better tease apart relation-ships at the base of the Habronattus tree and determineif the genus is monophyletic. We expect, however, thataddition of outgroups would not change the well-supported relationships within the major clade from H.paratus and geronimoi to the VCCR clade. The mtDNAtree (Fig. 3a), which included the additional outgroups(DFO) (DIL) (DFI) (DOL)0100020003000SNP Count + H. ophrys+ H. aestus+ H. ophrys+ H. aestus+ H. oregonensis+ H. decorus+ H. oregonensis+ H. decorusH. ophrys H. aestus H. oregonensis H. decorus* *0100020003000SNP Count + H. ophrys+ H. aestus+ H. ophrys+ H. aestus+ H. festus+ H. decorus+ H. festus+ H. decorusH. ophrys H. aestus H. festus H. decorus* * *0100020003000SNP Count + H. ophrys+ H. aestus+ H. ophrys+ H. aestus+ H. jucundus+ H. decorus+ H. jucundus+ H. decorusH. ophrys H. aestus H. jucundus H. decorus** *0100020003000SNP Count+ H. ophrys+ H. aestus+ H. ophrys+ H. aestus+ H. cambridgei+ H. decorus+ H. cambridgei+ H. decorusH. ophrys H. aestus H. cambridgei H. decorus*0100020003000SNP CountH. ophrys H. aestus H. zapotecanus+ H. ophrys+ H. aestus+ H. ophrys+ H. aestus+ H. zapotecanus+ H. decorus+ H. zapotecanus+ H. decorusH. decorus** *H. signatusH. aestusH. ophrysH. decorusH. zapotecanusH. cambridgeiH. oregonensisH. jucundusH. festusH. signatusH. aestusH. ophrysH. decorusH. zapotecanusH. cambridgeiH. oregonensisH. jucundusH. festusH. signatusH. aestusH. ophrysH. decorusH. zapotecanusH. cambridgeiH. oregonensisH. jucundusH. festusH. signatusH. aestusH. ophrysH. decorusH. zapotecanusH. cambridgeiH. oregonensisH. jucundusH. festusH. signatusH. aestusH. ophrysH. decorusH. zapotecanusH. cambridgeiH. oregonensisH. jucundusH. festusa bdfhjcegiFig. 7 Signals of introgression deeper in the phylogeny; DFOIL tests among the americanus group (H. ophrys and H. aestus), H. decorus, H.zapotecanus, H. cambridgei, H. oregonensis, H. jucundus, H. festus. On the left (a, c, e, g, i) are the total biallelic pattern counts for all DFOIL tests; *indicates significant difference at the 0.05 level after Bonferroni correction. On the right (b, d, f, h, j) are phylogenetic interpretations of the tests:no introgression was detected involving H. cambridgei (d), but introgression was detected between H. ophrys and H. decorus (a), and between theamericanus group and H. oregonensis, H. jucundus, and H. festus or their ancestors (f, h, j)Leduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 15 of 23Havaika, Harmochirus, Bianor, and others, is largelyconsistent with Fig. 1.The reconstructed phylogeny is generally consistentwith morphology [3] and the previous results from twogenes [5]. The major clade whose males have fringedfirst legs and modified third legs (the VCCR clade) is in-tact as monophyletic, as are the three contained groupsfirst recognized by morphology (viridipes, clypeatus, andcoecatus groups). In this regard, the phylogeny is moreconcordant with morphology than that of Maddison andHedin [5], whose analyses gave unexpected placementsfor H. jucundus (viridipes group) and split the VCCRclade. The americanus group, with distinctive genitaliaand relatively long first legs, is monophyletic in the tran-scriptome phylogeny. The agilis, amicus and tranquillusgroups hold together as the AAT clade, distinguished bycompact bodies, a relatively far-rotated bulb of the malepalp, and a tendency for dwelling above the ground inshrubs. As expected from the results of [5] the AATclade is sister to the remainder of the genus, withinwhich H. paratus (in the mitochondrial tree) is mostbasal with respect to the bulk of species, and H. geroni-moi next.The nuclear phylogeny provides new resolution ofmid-level relationships in Habronattus. The strong sup-port for the clade of H. decorus, H. altanus, H. zapoteca-nus and H. chamela indicates that their species groups(decorus, texanus, and banksi groups) form a clade, herecalled the DTB clade. The previously intractable H. hal-lani is strongly supported as sister to H. luminosus andH. pugillis. The placement of these three species as sisterto the VCCR clade is novel, as is their collective relation-ship with H. icenoglei, H. oregonensis, and H. cambridgei.For the first time there are well-supported clades withinthe coecatus group: H. festus, H. captiosus and H. borea-lis together (100% bootstrap support) will be referred toas the Northern clade (all specimens in this clade werecollected in Canada), while H. empyrus and H. pyrrithrixare sisters (100% bootstrap support), forming the South-ern clade (they are both found in the southern USAneighboring Mexico). H. virgulatus is sister to theNorthern coecatus clade with high support (97% boot-strap support). The relationships among the subgroupsof the VCCR clade are resolved with good bootstrapsupport: the viridipes, clypeatus, and coecatus groupsform a clade (nuclear bootstrap = 83%), with H. robertias sister to them; the clypeatus and coecatus group aresisters (100%); and each of the viridipes, clypeatus andcoecatus groups is monophyletic. H. jucundus groupswith H. calcaratus rather than with the oregonensisgroup (a poorly supported relationship found in [5]),confirming that the viridipes group is in fact monophy-letic. The internal relationships of the americanus groupalso have high bootstrap support.Most of these highly supported relationships are alsoreplicated in different nuclear subsets (see subset sup-port summaries at the nodes in Fig. 1a), suggesting thatthe phylogenetic signal for most branches is robust evenwith less data. However, there are a few nodes with highbootstrap support that are unstable across nuclear sub-sets. H. roberti, while well supported overall as sister tothe rest of the VCCR clade, also groups with the clypea-tus group in some nuclear subsets. H. virgulatus departsfrom its dominant concatenated nuclear position (boot-strap support = 97%) to group with either the Southerncoecatus clade or as a basal branch of the coecatus groupin some nuclear subsets. H. ophrys is positioned as thesister to H. sansoni and H. americanus in half of all datasubsets, despite being sister to H. tarsalis (bootstrapsupport = 100%) in the concatenated nuclear phylogeny.Even if Habronattus has a predominantly divergenthistory with a clear modal gene tree, there could still bea broad scatter of discordant gene trees through incom-plete lineage sorting or hybridization. Such discordanceis expected from the group’s youth (perhaps less than 5million years, [17]), and is seen in our data. For instance,the VCCR clade has very low dominant CFs from theBayesian Concordance Analysis (Fig. 3b), potentially in-dicating that the group is still in the early stages of diver-gence with widespread incomplete lineage sorting andpossibly also ongoing hybridization. The indications ofincomplete divergence are strongest in the coecatusgroup, where there are 196 secondary (conflicting) CFsthat are significant (greater than 0), though most arevery small and not clearly indicating ILS vs. introgres-sion. The americanus group, on the other hand, hasstronger genetic concordance (high dominant CF values)even though it is equally recently diverged ([12], fig. 4a).However, such differences in concordance could reflectdifferences between the groups in the density and loca-tions of sampling, and not necessarily a difference intheir evolutionary dynamics.Introgression in HabronattusMaddison & Hedin [5]’s evidence for introgression wasthe presence in two clypeatus group specimens of mito-chondrial 16SND1 more closely related to that of thecoecatus group, close enough to argue against ILS as thesource of discord. Our data from full mitochondrialtranscriptomes indicate another case of a clypeatus-group specimen with mitochondria falling with the coe-catus group, and hint of possible introgression elsewherein the phylogeny (in the americanus group, in the VCCRclade, and deeper in the tree).Our tests using BCA, DFOIL, and D-statistics indicatethat introgression extends to the nuclear genome as well.We detected nuclear introgression among closely relatedspecies, as also found by [21], ([4, 11]; pugillis group),Leduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 16 of 23and ([63]: H. americanus). However, our analyses alsofind clear signals of more distant nuclear introgression,among species groups: between H. roberti and the cly-peatus group, between the clypeatus and coecatusgroups, and apparently among even more distantly re-lated species. We will consider the evidence for intro-gression separately for different regions of thephylogeny.Introgression within the americanus groupThe nuclear phylogeny, by grouping H. americanus withH. sansoni, and H. tarsalis with H. ophrys, is in accordwith the informally-recognized distinction of the ameri-canus group into the americanus subgroup and the tar-salis subgroup. The former typically have black first legswithout lateral fringes and are usually found on groundcovered in rocks, sticks and litter (H. americanus, H.sansoni, H. waughi (Emerton, 1926), H. bulbipes(Chamberlin & Ivie, 1941), H. kubai (Griswold, 1979));the latter have yellow, green or brown first legs withlateral fringes, and are typically found on grassyground (H. tarsalis, H. kawini (Griswold, 1979), H.mustaciata (Chamberlin & Ivie, 1941), H. ophrys, H.gigas Griswold, 1987). The mitochondrial phylogeny(Fig. 2b), however, strongly supports a sister group re-lationship between H. ophrys and H. americanus tothe exclusion of H. sansoni. On its own, this resultcould be explained by either introgression or incom-plete lineage sorting. However, two introgressionevents in the americanus group were detected in thenuclear data.Introgression among H. ophrys, H. americanus and H.sansoni is indicated by both BCA and DFOIL (Fig. 4),linking H. ophrys to H. americanus and H. sansoni col-lectively, with no sign of a preferred link to either H.americanus or H. sansoni. However, H. ophrys and H.americanus in particular are resolved as sister taxa byseveral nuclear data subsets (Fig. 1a) and in the mito-chondrial phylogeny (Fig. 2b). Such a pattern could ariseby multiple introgression events involving H. ophrys witheither the ancestor of H. americanus and H. sansoni, orone of the two, or both separately, possibly combinedwith introgression between H. americanus and H. san-soni. H. ophrys is currently sympatric with H. ameri-canus along the Pacific Coast of British Columbia,Washington and Oregon, and so may have had more op-portunity to hybridize with it than with the allopatric H.sansoni. Introgression involving H. ophrys could explainthe distribution of one morphological trait that is dis-cordant with the nuclear phylogeny: tufts above the frontpair of eyes, male ornaments unique in Habronattus tothe americanus group, are found only in H. sansoni andclose relatives (H. kubai, some H. americanus popula-tions), and in H. ophrys.Also found in both the BCA and D statistic tests aresignals of introgression involving H. aestus and H. tarsa-lis (Fig. 4f ). The BCA found substantial, asymmetricsupport linking H. aestus with H. tarsalis only comparedto the conflicting clade H. aestus and H. ophrys (Fig. 4d).This pattern is also supported by ABBA-BABA tests,which detected an introgression signal between H. aestusand H. tarsalis (Fig. 4e).Results from the BCA provide some insights into theextent of nuclear introgression. The genome-wide con-cordance factor of a clade is an estimate of the propor-tion of the genome for which the clade is true [34]. Theprimary CF and species tree put H. ophrys and H. tarsa-lis together, but the CF for the discordant H. ophrys/H.americanus/H. sansoni clade is 0.305, indicating about30% of the genome has those conflicting relationships.Of this 30%, 9% can be explained as due to ILS, giventhat the alternative H. tarsalis/H. americanus/H. sansoniclade has a CF of 0.092. This leaves a difference of 21%that cannot be explained by ILS, suggesting that as muchas 21% of the genome could have been introgressed be-tween H. ophrys and H. americanus/H. sansoni. A simi-lar argument suggests up to 10% introgression betweenH. aestus and H. tarsalis.Introgression between H. roberti and other VCCR groupsAmbiguity in the placement of H. roberti in the phylo-genetic analyses hints to possible discordance. Theconcatenated nuclear tree and ASTRAL place H. robertisister to the remainder of the VCCR clade (Fig. 1), butsome partitions and SVDQuartets give other relation-ships, e.g. as sister to the clypeatus group, sister to theclypeatus plus coecatus groups, sister to the viridipesgroup, or sister to the coecatus group. The primary con-cordance tree from BCA places H. roberti as sister to theremainder of the VCCR clade, but with a weak concord-ance factor (CF = 0.267; Fig. 5a).Both BCA and D statistics suggest that the ambiguityin the placement of H. roberti is not due solely to shortbranches or incomplete lineage sorting, but involves astrong component of discordance from introgression.Discordant patterns show a strong asymmetry toward H.roberti plus the clypeatus group (BCA, Fig. 5a; ABBA-BABA, Figs. 5b-g) rather than the alternatives of H.roberti plus the viridipes group or the coecatus group.When the clypeatus group is excluded from tests, a sig-nal for introgression is seen between H. roberti and thecoecatus group (Figs. 5f and g), but this is likely as ghostlineage effect ([32, 63]), as this signal shifts to the cly-peatus group when included (Figs. 5d and e). This ghostlineage effect suggests a direction of the introgression:from the clypeatus group into H. roberti [59]. By Figs. 5dand e, the strong signal is between H. roberti and the cly-peatus group. If H. roberti were the donor speciesLeduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 17 of 23instead, into the clypeatus group, then we would expectto see sharing of alleles between just those two, failing topredict the ghost effect of sharing between H. robertiand the coecatus group (Figs. 5f and g). Following a simi-lar argument to that given for the americanus group, theexcess CF difference of 0.14 for the discordant clade H.roberti plus the clypeatus group compared to contrastingdiscordant patterns (CF 0.188 compared to 0.04-0.05)could suggest that 14% of the genome of H. roberti isintrogressed from the clypeatus group.The muddled nature of genetic relationships of H.roberti is echoed in its phenotype. The species was con-sidered part of the viridipes group by Maddison and He-din [5] based on its sharing the latter’s synapomorphy ofa raised ridge of setae between the male’s posterior eyes[35]. However, H. roberti has a clearly visible checkeredpattern of pigment in the male anterior median eyes,otherwise known only from the clypeatus group [35]. Italso has a dark medial ventral abdominal stripe, as in theclypeatus group. Males from some populations of H.roberti have red-purple bumps on their third leg’s patella[35], a trait found otherwise only in clypeatus group spe-cies H. cf. arcalorus “CHIH”, H. formosus (Banks, 1906),and H. velivolus Griswold, 1987 (a species sympatric withH. roberti). If H. roberti is sister to the VCCR clade as perthe nuclear phylogeny, those traits in which it resemblesthe viridipes group could be ancestral to the VCCR clade,with its clypeatus-group traits acquired by introgression.Introgression between the coecatus and clypeatus groupOur finding of a third case of mitochondrial sequencesin the clypeatus group falling in the coecatus group (Fig.3a) adds to the previous results of Maddison and Hedin[5]. Our H. clypeatus specimen has a mitochondrial gen-ome nestled well within the coecatus group, stronglysupported by the whole mitochondrial transcriptome(Fig. 2b). Maddison and Hedin [5] found specimens ofH. velivolus and H. cf. arcalorus “CHIH” with 16SND1sequences closely resembling those of the coecatusgroup. In each case, these species are polymorphic, withspecimens from the same or other locations showingtypical clypeatus-group mitochondrial genes. ILS is for-mally a possible explanation, requiring that the mito-chondrial lineage dominant in the clypeatus group wentextinct or unsampled in the coecatus group. However, ifILS were the cause of discord, the mitochondrial line-ages shared between unusual clypeatus group specimensand the coecatus group would have to extend deeperthan the common ancestor of the 35 described speciesin the clade containing the clypeatus and coecatusgroups. However, the sequence divergences are small —for example, the divergence between H. velivolus HA659(clypeatus group) and H. pyrrithrix HA010 (coecatusgroup) is about the same as that among the three H.virgulatus specimens from Arizona. If divergence be-tween HA659 and HA010 mitochondria was prior to thesplit of the two species groups, this would suggest thateither the three H. virgulatus mitochondria also extendthat deep and just by luck managed to sort themselvesinto three specimens of the same morphospecies, or therate of sequence divergence in H. virgulatus has drastic-ally sped up. Recent mitochondrial introgression is asimpler explanation.A signal for nuclear introgression between the coeca-tus and clypeatus groups, possibly many events, is foundvia DFOIL tests. Introgression is detected between H. pyr-rithrix and H. clypeatus (Figs. 6g–j), and between thecommon ancestor of H. aztecanus and H. clypeatus andH. pyrrithrix (Figs. 6e and f ). However, these signals arelost whenever H. empyrus, the sister taxon of H. pyrrith-rix, is included in the test, ruling out H. pyrrithrix-spe-cific introgression. The signal is also partially lost,producing DFOIL signatures inconsistent with a singleintrogression event, when H. mexicanus (Fig. 6e) is re-placed by H. borealis (Fig. 6 and k) as the coecatus spe-cies used for comparison. The complex pattern of DFOILsignatures that vary depending on the species includedfor comparison could be explained by multiple intro-gression events. It does not indicate a clear direction(coecatus group into clypeatus group or vice versa).However, the lack of an introgression signal in Figs. 6a–dhelps pinpoint the origin of the introgression signal in thecoecatus group to an unidentified ghost taxon from theSouthern clade that is closely related to H. pyrrithrix.Despite the clear signal of introgression from DFOIL,the BCA does not show any concordance factors ofgreater than 0.05 intermixing the clypeatus and coecatusgroups. Although there are 35 significant (not overlap-ping 0) CFs that link H. clypeatus with coecatus groupspecies, these are so small (averaging 0.01) that they can-not be attributed to introgression — low CF values arelikely to be overestimated [54], and are just as likely tobe a result of ILS or noisy gene trees [33].The mitochondrial and nuclear results could be an in-dication of frequent introgression between the clypeatusgroup (containing 10 species) and coecatus group (24species). This is unexpected, as these clades differ sub-stantially and consistently in morphology and courtshipbehaviour [5, 6, 64]. The male courtship display is highlycomplex in both groups, but differs in morphological or-naments, behaviours, and acoustic signals [6, 64]. Geni-tal morphology is also quite different in males andfemales of the two groups [3]. Species of the coecatusgroup lack an elbow on the male palp’s terminal apophy-sis, present in all clypeatus group species [3]. The cly-peatus group has a male tibial apophysis that is thinnerand hook-like; the coecatus group thicker and more tri-angular. The possibility of hybridization despite theseLeduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 18 of 23substantial differences related to mating is worth noting.Because of a predicted evolutionary lag between malesexually selected traits and female preferences [65],species-specific differences in male mating traits do notalways result in reproductive isolation [66]. However,given the considerable phylogenetic distance and differ-ences between clypeatus and coecatus groups, isolationwould have been expected.Introgression signal as artifact?Our inference of nuclear introgression is based on adataset sorted into proposed orthologs and assumed tobe free of contamination. Were some of our “loci” ablend of different paralogs, some resulting gene treesmay appear (falsely) to reveal introgression. However,these “loci” discordant with the primary species treewould not be expected to have a consistent direction ofdiscord, instead adding non-directional noise, unless wesuppose an elaborate scenario (e.g. chromosome duplica-tion in a shared ancestor, with the alternatives markedand later deleted wholesale in different descendants).There is no cytogenetic reason to suspect large-scale du-plications: except for H. banksi (Peckham & Peckham,1901) and H. zapotecanus, all Habronattus have the samenumber of chromosome arms and approximately thesame size of chromosomes (71 species studied in Maddi-son & Leduc-Robert [12], including most reported here).While paralogs might have survived our filters, it is un-likely that they generated our signals of introgression.A more serious concern would be cross-sample con-tamination [67], which could mimic introgressionclosely. Ballenghien et al. [67] noted that contaminationwas most likely between samples sharing a batch at a se-quencing facility. Our samples were run in five batchesseparated by many months or in different facilities: (1)H. festus and ophrys; (2) H. aztecanus, cambridgei, cap-tiosus, chamela, hirsutus, mexicanus, paratus, roberti,and zapotecanus; (3) H. signatus and ustulatus; (4) H.gilaensis; and (5) all other taxa. Of our four primary con-clusions of introgression (Figs. 1 b, c and d), one can beeasily defended: H. ophrys was in a different sequencingbatch from H. americanus and H. sansoni. The othertests of introgression could have been compromised bywithin-batch contamination. For example, H. robertishared a sequencing batch with H. aztecanus, a clypeatusgroup member sympatric with H. roberti, and thus apossible source of the patterns of Fig. 5 through eithercontamination or introgression. However, although con-tamination should be acknowledged as a possibility, it isunlikely to be the source of introgression signals in H.roberti. Most critically, contamination from H. aztecanuswould predict a pattern of biased discord similar in nu-clear and mitochondrial genomes, while introgressionwould predict a different pattern in the mtDNA, with asingle clear tree and no biased discord among sites. Thelatter is what we see. The nuclear ABBA-BABA test ofFigs. 5d for H. festus - gilaensis - roberti - ophrys showsa strong bias of 1051 for ABBA versus 476 for BABA(Additional file 7:Table S4), but the corresponding teston mtDNA (reordered H. festus - roberti - gilaensis -ophrys to match the inferred mitochondrial tree) showsno bias, with 89 sites ABBA and 88 BABA (and 215BBAA). There is no sign of a minority signal in themtDNA, against the predictions of contamination. Inaddition, our protocol would have scored as ambiguous,and thus excluded from D statistics tests, any nucleotideunless it reached more than 70% prevalence amongreads at a site in our protocol. Ballenghein et al.’s [67]data suggest it rare for a contaminant to achieve suchhigh prevalence: just 6 of their 446 cox1 samples exceeda 7:3 ratio of unexpected:expected reads (their fig. 2b). Afrequency of 6/446 is not high to yield the biases we ob-served (Additional file 7: Table S4). Finally, H. roberti isnot an outlier in percentage of sites with ambiguouscalls. Similar arguments can be given for the other con-clusions of introgression in Fig. 1, and thus contamin-ation is not a likely alternative explanation.Signals of broader introgressionAncient introgression deeper in Habronattus phylogenyis suggested by both the mitochondrial phylogeny and Dstatistics. Mitochondrial introgression between theamericanus group and the DTB clade could explain whythese groups form a clade in the mitochondrial tran-scriptome phylogeny (Fig. 2b) but are separate and dis-tant clades in the nuclear phylogeny (Fig. 2a). When weexplored this further in the nuclear data with D-statistics, we found introgression signals shared not onlybetween the americanus group and H. decorus, but alsowith many more distant species (Fig. 7). These results aredifficult to interpret because the phylogenetic distance ob-scures the donor and recipient lineages, and multipleintrogression events could create conflicting signals.Potential drivers of hybridizationOur results indicate a clear genomic signature of intro-gression in Habronattus among species groups withstrong phenotypic differences in courtship and genitalia,differences which might have been expected to providereproductive isolation. We found introgression not onlybetween the clypeatus and coecatus groups (as in [5]),but also between H. roberti and the clypeatus group, andpossibly other more distant species.Evidence for directional introgression involving thecoecatus group may hint at what processes have drivenhybridization despite divergent morphology and behav-iour. Courtship of the coecatus group is strikingly com-plex [6], more complex than any other known inLeduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 19 of 23Habronattus. Mitochondrial data suggest introgressionfrom the coecatus group into the clypeatus group(Figs. 2b, 3a), but nuclear data show hints of multiple in-trogressions with no clear directionality (Fig. 6). Direc-tionality of mitochondrial introgression but not nuclearintrogression could be the result of sex-biasedhybridization, driven by a difference in female discrimin-ation [68] causing coecatus females to sometimeshybridize with clypeatus group males but not the otherway around. In some Habronattus (outside of the coeca-tus and clypeatus groups), females have been found toprefer foreign males from divergent populations (H. pug-illis; [7]). Such xenophilia could be a consequence of anarms race between the sexes [27] involving sensory ex-ploitation: males would evolve novel exploitative signals;females evolve resistance to their own males, but not tonew signals evolved in foreign lineages to which the fe-males have not been exposed, but which share commonsensory biases [6, 7]. If females of a species group (e.g.the coecatus group) are particularly susceptible to novelsignals, then this could lead both to higher courtshipcomplexity in that group, and to their females’ tendencyto donate mitochondria to distant species. Thus, the ten-dency to hybridize might correlate positively with morecomplex courtship behaviour. However, other factorscould promote asymmetric introgression. Hedin and Low-der [21] suggested asymmetrical overlap in body size dis-tributions could explain the directional introgression theyfound in the amicus group of Habronattus. Demographiccharacteristics of hybridizing species [23], such as differ-ences in population size [69] and dispersal behaviours[70], could also cause asymmetric introgression.Evolutionary consequences of introgressionIntrogression in Habronattus may have been frequentenough to result in detectable nuclear admixture, but in-frequent enough so as not to be a homogenizing forcereducing genetic and morphological diversity [24]. Forinstance, Blackburn and Maddison [71] found evidencefor gene flow among parapatric subpopulations of H.americanus, and yet they showed consistent courtshipdisplay differences. Both substantial introgression andrapid diversification are found in the species-rich ameri-canus group (12 species) and the VCCR clade (42 spe-cies), both of which have remarkable diversities of maleornaments. In the VCCR group, as much as 14% of de-rived alleles are shared between distantly related clypea-tus/coecatus groups and H. roberti due to introgression,while closely related americanus group species H. ophrysand H. americanus/H. sansoni share up to 21% of de-rived alleles as a result of introgression.It is unclear whether introgression has played a cre-ative evolutionary role in Habronattus, either by pro-moting diversification or by influencing the phylogeneticdistribution of traits. Distant introgression, which we de-tected in several clades, is more likely to have adaptiveeffects on lineages because novel and potentially adap-tive genetic combinations are more likely to form as aresult of introgression when there is more time to accu-mulate genetic differences [72, 73]. Introgression couldalso create adaptive potential by increasing the standingvariation of hybridizing lineages, which could facilitatesubsequent diversification [73]. There are an increasingnumber of documented cases of adaptive introgression[74–76] and introgression-facilitated diversification inanimals [77, 78]. Given the strength of sexual selectionin Habronattus [1, 2, 11], loci implicated in sexually se-lected traits could be under strong selection if they wereexchanged between hybridizing species. We have sug-gested two courtship traits in Habronattus whose phylo-genetic distributions could be explained byintrogression: the red-purple bumps on the third legs ofmales of H. roberti and members of the clypeatus group,and the distinctive “eyebrows” of H. ophrys and H. san-soni. A denser phylogenetic sample of Habronattus spe-cies will be needed to better resolve whether suchsimilarities discordant with the species tree are best ex-plained by distant introgression.ConclusionsWe produced a highly resolved phylogeny of Habronattusand determined the contributions of hybridization and in-complete lineage sorting to genetic discordance in thegroup. We found that hybridization has been common inHabronattus phylogeny, and has resulted in considerablenuclear introgression in some instances (e.g., the ameri-canus group, H. roberti) and lesser nuclear introgressionaccompanied by strong mitochondrial introgression inothers (i.e., among the coecatus and clypeatus groups).However, we were unable to detect specific lineages anddirection of introgression in many cases, for which adenser sample of species will be important. Widespreadintrogression between both distant and closely related spe-cies indicates that only partial reproductive isolation hasevolved across much of the Habronattus phylogeny. Intro-gression has occurred between Habronattus speciesgroups with divergent male ornaments and courtship be-haviours. Further research should focus on determiningthe extent of introgression’s contribution to adaptationand diversification. In particular, widespread introgressionin the highly diverse VCCR clade could be indicative of acorrelation between introgression and diversification.Additional filesAdditional file 1: Table S1. Summary of RNA-sequencing and transcrip-tome assemblies for all species (except DNA sequencing of H. paratus).(PDF 616 kb)Leduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 20 of 23Additional file 2: NEXUS files of aligned matrices (five text files, NEXUSformat, totaling 131 mb, compressed to a single zip file of 12 mb) (ZIP15390 kb)Additional file 3: NEXUS file of the phylogenetic trees from variousanalyses (text, NEXUS format 250 kb) (TXT 244 kb)Additional file 4: Concordance file for americanus group. (TXT 74 kb)Additional file 5: Concordance file for VCCR clade (text 4.3 mb) (TXT4182 kb)Additional file 6: Table S2. Counts of alleles shared amongHabronattus species used for DFOIL tests. (PDF 405 kb)Additional file 7: Table S4.. Patterson’s D statistic results and counts ofallele patterns. (PDF 455 kb)Additional file 8: Table S3. DFOIL results. (PDF 476 kb)AbbreviationsAAT: The clade of species including the agilis, amicus, tranquillus groups.;BCA: Bayesian Concordance Analysis.; bp: Base pair (1 nucleotide).;CF: Concordance factor.; CI: Credibility interval.; DTB: The clade of speciesincluding the decorus, texanus, and banksi groups.; ILS: Incomplete lineagesorting.; kb: Kilobase (1 thousand nucleotides).; Mb: Megabase (1 millionnucleotides).; ML: Maximum likelihood.; MS: “Missing Species”, a subset of thedata.; NL: “Noncoding Loci”, a subset of the data.; VCCR: the clade of speciesincluding the viridipes, coecatus, clypeatus, and roberti groups.AcknowledgmentsWe thank Heather Proctor, Sam Evans, Abraham Meza López, and Tila MaríaPerez for help collecting specimens for this project in 2013 and 2014.Marshal Hedin and Megan Porter graciously contributed sequencing data forthree of the species. Allyson MisCampbell and Carol Ritland from theGenomic Data Centre at UBC, and Jeff Richards and his students, generouslyassisted with laboratory resources and advice. Anastasia Kuzmin facilitatedthe sequencing process and gave valuable advice. Saemundur Sveinssongave helpful bioinformatics assistance. Westgrid (CanadaComputes) suppliedaccess to their computer clusters. Marshal Hedin gave helpful comments onan earlier version of the manuscript. We thank also Sally Otto, LorenRieseberg, Darren Irwin, Edy Piascik, Junxia Zhang, Gwylim Blackburn, SamEvans, Hervé Philippe, and three anonymous reviewers for their helpfulfeedback and ideas.FundingThis work was supported by an NSERC Discovery grant to WPM. The fundershad no role in study design, data collection and analysis, decision to publish,or preparation of the manuscript.Availability of data and materialsThe sequencing reads of transcriptomes are archived in the Sequence ReadArchive [79]. Accession numbers are indicated in Additional file 1: Table S1.NEXUS files containing the aligned matrices of the 1877 loci in the primarydataset, and files detailing the results of analyses supporting the conclusionsof this article are included within its Additional files. The NEXUS files are alsoavailable in the Dryad repository [38].Authors’ contributionsWPM and GLR co-designed the study. WPM led the collecting of specimens,with contributions from GLR and collaborators (listed in Acknowledgments).GLR performed the lab work and analyses of data, including writing of scriptsfor processing data. WPM wrote custom modules in Mesquite for data pro-cessing. GLR wrote the first draft of the manuscript, which was subsequentlyco-written by GLR and WPM. Both authors have read and approved themanuscript.Ethics approval and consent to participateNot applicable.Consent for publicationNot applicable.Competing interestsThe authors declare that they have no competing interests.Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims inpublished maps and institutional affiliations.Received: 5 April 2017 Accepted: 15 February 2018References1. Peckham GW, Peckham EG. Observations on sexual selection in spiders ofthe family Attidae. Occas Pap Nat Hist Soc Wisconsin. 1889;1:117–51.2. Peckham GW, Peckham EG. Additional observations on sexual selection inspiders of the family Attidae: with some remarks on Mr. Wallace’s theory ofsexual ornamentation. Occas. Pap. Nat. Hist. Soc. Wisconsin. 1.3. NationalHistory Society of Wisconsin. 1890;1:3–60.3. Griswold CE. A revision of the jumping spider genus Habronattus FOP-Cambridge (Araneae; Salticidae), with phenetic and cladistic analyses. UnivCalif Publ Entomol Univ of California Press. 1987;107:1–344.4. Maddison WP, McMahon MM. Divergence and reticulation among montanepopulations of a jumping spider (Habronattus pugillis Griswold). Syst Biol.2000;49:400–21.5. Maddison WP, Hedin MC. Phylogeny of Habronattus jumping spiders(Araneae: Salticidae), with consideration of genitalic and courtshipevolution. Syst Entomol. 2003;28:1–22.6. Elias DO, Maddison WP, Peckmezian C, Girard MB, Mason AC. Orchestratingthe score: complex multimodal courtship in the Habronattus coecatus groupof Habronattus jumping spiders (Araneae: Salticidae). Biol J Linn Soc. 2012;105:522–47.7. Hebets EA, Maddison WP. Xenophilic mating preferences amongpopulations of the jumping spider Habronattus pugillis Griswold. Behav Ecol.2005;16:981–8.8. Scheidemantel DD. Behavioral and natural history studies of the jumpingspider Habronattus oregonensis and inquiry based secondary laboratorylesson development stemming from university research, MSc thesis. Tucson:The University of Arizona; 1997.9. Elias DO. Female preference for complex/novel signals in a spider. BehavEcol. 2006;17:765–71.10. Blackburn GS, Maddison WP. Stark sexual display divergence amongjumping spider populations in the face of gene flow. Mol Ecol. 2014;23:5208–23.11. Masta SE, Maddison WP. Sexual selection driving diversification in jumpingspiders. Proc Natl Acad Sci U S A. 2002;99:4442–7.12. Maddison WP, Leduc-Robert G. Multiple origins of sex chromosome fusionscorrelated with chiasma localization in Habronattus jumping spiders(Araneae: Salticidae). Evolution. 2013;67:2258–72.13. Maddison WP. XXXY sex chromosomes in males of the jumping spidergenus Pellenes (Araneae: Salticidae). Chromosoma. 1982;85:23–37.14. Zurek DB, Cronin TW, Taylor LA, Byrne K, Sullivan MLG, Morehouse NI.Spectral filtering enables trichromatic vision in colorful jumping spiders.Curr Biol. 2015;25:R403–4.15. Taylor LA. Color and communication in Habronattus jumping spiders: testsof sexual and ecological selection. PhD dissertation. Tempe: Arizona StateUniversity; 2012.16. Taylor LA, Amin Z, Maier EB, Byrne KJ, Morehouse NI. Flexible color learningin an invertebrate predator: Habronattus jumping spiders can learn to preferor avoid red during foraging. Behav Ecol. 2016;27:520–9.17. Bodner MR, Maddison WP. The biogeography and age of salticid spiderradiations (Araneae: Salticidae). Mol Phylogenet Evol. 2012;65:213–40.18. Degnan J, Rosenberg N. Gene tree discordance, phylogenetic inference andthe multispecies coalescent. Trends Ecol Evol. 2009;24:332–40.19. Grant PR, Grant BR, Petren K. Hybridization in the recent past. Am Nat. 2005;166:56–67.20. Donnelly MJ, Pinto J, Girod R, Besansky NJ, Lehmann T. Revisiting the role ofintrogression vs shared ancestral polymorphisms as key processes shapinggenetic diversity in the recently separated sibling species of the Anophelesgambiae Complex. Heredity (Edinb). 2004;92:61–8.21. Hedin MC, Lowder MC. Phylogeography of the Habronattus amicus speciescomplex (Araneae: Salticidae) of western North America, with evidence forlocalized asymmetrical mitochondrial. Zootaxa. 2009;2307:39–60.Leduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 21 of 2322. Pease J, Hahn M. Detection and polarization of introgression in a five-taxonphylogeny. Syst Biol. 2015;64:651–62.23. Toews DPL, Brelsford A. The biogeography of mitochondrial and nucleardiscordance in animals. Mol Ecol. 2012;21:3907–30.24. Slatkin M. Gene flow in natural populations. Annu Rev Ecol Syst. 1985;16:393–430.25. Twyford AD, Ennos RA. Next-generation hybridization and introgression.Heredity (Edinb). 2012;108:179–89.26. Rieseberg L, Raymond O, Rosenthal D, Lai Z. Major ecological transitions inwild sunflowers facilitated by hybridization. Science. 2003;301:1211–6.27. Holland B, Rice W. Perspective: chase-away sexual selection: antagonisticseduction versus resistance. Evolution. 1998;52:1–7.28. Ryan MJ, Keddy-Hector A. Directional patterns of female mate choice andthe role of sensory biases. Am Nat. 1992;139:S4–35.29. Schwenk K, Brede N, Streit B. Introduction. Extent, processes andevolutionary impact of interspecific hybridization in animals. Philos. Trans. R.Soc. Lond. B. Biol. Sci. 2008;363:2805–11.30. Yu Y, Than C, Degnan JH, Nakhleh L. Coalescent histories on phylogeneticnetworks and detection of hybridization despite incomplete lineage sorting.Syst Biol. 2011;60:138–49.31. Kronforst M, Young L, Blume L, Gilbert L. Multilocus analyses of admixtureand introgression among hybridizing Heliconius butterflies. Evolution. 2006;60:1254–68.32. Durand EY, Patterson N, Reich D, Slatkin M. Testing for ancient admixturebetween closely related populations. Mol Biol Evol. 2011;28:2239–52.33. Baum DA. Concordance trees, concordance factors, and the exploration ofreticulate genealogy. Taxon. 2007;56:417–26.34. Ané C, Larget B, Baum D, Smith SD, Rokas A. Bayesian estimation ofconcordance among gene trees. Mol Biol Evol. 2007;24:412–26.35. Maddison WP. New species of Habronattus and Pellenes jumping spiders(Araneae: Salticidae: Harmochirina). ZooKeys. 2017;646:45–72.36. Maddison WP, Maddison DR. Two new jumping spider species of theHabronattus clypeatus group (Araneae, Salticidae, Harmochirina). ZooKeys.2016;625:1–10.37. Maddison WP. A phylogenetic classification of jumping spiders (Araneae:Salticidae). J Arachnol. 2015;43:231–92.38. Leduc-Robert G, Maddison WP. Aligned data matrices and resulting trees forHabronattus transcriptome phylogeny. 2017. https://doi.org/10.5061/dryad.5kg33.39. Andrews S. FastQC: a quality control tool for high throughput sequencedata. 2010. http://www.bioinformatics.babraham.ac.uk/projects/fastqc.Accessed 8 July 2015.40. Aronesty E. ea-utils: command-line tools for processing biologicalsequencing data. 2011. https://expressionanalysis.github.io/ea-utils/.Accessed 20 Feb 2018.41. Schmieder R, Edwards R. Quality control and preprocessing ofmetagenomic datasets. Bioinformatics. 2011;27:863–4.42. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J,et al. De novo transcript sequence reconstruction from RNA-seq usingthe trinity platform for reference generation and analysis. Nat Protoc.2013;8:1494–512.43. Li B, Dewey C. RSEM: Accurate transcript quantification from RNA-Seq datawith or without a reference genome. BMC Bioinformatics. 2011;12:323.44. Fu L, Niu B, Zhu Z, Wu S, Li W. CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics. 2012;28:3150–2.45. Masta SE, Boore JL. The complete mitochondrial genome sequence of thespider Habronattus oregonensis reveals rearranged and extremely truncatedtRNAs. Mol Biol Evol. 2004;21:893–902.46. Katoh K, Standley D. MAFFT multiple sequence alignment software version 7:improvements in performance and usability. Mol Biol Evol. 2013;30:3059–66.47. Maddison WP, Maddison DR. Mesquite: a modular system for evolutionaryanalysis. Version 3.02. 2011. Accessed 7 Jan 2015.48. Stamatakis A. RAxML-VI-HPC: maximum likelihood-based phylogeneticanalyses with thousands of taxa and mixed models. Bioinformatics.2006;22:2688–90.49. Lanfear R, Calcott B, Ho S, Guindon S. PartitionFinder: combined selection ofpartitioning schemes and substitution models for phylogenetic analyses.Mol Biol Evol. 2012;29:1695–701.50. Mirarab S, Reaz R, Bayzid MS, Zimmermann T, Swenson MS, Warnow T.ASTRAL: genome-scale coalescent-based species tree estimation.Bioinformatics. 2014;30:i541–8.51. Chifman J, Kubatko L. Identifiability of the unrooted species tree topologyunder the coalescent model with time-reversible substitution processes,site-specific rate variation, and invariable sites. J Theor Biol. 2015;374:35–47.52. Mirarab S, Warnow T. ASTRAL-II: coalescent-based species tree estimationwith many hundreds of taxa and thousands of genes. Bioinformatics. 2015;31:i44–52.53. Swofford DL. PAUP*. Phylogenetic analysis using parsimony (* and othermethods). Version 4. Sunderland, MA: Sinauer Associates; 2002.54. Ané C. Reconstructing concordance trees and testing the coalescent model fromgenome-wide data sets. In: Knowles LL, Kubatko LS, editors. Estimating species trees:practical and theoretical aspects. Hoboken, NJ: Wiley-Blackwell; 2010. p. 35–52.55. Ronquist F, Huelsenbeck J. MrBayes 3: Bayesian phylogenetic inferenceunder mixed models. Bioinformatics. 2003;19:1572–4.56. Ronquist F, Huelsenbeck J, Teslenko M. MrBayes 3.2 Manual. 2011. http://mrbayes.sourceforge.net/mb3.2_manual.pdf. Accessed 30 Sept 2015.57. Larget B, Kotha S, Dewey C, Ané C. BUCKy: gene tree/species treereconciliation with Bayesian concordance analysis. Bioinformatics. 2010;26:2910–1.58. Durand E, Patterson N, Reich D, Slatkin M. Testing for ancient admixturebetween closely related populations. Mol Biol Evol. 2011;28:2239–52.59. Eaton D, Hipp A, González-Rodríguez A, Cavender-Bares J. Historicalintrogression among the American live oaks and the comparative nature oftests for introgression. Evolution. 2015;69:2587–601.60. Miller JR, Koren S, Sutton G. Assembly algorithms for next-generationsequencing data. Genomics. 2010;95:315–27.61. Logunov DV, Marusik YM, Rakov SY. A review of the genus Pellenes in thefauna of Central Asia and the Caucasus (Araneae, Salticidae). J Nat Hist.1999;33:89–148.62. Metzner H. Die Springspinnen (Araneae, Salticidae) Griechenlands. Andrias.1999;14:1–279.63. Eaton DAR, Ree RH. Inferring phylogeny and introgression using RADseqdata: an example from flowering plants (Pedicularis: Orobanchaceae). SystBiol. 2013;62:689–706.64. Elias DO, Mason AC, Maddison WP, Hoy RR. Seismic signals in a courtingmale jumping spider (Araneae: Salticidae). J Exp Biol. 2003;206:4029–39.65. Schluter D, Price T. Honesty, perception and population divergence insexually selected traits. Proc R Soc London B Biol Sci. 1993;253:117–22.66. Panhuis TM, Butlin R, Zuk M, Tregenza T. Sexual selection and speciation.Trends Ecol Evol. 2001;16:364–71.67. Ballenghien M, Faivre N, Galtier N. Patterns of cross-contamination in amultispecies population genomic project: detection, quantification, impact,and solutions. BMC Biol. 2017;15:25.68. Avise J. Cytonuclear genetic signatures of hybridization phenomena:rationale, utility, and empirical examples from fishes and other aquaticanimals. Rev Fish Biol Fish. 2000;10:253–63.69. Lepais O, Petit RJ, Guichoux E, Lavabre JE, Alberto F, Kremer A, Gerber S.Species relative abundance and direction of introgression in oaks. Mol Ecol.2009;18:2228–42.70. Funk DJ, Omland KE. Species-level paraphyly and polyphyly: frequency,causes, and consequences, with insights from animal mitochondrial DNA.Annu Rev Ecol Evol Syst. 2003;34:397–423.71. Maddison WP, Blackburn GS. Insights to the mating strategies ofHabronattus americanus jumping spiders from natural behaviour and stagedinteractions in the wild. Behaviour. 2015;152:1169–86.72. Stelkens R, Seehausen O. Genetic distance between species predicts noveltrait expression in their hybrids. Evolution. 2009;63:884–97.73. Seehausen O. Conditions when hybridization might predispose populationsfor adaptive radiation. J Evol Biol. 2013;26:279–81.74. Huerta-Sánchez E, Jin X, Asan, Bianba Z, Peter BM, Vinckenbosch N, et al.Altitude adaptation in Tibetans caused by introgression of Denisovan-likeDNA. Nature. 2014;512:194–7.75. Norris LC, Main BJ, Lee Y, Collier TC, Fofana A, Cornel AJ, Lanzaro GC.Adaptive introgression in an African malaria mosquito coincident with theincreased usage of insecticide-treated bed nets. Proc Natl Acad Sci U S A.2015;112:815–20.76. Heliconius Genome Consortium. Butterfly genome reveals promiscuousexchange of mimicry adaptations among species. Nature. 2012;487:94–8.77. Streicher JW, Devitt TJ, Goldberg CS, Malone JH, Blackmon H, Fujita MK.Diversification and asymmetrical gene flow across time and space:lineage sorting and hybridization in polytypic barking frogs. Mol Ecol.2014;23:3273–91.Leduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 22 of 2378. Seehausen O. Hybridization and adaptive radiation. Trends Ecol Evol. 2004;19:198–207.79. Leduc-Robert G, Maddison WP. Sequence reads for Habronattustranscriptomes. 2017. https://www.ncbi.nlm.nih.gov/bioproject/PRJNA421554.• We accept pre-submission inquiries • Our selector tool helps you to find the most relevant journal• We provide round the clock customer support • Convenient online submission• Thorough peer review• Inclusion in PubMed and all major indexing services • Maximum visibility for your researchSubmit your manuscript atwww.biomedcentral.com/submitSubmit your next manuscript to BioMed Central and we will help you at every step:Leduc-Robert and Maddison BMC Evolutionary Biology (2018) 18:24 Page 23 of 23