The grass subtribe Loliinae has great ecological and economic importance, as it includes community-dominant species of mountain grasslands and the most extensively cultivated pasture, fodder, and turf grasses (fescues, ryegrasses). Resolving the phylogeny of recently evolved Loliinae lineages has proven challenging due to frequent introgressions and polyploidizations that occurred throughout their history. Here we present the first target-capture phylogeny of Loliinae using 270 orthologous single-copy nuclear coding loci for a large sample of 132 representative taxa, covering all its 29 evolutionary lineages. Additionally, we assembled plastome sequences to complement the inferred hybrid speciation history of the Loliinae. Concatenated maximum likelihood and multispecies coalescent trees of ortho-homeolog single-copy genes showed well-supported relationships for major lineages, which were generally consistent across analyses and genomes, and with previous taxonomic and phylogenetic findings. However, they also revealed high levels of both nuclear and cytonuclear discordances estimated to be caused by hybridization and incomplete lineage sorting. We complemented this with a phylogenetic network analysis with representative samples from the main clades to infer reticulation events in the evolution of these grasses.  Furthermore, we performed gene tree – species tree reconciliation methods using gene duplication and loss models and multi-labeled trees for polyploidy analysis to estimate the proportion of duplicated genes and the nature of polyploidy of the major Loliinae lineages. These analyses agreed that the Fine-Leaved (FL) Loliinae clade could have evolved from hybridization between more ancestral Broad-Leaved (BL) Loliinae lineages and that both groups underwent ancestral and recent hybridizations. The time-calibrated phylogeny of for the main Loliinae clades supports an early Miocene origin for Loliinae and Mid-Late Miocene splits for its main BL and FL lineages, while current species-rich groups radiated in the Late Miocene. Hybridization tests of nuclear data and topological incongruence assays between the nuclear single-copy genes and plastome-based trees using various approaches and different sampling subsets confirmed the rampant hybridization experienced by Loliinae at deep and shallow nodes. However, hybridization rates differed from lineage to lineage within the major clades and were not correlated with time or ploidy level, but rather depended on their different propensities to hybridize with species within and/or outside their own clade. Our analyses detected high hybridization rates in four broad-leaved (Subulatae-Hawaiian, Tropical-South African, Mexico-Central-South American (MCSA I-II), and Leucopoa p. p.) and five fine-leaved Loliinae lineages (American II, Aulaxyper, Afroalpine, American-Neozeylandic, Australia-Tasmania) containing rogue species that probably originated from trans-clade crosses and are more likely to hybridize greatly. In contrast, they recovered low hybridization rates in four broad-leaved (Schedonorus-Lolium, Subbulbosae, Drymanthele-Pseudoscariosae-Lojaconoa, Leucopoa p. p.) and six fine-leaved lineages (Festuca, Psilurus-Vulpia, American I, Exaratae-Loretia, American-Vulpia-Pampas, Eskia), with species derived from single ancestors that hybridize only with close congeners. Inferences from gene duplications and allopolyploidizations, along with inheritance probabilities from the phylogenetic network, point to the BL Loliinae lineages as ancient hybrids or paleo-allopolyploids, while the FL lineages, especially those of the core FL clade, correspond to more recent meso- or neo-allopolyploids.

Sampling

Samples of 132 species from eight genera of Loliinae (Festuca, 118; Lolium, 5; Vulpia, 4; Micropyropsis, 1; Megalachne, 1; Wangenheimia, 1; Psilurus, 1; Hellerochloa 1) and three related outgroup grasses (Oryza sativa, Hordeum vulgare, Alopecurus aequalis) were included in the study (Supplementary Table S1; Supplementary Fig. S1). The 132 Loliinae samples were first analyzed to capture single-copy nuclear gene targets; 79 of them were analyzed for plastome data for the first time, and 30 were not previously included in phylogenetic studies (Supplementary Table S1). Sampling was carried out from fresh materials collected in the field, collections of silica-dried leaves, and herbarium vouchers (Supplementary Table S1). Plastome data retrieved from 49 previously studied taxa (Moreno-Aguilar et al. 2022a; Moreno-Aguilar et al. 2022b) were incorporated into the analysis (Supplementary Table S1). The selected taxa represent the 28 evolutionary lineages currently recognized within Loliinae (Minaya et al., 2017; Moreno-Aguilar et al., 2022b).

