Evolving view on phylogenetic networks

Shifting into the era of network thinking

Science is driven forward by paradigm shifts [1]: moments when scientists recognize that anomalies are not just margin points, but aspects from reality ignored by our current paradigm. The acknowledgment of anomalies is a difficult task because it involves thinking outside the current scientific model, but it is only through such a change in perspective that knowledge can expand and evolve.

The tree paradigm fails to explain certain patterns (anomalies under the current paradigm) just like any scientific paradigm. Gene flow is such an anomaly. It represents a real evolutionary process that violates the key assumption of tree thinking: vertical inheritance of genetic material. Historically, gene flow in the form of horizontal gene transfer or hybridization has often been set aside to simplify the reconstruction of the evolutionary histories of species. Hybrid species or genes with strong conflicting histories tend to be removed from analysis as noise that does not fit the overall tree-like pattern. Gene flow, however, is not rare. One could argue it is even the norm in some taxonomic groups. Gene flow appears at all levels in the tree of life as shown by a large body of evidence [2, 3, 4, 5, 6, 7, 8, 9, 10, 11 e.g. in recent years], among which we can highlight pervasive hybridization during rapid plant radiations (e.g., Tillandsia [12]) and extensive introgression across both deep and shallow phylogenetic timescales in animal systems such as Heliconius butterflies [13].

The discovery of such reticulate evolutionary histories is only possible through the development of hybrid detection and network inference methods, as traditional tree-based methods are unable to capture these processes. Recently, gene flow has begun to appear as an integral part of phylogenomic analyses, effectively shifting the paradigm from tree thinking to network thinking.

The era of network thinking brings several challenges that will be described in this perspective. Even the concept of a phylogenetic network is not straightforward, as different classes of networks reflect distinct biological interpretations. Furthermore, estimating phylogenetic networks is computationally more demanding than inferring trees and necessitates larger quantities of high-quality data. In particular, statistical identifiability and computational scalability remain major challenges due to the complexity of models that must simultaneously account for multiple biological processes like incomplete lineage sorting (ILS) and gene flow. Furthermore, optimization within the expansive network space remains operationally difficult.

Despite these challenges, the continued advancement of methodological approaches is increasingly allowing evolutionary biologists to infer reticulate events across diverse lineages. The more we untangle the truth about evolution and species diversity, the more complicated it appears, but it is only because we are seeing it closely. The reality is more complex than we originally envisioned. So, as we zoom in, pieces will certainly not fit in the way we initially imagined. As we expand our understanding and inference capabilities of phylogenetic networks, we are able to include more reticulate organisms in the tree – or rather network – of life, effectively providing a broader and more complete picture of the evolutionary history of biodiversity on Earth.

A web of definitions for networks

Mathematically, trees are structurally simple. Trees are directed acyclic graphs composed of external nodes (or leaves) which represent extant taxa and internal nodes that represent speciation events of ancestral populations [14]. Edges connecting nodes represent ancestry, and along the edges, genetic mutations are expected to accumulate. Thus, a tree is a mathematical abstraction used to represent the evolutionary process (sequence of speciation events) by which an ancestral population gives rise to a collection of extant taxa.

At first glance, networks appear to be a straightforward extension of trees by the simple inclusion of a third type of node: hybrid nodes with two parental incoming edges and one outgoing child edge (Figure 1 right). Hybrid nodes represent reticulate processes such as hybrid speciation, introgression or horizontal gene transfer. However, one immediate reality complicates our understanding of networks: there is in fact no universal definition of what internal nodes (or edges) represent in what is denoted (quite loosely) a phylogenetic network. Indeed, scientists in different domains have used very different definitions of phylogenetic networks.

One such example is haplotype networks [15, 16, 17] and their median [18, 19] or median-joining [20] variants which represent the genealogical relationships among distinct haplotypes (unique combinations of linked alleles) observed within a population or species, typically inferred from non-recombining genetic loci. Nodes represent unique haplotypes, edges represent mutations between haplotypes, and internal nodes represent inferred ancestral haplotypes. These graphs can accommodate reticulations in the form of cycles, and thus, they extend tree models to account for recombination. Despite not representing evolutionary relationships, haplotype networks continue to be used to represent mutational relationships, especially within recently diverged species or populations [21, 22, 23] such as SARS-CoV-2 [24, 25, 26, 27].

Another such example is implicit or split networks [28, 29, 30], also called abstract or data-display networks. A split is any bipartition of the taxa into two non-empty, disjoint and exhaustive taxon sets, and a split network represents every split as an array of parallel edges (Figure 1 left). This non-tree-like representation is well-suited for visualizing reticulate evolution, but it cannot distinguish whether observed discordance is due to reticulation or to other processes such as ILS. Indeed, internal nodes in split networks are devoid of biological meaning, and thus, they are better suited as data-display objects rather than as a representation of true evolutionary histories. Extensive software have been implemented to reconstruct split networks such as SplitsTree [31, 32], NeighborNet [33, 34], QNet [35] and Quartet-Net [36], and much theory has been developed in the reconstructability of split networks (i.e., conditions under which a network can be uniquely recovered) from different types of data [37, 38, 39, 40, 41, 42 e.g.].

