Addressing incomplete lineage sorting and paralogy in the inference of uncertain salmonid phylogenetic relationships

Sequence data acquisition and processing

We obtained sequence data from representatives of all major salmonid genera (n = 15), as well as from representatives of the two extant esociform families Esocidae and Umbridae, which are the closest diploid relatives of salmonids (Ishiguro, Miya & Nishida, 2003; López, Chen & Ortí, 2004). We generated genetic sequence data for phylogenetic analysis by sequence capture of 500 UCE loci following the protocol outlined for the 500 UCE actinopterygiian probe set (Faircloth et al., 2013) modified according to experimental conditions described in Campbell et al. (2017). We sequenced these loci for 11 taxa and incorporated data for six additional taxa from previous studies described in Table S1 (Campbell et al., 2017; Faircloth et al., 2013).

Demultiplexed reads were trimmed of low-quality bases and filtered for minimum length with Trimmomatic version 0.32 using the following command line specifications and adapter fasta file available for download with the program: ILLUMINACLIP: TruSeq3-PE-splitAdapter.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 (Bolger, Lohse & Usadel, 2014). Sequences from a subset of samples were assembled with Velvet version 1.2.10 (Zerbino & Birney, 2008) to establish a range of suitable k-mer lengths for sequence assembly. We initially attempted optimal assemblies using various criteria with VelvetOptimiser version 2.2.5 (Gladman & Seemann, 2012). However, because some regions of the salmonid genome are tetrasomic, and thus do not meet the ploidy assumptions made in VelvetOptimiser, we established a range of k-mer lengths, 57–83 sites representing between 1/2 and 2/3 of total read length, and assembled salmonid sequences from all odd k-mer values within that range. We evaluated each of the resulting assemblies for the presence of UCEs using the PHYLUCE pipeline (Faircloth, 2015) and retained only the assembly that contained the maximum number of UCE loci, following the approach described in Campbell et al. (2017).

We identified homologous UCE loci and prepared them for alignment using the PHYLUCE pipeline. The PHYLUCE package screens out reciprocally duplicate enriched loci (i.e., potential paralogs) and performs a multiple sequence alignment of each of the remaining loci using MAFFT version 7.130b (Katoh et al., 2002; Katoh & Standley, 2013; Katoh & Toh, 2008). We screened tetrasomic loci in our dataset by searching previously published assemblies of Atlantic salmon (Salmo salar) (Lien et al., 2016) and rainbow trout (Pearse et al., 2019) for the UCEs sequenced in our study. Within Atlantic salmon, nine pairs of chromosome arms were indicated by Lien et al. (2016) that may be tetrasomic, seven of which were identified LORe regions by Robertson et al. (2017). From rainbow trout, seven chromosome pairs that were clearly tetrasomic in Pearse et al. (2019) were screened constituting homeologs to the seven LORe pairs screened in Atlantic salmon (Blumstein et al., 2019). BLASTN (Altschul et al., 1990) alignments of assembled UCEs that presented two alignments of >95% overlap and >97% similarity in the tetrasomic regions of either Atlantic salmon or rainbow trout were considered tetrasomic. UCEs that were not assembled and found to be single copy for all targeted species were also removed to create a “100% coverage alignment.”

To understand the effects of relaxing strict filtering of loci on phylogenetic inference, we generated four additional datasets: (1) a “95% coverage alignment” allowing one missing taxon per locus, (2) a “93% coverage alignment” allowing up to two missing taxa per locus, (3) a “75% coverage alignment” allowing up to five missing taxa per locus, (4) a “tetrasomic loci alignment” including only tetrasomically inherited loci permitting 5 missing taxa per locus.

Phylogenetic inference with the multispecies coalescent

We conducted joint GT–ST estimation of each of the five datasets with the StarBEAST2 Bayesian framework (Ogilvie, Bouckaert & Drummond, 2017) as implemented in BEAST version 2.4.8 (Bouckaert et al., 2014). For analysis in StarBEAST2, our data sets were partitioned by each UCE locus and the Akaike Information Criterion (AIC) was applied to select the best-fitting model of nucleotide substitution for each locus in each of our datasets using ModelGenerator version 0.85 (Keane et al., 2006). We ran several iterations of BEAST with varying length and sampling frequency and evaluated convergence of the runs and effective sample sizes (ESS) with Tracer version 1.6 (Rambaut et al., 2018) to determine proper sampling frequency and MCMC chain run length. We required that ESS values be 200 or greater to assure an adequate sample size for our posterior parameter estimations. Analysis of our initial BEAST runs showed that the estimates of parameter values failed to converge for datasets with more than 13 UCEs even after 300 million to 2 billion MCMC generations and showed evidence of over-parameterization, such as a lack of convergence in each chain and low consistency in LnL/prior values. To account for this lack of convergence in the larger datasets and test the effects of using a simpler nucleotide substitution model, we conducted an additional set of phylogenetic inferences in BEAST of all of our datasets and partition schemes thereof but specified the HKY+Γ nucleotide substitution model with four categories of gamma distributed rate variation and empirical base frequencies for each partition. For each of our analyses, we combined the independent chains and applied a burnin with LogCombiner version 2.1.3, and generated a maximum clade credibility tree using TreeAnnotator version 2.1.2, specifying a posterior probability (PP) limit of 0.50 and median node heights.

