Incomplete lineage sorting and introgression among genera and species of Liliaceae tribe Tulipeae: insights from phylogenomics

Abstract

Background: Phylogenetic research in Tulipa (Liliaceae), a genus of significant economic and horticultural value, has relied on limited nuclear (mostly nuclear ribosomal internal transcribed spacer, nrITS) and plastid DNA sequences, resulting in low-resolution phylogenetic trees and uncertain intrageneric classifications. The genus, noted for its large genome, presents discordant relationships among Amana, Erythronium, and Tulipa, likely due to incomplete lineage sorting (ILS) and/or reticulate evolution. Thus, phylogenomic approaches are needed to clarify these relationships and the conflicting signals within the tribe Tulipeae.

Results: We newly sequenced 50 transcriptomes of 46 species of tribe Tulipeae (including multiple accessions of all four genera) and one outgroup species of the sister tribe Lilieae (Notholirion campanulatum), and downloaded 15 previously published transcriptomes of tribe Tulipeae to supplement the sampling. One plastid dataset (74 plastid protein-coding genes, PCGs) and one nuclear dataset (2594 nuclear orthologous genes, OGs) were constructed, with the latter used for species tree inference based on maximum likelihood (ML) and multi-species coalescent (MSC) methods. To investigate causes of gene tree discordance, "site con/discordance factors" (sCF and sDF1/sDF2) were calculated first, after which phylogenetic nodes displaying high or imbalanced sDF1/2 were selected for phylogenetic network analyses and polytomy tests to determine whether ILS or reticulate evolution best explain incongruence. Key relationships not resolved by this technique, especially those among Amana, Erythronium, and Tulipa, were further investigated by applying D-statistics and QuIBL.

Conclusions: We failed to reconstruct a reliable and unambiguous evolutionary history among Amana, Erythronium, and Tulipa due to especially pervasive ILS and reticulate evolution, likely caused either by obscured minority phylogenetic signal or differing signals among genomic compartments. However, within Tulipa we confirmed the monophyly of most subgenera, with the exception of two species in the small subgenus Orithyia, of which Tulipa heterophylla was recovered as sister to the remainder of the genus, whereas T. sinkiangensis clustered within subgenus Tulipa. In contrast, most traditional sections of Tulipa were found to be non-monophyletic.

Keywords: Incomplete lineage sorting (ILS); Introgression; Phylogeny; Reticulate evolution; Transcriptome sequencing; Tulipeae.

Fig. 1

Representative members of Liliaceae tribe Tulipeae. A Gagea granulosa. B G. serotina (= Lloydia serotina). C Amana edulis. D Erythronium sibiricum. E Tulipa heterophylla. F, T. heteropetala. G T. uniflora. H T. clusiana (subgenus Clusianae). I T. montana (subgenus Clusianae). J T. dasystemon (subgenus Eriostemones). K T. saxatilis (subgenus Eriostemones). L T. turkestanica (subgenus Eriostemones). M T. sinkiangensis. N T. kaufmanniana (subgenus Tulipa). O T. thianschanica (subgenus Tulipa). P T. altaica (subgenus Tulipa). Photograph H by Xiang-Yu Zhou, all others by Xin Zhong

Fig. 2

ML tree based on 74 PCGs. Support values are marked in the form of SH-aLRT/bootstrap value. The support values from low to high corresponds to the gradient color from red to blue

Fig. 3

ML tree based on different datasets. Support values are marked in the form of SH-aLRT/bootstrap value. Support values from low to high corresponds to the gradient color from red to blue. E1, E2, T1, T2, alt, acu, and eic are shorthands for specific clades. The phylogenetic relationships different from that of original tree (2594 OGs) are marked with red dotted lines. a1 Phylogenetic relationships reconstructed from concatenated 2594 OGs, among them, 2037 OGs were considered suitable for phylogenetic analysis by IQ-TREE2, thus a2 was also generated. b Phylogenetic relationship reconstructed by 1594 OGs with relatively low copy number (no topology change). c Phylogenetic relationship reconstructed by 1000 OGs with relatively high copy number

Fig. 4

Coalescent trees based on different datasets. ASTRAL support values (local posterior probability, LPP) of each node are marked. The values from low to high corresponds to the gradient color from red to blue. E1, E2, T1, T2, alt, acu, and eic are shorthands for specific clades. h, sc, T1, T2, t1, and t2 are shorthands for specific clades. The phylogenetic relationships different from that of original tree (2594 OGs, BS > 10) are marked with red dotted lines. a1–a5 Phylogenetic relationships reconstructed from 2594 gene trees with five different threshold bootstrap values ranging from 10 to 50. b Phylogenetic relationships reconstructed by 1594 OGs with relatively low copy number. c Phylogenetic relationship reconstructed by 1000 OGs with relatively high copy number (no topology change)

Fig. 5

ML trees based on 2594 OGs. All nodes whose sDF1 and sDF2 are imbalanced or whose sCF is lower than 40 are marked the scale of site (dis)concordance factors as pie charts. The nodes with AGTs are denoted by red letters

Fig. 6

Coalescent trees based on 1594 OGs (BS > 10). All nodes whose sDF1 and sDF2 are imbalanced or whose sCF is lower than 40 are marked the scale of site (dis)concordance factors as pie charts; One node with AGTs is denoted by red letters

Fig. 7

Four nodes with P-value > 0.01 in polytomy tests

Fig. 8

Best networks predicted byPHYLONETWORKS. The gene flows involving problematic taxa are the focuses, which are shown in red. The unimportant gene flows are represented in blue and not highlighted. The gene flows appeared many times in suboptimal networks are also displayed in form of green arc dotted line. Those values next to each hybrid edge indicate inheritance probabilities