Loricaria is a genus restricted to the high elevation habitats of the tropical Andes. The genus comprises small dioecious shrubs with scale-like leaves, which is a morphological convergence to certain gymnosperm genera (Cuatrecasas, 1954). Currently, 19 species are accepted, 17 as a result of the synopsis of the genus by Cuatrecasas (1954), and two species described by Dillon and Sagastegui Alva (1986), and Hind (2004). Morphological investigations of herbarium specimens plus signatures of potential cryptic speciation within L. thuyoides (Lam.) Sch. Bip. (Kolář et al., 2016) suggested that there is a substantial degree of taxonomic uncertainty and that there potentially exist more species than are described to date.

The genus has been divided into three different sections (Table 1), primarily based on the position of the flower heads, which is axillary in sect. Thyopsis and sect. Graveoleum and terminal in sect. Terminalia. Section Graveoleum is differentiated from the other two sections by a glandulose-pilose ovary and glandulose-pubescent leaves. Loricaria graveolens (Sch. Bip.) Wedd. is the only member of sect. Graveoleum. Distribution information and section assignation are taken from Cuatrecasas (1954), Dillon and Sagastegui Alva (1986), and Hind (2004).

The genus belongs to the tribe Gnaphalieae, which has its species richness concentrated in the southern hemisphere (Nie et al., 2016). Hybridization has been inferred for this tribe (Galbany-Casals et al., 2014; Barker et al., 2016), thus phylogenetic discordance can be expected. Further, within the tribe Gnaphalieae the most recent common ancestor (mrca) of the FLAG-clade, which includes Loricaria [defined in Galbany-Casals et al. (2010), the acronym stands for the species-rich genera within this clade: Filago L., Leontopodium R. Br. ex Cass., Antennaria Gaertn., and Gamochaeta Wedd.], likely underwent a hybridization plus WGD event (Smissen et al., 2011). Further, the Gnaphalieae experienced a WGD event about 10 Ma ago (Huang et al., 2016; Zhang et al., 2020b). These WGD events are expected to add phylogenetic discordance.

Sampling was carried out over a period of 12 years (2006–2018) by one of the authors. Leaves were dried on silica for DNA extraction. Herbarium specimens are deposited in PRC, and duplicates stored in QCA, QCNE, and AAU. Further, we got additional material from AAU, B, MA, and BONN. We sampled 15 out of the 19 accepted Loricaria species and, in addition, four new morphological groups representing potentially new species (nine individuals; Supplementary Table 1).

For Hyb-Seq, we included between 1 and 13 samples per species to test for their monophyly, with a stronger focus on sect. Thyopsis, as it includes more widespread species than the other sections. The outgroup was complemented with sequences from Mandel et al. (2019), who used the same probe set for target enrichment (see below). In total, 13 species from seven genera of the Gnaphalieae were sampled as outgroup taxa. Overall, we sampled 63 individuals for this study, including 13 outgroup samples.

DNA extraction, genomic library preparation, and bait hybridization followed Gizaw et al. (2021). We used the Compositae1061 probe set (Mandel et al., 2014), implemented in the myBaits Expert Compositae1061 target capture kit (Arbor Biosciences, Ann Arbor, MI, United States). Enriched libraries were mixed with unenriched libraries in the ratio 2: 1 (run2), 1.5: 1 (run3), and 1: 1 (run5), respectively. Samples were sequenced on different sequencers, either an Illumina (San Diego, CA, United States) NextSeq at the Genomics Core Facility of CEITEC (Brno, Czechia) or an Illumina NovaSeq at IAB (Olomouc, Czechia); in all cases 150 base pairs (bp) paired-end reads were obtained. Raw reads are available under NCBI SRA BioProject PRJNA777419.