In addition to the joint GT–ST estimation conducted in BEAST, we also estimated the ML tree of each UCE locus and used summary coalescence analyses to estimate a species tree from the resulting set of gene trees a posteriori. We conducted this analysis in addition to the joint GT and ST analysis conducted in StarBEAST2 because this approach can show high accuracy and is computationally more tractable than a full probabilistic approach (Mallo & Posada, 2016). For each of our datasets, we estimated the ML tree of each UCE locus using RAxML version 8.0.19, specifying the GTR+Γ model of sequence evolution. We passed these trees to ASTRAL version 5.5.9 and MP-EST version 2.2.0, which each estimated a species tree from the collection of UCE loci.

Concatenated analyses

A concatenated maximum likelihood (ML) analysis was conducted with each data set (100%, 95%, 93%, 75% coverage alignments and the tetrasomic loci alignment). We conducted both partitioned (i.e., by gene, or objectively defined) and un-partitioned analyses under the GTR+Γ model of sequence evolution using RAxML version 8.2.10 (Stamatakis, 2014) with rapid bootstrap stopping. We used the program PartitionFinder version 1.1.1 (Lanfear et al., 2012) to generate an objective optimal partitioning strategy with partitions of individual UCE loci specified with the following options: the GTR+Γ model, greedy algorithm (Lanfear et al., 2014) and likelihoods estimated by RAxML (–raxml option). If the objective partition strategy differed from either unpartitioned or partitioned strategies, a third analyses of that data set was undertaken with RAxML.

Salmonid subfamily relationships

Though the composition and monophyly of salmonid subfamilies is well-established (Norden, 1961; Sanford, 1990), the relationships between them remain contentious. However, if we assume that each subfamily is truly monophyletic, and that their relationships take the form of a bifurcating tree, then there are only three possible arrangements of the salmonid subfamilies. This tractable number of potential combinations allows us to consider each possible arrangement of the taxa explicitly and test for the most preferred arrangement given our data. To avoid test-specific biases, we evaluated the support for each possible relationship of the subfamilies using four different approaches: (1) triplet analysis of independently estimated gene trees, (2), Bayes factors computed from the posterior trees of StarBEAST2 analyses, (3), a test under the multi-species coalescent with the program MP-EST, and, (4), approximately unbiased (AU) tests of concatenated trees.

Sequence alignment characteristics

Of the 500 UCE loci that we targeted, the number of total captured loci and the number of putatively single copy loci sequenced from each taxon varied (Table S1). For diploid outgroups, Esox and Umbra, 415 and 409 of the captured UCEs, respectively, were determined to be single-copy. Within salmonids, 166–275 of the captured loci were indicated to be single copy. Overall, we retained only eight single-copy UCE loci across all sampled taxa. Two of these loci were identified as coming from tetrasomically-inherited regions of the salmonid genome, and so were removed. Thus only six out of 500 candidate UCE loci passed the most stringent screening process and were present in all targeted taxa for a total of 2,331 aligned nucleotide sites (Table 1). The datasets that we generated when allowing for missing data yielded more retained UCE loci, ranging between 13 and 105, with aligned lengths between 5,297 and 36,835 aligned nucleotide sites (Table 1). We identified 33 tetrasomically inherited loci for analysis in the “tetrasomic loci alignment” (see Table 1) and this dataset contains a total of 12,758 aligned nucleotide sites. Across all datasets, the proportion of gaps and other ambiguously alignment regions ranged from 15.63% to 33.83%.

Model selection and MCMC convergence