Target gene capture and genome skimming sequencing

Genomic library preparation and targeted enrichment of Loliinae single-copy nuclear genes were performed with the Angiosperms353 target capture kit at Arbor Biosciences (Michigan, USA), following the protocols outlined in Johnson et al. (2019). Targeted sequence capture was performed in pools of 12 libraries following the myBaits v5 manual (https://arborbiosci.com/mybaits-manual/). Capture reactions were pooled in equimolar ratios to form a sequencing pool, which was sequenced on a partial lane of the Illumina NovaSeq 6000 platform in paired-end (PE) mode (2 x 150 bp). Output reads from target enrichment sequencing were checked with FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and trimmed with Trimmomatic (Bolger et al. 2014), preserving reads with a minimum length of 50bp (minlen:50) to reduce the risk of potential misalignments of short reads to genes and trimming bases at the end of the reads with low quality (trailing:30).

Genome skimming sequencing was performed on 79 Loliinae DNA samples at the Centro Nacional de Análisis Genómico (CNAG) and Macrogen (Spain) to assemble their plastome sequence (Supplementary Table S1). Genomic sequencing of a multiplexed pool of PCR-free KAPA libraries was performed on a HiSeq4000 or HiSeq 2500 (TruSeq SBS Kit v4, Illumina, Inc) in PE mode (2 x 101 bp) as described in Moreno-Aguilar et al. (2020). Illumina PE reads were quality-checked using FastQC, and adapters and low-quality sequences were trimmed and removed with Trimmomatic. Loliinae genomic samples used in subsequent analysis contained between 1.5 – 35.0 million reads (average 16.0 million reads) with insert sizes ranging between 69 – 260 bp (Supplementary Table S1).

Sequence assemblies of target capture data, gene orthology assessments, and hybrid-class detection (AD and LH metrics)

Loliinae single-copy nuclear genes were recovered from the target-enriched PE reads using the HybPiper v.2.3.2 pipeline (Johnson et al. 2016). This was done by mapping the filtered reads against the template sequences of the 353 low-copy nuclear genes (available at https://github.com/mossmatters/Angiosperms353https://github.com/mossmatters/Angiosperms353) using the Burrows-Wheeler Alignment (BWA v.0.7) (Li et al. 2009), and then through de novo assembly of mapped reads for each gene separately using SPAdes v. 3.13 (Bankevich et al. 2012), with the default minimum coverage threshold of 8× (Supplementary Fig. S2). The HybPiper script paralog_retriever was used to obtain a data set from Loliinae, made of nuclear single-copy genes. Following this initial HybPiper assembly step, a total of 352 genes were recovered across the 132 Loliinae samples studied. This unfiltered dataset (scg-raw dataset; Appendix 1), which contains all genetic and allelic diversity, was used for gene duplication analyses and to assess the nature of polyploidy in Loliinae lineages. The homology of these 352 nuclear genes was circumscribed using Orthofinder v. 3.0.1 (Emms and Kelly 2015) to incorporate the coding sequence from the reference genomes of the diploid outgroups Oryza sativa (NCBI accession GCF_034140825.1), Hordeum vulgare (subsp. vulgare) (GCA_904849725.1), and Alopecurus aequalis (GCA_964340505.1), and the diploid ingroup Lolium perenne (P226/135/16) into the Loliinae gene capture dataset. This initial filtering of homologous genes was necessary to select appropriate outgroup reference sequences for the target genes studied in the Loliinae dataset and did not affect subsequent tree-based orthology filtering approaches, which focused primarily on ingroup sequences. We then used ASTRAL-Pro (Zhang et al., 2020) to construct a species tree using 325 of the 352 total loci that meet the constraint of fewer than 27 (20% of total taxa) missing taxa. We filtered the 352 groups based on this criterion and the constraint of the maximum number of copies of the same taxon less than 5, resulting in a total of 302 homologous clusters. These clusters (scg-full dataset; Appendix 2) were used for gene tree - species tree reconciliation analysis using Notung and GRAMPA (see below). They were further analyzed using the Yang and Smith (2014) pipeline to obtain orthologs for all Loliinae samples under study. We initially tested all three alternative tree-based orthology inference approaches, including maximum inclusion (MI), monophyletic outgroups (MO), and one-to-one (1-to-1) methods (Yang and Smith 2014). For all these analyses, the minimum number of taxa was 108 (allowing a maximum of 27 missing taxa), and for the MI method, the long tip cutoff was set at “0.2”. We selected the 270 genes filtered by MI as a strict dataset of reliable orthologous single-copy genes (Supplementary Table S2). This was because the 270 gene trees retrieved from MI (scg-strict and inclusive datasets; Appendix 3) yielded a consensus phylogeny regarding the placement of major clades (concatenated IQ-TREE2 and ASTRAL trees) similar to those obtained from the 184 gene trees retrieved from MO and the 175 genes retrieved from 1-to-1 (see Results and Appendices 3 and 4), but with the highest average bootstrap support percentages or posterior probability among the three datasets.  Furthermore, since the MO process retrieves MO genes and 1-to-1 genes simultaneously, and of the 184 genes retrieved by the MO process, 175 overlap with the 1-to-1 process, we will only show the result of the MO analysis in the supplementary files (see below).

For the HybPiper 270 MI-filtered orthologous single-copy-genes data set (Supplementary Table S2), individual genes were aligned with MAFFT v.7.490 (Katoh and Standley 2013) using the iterative refinement method --maxiterate 1000. Empty genes for any Loliinae sample were removed from downstream analysis. trimAl v. 1.2rev59 (Capella-Gutiérrez et al. 2009) was used to remove sequences with insufficient coverage (-automated1) in well-occupied columns of each gene alignment.  Gene alignments were visually inspected with Geneious Prime to detect potentially misaligned sequences. Sequences with more than 60% of the total alignment length consisting of gaps were removed to enhance the quality of the data sets, filtering out species with potentially low phylogenetic information. This process generated the respective multiple sequence alignments (MSAs). To check alignment quality, we estimated summary statistics of gene alignments using AMAS v.0.98 (Borowiec, 2016). These statistics included alignment length, missing data, and the number of parsimony informative sites.

We used HybPhaser v2.0 (Nauheimer et al. 2021) to detect potential hybrids in our 270 orthogroups data set of Loliinae. This tool maps raw sequence data to the HybPiper contigs, taking into account SNP variation using nucleotide ambiguities and quantifying divergence between gene variants. The gene variants correspond either to paralogs or to hybrids (“homeologs” in allopolyploids; Sancho et al. 2022) (Supplementary Table S3). Samples showing high SNP content in genes are likely hybrids or polyploids and are phased for the different copy variants (alleles) into phased haplotypes (Nauheimer et al. 2021; Hendriks et al. 2023). We ran HybPhaser to detect potential hybrid/polyploid samples without ruling out putative “paralogs” and generated a data set of phased accessions of Loliinae (scg-inclusive data set) (Supplementary Table S3, Appendix 3). We used R scripts to compute two metrics indicative of hybridization, sample allele divergence (AD, percentage of SNPs in all genes) and locus heterozygosity (LH, percentage of genes with SNPs). We classified the Loliinae samples into three LH-vs-AD hybrid-type classes [1) low-hybrid level: low-to-medium LH (≥10 - ≤82.4) and low AD (>0 - ≤1.7); 2) medium-hybrid level: medium LH (≥70 - < 90) and medium AD (>1.71 - ≤2.8); 3) high-hybrid level: high LH (≥80) and high AD (>2.8)], which are biologically meaninful and reflect their hybridization history (Nauheimer et al. 2021; Hendriks et al. 2023). 