In contrast to split networks, the nodes in explicit networks represent a specific biological process: tree nodes (internal nodes with one parental edge and two children edges) represent speciation events and hybrid nodes (internal nodes with two parental edges and one child edge) represents reticulation events (Figure 1 right). One could argue that the term phylogenetic networks should be restricted to the case of explicit networks as those are the ones that have an evolutionary interpretation. For the remaining of this work, we will refer to explicit phylogenetic networks whenever we use the term phylogenetic networks.

This perspective centers on explicit phylogenetic networks. Yet, even when restricting to explicit networks, the definition is not free of entanglements. Unlike species trees [43] where internal nodes have one clear biological interpretation (speciation), hybrid nodes in networks do not have a single biological meaning. Indeed, when inferring a phylogenetic network, hybrid nodes present a non-tree-like reticulate process, but no existing method is able to distinguish if the reticulate signal comes from hybridization or introgression or other gene flow process. In fact, depending on the method used to infer the phylogenetic network, we cannot even know if the signal is coming from other biological sources like duplications or losses, or even data errors or systematic noise. Furthermore, unlike trees where there is not a “wrong placement” of internal (speciation) nodes, hybrid nodes cannot be freely placed anywhere in the phylogeny due to risks of internal incongruences such as ેણ flow of genetic material contradicting temporal constraints imposed by outgroups or the root. To further complicate matters, phylogenetic networks are classified into different groups by the number of hybrid nodes and their placement in the phylogeny (Figure 2). These groups differ widely in their ability to be estimated from data and their biological interpretability, as we will see in the next sections.

Exploring the untamed space of phylogenetic networks

Just like a jungle on an alien planet, the space of phylogenetic networks is filled with unknown graphical creatures for which we know very little about their statistical detectability (can we estimate them accurately with the data we have?), their identifiability (could our data have been equally generated under two different networks?), or their biological interpretability (can this network represent an actual evolutionary history or it is just a mathematical object disconnected from reality?).

Early attempts to make sense of this landscape took an almost taxonomic approach, defining rigorous classes of phylogenetic networks. In a foundational textbook, Huson et al. [30] cataloged a variety of network types such as tree-based, tree-child, galled, and level-$k$ networks, each imposing different structural constraints (Figure 2). This figure intends to simply emphasize the richness and fragmentation of network space rather than to serve as a classification guide, and thus, to learn more about network classes, please refer to Kong et al. [44] who not only present a comprehensive list of network classes, but also show how some of these classes correspond to particular evolutionary mechanisms. This growing taxonomy, however, only partially tames the jungle: many network classes lack a clear biological interpretation, and for most, fundamental questions of statistical identifiability remain unanswered as we will see in future sections.

Most methods for inferring networks have specialized in level-1 networks. Technically, the two incoming parental edges on a hybrid node create what is denoted as a hybridization cycle (Figure 1 right), and in a level-1 network, the hybridization cycles do not share edges. Level-1 networks have become the model organism of network space as the class on which most theoretical results on identifiability and reconstructability have been established (see sections below). Even methods that are capable of inferring more complex level-$k$ ($k > 1$) networks are routinely evaluated only on level-1 networks. So, ultimately, most of what we know about phylogenetic network inference is from the class of level-1 networks, and very few studies have tested how methods perform outside of this class [46].

The network of life is unlikely to be level-1. As a field, we are used to making simplifications when representing biological processes with mathematical graphs. For example, the tree assumption of instantaneous isolation after speciation (hence represented by one lineage splitting into two at an internal node) is a well-recognized biological idealization, for which alternative models exist but are not routinely implemented. Restricting phylogenetic networks to level-1, however, is largely driven by current theoretical and computational constraints, rather than by a comparable biological rationale. Thus, expanding our knowledge about identifiability and computational scalability to other classes of networks will allow us to include broader, more realistic evolutionary relationships of highly reticulate organisms in the network of life.

Inferring the Network of Life

Networks as trees with added hybrid arcs

The evolution of network methods has followed an intuitive path. Given our field’s success in inferring trees, one of the initial strategies to infer a network was to first infer a tree, and then test subsets of taxa for signal of reticulations. These so-called “hybrid detection methods” such as the D-statistics (ABBA-BABA) [47] and similar D tests [48, 49, 50, 51] became universally embraced to identify taxa of hybrid origin. Unfortunately, recent work [52, 53, 54, 55, 56] has shown that D tests often exhibit elevated false discovery rates, erroneously classifying non-hybrid taxa as having hybrid origin, especially in cases of ancestral gene flow, ghost lineages [57, 58] or rate variability [54].