For each locus in the “100% coverage” alignment (see Table 1), we initially specified the model of nucleotide substitution that received the lowest AIC score from ModelGenerator (Table S2). We obtained sufficient ESS and convergence of the parameter estimates with three independent MCMC chains of 400 million generations, sampled every 40,000 generations, with a 10% burn-in. For the “95% coverage” dataset, we specified nucleotide substitution models for each locus which received the lowest AIC score in our analysis using ModelGenerator. We obtained sufficient ESS and convergence with two chains of 50 million generations, sampled every 5,000 generations, with a 10% burn-in. For all remaining datasets, analyses with models of nucleotide substitution specified by ModelGenerator failed to generate sufficient ESS and convergence of parameter estimates in our MCMC chains. In response, we broadly reduced parameterization by simplifying nucleotide substitution models.

Analyses of the “93% coverage” alignment, the “75% coverage”, and “tetrasomic loci” alignments generated sufficient ESS and convergence of the parameter estimates with two independent MCMC chains of 200 million generations sampled every 20,000 generations, three chains of 800 million generations sampled every 80,000 generations, and four chains of 200 million generations sampled every 20,000 generations, respectively, and in each case we applied a 10% burn-in.

Phylogenetic inference

The results of joint GT-ST analyses are shown in Figs. 4 and 5 and summarized in Table 2. The clade Coregoninae + Thymallinae (C + T) is found in our 100% coverage alignment with ModelGenerator-specified nucleotide substitution models with strong support (0.98 PP) and in all analyses with reduced coverage datasets with very weak support (0.34–0.52 PP; see Table 2). Analysis of the tetrasomic loci alignment strongly supports Coregoninae + Salmoninae (C + S, 0.99 PP). The monophyly of Coregonus is not supported in the results of the 100% coverage species-tree, but is weakly supported (0.43–0.45 PP) in two of three of our reduced coverage analyses and strongly supported in our tetrasomic loci species-tree (1.00 PP). Salmo is inferred as the sister group of Parahucho in the 100% coverage species-tree with strong support (0.99 PP) and in all the reduced coverage analyses (0.98–0.99 PP). In contrast, the species tree inferred from the tetrasomic loci alignment places Oncorhynchus and Parahucho as sister lineages (0.38 PP).

The results of the summary coalescent analyses conducted in ASTRAL support a sister group relationship between Salmoninae and Thymallinae (4/5 analyses), but with low support (PP ranges from 0.38 to 0.57), and, from tetrasomic loci, Coregoninae + Salmoninae (0.91 PP, see Table 2). Monophyly of Coregonus is supported in the results of 3/5 of the analyses, though with only weak support (PP ranges from 0.68 to 0.74). Within Salmoninae, Parahucho is variously supported as sister to several different taxa with weak support (BPP < 0.70), and to a Oncorhynchus + Salvelinus clade with moderate support from tetrasomic loci (0.85 BPP, see Table 2). The results of the summary coalescent analyses conducted in MP-EST are inconsistent and reported in Table 2.

Concatenated analyses

Concatenated ML analyses supported either Coregoninae + Salmoninae or Salmoninae + Thymallinae relationships depending on the number of loci and partitioning strategy (Table S3, trees provided in Data Supplement). Partitioning strategy of alignments only altered inferred subfamily relationships in one case (95% coverage alignment), but caused substantial changes in bootstrap support (BS) values. Monophyly of Coregonus was only supported in three of thirteen analyses (BS range 70–100%), and the sister of Parahucho varied from Salmo (5/13 analyses, BS range 24–100%), to Oncorhynchus (5/13 analyses, BS range 50–99%) to Salvelinus (3/13 analyses, BS range 54–67%).

Salmonid subfamily relationships

A majority (50.79%) of trees generated in the triplet analyses support a Coregoninae + Thymallinae clade. The Salmoninae + Coregoninae (31.31%) and Salmoninae + Thymallinae (17.90%) clades are found in smaller but substantial proportions of the remaining trees. Bayes factors from the 100% coverage alignment indicate that a sister group relationship of Coregoninae and Thymallinae should be strongly preferred. Bayes factors for all other datasets when comparing C + T to alternatives do not provide strong support for C + T (Table 3). MP-EST analysis indicated that a sister group relationship for C + S was significantly supported via likelihood ratio tests for the tetratomic loci; no other data sets could distinguish among alternative topologies. The concatenated AU test results were not significant with all 100% and 95% coverage alignment tests, indicating that the alternative C + T hypothesis could not be rejected in analyses of the reduced data sets (Table S3).

Salmonid subfamily relationships