Sequence assemblies of plastomes

Whole-plastome sequences for 79 new Loliinae samples were assembled from their respective genome skimming PE reads, following the procedures indicated in Moreno-Aguilar et al. (2020). Plastome assembly was performed with Novoplasty v.2.7.1 (Dierckxsens & al., 2017) using the plastome of Festuca pratensis (JX871941) as reference** sequence and standardized parameters (k-mer: 30-39, insert size: ∼69-200 bp, genome range: 120,000–140,000 bp, and PE reads: 101-150 bp). Furthermore, to recover plastome sequences in data with a low number and quality of total PE reads, plastome assembly was performed using a read-mapping strategy to, respectively, closely related Festuca plastomes using Geneious Prime (Supplementary Table S1). The data for 49 representative Loliinae taxa from previous studies (Moreno-Aguilar et al., 2020; Moreno-Aguilar et al., 2022a, 2022b) were incorporated into the existing plastome dataset (Supplementary table S1, Appendix 5). Whole plastomes sequences were aligned separately with MAFFT v.7.031b, and trimAl v. 1.2rev59 was used to remove poor quality regions from the MSA with the -automated1 parameter.

Phylogenomic reconstructions and intragenomic nuclear discordance

We used three different approaches to analyze our genomic datasets and reconstruct Loliinae phylogenies. First, trimmed MSAs from the Loliinae scg-strict dataset (270 MI-filtered genes; Supplementary Table S2, Appendix 3) and different numbers of MO- (184) filtered genes (Appendix 4) were concatenated into their respective supermatrices. Subsequently, they were used for maximum likelihood (ML) phylogenetic reconstruction using IQ-TREE2 v. 2.2.2.6 (Nguyen et al., 2015; Minh et al., 2020a), applying the “proportional to edges” partitioning model with individual loci as partitions and model selection implemented using ModelFinder (Kalyaanamoorthy et al., 2017) for each partition. Branch support was calculated from 1000 UltraFast Bootstrap replicates (Hoang et al., 2018) and gene concordance factors (gCF) and site concordance factors (sCF, parameter –scf 1000; Minh et al., 2020b) for all nodes. Second, individual ML trees of each single-copy nuclear gene, as well as ML trees of the entire plastome MSA, were calculated separately using the same IQ-TREE2 procedure. Oryza sativa, Hordeum vulgare, and Alopecurus aequalis were used to root all trees. Third, to consider potential ILS events between closely related Loliinae lineages, a species tree was inferred by analyzing the Loliinae 270 scg-strict dataset under the multispecies coalescence (MSC) using the ASTER package v.1-3.5 ASTRAL (Zhang et al., 2025), which was fed with the individual ML gene trees from IQ-TREE2. All gene trees were rooted using Oryza sativa, Hordeum vulgare, and Alopecurus aequalis, with the pxrr function in Phyx (Brown & al., 2017), except for four trees that were rooted using the early divergent lineages of Festuca lasto and F. drymeja as the outgroup. Branches with bootstrap support (BS) values <30% in the gene trees were collapsed using nw_ed from Newick Utilities 1.6.0 (Junier and Zdobnov, 2010). To estimate intragenomic discordance in the nuclear dataset, we analyzed in R the normalized quartet scores generated for the main topology and the first and second alternative topologies when inferring the MSC species tree. An ASTRAL topology was also generated for the MO-(184) dataset. While the MSC model is only consistent when ILS is the sole source of gene tree discordance, comparison with the supermatrix ML tree can identify highly supported clades (Nauheimer et al. 2021).