Not all hybrid detection methods are inspired by the D-statistics. HyDe [59] and PhyloDiamond [60] use phylogenetic invariants under the multispecies coalescent (MSC) model, on site pattern frequencies and concordance factors respectively, to detect signal of reticulation on subsets of taxa, and MSCQuartets [61, 62] tests deviations from expectations under the MSC on a 4-taxon species tree. Likelihood ratio tests [63] have been introduced as more powerful alternatives to test for reticulation signal in site pattern frequencies. Some of these tests, however, still show high false negative rates when involving extinct or unsampled “ghost” lineages and under multiple hybridization events, implying that complex hybrid histories can attenuate rather than confuse the signal of gene flow [55].

Furthermore, all hybrid detection methods are affected by biases in tree inference under gene flow. Specifically, gene flow can generate anomalous gene trees that are overrepresented in the data, leading to inconsistent species tree estimation [64, 65, 66, 67]. Gene flow can also lead to overestimated population sizes and underestimated species divergence times, and can result in lower posterior support for true clades, potentially yielding high confidence in an incorrect topology [68]. Even low levels of gene flow can drastically distort species tree estimates, especially after rapid radiations with consensus and distance-based methods outperforming MSC likelihood approaches that ignore gene flow [69].

Machine learning has recently been introduced to distinguish histories of speciations from histories of hybridization with high accuracy [70]. Once trained, machine-learning models can be rapidly applied to new datasets, yet they require labeled training datasets which, in phylogenetics, must be generated by simulation. Despite efforts to emulate realistic evolutionary scenarios across broad parameter ranges, such simulations inevitably introduce bias [71].

Networks as amalgamation of (sub)trees

Early methods to reconstruct phylogenetic networks started as heuristic or distance-based with limited incorporation of explicit evolutionary models [72]. Tools like T-REX [73, 74, 75, 76] were among the first to explicitly model reticulate evolutionary patterns using a distance-based approach. However, because these algorithms are heuristic in nature, reticulation edges lack biological interpretation and may confound multiple sources of conflict, such as ILS or measurement error.

Developments on the mathematical properties of splits, clusters and subgraph patterns within explicit networks gave rise to the formalization of theory and algorithms to build explicit networks from integral subparts. For example, we can build explicit networks from triplets [3-taxon trees; [77, 78, 79, 80, 81]], quartets [4-taxon trees; [82, 83]], binets [2-taxon subnetworks; [84, 85]], trinets [3-taxon subnetworks; [85, 86, 87, 88, 89, 90]], quarnets [4-taxon subnetworks; [91, 92, 93, 94]], clusters [95, 96, 97], and even ancestral profiles [98].

Substantial theoretical literature has characterized level-$k$ networks showing, for example, their uniqueness from triplets despite NP-hard reconstruction [99], their decomposition into level-$k$ generators [100], and the possibility of polynomial-time reconstruction from clusters [101]. Related work characterizes which level-1 networks can be encoded from trees, clusters, and triplets [102], and shows how classes such as tree-child networks can be counted and enumerated for a fixed number of leaves and reticulations [103]. These results directly informed algorithms and software including Dendroscope [104], methods reconstructing level-1 networks from triplets [Lev1athan; 79], trinets [TriLoNet; 88], and quarnets [Squirrel; 105], distance-based extensions of neighbor-joining [106, 107, 108], and gene tree approaches such as RF-Net [109].

All these algorithms are built upon solid mathematical and computational theory, yet they do not rely on explicit probabilistic models of data generation. By amalgamating subgraphs into networks, these methods often fail to account for different sources of gene tree discordance or gene tree estimation error. Their heuristic setting makes it challenging to extend to different biological processes beyond reticulations, and thus, while these algorithms do build explicit networks that have an evolutionary interpretation, they are sometimes better suited as data-display objects just as split networks are.

Networks as parameters in probabilistic models

A shared dream of statistics and biology is to define a unified probabilistic model that accounts for all possible biological processes involved in the generation of the observed data. Under the assumption of no recombination within loci, extending our tree probabilistic framework to networks only requires replacing the stochastic process that generates gene tree variability within a species tree with the corresponding process on a species network, given that the mutational process on gene trees that gives rise to sequences remains unchanged.

The multispecies coalescent (MSC) model on trees [43, 110, 111] tracks the genealogical history of gene sequences going backwards in time and allows the estimation of species tree topology, divergence times, and population sizes. [112] extended the MSC model to include hybridization events by incorporating a new numerical parameter: the inheritance probability $\gamma$ that models the probability of a single gene from the hybrid species to follow a specific path towards the most recent common ancestor. That is, in this model, if $B$ is the hybrid species with two sister species $A$ and $C$ (Figure 1 right), when a single gene from the hybrid species is selected for analysis, its most recent common ancestor occurs with taxon $C$ with probability $\gamma$ and with taxon $A$ with probability $1-\gamma$. Upon selecting a species tree for a given gene, the coalescent process is applied to that tree, yielding probabilities for the various gene tree topologies.