A sister group relationship of Coregoninae and Thymallinae (C + T) is broadly preferred to alternative topologies in our analyses (e.g., Fig. 5). This relationship is recovered in the joint GT–ST estimations of all data sets that exclude tetrasomically-inherited loci (PP = 0.34–0.98) and is also the preferred hypothesis indicated by Bayes factors and the triplet analysis. While no concatenated results support C + T, this result is obtained with as few as six UCE loci in joint GT-ST estimation with declining support as the number of loci increased. A Coregoninae + Thymallinae clade has been identified in previous molecular studies that either implicitly or explicitly account for potential polyploidy in their data. The C + T clade was recovered by analyses of mitochondrial genomes, which are not affected by issues of paralogy and can be expected to be less prone to ILS due to their haploid mode of transmission (Campbell et al., 2013; Horreo, 2017). Likewise, constructing phylogenies with both copies of Ss4R duplicated genes, an explicit approach, leads to support of a sister Coregoninae + Thymallinae (Macqueen & Johnston, 2014; Robertson et al., 2017).

Our results do not support the conclusions of previous studies that examined morphological characters of the salmonids and unanimously placed Thymallinae as the sister group to Salmoninae (Kendall & Behnke, 1984; Sanford, 1990; Stearley & Smith, 1993; Wilson & Li, 1999). Furthermore, our results combined with mitogenomic and other studies explicitly addressing polyploidy strongly support C + T over the two alternative hypotheses apparent in molecular phylogenetic studies of salmonids (Alexandrou et al., 2013; Crête-Lafrenière, Weir & Bernatchez, 2012; Shedko, Miroshnichenko & Nemkova, 2012). In light of growing molecular evidence conflicting with the morphology-based hypothesis, a critical reassessment of the morphological evidence that supports that relationship is needed but is beyond the scope of the present study.

Though C + T is the preferred grouping in our most stringent analyses, not all of the inference strategies implemented here supported this relationship. Analysis of the concatenated dataset and GT-ST estimation using ASTRAL produced inconsistent results and generally did not support the C + T clade. The inconsistent results between summary coalescence approach and joint estimation of the species tree and gene trees may indicate loss of a weak phylogenetic signal in the summarization step. Alternatively, the phylogenetic signal may have been blurred as a result of treating all UCE loci as a single partition evolving under a common substitution mode. The triplet analysis did not meet our prediction that, given three possible topologies, the true relationship should predominate and the alternative, incorrect arrangements should receive approximately equal portions of support. The triplet analyses yielded the putatively incorrect topologies at unequal frequencies (31.30% and 17.90%). This could be a result of the relatively small size of our “100% coverage” alignment or be a product of resampling selecting certain individuals frequently due to few representatives in some subfamilies (two thymallinine species are in the dataset).

Monophyly of Coregonus

Our results inconsistently support the monophyly of Coregonus (e.g., Figs. 4 and 5). Though the joint GT–ST analysis of tetrasomic loci and the ML analysis of the concatenated 100% coverage data set shows strong support for a monophyletic Coregonus, all other analyses that recovered the clade did so with weak support (PP < 0.75, BS = 71%, see Table 2). Even within a given analytical framework, support for monophyly of this genus varied. For example, MP-EST analyses showed a monophyletic Coregonus for the two most stringent datasets, but not in the analysis of any others. The uncertainty in this relationship within our data set may be derived from the relatively young age (i.e., ~10 mya, Horreo, 2017) of the clade that includes Stenodus and Coregonus resulting in low variation in the conserved UCE loci that were indentified to be single copy across a broad range of salmonid species. Additionally, during the analytical study design of this manuscript, a chromosomal-level assembly of a coregonine was unavailable leading to uncertainty in the tetrasomic status of loci in this subfamily based of off salmonine genome assemblies.

Previous studies examining the monophyly of Coregonus with respect to Stenodus have yielded conflicting results (Bernatchez, Colombani & Dodson, 1991; Bodaly et al., 1991; Crête-Lafrenière, Weir & Bernatchez, 2012; Horreo, 2017; Sajdak & Phillips, 1997; Vuorinen et al., 1998). Most recently, analysis of mitogenomic data has shown Coregonus to be monophyletic, with Stenodus as its sister taxon (Horreo, 2017). However, the phylogenetic placement of the enigmatic Coregonus huntsmani, a divergent member of the genus, remains untested in any mitogenomic or phylogenomic study.