To further characterize the impact of ILS versus introgression/ hybridization (IH) in Loliinae lineages, we estimated the ILS and IH indices using Phytop v.1.0 (Shang et al. 2025), a pipeline that quantitatively assesses ILS/IH magnitudes based on the nodal quartets’ supports (q1, q2, q3) derived from an ASTRAL nuclear species tree inferred from 270-MI individual ML gene trees. The inferred species tree served as the reference topology for Phytop comparative analysis, allowing for the quantification and visualization of Loliinae lineage-specific ILS/IH values using topology matching algorithms (Shang et al. 2025; Chen et al. 2025).  

Species network analysis, whole genome duplications, and polyploidy events

To assess the impact of reticulated evolution on the observed conflicting nodes of the Loliinae nuclear coalescent tree, we applied a phylogenetic network approach based on a maximum pseudo-likelihood method (Yu and Nakhleh, 2015) implemented in PhyloNet v. 3.8.0 (Than et al., 2008; Wen et al., 2018) using the 270 scg-strict dataset. This approach considers ILS in the coalescent model and hybridization at the reticulating nodes of the network (Keuler et al., 2020). Due to the computational burden, we restricted our searches to 20 species representative of major Loliinae lineages recovered in the IQ-TREE2 and ASTRAL trees (BL: Leucopoa, Tropical-South African, MCSA I, MCSA II, Schedonorus-Lolium; FL: Eskia, American-Neozeylandic, American I, American II, Subulatae-Hawaiian, Afroalpine, Exaratae-Loretia, Aulaxyper, Festuca) and Oryza sativa as outgroup. We performed network searches based on the reduced 270 scg-strict gene trees. Missing taxa were also allowed for some orthologous gene trees in this 20-taxa dataset. Network searches were performed with the allowed networks (1 for the first search, 2 for the second, and subsequently up to 10 for the tenth). For each search, 10 runs were performed with the parameter "-x 10". The best-fitting networks obtained for each allowed reticulation model were further optimized for branch lengths and inheritance probabilities using a likelihood framework, and the best-fitting network was selected using the bias-corrected Akaike Information Criterion (AICc, Sugiura 1978) and the Bayesian information criterion (BIC, Schwarz 1978) by applying the number of parameters as the number of branches plus the number of inheritance probabilities [(2n-3)+(2h), where h is the number of hybridizations in the network], following Morales-Briones et al. (2021). Complementarily, we applied a gene tree reconciliation method to detect gene duplication events in the ASTRAL phylogeny of Loliinae, applying the gene duplication and loss model implemented in Notung v. 2.6 (Stolzer et al., 2012). The 302 individual ML gene trees (scg-full dataset; Appendix 2) were used as input with a duplication cost of 1.5, following Koenen et al. (2021). To account for gene tree estimation error, we also filtered out genes whose outgroup did not form a monophyletic clade and whose average bootstrap was less than 60, and rerun Notung with these remaining 167 trees with the same parameters. The 167 filtered individual ML genetic trees were used as input for Notung, following the same procedure for those 302 individual gene trees. The proportion of duplicated genes for each branch was tallied from 302 homologous gene trees and 167 homologous gene trees, respectively. In addition, we assessed whether branches with more than 10% duplicated genes could have undergone allopolyploidization, autopolyploidization, or no-polyploidization, and the identification of the possible parental lineages that formed the allopolyploids, using the species tree reconciliation approach from gene trees to multi-labeled trees with GRAMPA v. 1.4.0 (Thomas et al. 2017). GRAMPA performed the least common ancestor reconciliation of all gene family trees against both the singly-labeled species tree and all possible multi-labeled (MUL). The program assigned a reconciliation score for each MUL considered and the origin of homologs, scoring polyploid lineages (H1 nodes) together with the possible locations of the second parental lineage (H2 node) to infer genetic origin (Yang et al., 2023). We tested possible WGD scenarios with the lowest parsimony scores by comparing the scores obtained for the singly-labelled species tree (not WGD) with those of the MUL tree (the node subject to WGD forms sister branches with descendant paralogs (autopolyploidization) or not (allopolyploidization)) to decide whether WGD was supported by GRAMPA and what type of polyploidization was inferred. In our case, we tested the potential polyploid nodes identified by Notung with more than 10% of gene duplications with all the nodes along the phylogeny, and to incorporate more gene copy variation, we used the parameter “-c 14” when running GRAMPA.