We developed a data analysis workflow that implements data filtering, paralog detection and utilization for phylogenetic reconstruction, and investigation of ILS and hybridization to robustly infer the phylogeny of a young radiating group (Figure 1). Most scripts used herein are part of HybPhyloMaker (Fér and Schmickl, 2018; scripts are available at https://github.com/tomas-fer/HybPhyloMaker), which we indicate as ‘‘HPM’’ followed by the number of the respective script. Customizing the reference, paralog detection, and orthologous alignment building was performed using ParalogWizard¹ (scripts are available at https://github.com/rufimov/ParalogWizard), to which we refer to as “PW” followed by the number of the respective script. In the case of running scripts outside of these two bioinformatic pipelines we refer to the scripts directly in the respective methodological section, and if steps need to be done by the user manually, we denote this as “manual.” All steps and scripts are also summarized in the Supplementary Table 3.

As missing data have a substantial impact on correct species tree estimation, especially under a high degree of ILS (Nute et al., 2018), we chose a relatively stringent tradeoff between the number of assembled loci per sample and the number of samples removed from the dataset due to low quality. We excluded samples with a recovery of less than 80% across all loci from both the nuclear and plastome datasets, which resulted in the removal of 26 and 28 samples for the nuclear and plastome dataset, respectively. In total, we used three outgroup samples and 34 samples from the ingroup for the nuclear dataset, reducing the dataset from 63 to 37 samples. For the plastome dataset, we used 10 outgroup samples and 25 samples from the ingroup, totaling 35 samples. Alignments were built for each dataset separately, and after sample filtering resulted in 13 species plus our additional four new morphological groups for the nuclear alignments, and 13 species plus additionally three of the new morphological groups of Loricaria for the plastome alignments.

After trimming and de-duplicating the raw reads, between 0.9 and 14 million reads per sample were retained. The nuclear dataset 1, from which the initial ASTRAL species tree was reconstructed, consisted of 973–1150 loci per sample. To avoid data loss by removing paralogous loci, we differentiated between orthologous and paralogous copies and built separate alignments from those, which increased the number of loci for phylogenetic reconstruction by about 25%. We did not find a higher degree of paralogy in particular species or clades (Supplementary Table 2).

Without employing our stringent sample filtering by removing samples with less than 80% of the exons recovered, the support for clades within Loricaria was low [<0.95 local posterior probability (LPP); data not shown]. After removing those samples (HPM 5), the number of loci used during alignment building increased, as these low-quality samples did not longer have a dominating effect on the removal of loci that had less than 80% of all samples present (Supplementary Figure 1A). Further, gene tree discordance decreased (Supplementary Figure 1A), and the support for the major clades (those discussed later) increased, from only four nodes within Loricaria having a support greater than 0.95 LPP to eight such nodes (Figure 2A; not counting supported nodes within species-complexes).

The number of mapped plastid reads ranged from 5027 to 199,001, with an average proportion of missing data of 3.3% (min-max: 0.06–17.6%). The length of the concatenated plastome dataset after removing samples with more than 20% missing data was 261,701 bp.

The monophyly of Loricaria was strongly supported based on the nuclear data, both in the ASTRAL species tree (1 LPP) and the concatenated tree (100% bootstrap support [BS]; Figure 2; dataset 1). The earliest diverging taxon to the outgroup, L. graveolens, the only representative of the sect. Graveoleum, was the only species with support for monophyly (ASTRAL: 1 LPP; concatenation: 100% BS). Based on the ASTRAL and concatenated tree the two main clades retrieved grouped species of high morphological similarity (for dataset 1: Figure 2; for dataset 2: Supplementary Figure 6), and we hereafter refer to these main clades as the Terminal- (1 LPP; 100% BS for dataset 1 and 2) and Axilliary-clade (0.88 LPP/0.96 LPP for dataset 1/dataset 2 and 96%/76% BS, respectively). The phylogenetic placement of most samples followed the sectional classification of the genus, with the exception of L. unduaviensis Cuatrec. and one new morphological group, L. “leptothamna-like,” which were placed with samples from the sect. Terminalia. Loricaria unduaviensis and L. leptothamna were placed in the sect. Thyopsis by Cuatrecasas (1954), irrespective of the position of their capitulum, which is terminal, not axilliary. Similarly, Loricaria colombiana has axilliary capitulas, but is part of the Terminal-clade and assigned to sect. Terminalia. The samples identified as L. “ollgaardii-like,” have high morphological similarity to L. ollgaardii except for the position of their capitulas, which is axilliary in our samples but terminal in the original species description. Accordingly, our samples were found in the Axilliary-clade. In addition to the L. “ollgaardii-like” samples, the Axilliary-clade consisted of two subclades, the scolopendra-complanata-clade (1 LPP; 100% BS for dataset 1 and 2) and the thuyoides-complex, which comprised samples from L. pauciflora Cuatrec., L. azuayensis Cuatrec., and L. thuyoides as well as three of the four new morphological groups. Within the thuyoides-complex, samples of the species L. thuyoides were found to be non-monophyletic: samples of L. pauciflora and L. azuayensis were nested among L. thuyoides samples. Further, four L. thuyoides samples formed a sister clade relationship to the remaining samples of the Axilliary-clade (1 LPP; 100% BS for dataset 1 and 2), denoted as thuyoides-clade2. The Terminal-clade comprised the ilinissae-clade (1 LPP; 100% BS for dataset 1 and 2), including L. ilinissae (Benth.) Cuatrec., L. puracensis Cuatrec., L. antisanensis Cuatrec., and L. colombiana Cuatrec., and the leptothamna-clade (dataset 1: 0.99 LPP; 100% BS; dataset 2: 0.97 LPP; 100% BS), composed of L. “leptothamna-like,” L. lycopodinae Cuatrec., L. ferruginea (Ruiz & Pav.) Wedd., and L. unduaviensis.

We recovered the same clades for the concatenated phylogeny as we did based on the ASTRAL species tree, with only a few samples showing different placements within these clades (Figure 2). Analysis of gene tree discordance using PhyParts showed that each node is highly discordant, with only a small proportion of gene trees supporting the species tree (Supplementary Figure 2). When investigating the frequency of the different topologies retrieved among the gene trees using DiscoVista (method 5), for the node leading to the genus, one topology was most frequent (>50%, Figure 3A). After removing samples that showed gene flow (dataset 2), the split separating the early diverging species L. graveolens from the remaining clades was found in more than 50% of the topologies (Figure 3B). The remaining nodes showed almost equal frequencies for all three topologies (highest frequency of maximally 40%; Figure 3A). The percentage of alternative quartets according to ASTRAL gave similar results, mostly having similar frequencies for all three quartets (Supplementary Figure 3). Only the node indicating genus monophyly and the species L. graveolens were well supported with more than 50% of the quartets showing the same topology. Clades supported by more than 40% of the quartets were the ilinissae-clade, the scolopendra-complatana-clade, and the thuyoides-clade2, as well as the mrca of the Axilliary- and Terminal-clade (Supplementary Figure 3).

Cytonuclear discordance was strong in our dataset (Supplementary Figure 4). The plastome tree was reconstructed using a slightly different selection of samples than for the nuclear trees, due to our missing data filter approach. While most clades recovered in the nuclear dataset were present in the plastome tree, there were major differences. First, the genus was non-monophyletic: Belloa schultzii (Wedd.) Cabrera formed a clade with the earliest diverging species L. graveolens (100% BS). Second, the only sample available from the scolopendra-complanata-clade (14281) for the plastome dataset plus two samples of L. thuyoides (11523, 12246) formed a clade together with the samples from the ilinissae-clade (100% BS), while other members of the thuyoides-clade1 did not (Supplementary Figure 4). Additionally, the ilinissae-clade is part of the Axilliary-clade in the plastome tree (100% BS).

According to the distance-based network, we identified a different set of misplaced samples compared to the ones identified in the plastome tree as well as those that showed signs of ILS: Samples from two of the new morphological groups, i.e., L. “ollgaardii-like” samples (12066, LR_011) and the L. “scolopendra-like 1” sample (12002; Supplementary Figure 5) grouped with different samples in the network than in the phylogenetic reconstructions.

Calling SNPs from our Hyb-Seq data resulted in an initial set of 16,000 SNPs before filtering according to our minimum threshold for SNP presence across samples and accounting for linkage. Approximately 50% of the SNPs were removed due to their presence in less than 80% of the samples, and restricting the SNP dataset to unlinked SNPs resulted in a further reduction to final 1515 SNPs. Based on the Dsuite analyses, we identified seven samples showing introgression greater than 5% (Table 2). All these samples belonged to the thuyoides-clade1: three out of the four new morphological groups (L. “ollgaardii-like”: 12066, LR_011; L. “scolopendra-like 1”: 12002; L. “scolopendra-like 2”: 11160), two samples of L. thuyoides (12444, 11573), and L. azuayensis (11122).

While the samples identified as discordant varied between the different methods, they all belonged to the Axilliary-clade (Table 2). The Dsuite analyses indicated that most samples showed introgression with members of the ilinissae-clade. Removing those samples (dataset 2) reduced gene tree discordance for the mrca of the Terminal- and Axilliary-clade (Figure 3B) and increased support for the phylogenetic backbone (Supplementary Figure 6).

Only after removing those samples showing introgression with samples from different clades (dataset 2) were the PhyloNet analyses able to finish. Thus, we used PhyloNet to test if all potential hybridization events were detected using Dsuite. Allowing for a maximum of three hybridization events resulted in the best model according to the AIC in PhyloNet (Supplementary Table 4). The structure of the five best networks within this analysis was mostly consistent, indicating that within Loricaria there were between one to three reticulation events, the γ-parameter ranging from 0.9 to 0.49, resulting in a hybridigenous origin of either the thuyoides-clade1 or -clade2 and/or the complanata-scolopendra-clade with an ancestor of the Terminal-clade (Figure 4). In four out of the five best networks, six out of 11 introgression events showed γ ≥ 0.3 within at least one event per network (Figure 4); the remaining events showed lower levels of gene flow.

Target enrichment protocols employing probe sets that may be either customized (often genus- or tribe-specific; designed to match exons across a relatively small number of species) or universal (family- or order-specific; designed to match exons across larger ranges of the tree of life) are widely used to generate hundreds of nuclear loci for samples from evolutionary lineages ranging from deep to shallow phylogenetic scales (e.g., Carlsen et al., 2018; Johnson et al., 2019; Bagley et al., 2020). Although custom and universal probes often show a similar relative performance (Larridon et al., 2020; Shah et al., 2021; Ufimov et al., 2021), universal probes are frequently preferred for non-model organisms, for which the genomic resources, which are a prerequisite for the probe design, are often not available and too costly to generate.

In this study, we used the universal probe set Compositae1061 (Mandel et al., 2014), but customized the reference for “read fishing” before de novo assembly, which likely increases locus recovery and the length of recovered exons, as was shown in McLay et al. (2021) and Ufimov et al. (2021). (Universal) target enrichment kits target conserved exons that are present across a wide range of taxa. Hence, the phylogenetic signal is reduced, especially compared to introns (Folk et al., 2015; Bagley et al., 2020; Gardner et al., 2021). Nevertheless, the utility of the Compositae1061 probe set to resolve young phylogenies was shown before (Gizaw et al., 2021). By using ParalogWizard to detect paralogous loci and utilize them for phylogenetic reconstruction, we restricted our data analysis to exons, as ParalogWizard was developed for exons only. In contrast, HybPiper, the standard data analysis pipeline for Hyb-Seq data, permits to use flanking introns or supercontigs (exons plus flanking introns), potentially providing more phylogenetic informative characters (Jones et al., 2019; Ogutcen et al., 2021; Ufimov et al., 2021). However, an earlier study showed that for Antennaria, a close relative of Loricaria, very few supercontigs remained after alignment trimming and the remaining ones had low degrees of informative characters, a pattern generally prominent for Gnaphalieae (Jones et al., 2019).

While nuclear Hyb-Seq data are often able to resolve phylogenetic relationships with high support, gene tree discordance is usually high in these datasets (Jones et al., 2019; Smith et al., 2020), and the species tree topology can even be represented by only a minority of gene trees (so-called anomaly zone; Degnan and Rosenberg, 2009; Liu and Edwards, 2009; Roch and Steel, 2015). The degree of gene tree discordance tends to be lower for custom probe sets (Bagley et al., 2020; Siniscalchi et al., 2021), especially if the custom probes target longer loci compared to the universal probes (Ufimov et al., 2021). In the case of Loricaria, we found a high degree of gene tree discordance also in comparison to other young Asteraceae groups. However, we were able to show that using stringent filters for missing data and an elaborate analysis workflow can reduce gene tree discordance, and at it least partly explains its underlying biological processes.

Repeated rounds of WGD are common for the angiosperms (Wendel, 2015), also for the tribe Gnaphalieae and family Asteraceae (Smissen et al., 2011; Barker et al., 2016; Huang et al., 2016; Zhang et al., 2020b). We, thus, accounted for paralogy during alignment building to remove this source of gene tree discordance before addressing the effect of ILS and hybridization on discordance. Even though the Compositae1061 probe set was designed to comprise exclusively single-copy loci, a certain proportion of the loci were flagged as paralogous using HybPiper in recent works (Jones et al., 2019; Siniscalchi et al., 2019). We inferred that about 25% of the loci included paralogous copies in our dataset. The high number of paralogous loci can likely be attributed to the multiple WGD events within the family Asteraceae and the tribe Gnaphalieae in particular. As previous genome size estimates indicated that the genus lacks neopolyploids (Kolář et al., 2016), the duplicated loci are likely the result of WGDs in the tribe.

It should be noted that we do not address certain methodological artifacts as sources of gene tree discordance, namely the effect of (a) collapsing weakly supported nodes in gene trees (HPM 10), (b) removing gappy regions in the alignment (HPM 4a3), (c) selecting the most parsimony informative alignments (run manual) or (d) excluding loci showing signs of recombination (PhiPack; Bruen et al., 2006; run manual). Initial analyses showed that the effect of these artifacts on gene tree discordance was weak for our datasets (data not shown). While gene tree discordance tends to decrease with increasing data completeness (Siniscalchi et al., 2019), our stringent removal of low-quality samples and alignment filtering likely reduced the possible effect of methodological artifacts on gene tree discordance to a minimum (Supplementary Figure 1). Although several highly interesting samples were removed in the process of filtering for missing data, it was the only possibility to reduce gene tree estimation errors and, thus, be able to focus on the biological processes as sources of gene tree discordance (Supplementary Figure 1).

The genus Loricaria is monophyletic according to our nuclear phylogeny (Figure 2), while Belloa is nested within the earliest diverging lineage of Loricaria according to the plastome phylogeny (Supplementary Figure 4). Employing a less stringent sample filter allowed to include Belloa in the nuclear dataset, and in this case Belloa did not belong to the genus Loricaria in both the ASTRAL species tree and concatenated tree (data not shown). The placement of Belloa within Loricaria based on the plastome dataset provides evidence for ILS or hybridization between lineages across genera in the Gnaphalieae. This highlights that a well-sampled outgroup and good knowledge about sister lineages through broader sampling of the tribe is required to gain a better understanding of the evolution of the genus and the tribe.

The phylogeny is split into three major clades, reflecting mainly the three different sections within the genus. Strong support for species monophyly was only found for L. graveolens, whereas all other species were not supported to be monophyletic. Whether this is due to gene flow between species, overdescription of taxonomic species or limited phylogenetic signal for young high-altitude Andean groups for the universal Compositae1061 loci needs to be evaluated. The widespread species L. thuyoides is highly polyphyletic, as indicated in Kolář et al. (2016). We detected introgression between seven members of the thuyoides-clade1 and members of the ilinissae-clade using Dsuite, and another three samples showed strong cytonuclear conflict (Supplementary Figure 4). These results were confirmed by the PhyloNet analyses. The high degree of gene tree discordance (Supplementary Figure 1), while accounting for paralogy due to WGD during alignment building, suggests that ILS and hybridization played an important role in the evolution of Loricaria.

Although the initial species tree after sample and alignment filtering (dataset 1) supported most of the major clades in the phylogeny with support > 0.95 LPP, three different gene tree topologies were almost equally likely for most nodes, indicating a substantial degree of ILS (Figure 2A). It needs to be noted that a comparison between the ASTRAL species tree and the concatenated tree resulted in the same set of supported clades, although within these clades some samples showed different placements (Figure 2). This low discordance indicates a moderate degree of ILS, as RAxML is more sensitive to ILS than ASTRAL (Mirarab et al., 2016). Collapsing nodes with low support (i.e., 10 or 30% BS) before species tree calculation did not decrease the degree of discordance (data not shown), indicating that the discordance is not due to low phylogenetic signal, but rather due to ILS.

In rapidly radiating lineages, the degree of ILS is expected to be high (Whitfield and Lockhart, 2007), due to insufficient time for alleles to coalesce. Earlier molecular dating efforts estimated a crown age of 4 Ma for Loricaria (Nie et al., 2016), which resulted in an approximate net diversification rate of 0.74 species per Million years [ln(N)/t; N: number of species, t: crown age; Magallon and Sanderson, 2001]. This is comparable to the high diversification rates for other plant groups from tropical high elevations in South America (Madriñán et al., 2013), as well as rates that were found in the biodiversity hotspot of the Cape Floristic Region (Pirie et al., 2016).

The degree of ILS increases with large population sizes (Slatkin and Pollack, 2006). Loricaria evolved during the uplift of the Andes in the northern parts of South America, where during Pleistocene glaciations the alpine belt shifted downwards (Van der Hammen, 1985; Hooghiemstra and Van der Hammen, 2004), resulting in larger potentially suitable habitats. This might have resulted in larger effective population sizes during these intervals throughout the evolution of Loricaria. In a future study, we will investigate the demography of Loricaria species using a population genomics approach.

The sum of our results (cytonuclear discordance, the NeighborNet as well as the Dsuite and PhyloNet analyses) indicate that ILS alone cannot explain the high degree of discordance that we observed within Loricaria (Figure 2, Table 2, and Supplementary Figure 2). The removal of samples with introgression according to Dsuite (dataset 2) increased the support for the major split between the two main subclades in the genus from about 40% of all gene trees to above 50% (Figure 3 and Supplementary Figure 2). Nevertheless, standard methods to illustrate genomic discordance did not show major improvements (based on node support and PhyParts; Figure 2 and Supplementary Figures 1B, 2) between dataset 1 and dataset 2. This highlights the importance to investigate species tree hypotheses beyond support values and standard methods of measuring discordance and to thoroughly test for all potential sources of discordance, especially in highly understudied lineages.

Hybrids are unknown for the genus Loricaria based on morphological evidence (Cuatrecasas, 1954). In addition, previous genome size estimations of several Loricaria species revealed only relatively small differences (8.76–11.69 pg DNA; measurements available for L. ilinissae, L. scolopendra, L. thuyoides, and L. complanata; Kolář et al., 2016), suggesting the absence of hybridization, under the assumption that hybrids are frequently stabilized by polyploidization. Using a diverse spectrum of methods, we detected multiple hybridization events within Loricaria. This was not surprising given that the genus belongs to the tribe Gnaphalieae, for which many hybridization events have been reported (e.g., Smissen et al., 2011; Galbany-Casals et al., 2014; Barker et al., 2016; Huang et al., 2016; Vargas et al., 2017; Watson et al., 2020; Zhang et al., 2020b). Further, several radiating plant groups in the tropical high altitude areas of South America show hybridization (Lachemilla: Morales-Briones et al., 2018; Lupinus: Nevado et al., 2018; Diplostephium: Vargas et al., 2017; Espeletiinae: Cortés et al., 2018). The dynamic nature of this ecosystem with multiple range expansions and contractions during the Pleistocene (Flantua et al., 2019) may have facilitated the contact between geographically isolated species that probably did not yet exhibit strong barriers to gene flow. Using Dsuite, we identified a total of seven samples showing introgression with samples of other clades in the species tree. After removing those samples, we were able to detect one to three clades within the genus that are of hybridigenous origin according to PhyloNet (Figure 4). Unfortunately, PhyloNet and related methods (SNAQ; Solís-Lemus et al., 2017) are difficult to use for large datasets with hundreds of samples and a high number of hybridization events. We ran PhyloNet using the MPL algorithm, after unsuccessfully attempting to utilize the “divide and conquer” method (Zhu et al., 2019) using a maximum likelihood implementation, which did not finish (<14 days and 14 CPUs of computation) for a subset of quartets, even though the method is intended for large datasets. The subset of quartets that did not finish included to a large extent those that were subject to hybridization events based on MPL (observations during trials). While PhyloNet accounts for gene flow and ILS as a source of discordance, the type of gene flow can be the result of hybridization, introgression or horizontal gene transfer. These processes are biologically very similar and cannot be differentiated methodologically by this method. However, different degrees of gene flow between species or lineages, depicted by the γ-parameter, may hind-cast the different processes (Solís-Lemus et al., 2017), in our case supporting hybridization between early lineages within Loricaria (six out of 11 introgression events showed γ ≥ 0.3, Figure 4). Nevertheless, the different degree of gene flow detected, ranging from γ = 0.09–0.49, suggests pure hybrids as well as hybridization with extensive backcrossing.

The polyphyletic nature of L. thuyoides, with some samples showing cytonuclear discordance, and others exhibiting introgression with members of the ilinissae-clade according to Dsuite analyses, indicates that L. thuyoides was subject to chloroplast capture and hybridization early in its history. While L. thuyoides is described to be morphologically variable (Cuatrecasas, 1954), we could not find any morphological characters that enable samples showing introgression to be distinguished from pure samples. Chloroplast capture is the result of two species hybridizing with extensive backcrossing to one of the parents (Rieseberg and Soltis, 1991). Due to the extensive backcrossing, the nuclear signal of the hybridization event is swamped out, but the novel plastid from the hybridization event remains. The two L. thuyoides samples showing indication of chloroplast capture group in the plastome phylogeny with ilinissae-clade samples from close geographic proximity, a pattern common for chloroplast capture (Acosta and Premoli, 2010; Liu et al., 2020). The clades we identified to have a potentially hybridigenous origin, the scolopendra-complanata-clade and both thuyoides-clades (Figure 4), overlap geographically with members of the ilinissae-clade, the potential hybridization partner.

The samples of the Axilliary-clade, which showed signatures of introgression according to the Dsuite, are predominantly found in southern Ecuador, close to the Huancabamba Depression in northern Peru, which exhibits a partial interruption of the Andes by low-elevation river systems. Some works suggested that this area poses a barrier to gene flow for high altitude species (Cosacov et al., 2009; Richter et al., 2009), a pattern that cannot be confirmed by our study. Members of the two main clades, the Terminal- and Axilliary-clade, are found on both sides of this depression. Further, two of the new morphological groups were found in close proximity to this area, L. “vermicularis” and L. “ollgaardii-like,” and the latter was found to be of hybridigenous origin. Due to the shifts of the alpine belt during the Pleistocene glaciation and deglaciation cycles (Flantua et al., 2019), populations of the different clades likely came into contact in the Huancabamba Depression, which might have facilitated hybridization (L. “ollgaardii-like,” L. “scolopendra-like” 1 and 2) as well as speciation (L. “vermicularis”). As such, for species of Loricaria that are exclusively found in the páramo ecosystem, the Huancabamba area does not seem to be a barrier to gene flow. The Huancabamba Depression has also been identified as a center of diversity for montane species (Mutke et al., 2014; Quintana et al., 2017). We cannot confirm, however, if this secondary contact was also facilitated in the southern part of the Huancabamba Depression, as we lack good sampling from the northern parts of Peru.

Despite evidence for hybridization within Loricaria, the exact parents and clades subject to the hybridization events could not be determined. Due to our simplification in the PhyloNet analysis, by mapping samples to supported clades, we could not differentiate if only some of the species in the clade or the clades as a whole were subject to the hybridization events. To elucidate which species are of hybrid origin and which parental species gave rise to these hybrids, further sampling in the area is needed as well as population-level analyses.

Hundreds of loci and thorough testing of potential causes of discordance provided a better understanding of the evolution of the genus. And yet, while nowadays detecting hybrids using genomic data is easier than during Sanger sequencing times, the lack of knowledge about lineages and missing taxonomic expertise in young radiations complicate our understanding of their evolution.