Crête-Lafrenière, Weir & Bernatchez (2012) found C. huntsmani to be the sister lineage to Stenodus and all other Coregonus species and earlier isozyme and mitochondrial restriction data show C. huntsmani to be distinct from the other main Coregonus subgenera (Coregonus and Leucichthys) but did not include Stenodus in the sampling (Bernatchez et al., 1991). We examined publicly available COI data from coregonines including two haplotypes from C. huntsmani and found that species to be sister group to a monophyletic Coregonus with low support (BS = 58%, alignment tree and methods are supplied in the Data Supplement). At present molecular data from C. huntsmani are very limited, however including this taxon in future datasets will be critical for conclusively testing the monophyly of Coregonus and its relationship with Stenodus. Further investigations of phylogenetic relationships of Coregonus and Stenodus should also consider the role and prevalence of hybridization within Coregonus and between Stenodus and Coregonus, (McPhail, 2007) as the signatures of hybridization and ILS are difficult to disentangle and our analyses assumed such signatures derived from ILS (Yu et al., 2011).

Placement of Parahucho

The Sakhalin taimen Parahucho perryi is a rare species and the sole member of its genus (Rand, 2006; Vladykov & Gruchy, 1972). This critically endangered fish occupies a limited geographic range and has a narrow range of suitable habitats in Northern Japan and the Russian Far East (Fukushima et al., 2011; Kimura, 1966; Rand, 2006). The placement of Parahucho as sister to Salmo is strongly supported by all joint GT–ST analyses of non-tetrasomic loci (4 of 4, 0.98–0.99 PP) and by concatenated analysis of our 100% and 95% coverage alignment (BS = 66–100%). The original description of Sakhalin taimen placed it as a congener of Atlantic salmon (Salmo salar), as Salmo perryi (Brevoort, 1856). Subsequently, the similarities of Sakhalin taimen to the European huchen Hucho hucho were noted, and these two species were placed in the same genus, Hucho (Günther, 1868; Jordan & Snyder, 1902). Hucho and Parahucho are both large bodied piscivores, but differ in that the latter has 10 less vertebrae, basibranchial teeth, a median set of teeth on the supralingual, and 70–100 less scales in the lateral line (Vladykov, 1963; Vladykov & Gruchy, 1972). Indeed, the two genera are not closely allied in any of the results presented here, nor has a close relationship of the two genera been reported in any other molecular phylogenetic studies (Campbell et al., 2013; Crespi & Fulton, 2004; Matveev, Nishihara & Okada, 2007; Oakley & Phillips, 1999). Though there has been considerable disagreement in the placement of Parahucho, the emerging consensus from nuclear gene sequences supports a Parahucho + Salmo sister relationship (Lecaudey et al., 2018), as is recovered here.

Phylogenetic results from disomic and tetrasomic loci

Resolving subfamily relationships was a primary aim of this study, but our results repeatedly showed a conflict in this segment of the phylogeny between analyses of the most stringently assembled dataset (i.e., 100% coverage) and of the remaining alignments, and particularly of the tetrasomic loci. For example, in GT–ST analyses conducted with StarBEAST2 of the 100%, 95%, 93% and 75% coverage alignments, Coregoninae and Thymallinae are supported as sister lineages. Analyses of the 100% coverage data set yield very strong (0.98 PP) support for that grouping, which declines substantially in analyses of the other datasets (0.34–0.52 PP). On the other hand, the 33 tetrasomic loci we removed strongly supported a Coregoninae + Salmoninae clade (0.99 PP). Likewise, we observed strong support for a sister relationship between Salmo and Parahucho in the more reduced datasets (0.98–1.00 PP), and less support in the more extensive 75% alignment (0.94 PP). The topology of the species tree inferred from StarBEAST2 analysis of tetrasomic loci indicates that Salmo and Parahucho are not closely related (Fig. 4D). In summary, two of three sets of alternative placements we considered in detail (i.e., subfamily relationships, Coregonus monophyly, and placement of Parahucho) in this study yield conflicting results when comparing datasets of loci with differing ploidy levels. We also find that support values for the C + T clade decline with greater proportions of missing data in the underlying dataset.

We attempted to identify tetrasomic loci in the sequenced dataset through the assembly of a single contig matching a UCE locus that was placed within genomic locations of known tetrasomy in rainbow trout and Atlantic salmon. While there is evidence for broad conservation of homologous tetrasomic genomic regions in salmonids (Brieuc et al., 2014; Kodama et al., 2014; Lien et al., 2011; Sutherland et al., 2016), there is also evidence of substantial lineage-specific rediploidization in salmonids and polyploids in general taxa (Robertson et al., 2017; Spoelhof, Soltis & Soltis, 2017). Our definition of tetrasomic regions from two salmonines is taxonomically restrictive and may have allowed the inclusion of regions of tetrasomic inheritance in other species sequenced. Consequently, increasing the proportion of missing data may have introduced paralogous loci into our analyses due to lineage-specific rediploidization processes (e.g., Fig. 3B).