Evaluation of cytonuclear topological discordance and dating analysis.

To assess the degree to which the phylogeny retrieved by the nuclear coding genome agrees with that of the organelle (plastome) genome, we also applied the Procrustean Approach to Cophylogeny (PACo) procedure implemented in R (Balbuena et al. 2013; Hutchinson et al., 2017) to the nuclear 270 scg-strict supermatrix ML trees versus the plastome supermatrix ML trees (bootstrap trees) of 128 common Loliinae species. PACo is a global-fit method used to detect coevolutionary associations between different organisms or different nuclear and (endosymbiotic) organellar genome trees of plants (Balbuena et al. 2013; Pérez-Escobar et al. 2016). It assesses phylogenetic congruence with the explicit goal of testing the dependence of one phylogeny on another, accommodating large-scale data and allowing for the examination of phylogenetic congruence between highly enriched clades (Hutchinson et al., 2017). Due to the high degree of intragenomic congruence of the Loliinae plastome-coding genes (Minaya et al. 2017; Moreno-Aguilar et al. 2020), whole-plastome alignment was employed to generate 1,000 bootstrap plastome trees. To do this, a phylogenetic reconstruction was carried out using IQ-TREE2, providing the program with partition information for each nuclear and plastome alignment, executing 1,000 ultrafast bootstrap replicates, and saving the respective bootstrap trees (ufboot). This procedure evaluates the similarities between any pair of topologies by comparing the Euclidean distances that separate the terminals in both trees through the Procrustean superposition (Balbuena et al., 2013). The sum of the squared residuals (the disparity between an observed value and a fitted value derived from a model) for each association and each pair of topologies evaluated can be interpreted as a concordance score because it is directly proportional to the magnitude of the topological conflict for the pair of terminals considered (Pérez-Escobar et al. 2016). Differences in terminal position between nuclear and plastome ML trees were summarized in bar plots using the R package ggplot2 (Wickham 2016). The sum of squared residuals for each pair of nuclear and plastome gene terminals was classified into quartiles, and the magnitude of discordance was assessed by the proportion of genes binned in quartiles 3 and 4 (50% and 75%) in each terminal. The higher the number of genes binned in these quartiles, the greater the discordance is (Pérez-Escobar et al. 2021). Furthermore, we examined the congruence of nuclear vs. plastome topologies using generalized Robinson-Fould (gRF) distances, a metric that allows matching of similar splits in the compared trees, giving similarity scores to each pair of splits and the overall similarity measure. gRF distances were calculated for the Loliinae tree and for the major BL and FL lineages using the R package TreeDist, following Smith (2020) and Hendriks et al. (2023). 