This extended MSC model had one limitation: if there was more than one hybrid species (or multiple individuals sampled from the one hybrid species), all lineages below the hybrid node were forced to follow the same species tree. [113] extended this model to allow any hybrid lineage to follow any path, independently of other lineages, and this was further parameterized by [114] to include a $\rho$ parameter to determine the extent of correlation among lineages with $\rho= 0$ corresponding to independent inheritance as in [113] and $\rho= 1$ corresponding to common inheritance as in [112].

This new probabilistic model that combined ILS and hybridization was denoted indistinguishably the multispecies network coalescent (MSNC) or the network multispecies coalescent (NMSC) and allows us to define a likelihood function to implement a maximum likelihood procedure that will traverse network space in search of the optimum network to fit the data. This methodology, originally implemented in PhyloNet [115, 116], allows multilocus sequences or estimated gene trees as input data, later extended to Bayesian inference in BEAST2 [117] and in PhyloNet [118].

The shortage of sufficiently long genes to infer accurate gene trees in some taxonomic groups motivated the development of methods based on biallelic markers [119]. Algorithms such as PhyloNet [120, 121], SnappNet [122], and the pioneering TreeMix [123] infer networks from SNP data, although TreeMix is rarely considered within the realm of phylogenetic network methods given its focus on “admixture graphs" (see Discussion).

Crucially, by defining a likelihood function (whose exact form depends on the data type and modeling assumptions), probabilistic approaches enable the inference of numerical parameters such as branch lengths and inheritance probabilities alongside network topology. Focusing on numerical parameter inference, BPP [124] implements two models of gene flow under the MSC: 1) the MSC with introgression (MSC-I or MSci) [125, 126], which is the same as the MSNC, assumes episodic gene flow at discrete points in time, and 2) the MSC with migration [126, 127] (MSC-M, also known as the isolation-with-migration or IM model [128]) which assumes continuous gene flow at a fixed rate after divergence. While BPP does not currently infer the network topology, it can test for gene flow and estimate numerical parameters such as migration or introgression rates, as well as species divergence times and population sizes.

The rise of likelihood-based approaches required not only the development of probabilistic models to define likelihoods (and priors), but also of computational tools to explore the vast space of networks efficiently. Several studies have formalized network rearrangement moves that guarantee connectedness of network space. For example, [129] generalized nearest-neighbor interchange (NNI) operations to level-1 networks, and [130] extended NNI and subtree pruning and regrafting (SPR) moves to include “horizontal” moves connecting networks of the same complexity and “vertical” moves connecting networks of different complexities. Other work has shown that connectedness can also be achieved using only tail moves [131] or head moves [132], and that tree-child networks are connected under a generalized SPR operation for a fixed number of hybridizations [133]. Because the root is generally unidentifiable under likelihood-based models (see Statistical Identifiability section), many methods infer semi-directed networks, in which hybrid edges and their descendants are directed while all other edges remain undirected. In this context, [134] proposed a rearrangement move for semi-directed networks (cut edge transfer, CET) and showed that the space of semi-directed level-1 networks with a fixed number of reticulations is connected.

Although probabilistic models now routinely accommodate ILS and gene flow, extensions beyond these processes are rare. The dream of integrating several (all) biological processes into one probabilistic framework is complicated for statistical and computational reasons, as we will see in the next section. Network models described so far account for ILS and gene flow; none can tackle allopolyploidization. Polyploids can arise through hybridization between distinct species, resulting in distinct subgenomes that are allopolyploids, but this biological process is ignored by most existing network models that assume species are diploid. One notable exception is AlloNet [135, 136], which infers allopolyploid phylogenetic networks in the presence of ILS, but is restricted to a single hybridization event between 2 diploid species resulting in a genomic allotetraploid. Further efforts are required to continue to develop probabilistic models that encompass a broader array of biological processes, while rigorously addressing the inevitable statistical and computational challenges, some of which we describe next.

The challenges of network inference

Built upon biologically-informed probability models, likelihood-based network inference methods have established as accurate methodologies to infer reticulate evolutionary histories. These methods, however, have been confronted with two primary challenges. First, the prohibitive computational cost of full likelihood evaluation on phylogenetic networks rendered these methods impractical for large empirical datasets (lack of scalability). Second, as the complexity of the probability models increases, it becomes more likely that different parameter settings produce similar patterns in the data such as gene trees or site frequencies (lack of identifiability). For example, for the case of network topologies, as we increase the network complexity (number and placement of hybridizations), distinct network topologies can generate identical data (e.g. gene tree or site pattern distributions), rendering them statistically indistinguishable.

Computational scalability

Dividing and conquering. Likelihood computation can be accelerated by decomposition into smaller, more tractable components. SNaQ [137] pioneered a composite (pseudo) likelihood framework by decomposing the full likelihood into quarnets, a strategy also implemented in PhyloNet [138] and in PhyNEST [139]. However, composite likelihood strategies alone are not likely to achieve the scalability of tree methods that can routinely handle thousands of taxa [140 e.g.]. Even when likelihood optimization is indeed accelerated, the optimization must still traverse network space which grows super-exponentially with the number of taxa.

FastNet [polynomial complexity; [141]] and later InPhyNet [linear complexity; [142]] provide a decomposition strategy to split the full problem (network inference on the whole set of taxa) into smaller subproblems (network inferences on a subset of taxa), and then merge solutions of the subproblems into a solution for the full problem. InPhyNet is the first network reconstruction method to achieve the same levels of number of taxa as trees by building the phylogenetic network for the 10,000 plant genomes [140].

Other divide-and-conquer strategies combine probability models and heuristic approaches. For example, NANUQ [143] infers level-1 networks by analyzing quartets displayed on gene trees. Its updated version NANUQ+ [144] first infers the tree of blobs [145]: a tree where hybridization biconnected components (blobs) are collapsed into multifurcating nodes, and then resolves multifurcations as hybridization events. TINNiK [146], used within NANUQ+, infers trees of blobs from gene quartet distributions under the MSNC. Similarly, Squirrel [105] infers candidate trees of blobs and then resolves them into level-1 networks. Some other methods simplify computation by ignoring biological processes, such as NetRax [147] that performs maximum likelihood network inference without modeling ILS.

Bypassing optimization. Algebraic geometry provides a theoretical framework that enables faster and more efficient phylogenetic network inference through the use of algebraic invariants: polynomial equations that model parameters have to satisfy for internal consistency. For example, concordance factors [148] represent the proportion of genes that support each of the three possible quartet representations for a 4-taxon subset: $ab|cd$, $ac|bd$, $ad|bc$, and these three concordance factors need to sum up to one. Thus, the equation $CF_{ab|cd}+CF_{ac|bd}+CF_{ad|bc}=1$ needs to be satisfied, and thus, the polynomial equation $f(x,y,z)=x+y+z-1$ for $x=CF_{ab|cd}$, $y=CF_{ac|bd}$, $z=CF_{ad|bc}$ is an algebraic (phylogenetic) invariant.

Phylogenetic invariants have been used to infer 4-taxon networks under a variety of models of evolution such as Cavender–Farris–Neyman [149] and Jukes–Cantor [150] for level-1 networks, with subsequent extensions to level-2 networks under Jukes–Cantor [151, 152]. [153] extends to any 4-taxon semi-directed network under the Jukes–Cantor and the Kimura 2-parameter models, and [154] under the random walk 4-state Markov model. These invariants account for the mutation process of sequences under networks, but not for gene tree variability under the coalescent model. In contrast, Wu and Solís-Lemus [60] use phylogenetic invariants defined on concordance factors under the MSNC, which enables inference of level-1 networks for $4<n<9$ taxa and serves as a hybrid detection tool for larger $n$.

Because these methods evaluate polynomial equations instead of exhaustively traversing network space, they achieve greater computational efficiency than traditional likelihood-based approaches. However, defining invariants becomes increasingly challenging as network size grows, limiting their applicability primarily to smaller networks where likelihood-based methods remain tractable. Future research could focus on integrating these strategies, for example, by using invariant-based approaches to preselect candidate networks for subsequent likelihood optimization and parameter estimation.

Statistical identifiability

Intuitively, two networks that generate the same data patterns would not be distinguishable from that data alone. Therefore, it is essential to determine the model assumptions for which a specific data type uniquely identify a single network topology and its numerical parameters. That is, we need to assess statistical identifiability of models on networks. These results depend on both the form of the data (gene trees, quartet concordance factors, or sequence data), and the underlying modeling assumptions such as substitution models and whether gene tree variation is explicitly modeled under the MSNC. Statistical identifiability transcends purely theoretical considerations since understanding which local network features are most readily identifiable (and with which type of data) can help guide the design of sampling strategies that maximize detectability of reticulations.

Building intuition: early studies on identifiability. [137] began the investigations on identifiability (or more precisely, detectability) of hybridization cycles on level-1 networks under the MSNC from quartet concordance factors (CFs). They found that the root is not identifiable which motivated further development of theory and methods on semi-directed networks. This early study focused on detectability (that is, whether hybridization events could be statistically detected from quartet CFs) rather than on uniquely identifying the underlying network, and showed that our ability to detect hybridization cycles (from quartet concordance factors under MSNC) depends on the number of nodes in that cycle, with 2-cycles (cycles with 2 nodes) never detectable; 5-cycles (or bigger) always detectable, and 3- and 4-cycles detectable under certain taxon sampling conditions.