Due to the relatively high cytonuclear concordance detected for the main backbone Loliinae topologies, we dated the origins of the Loliinae lineages using the concatenated IQ-TREE2 ML phylogeny with branch length information and the treePL software (Smith and O’Meara 2012). We constrained four nodes of the tree using dates inferred from our previous studies. Three nodes were calibrated using secondary age constrains for the crown nodes of the BOP (Oryza + Pooideae) clade (normal prior maximum = 55.08, minimum = 47.76 Ma), Loliinae clade (normal prior maximum = 23.13, minimum = 16.02 Ma), BL Loliinae clade (normal prior maximum = 20.45, minimum = 12.15 Ma), and a fourth node was calibrated using a minimum age constrain for the crown node of fine-leaved Loliinae (lognormal prior maximum = 19.75, minimum = 14.13 Ma) based on a Festuca sect. Festuca leaf macrofossil dated to the Early Miocene (Moreno-Aguilar et al., 2020). The analyses were performed using smoothing values of 0.1 to estimate the best optimization parameter values (3-3-5), which were then used in a subsequent analysis to calculate divergence times. Confidence intervals of nodal ages estimates were calculated from 1000 bootstrap replicates, summarized in TreeAnnotator 2.7.7 (Bouckaert et al., 2014) (Appendix 7).