At the same time, [155] were examining network identifiability under a no-ILS assumption, which allows the full network likelihood to be decomposed into the likelihoods of its displayed trees (each obtained by removing one hybrid edge per hybridization). They showed that canonical networks, those in which no node has multiple incoming edges and a single outgoing edge, are generally identifiable, whereas certain configurations with stacked hybridizations (hybrids giving rise to hybrids) remain indistinguishable. [156] later showed that networks deemed not identifiable under no ILS by [155] are actually distinguishable under the MSNC (which accounts for ILS) when more than two lineages are sampled below the hybrid. A central intuition emerging from these initial investigations is that taxon sampling incorporating multiple alleles per species, and gene tree variability due to ILS can help network distinguishability. Moreover, even when full networks remained indistinguishable, certain features such as their canonical, semi-directed, 2-cycle free or triangle (3-cycle) free forms may still be identifiable.

Level-1 networks are identifiable across models and data types. [157] proved full (generic) identifiability of level-1 semi-directed networks under the MSNC from gene quartet distributions with cycles of size at least 4, later extended to 3-cycles by [158]. Generic identifiability means for all values of numerical parameters except some in a set of measure zero (almost all values of numerical parameters). These works allow us to conclude that level-1 semi-directed networks (topology and numerical parameters) are generically identifiable under the MSNC from gene quartet distributions except in the presence of 2-cycles and, in certain cases, hybrid nodes and adjacent edges within 3-cycles.

Similar results have been proven under other models and data types. Level-1 networks with cycles of size at least 4 were proven identifiable under the Jukes–Cantor substitution model for sequence data [159], later extended to the Kimura 2- and 3-parameter models [160], with subsequent work demonstrating that even the number of reticulations in a level-1 network is identifiable under Jukes–Cantor [152]. Level-1 ultrametric networks, under MSNC combined with a mixture of General Time-Reversible evolution models for sequence data, are identifiable from logdet distances when cycles have size at least 4 [161]. [162] showed that level-1 networks with positive internal edge lengths are identifiable from average pairwise distances, except for local ambiguities in small (4-) cycles where hybrid edge orientation and adjacent branch lengths and inheritance parameters are not identifiable.

Some networks are identifiable beyond level-1. Stepping outside the class of level-1 networks requires strong assumptions on the network structure to guarantee identifiability. Specifically, it seems that “hybridizations on top of hybridizations" can be quite problematic. Scenarios where hybrid nodes have hybrid nodes as children or there is a path from the root to the tips that traverses multiple hybridizations cannot be easily detected from different data types. Once we restrict these scenarios, we can have identifiability results outside the level-1 class. Namely, [163] proved identifiability of more general networks from quartet concordance factors under MSNC. Ruling out small (2- or 3-) cycles, a network is identifiable if it is galled and all hybridization blobs are tree-child, meaning that hybridization events are isolated rather than layered on top of one another. Intuitively, galled networks restrict each hybridization to a single undirected (tree) cycle formed by two disjoint paths, while the tree-child condition ensures that every internal node has at least one descendant that follows a tree-like lineage, thereby ruling out stacked hybridizations in which a hybrid node gives rise to another hybrid node. Their study also corroborated what we knew from early identifiability studies: multiple samples per taxon aid in network identifiability compared to single samples.

[152] showed that level-2, triangle (3-cycle) free and strongly tree-child networks are identifiable from sequence data under the Jukes–Cantor model, assumptions that again exclude stacked hybridizations. The recurring exclusion of stacked hybridizations in identifiability results may indicate that the canonical network classes (which does not allow stacked hybridizations) may capture the identifiable structure even in complex reticulation scenarios.

General network features emerge as identifiable across models and data types. Beyond topology and numerical parameters, several network features are also identifiable under different models and data types. For example, [145] proved that the unrooted tree of blobs can be determined from the distribution of gene quartet topologies under the MSNC and is also identifiable under Jukes–Cantor and the Kimura-2-Parameter substitution models from sequence data [152]. A refinement of this tree, the distance split tree, is identifiable from pairwise distances [162]. As mentioned before, the identifiability of the tree of blobs has opened the door to multiple inference network methodologies that first infer a tree of blobs to later identify the hybridization blobs [105, 144]. Other features, such as the circular order of subnetwork attachments to a blob, are also identifiable [152, 164].

Proven cases of unidentifiability. Because the root is not identifiable across models and data types [137, 157, 162], many methods these days infer semi-directed networks. Yet, as evolutionary histories inherently require a root for interpretation, these networks must be rooted post hoc, a process that can introduce interpretation challenges (see Discussion). Small cycles of length 2 or 3 tend to also be unidentifiable across models and data types [137, 152, 160, 162]. Some 3-cycles in particular can create anomalous networks [165] where an unrooted tree displayed by the network may have lower probability under the MSNC than a tree not displayed. Such anomalies complicate inference from gene tree data. Biologically, simply discarding small cycles from networks is not ideal. Highly prevalent gene flow events in the network of life are expected to appear in small cycles. For example, introgression or migration between sister extant taxa is by definition a 3-cycle. These reticulation events are common among reticulate organisms, yet most of our current theory rules them out as unidentifiable. Numerical parameters around hybridizations can also be hard to infer, especially around small cycles [158, 162].

In terms of network classes, it has been proven that level-2 networks are not identifiable from average distances, not even generically [162]. Looking ahead and borrowing ideas outside of probabilistic modeling, [93] recently showed that level-3 networks are not encoded by quarnets, and [152] demonstrated that any two binary semi-directed phylogenetic networks are only distinguishable if they do not display exactly the same quartets. Together, these results suggest that level-3 networks could be challenging to identify under some settings. Stacked hybridizations are also hard or impossible to detect [152]. This might motivate the focus on canonical (or somehow reduced) networks [155, 166, 167] which would prompt us to look beyond binary networks.

As we move farther from level-1-ness, more restrictions will need to be imposed on the networks to prove statistical identifiability, and this investigation will shape the type of network we can infer from different data types. Evolutionary paths from root to extant species that involve multiple hybridization events are hard to infer with our current models and data types. This suggests exploring new types of data, different sampling strategies, or simplified network structures that focus on the most identifiable features. Just as the field quickly moved into semi-directed networks as the estimable object of interest, as we move outside level-1-ness, we could also think of canonical or other forms of reduced networks as our identifiable parameter of interest.

Discussion

This perspective centers on explicit phylogenetic networks inferred under an assumed probabilistic model for data generation. This choice was subjective. As a statistician, I favor probabilistic approaches. As an evolutionary biologist, I favor biological interpretability. However, the field of phylogenetic networks is broad and rich, and any of the subtopics I tangentially discussed could be the focus of its own perspective piece.

This perspective has also focused on the methods to build phylogenetic networks rather than on the extensive body of research on reticulations across the diverse array of organisms comprising the network of life (a small subset of those studies mentioned in the Introduction). This perspective is also not meant to guide a novice into the field of phylogenetic networks. Other reviews on phylogenetic networks (better suited as introductions to the topic) have been published over the years [72, 139, 168, 169, 170, 171, 172, 173, 174, 175].

I conclude this Discussion by highlighting additional challenges, each also deserving a review of their own.

Interpretations. Interpreting trees, even unrooted trees, is relatively straightforward. A tree represents a sequence of speciation events and is rarely subject to ambiguous visualization. It is not the same for networks. The manner in which networks are plotted (many times a subjective decision by the plotting software) can impact the way in which they are interpreted. Among the plotting software for phylogenetic networks, we can highlight Dendroscope [104], PhyloPlots [176], tanggle [177], and ggret [178]. Each software has a default way to plot phylogenetic networks.

Figure 3 shows three ways to plot the same network topology, but with different biological interpretations. On the left, the network seems to represent introgression or horizontal gene transfer, that is, a reticulate process in which there is gene flow from one population (ancestor of $E$) to another population (ancestor of $A$). The center figure seems to indicate that taxon $A$ arises through a process of hybrid speciation. Unlike the left network, this center network does not have estimates for the inheritance probabilities to aid in the interpretation, and the manner in which it is plotted could suggest equal contribution on both hybrid edges. The figure on the right represents gene flow from an unsampled or extinct taxon $X$ into taxon $A$. It should be emphasized that, from a mathematical standpoint, the three network topologies are identical. Any perceived differences arise solely from artifacts of their graphical representation.

So, how can we interpret inferred networks? First, numerical parameters are important components of the inferred network. In trees, the branch lengths do not alter the pattern of speciation events, so a tree topology (even without branch lengths) can still be biologically interpreted. For networks, branch lengths and inheritance probabilities are necessary to provide biological meaning. Branch lengths in calendar time could provide information about the time points of the origin and the target of gene flow, and whether these two time points coincide. Estimated branch lengths of 0 for the hybrid edges are represented by perpendicular lines as in Figure 3 center. Given that gene flow should happen between contemporary populations, a hybrid edge with positive length could provide evidence of hybridization with an unsampled or extinct taxon (as in Figure 3 right). The inheritance probability can also provide information on whether the gene flow event is introgression (small $\gamma$) or hybrid speciation ($\gamma$ close to 0.5). We should emphasize, however, that the interpretation of inheritance probabilities depends on whether the model represents episodic or continuous gene flow. In a continuous gene flow model, it is possible to have values of $\gamma$ near 0.5 that do not invoke hybrid speciation [179], so caution is needed to interpret $\gamma$ as the sole source of evidence for hybrid speciation versus introgressive hybridization.

When interpreting networks, it is also important to understand the constraints under which the inference method operates. For example, a method that can only infer level-1 networks will have a limit on the number of hybridization events it can place on a phylogeny. Selecting $h = 5$ hybridization events to be estimated, but obtaining a network with only $h = 3$ does not mean that there is biological or even statistical evidence for only 3 hybridization events instead of 5. The method might have not been able to place more hybridization events given the level-1 constraint. Some methods can place constraints on hybrid edges as well, sometimes due to identifiability issues, forcing the lengths to be equal to 0.0, and thus, they cannot be used as evidence for contemporary gene flow or hybrid speciation as they are an artifact of the optimization constraints.

Robustness. Identifiability results clarify what is theoretically detectable from data, but they do not address how network inference methods behave when model assumptions are violated. Despite extensive work on the robustness of tree methods under gene flow [64, 65, 66, 67, 180 e.g.], analogous studies for phylogenetic networks remain scarce. Existing work shows that data errors and model misspecification can spuriously increase inferred network complexity: sequencing and gene-tree errors may introduce non-biological reticulations [181], while ignoring gene-tree uncertainty or assuming homogeneous substitution rates can have similar effects [182]. Practical evaluations of software such as PhyloNet have begun to translate these findings into guidance for users [46]. More studies are needed to assess the robustness of our network inference methodologies as these methods tend to have several topological and numerical assumptions that could be violated in empirical data.

Distances and consensus. Computational operations that are straightforward on trees such as calculating distances or dissimilarities, and finding consensus graphs become substantially more complex when extended to networks. Key challenges arise because distance functions that work well within one class of networks often fail to behave as proper distances in other classes, and commonly used units to build distances or consensus such as splits or clades do not map to unique edges in networks. More work is needed to find better distance functions that properly capture the notion of similarity among phylogenetic networks [30, 183, 184, 185, 186, 187, 188, 189 e.g.], as well as algorithms for computing consensus networks from collections of networks (not just trees) [190] to support Bayesian inference and other uncertainty quantification frameworks.

Traits. Phylogenetic comparative methods such as phylogenetic regression, phylogenetic ANOVA, ancestral trait reconstruction, and Pagel’s $\lambda$ test of phylogenetic signal study the distribution of quantitative traits among related species. Karimi et al. [4], Jhwueng and O’Meara [191], Bastide et al. [192], Teo et al. [193, 194] present comparative methods where the phylogenetic correlation structure of a trait is governed by a network. Existing models, implemented in PhyloTraits [195], already allow for network-specific tests such as transgressive evolution. However, substantial work remains to match the breadth of trait models available for trees.

A network by any other name. As noted earlier, the term ’phylogenetic network’ can have different meanings on different domains. It turns out that the converse is also true: there exist terminologies that refer to objects that mathematically are phylogenetic networks under different names. Admixture graphs [123, 196, 197, 198, 199, 200, 201, 202], highly popular in human population genetics, are graphical models that represent divergence and admixture events in a population’s history. Although mathematically equivalent, admixture graphs and phylogenetic networks have evolved in parallel fields, with largely non-overlapping theoretical foundations and methodologies.

A similar isolation exists between phylogenetic networks and ancestral recombination [203, 204, 205, 206, 207, 208] or reassortment [209, 210, 211, 212] graphs, both ambiguously referred to as ARGs. While these ARGs capture gene-level reticulate processes such as recombination or reassortment, their graphical structure is mathematically equivalent to that of a phylogenetic network.

While the theoretical and methodological paths of these network objects have diverged, perhaps it is time to recombine them by focusing on their shared characteristics, and to exploit the substantial opportunities for cross-fertilization that remain unexplored. Network theory aims to understand the mathematical objects that we denote networks, to study their statistical, mathematical and computational properties, and to develop efficient, accurate and scalable ways to infer them from available data. Although different domains rely on different empirical data and probabilistic models tied to specific biological processes, the inherent graphs are networks that share commonalities that warrant a unified and comprehensive investigation.

Just like species on Earth, our methodologies evolve. The more deeply we study reality, the more anomalies we find, and the bubble bursts: our current methodologies cannot handle some realities. A tree structure is not enough to represent the evolutionary history of all species on Earth. As we improve our network models and better ways to extract information about reticulate processes from empirical data, we can help transition to a network thinking era, a new paradigm that is not meant to overrule the idea of a tree of life, but rather to expand it into a broader picture of evolution.

Funding

This work was partially supported by the National Science Foundation (NSF DEB-2144367 to CSL).

Acknowledgements

I want to thank Chris Simon and Stacey Smith for inviting me to contribute to the 75th anniversary edition of the journal of Systematic Biology as well as Bob Thomson, Paul Bastide and Laura Kubatko for the thorough comments that greatly improved this work. It was an honor to acknowledge the contributions of past and present researchers who dared venture into the vast network space, laying the methodological foundations that make it possible to reveal the hidden patterns of reticulate evolution in data. I am also very grateful for the continuous support throughout my career of the Society of Systematic Biologists.

Disclosure statement

No conflicts of interest to disclose.

Data Availability Statement

No data associated with this perspective article.