Characterization of Genome-wide Phylogenetic Conflict Uncovers Evolutionary Modes of Carnivorous Fungi

20 Mass extinction has often paved the way for rapid evolutionary radiation, resulting in the 21 emergence of diverse taxa within specific lineages. While the emergence and diversification of 22 carnivorous nematode-trapping fungi (NTF) in Ascomycota has been linked to the 23 Permian-Triassic (PT) extinction, the processes underlying NTF radiation remain unclear. Here, 24 we conducted phylogenomic analyses using 23 genomes spanning three NTF lineages, each 25 employing distinct nematode traps — mechanical traps ( Drechslerella spp.), three-dimensional 26 (3-D) adhesive traps (Arthrobotrys spp.), and two-dimensional (2-D) adhesive traps ( Dactylellina 27 spp.), and one non-NTF species as the outgroup. This analysis revealed how diverse mechanisms 28 contributed to the tempo of NTF evolution and rapid radiation. The genome-scale species tree of 29 NTFs suggested that Drechslerella emerged earlier than Arthrobotrys and Dactylellina. Extensive 30 genome-wide phylogenetic discordance was observed, mainly due to incomplete lineage sorting 31 (ILS) between lineages (~81.3%). Modes of non-vertical evolution (i.e., introgression and 32 horizontal gene transfer) also contributed to phylogenetic discordance. The ILS genes that are 33 associated with hyphal growth and trap morphogenesis (e.g., those associated with the cell 34 membrane system and cellular polarity division) exhibited signs of positive selection. 35 36 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 3

37 Mass extinctions result in vacated ecological niches that can be occupied by novel species and 38 drive subsequent radiation events (Sepkoski 1998; Jablonski 2001). Mass extinction and 39 concomitant radiations have been documented in multiple lineages, including angiosperms 40 (Silvestro et al. 2015), planktic foraminifera (Lowery and Fraass 2019), snakes (Grundler and 41 Rabosky 2021), modern birds (Jarvis et al. 2014), and mushrooms (Varga et al., 2019). 42 Comparative genomics has facilitated systematic identification of candidate genetic changes 43 underlying speciation and adaptive radiation (Marques et al. 2019). 44 Carnivorous nematode-trapping fungi (NTF) emerged after the Permian-Triassic (PT) extinction 45 and radiated into multiple lineages that form distinct trapping devices to capture free-living 46 nematodes (Yang et al. 2012). The emergence of NTF from saprophytic fungal species is thought 47 to be driven by nematode proliferation that occurred in the wake of the PT extinction, an event 48 that resulted in a carbon-rich and nitrogen-poor environment (Barron, 1977; Gray, 1983; Liu et al. 49 2014; Fan et al. 2024). The ability to capture and consume nematodes allows NTF to obtain extra 50 nitrogen, likely conferring a competitive advantage over saprophytic fungi (Yang et al. 2012) and 51 driving diversification. NTF have radiated into three clades, each with unique genus designations 52 and trapping systems—Arthrobotrys spp. employ three-dimensional (3D) adhesive traps 53 (networks); Dactylellina spp. utilize two-dimensional (2D) adhesive traps (knob, column, 54 non-constricting ring); and Drechslerella spp. form mechanical traps (constricting ring) (Jiang et 55 al., 2017). 56 Phylogenomics has greatly advanced our understanding of the Tree of Life, mechanisms of gene 57 and genome evolution, and relationship between genomic and phenotypic divergence during 58 speciation (Jin et al. 2021; Steenwyk et al. 2023). Phylogenomics has also revealed how individual 59 genes have undergone evolutionary histories that are distinct from the phylogenetic history of 60 species carrying these genes (Steenwyk et al. 2019; Salichos and Rokas 2013). Theoretical and 61 empirical studies have shown that this discordance or incongruence can be caused by analytical 62 factors, such as errors in taxon sampling and gene tree estimation, and biological mechanisms, 63 such as incomplete lineage sorting (ILS), horizontal gene transfer (HGT), and introgressive 64 hybridization (Steenwyk et al. 2020a; Lopes et al. 2021; Shen et al. 2021; Feng et al. 2022; 65 Steenwyk et al. 2023). Other factors also influenced species diversification during radiation. For 66 example, adaptive evolution punctuated by positive selection occurs more frequently in radiating 67 lineages than in slowly diversifying ones (Nevado et al. 2019). While ILS, HGT, introgression, 68 and positive selection have been documented in several eukaryotic lineages, such as cichlids 69 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 4 (Brawand et al. 2014), wild tomatoes (Pease et al. 2016), honeybees (Fouks et al. 2021), big cats 70 (Figueiró et al. 2017), and Populus species (Wang et al. 2020), their impact on fungal radiation 71 events remains poorly understood. 72 Here, we characterized the genome-wide patterns and drivers of phylogenetic discordance in the 73 three NTF lineages. Patterns of genome-wide phylogenetic discordance showed that ILS between 74 lineages caused most of the observed discordances. In contrast, introgression and HGT contributed 75 less to the incongruence between the species tree and gene trees. Positive selection of ILS genes 76 associated with growth and trap morphogenesis were also observed. Similar to previous studies of 77 other lineages, our phylogenomic analyses revealed how diverse evolutionary mechanisms 78 contributed to the tempo of NTF evolution and rapid radiation. 79

80 Extensive phylogenomic discordance among NTF 81 To investigate the evolutionary history of NTF in Ascomycota, we analyzed 23 NTF genomes 82 (Supplementary Table S1). Our taxon sampling covered three major lineages that underwent 83 radiation and evolved distinct mechanisms of nematode trapping, including 3-D adhesive 84 networks (Arthrobotrys spp.), 2-D adhesive traps (Dactylellina spp.), and mechanical traps 85 (Drechslerella spp.), and Dactylella cylindrospora, a non-NTF closely related to NTF as the 86 outgroup. 87 Single-copy orthologous genes (2,944 in total; Supplementary Table S2) present in all species 88 were combined to construct maximum likelihood species tree using two alignment and trimming 89 strategies (Clustal-Omega + ClipKIT and MAFFT + Gblocks) (Castresana 2000; Katoh and 90 Standley 2013; Sievers and Higgins 2018; Steenwyk et al. 2020b). The species tree topologies 91 were consistent under both strategies (Figure 1, Figure S1), suggesting our analyses are robust to 92 some analytical sources of error associated with alignment and trimming strategy (Steenwyk et al. 93 2023). The genome-scale phylogeny was consistent with our previously published mu ltiple-gene 94 phylogeny (Yang et al. 2007) and strongly supported the placement of Arthrobotrys and 95 Dactylellina as sister genera (Figure 1). The species tree supported two notable evolutionary 96 events: the divergence of those forming mechanical traps (Drechslerella) and the lineage that 97 produces adhesive traps and the subsequent divergence of 2-D ( Dactylellina) and 3-D 98 (Arthrobotrys) adhesive traps. 99 Maximum likelihood trees of single-copy orthologous genes were also constructed using 100 Clustal-Omega + ClipKIT and MAFFT + Gblocks approaches. The resulting trees were largely 101 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 5 consistent between each method, suggesting analytical errors associated with software choice are 102 minimal (Steenwyk et al. 2023). Nonetheless, discordance between the single-gene trees and the 103 species tree was abundant (Figure 2, Supplementary Table S2). Densitree plots depicted numerous 104 topological conflicts among the gene trees (Figure 2a), and MDS analysis based on 105 Robinson-Foulds (RF) distances revealed differences between the gene trees and the species tree 106 (Figure 2b). Concordance analyses based on IQ-TREE showed that there was a high rate of 107 conflict between gene trees and the species tree at the divergence points between mechanical traps 108 (Drechslerella) and adhesive traps, as well as between 2-D (Dactylellina) and 3-D (Arthrobotrys) 109 adhesive traps (gene-concordance factors (gCF) < 60%, Figure 1). Furthermore, there were two 110 and one nodes with high conflict (gCF < 60%, Figure 1) within Arthrobotrys and Dactylellina, 111 respectively. 112 ILS is largely responsible for phylogenetic discordance 113 To further rule out analytical sources of error, we identified single gene trees that were consistent 114 between the two alignment and trimming strategies — Clustal-Omega + ClipKIT and 115 MAFFT+Gblocks (Castresana 2000; Katoh and Standley 2013; Sievers and Higgins 2018; 116 Steenwyk et al. 2020b). Among the 2,944 single-copy orthologous genes, 64 orthologous genes 117 yielded inconsistent gene trees with the two strategies (see Supplementary Table S2); inconsistent 118 genes, which are likely subject to analytical errors, were removed from subsequent analyses. 119 Among the remaining 2,880 gene trees, 496 exhibited average bootstrap support below 80% 120 (Supplementary Table S1), suggesting that errors in phylogenetic inference may have affected 121 these trees. Among the remaining 2,392 trees with high bootstrap support, 978 trees (40.9%; a 122 group designated as Tree1) supported the species tree, whereas 1,414 trees (59.1%) were 123 inconsistent with the species tree. 124 The Multispecies Coalescent (MSC) model was employed to investigate whether the observed 125 topologies of gene trees across sets of four lineages could be attributed to ILS. By employing 126 hypothesis testing on 1,414 gene trees, we assessed the concordance between observed gene tree 127 distributions and those predicted under the MSC model. The results showed that at the 0.0001 128 significance level, 81.3% of the four-lineage scenarios supported the hypothesis that ILS shaped 129 the topology. Whereas 18.7% of the scenarios rejected the hypothesis (Figure 3a), suggesting other 130 evolutionary modes, like introgression and HGT, may influence the history of these loci. 131 Examination of genome-wide D-statistics analysis (also known as the ABBA-BABA test; Hibbins 132 and Hahn, 2022; Bjornson et al. 2023), which test for introgression, revealed insignificant 133 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 6 amounts of introgression among the three NTF lineages (Figure 3b; D = -0.0075, Z = -0.114). 134 However, phylogenetic network analysis revealed two gene introgression events in Arthrobotrys 135 lineage and one in Dactylellina lineage (Figure 3c); notably, these nodes high degrees of conflict 136 among gene trees and the species tree in Figure 1. 137 Among the 1,414 genes displaying topological structures that conflict with the species tree, 36 138 appeared to have been acquired via HGT. These HGT genes predominantly originated from 139 bacteria, with Pseudomonadota being the main donor phylum. Additionally, some HGT events 140 from fungi, particularly from the sister phylum Basidiomycota, were also observed. 141 The remaining 1,378 trees were categorized into three groups (Figure 4): 7.0% (97) placed 142 Arthrobotrys and Drechslerella as sister groups (Tree2), 18.0% (245) clustered Drechslerella and 143 Dactylellina (Tree3), and 75.0% (1,036) did not align with their corresponding generic clades 144 (designated as Unclassified). 145 The branch lengths at the divergence nodes of the gene trees likely affected by ILS were longer 146 than those in the species tree, a significant signal supporting ILS (Song et al. 2023). We compared 147 the divergent branch lengths between the ancestral node (node1) and the next divergence node 148 (node2), which represents the duration of nematode-trapping device divergence in the three 149 different types of gene trees (Figure 4a). The mean divergent branch lengths for Tree2 and Tree3 150 (0.5746 and 0.5895, respectively) were significantly shorter than that for the species tree (0.4399, 151 p < 0.0001), supporting the contribution of ILS to the divergence of the three NTF lineages. 152 The phylogenetic conflicts in those categorized as "Unclassified" (1,036 trees) were likely caused 153 by ILS. The MSC analysis indicated that 84.44% of the conflicts in the four lineages could not 154 reject the hypothesis that they arose from ILS (Figure S2). ILS events involve random fixation of 155 ancestral sequences, leading to a plethora of topologies spanning the NTF lineages. A substantial 156 number of gene trees exhibiting inconsistency with the lineages may be due to the stochastic 157 nature of ILS. At the same time, the lack of correspondence between these gene trees and the 158 branches of the lineage suggest that these are more ancient ILS events, and the gene sorting may 159 have occurred before the lineage divergence. Compared to Tree1, Tree2, and Tree3, the 160 Unclassified type trees have significantly shorter cumulative branch lengths (Figure 4c), 161 suggesting lower evolutionary rates (Steenwyk et al. 2021) (Supplementary Table S3). 162 ILS genes under positive selection are broadly associated with growth and trap morphogenesis. 163 Natural selection during rapid evolutionary radiation frequently leads to accelerated gene 164 evolution and resulting phenotypic changes (Nevado et al. 2019; Hines and Rahman 2019). 165 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 7 Positive selection among ILS genes was detected using CodeML with the site model. Sixteen 166 single-copy orthologous genes exhibited signs of significant pos itive selection (Supplementary 167 Table S2) and were enriched in functions related to the cell membrane system and cellular polarity 168 division (Figure 5). For example, functions related to the cell membrane system include the 169 nuclear outer membrane (GO:0005640), plasma membrane (GO:0005886), outer membrane 170 (GO:0019867), and endoplasmic reticulum (GO:0005783). Meanwhile, functions related to 171 cellular polarity division include the cellular bud tip (GO:0005934) and neck (GO:0005935), 172 cellular bud (GO:0005933), and site of polarized growth (GO:0030427). Since the cell membrane 173 system and polarity division are cellular bases for morphological innovation, the positively 174 selected functions of these ILS genes are likely related to the morphogenesis of NTF trap 175 structures. 176 Among the conserved genes categorized as "Unclassified", 35 gene families showed significant 177 evidence of positive selection. The functions of these gene families are primarily enriched in 178 processes related to the RNA polymerase, cell nucleus, and transcription (Figure S3, 179 Supplementary Table S2). These functions are crucial for performing conserved cytological 180 processes. 181

182 We investigated the evolutionary history of carnivorous NTF in Ascomycota by analyzing the 183 genome-wide pattern of phylogenetic discordance and positive selection using the genomes of 21 184 species (23 strains) representing three NTF lineages. We generated the first genome-scale species 185 tree for these NTF. While the genome-scale species tree (Figure 1) was consistent with previously 186 published phylogenetic trees (Yang et al. 2012; Yang et al. 2007), we found extensive 187 phylogenetic discordance across the genome and the nodes of the species tree. The ILS between 188 lineages caused 81.3% of the phylogenetic discordance, while 18.7% were attributed to 189 post-speciation introgression within the lineage or HGT. The reticulate phylogenetic inference 190 indicates that introgression only led to differentiation within certain NTF genera. Although HGT 191 events caused some conflicts between gene trees and the species tree, they are not the primary 192 driver of the widespread phylogenetic inconsistencies. These results suggest that the PT extinction 193 led to the rapid stochastic fixation of ancestral polymorphisms and diverged along the lineages in 194 NTF. Subsequent positive selection accelerated the evolution of genes associated with carnivory. 195 During this process, a small number of HGT events might have contributed to genetic 196 polymorphism in carnivorous fungi. Moreover, gene flow between NTF lineages was restricted, 197 with only limited introgression occurring within each lineage. 198 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 8 The main sources of phylogenetic discordance between gene trees and species trees are ILS, 199 introgression, and HGT. Genome-wide signatures of ILS and introgression can be distinguished 200 because coalescence times for regions under ILS should be older than the speciation events, 201 whereas hybridization is the post-speciation events (Feng et al. 2022). The observation that the 202 branch lengths from the ancestral nodes to the lineage differentiation nodes in Tree2 and Tree3 are 203 longer than those in the species tree supports the hypothesis that ILS is the primary cause of the 204 observed phylogenetic discordance. ILS causes ancestral genetic polymorphisms to persist during 205 rapid speciation (Hibbins et al. 2020), and ILS events have been detected in many lineages, 206 including marsupials (Feng et al. 2022), peat moss (Meleshko et al. 2021), butterflies (Edelman et 207 al. 2019), and eared seals (Lopes et al. 2021). Our study indicates that the evolution of NTF 208 represents a new case of ILS-driven evolution. 209 We also observed signals of introgression within Arthrobotrys and Dactylellina, with the 210 occurrence of a reticulate phylogenetic relationship within each lineage. Consequently, the effect 211 of introgression was more pronounced among the closely related species within the generic 212 lineage. The role of gene introgression events in species evolution has garnered increasing 213 attention because numerous studies have highlighted their significant effects on ecological 214 adaptability and evolution in species such as primates, butterflies (Edelman et al. 2019), gray 215 snub-nosed monkey (Wu et al. 2023) and foxes (L Rocha et al. 2023) . Future studies should 216 explore the effects of introgression within each NTF lineage. 217 Some inconsistencies between the gene trees and species tree were caused by HGT. Most HGT 218 genes originated from bacteria, but some originated from fungi in the phylum Basidiomycota. 219 Although HGT events may not be the main factor driving the divergence of NTF lineages, they 220 typically introduce traits that play a crucial role in evolution (Li et al. 2022), which may also hold 221 true for carnivorous fungi. Functional characterization of such genes should be performed to 222 assess their significance in the evolution of NTF (Fan et al. 2024). 223 The most conflict-rich regions tend to be associated with the highest rates of phenotypic 224 innovation, which have been detected in six clades of vertebrates and plants (Parins-Fukuchi et al. 225 2021). The most conflict-rich nodes in this study also coincide with the differentiation nodes of 226 NTF nematode traps, which also implies that these genes undergoing ILS may be associated with 227 morphological innovation in NTF. Here, we found that some ILS genes underwent positive 228 selection, especially genes involved in cell membrane system have been shown to be involved in 229 trap morphogenesis (Bai et al., 2023; Chen et al., 2022) and inflation of the constricting ring 230 (Chen et al., 2023). The role of positive selection in the adaptive radiation of cichlids, wild 231 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 9 tomatoes, and Jaltomata has also been demonstrated, even though there are gene tree discordances 232 in their evolutionary process (Brawand et al. 2014; Pease et al. 2016; Wu et al. 2018). The use of 233 gene trees for each ILS gene instead of the species tree in our positive selection analysis helped to 234 reduce the risk of false positives. This underscores the significance of positive selection as an 235 evolutionary driver that accelerates the adaptive radiation of carnivorous fungi. 236 Gene tree discordance represents another source of substitution rate variation that can lead to false 237 inferences regarding positive selection (Mendes and Hahn 2016). Genes linked to adaptive traits 238 might not align with the species tree, causing changes in substitution rates and potentially 239 misleading conclusions about positive selection. Therefore, interpretation of positive selection and 240 adaptive radiation requires caution. Our study detected positive selection in the genes associated 241 with carnivorous traits. The use of gene trees for each ILS gene instead of the species tree in our 242 positive selection analysis helped to reduce the risk of false positives. This underscores the 243 significance of positive selection as an evolutionary driver that accelerates the adaptive radiation 244 of carnivorous fungi. 245 Additionally, many genes that did not align with their corresponding generic clades were detected 246 and are likely to have originated prior to the divergence of the three NTF lineages. ILS typically 247

in the random retention of ancestral sequences (Korstian et al. 2022; Rivas-González et al. 248 2023), and this stochastic process is responsible for the generation of gene trees that do not align 249 with the clades of the lineage. The significantly shorter cumulative branch lengths observed in 250 these gene trees (Figure 4c) suggest their ancient origin and conservation, indicating their role in 251 conserved functions related to basic life processes rather than those associated with carnivorous 252 lifestyle. Our findings highlight, for the first time, the importance of these genes. 253

254 Through phylogenomic analyses, the evolutionary history of NTF in Ascomycota, a phylum to 255 which most known carnivorous fungi belong, was investigated. Their evolution was facilitated by 256 the PT extinction, which led to rapid radiation driven by ILS, coupled with positive selection of 257 the genes associated with various carnivorous traits between generic lineages, and introgression 258 within each lineage of two genera that form adhesive traps. These analyses advanced our 259 understanding of the genetic mechanism underlying fungal adaptive radiation and evolution. 260

261 Genome mining 262 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 10 Genomes with published protein-coding gene predictions were obtained from the National Center 263 for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov/bioproject/791178). 264 Considering the frequent gene family expansion during fungal evolution, only single-copy genes 265 present in all species, with cutoffs of >70% identity and >90% coverage for their cDNAs, were 266 used in this study. In total, 2,944 gene groups were identified (Table S2, S3) using OrthoFinder v 267 2.5.6 (Emms and Kelly 2019). The nucleotide and protein sequences of these genes were matched. 268 Conserved protein domains were predicted using pfam-scan (Mistry et al. 2021). Gene Ontology 269 (GO) terms based on the functional domains were obtained using pfam2go 270 (http://geneontology.org/external2go/pfam2go). Detailed gene functions were predicted using 271 InterPro Scan (http://www.ebi.ac.uk/interpro/). 272 Phylogenetic analyses 273 To minimize the impact of phylogenetic inference errors on subsequent analyses, we employed 274 two methodologies for phylogenetic analysis, resulting in two sets of species and gene tree 275 datasets. 276 The first approach involved aligning the nucleotide sequences of all single-copy orthologous 277 genes using MAFFT v 7.520 and Gblock v 0.91b (Castresana 2000; Katoh and Standley 2013). 278 The combined sequences were used to construct species trees using IQ-TREE v 2.2.2.7 with 1,000 279 replicates (Minh et al. 2020). Individual gene trees based on nucleotide sequences were 280 constructed using IQ-TREE v 2.2.2.7 with 1,000 replicates. The species and gene trees were 281 rooted using the corresponding sequences of D. cylindrospora. 282 The second approach involved aligning the nucleotide sequences of all single-copy orthologous 283 genes using Clustal-Omega v 1.2.4 and ClipKIT v 2.2.2 (Sievers and Higgins 2018; Steenwyk et 284 al. 2020b). The combined sequences were used to construct species trees using IQ-TREE v 2.2.2.6 285 with 1,000 replicates (Minh et al. 2020). Individual gene trees based on nucleotide sequences were 286 constructed using IQ-TREE with 1,000 replicates. The species and gene trees were rooted using 287 the corresponding sequences of D. cylindrospora. 288 Tree types were identified using classify_tree.py, a tool available on GitHub 289 (https://github.com/dengweihx/classifytree). By comparing the classification results of the two 290 datasets, only gene trees consistent across both datasets were used for further analysis. 291 The incongruence coefficients of gene trees at each branch node of the species tree were 292 calculated using IQ-TREE v 2.2.2.6, and the species tree was presented using Interactive Tree Of 293 Life (iTOL) v5 (Letunic and Bork, 2021). Densitree plots of conflicting gene tree topologies were 294 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 11 drawn using DensiTree v 3.0.2. 295 (https://www.cs.auckland.ac.nz/~remco/DensiTree/download.html). Pairwise Robinson-Foulds 296 (RF) distances between gene trees were calculated using the ape package for R 4.1.3, and the RF 297 distances were then analyzed and plotted by multidimensional scaling (MDS) (Duchene et al. 298 2018; R Core Team 2023). The evolutionary rate of each gene tree was calculated using PhyKIT 299 (Steenwyk et al. 2021). 300 Incomplete lineage sorting analysis 301 ILS signals were detected by calculating the branch lengths of the differentiated nodes of the gene 302 trees using the Internal Branch Statistics feature of the PhyKIT toolkit (Steenwyk et al. 2021). 303 Differences in branch length were determined using the t-test. D values were detected by z-test 304 against whole genome backgrounds (see https://github.com/simonhmartin/ 305 tutorials/tree/master/ABBA_BABA_whole_genome for D statistics). ILS analyses based on the 306 four-taxon branch length chi-square test were performed and plotted using the MSCquartets 307 package R 4.1.3 (Rhodes et al. 2021). Reticulated phylogenetic inference based on the 308 InferNetwork_MP model was performed using PhyloNet v 3.8.2 309 (https://phylogenomics.rice.edu/html/tutorials.html). Detection of genome-wide HGT events was 310 performed using HGTector2 (https://github.com/qiyunlab/HGTector). 311 Analysis of positive selection 312 Positive selection on 2,944 single-copy orthologous genes was evaluated using CodeML and 313 PAML (Yang 2007) based on the GWideCodeML package for Python 3.10.12 314 (https://github.com/lauguma/gwidecodeml). The dn/ds values for each clade were calculated using 315 the site model. To correct for errors in substitution rate estimation due to ILS, we performed 316 branch site model calculations for the genes subjected to ILS based on their gene trees. Results 317 from the GO term enrichment analysis were presented using the clusterProfiler package for R 318 4.3.2 (R Core Team 2023). 319

320 We are deeply grateful to Dr. Antonis Rokas, Department of Biological Sciences at Vanderbilt 321 University, and Dr. Yafei Mao, Bio-X institutes, Shanghai Jiaotong University, for providing 322 insightful advice. This work was supported by a Major International Joint Research Project grant 323 from the Natural Scientific Foundation of China (Grant no. 32020103001) and the Startup Fund 324 from the Nankai University to XZL. SK acknowledges support from the USDA-NIFA and Hatch 325 Appropriation (PEN4839). JLS is a Howard Hughes Medical Institute Awardee of the Life 326 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 12 Sciences Research Foundation. 327 Competing interest 328 JLS is an advisor for ForensisGroup Inc. The authors have no relevant financial or non-financial 329 interests to disclose. 330

331 Bai N, Xie M, Liu Q, Wang W, Liu Y , Yang J. AoSte12 Is Required for Mycelial Development, 332 Conidiation, Trap Morphogenesis, and Secondary Metabolism by Regulating Hyphal Fusion in 333 Nematode-Trapping Fungus Arthrobotrys oligospora. Microbiol Spectr. 2023:11(2):e0395722. 334 https://doi.org/10.1128/SPECTRUM.03957-22 335 Barron GL .1977. The nematode-destroying fungi. In: Topics in mycobiology, No. 1. Canadian 336 Biological Publications Ltd., Guelph, ON, Canada 337 Bjornson S, Upham N, Verbruggen H, Steenwyk J. Phylogenomic Inference, Divergence-Time 338 Calibration, and Methods for Characterizing Reticulate Evolution. Preprints. 2023, 2023090905. 339 https://doi.org/10.20944/preprints202309.0905.v1 340 Brawand D, Wagner CE, Li YI, Malinsky M, Keller I, Fan S, Simakov O, Ng AY, Lim ZW, 341 Bezault E, et al. The genomic substrate for adaptive radiation in African cichlid fish. Nature. 342 2014:513(7518):375-381. https://doi.org/10.1038/NATURE13726 343 Castresana J. Selection of conserved blocks from multiple alignments for their use in phylogenetic 344 analysis. Mol Biol Evol. 2000:17(4):540-552. 345 https://doi.org/10.1093/OXFORDJOURNALS.MOLBEV .A026334 346 Chen Y, Liu J, Fan Y, Xiang M, Kang S, Wei D, Liu X. SNARE Protein DdVam7 of the 347 Nematode-Trapping Fungus Drechslerella dactyloides Regulates Vegetative Growth, Conidiation, 348 and the Predatory Process via Vacuole Assembly. Microbiol Spectr. 2022:10(6):e0187222. 349 https://doi.org/10.1128/SPECTRUM.01872-22 350 Chen Y, Liu J, Kang S Wei D, Fan Y, Xiang M, Liu X. A palisade-shaped membrane reservoir is 351 required for rapid ring cell inflation in Drechslerella dactyloides. Nat Commun. 2023:14(1):7376. 352 https://doi.org/10.1038/S41467-023-43235-W 353 Duchene DA, Bragg JG, Duchene S, Jofre GI, D’Agostino ERR, Mavengere H, Tate AD, Matute 354 DR. Analysis of Phylogenomic Tree Space Resolves Relationships Among Marsupial Families. 355 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 13 Syst Biol. 2018:67(3):400-412. https://doi.org/10.1093/SYSBIO/SYX076 356 Edelman NB, Frandsen PB, Miyagi M, Clavijo B, Davey J, Dikow RB, García-Accinelli G, Van 357 Belleghem SM, Patterson N, Neafsey DE, et al. Genomic architecture and introgression shape a 358 butterfly radiation. Science. 2019:366(6465):594-599. 359 https://doi.org/10.1126/SCIENCE.AAW2090 360 Emms DM, Kelly S. OrthoFinder: phylogenetic orthology inference for comparative genomics. 361 Genome Biol. 2019:20(1):238. https://doi.org/10.1186/S13059-019-1832-Y 362 Fan, YN, Du MH, Zhang WW, Deng W, Yang EC, Wamg SX, Zhang L, Kang SC, Steenwyk JL, 363 An ZQ, et al. Genomic changes underlying fungal carnivorism emerged after the Permian-Triassic 364 Extinction. Mol Biol Evol. 2024, forthcoming. 365 Feng S, Bai M, Rivas-González I, Li C, Liu S, Tong Y, Yang H, Chen G, Xie D, Sears KE, et al. 366 Incomplete lineage sorting and phenotypic evolution in marsupials. Cell. 367 2022:185(10):1646-1660.e18. https://doi.org/10.1016/J.CELL.2022.03.034 368 Figueiró H V., Li G, Trindade FJ, Assis J, Pais F, Fernandes G, Santos SHD, Hughes GM, 369 Komissarov A, Antunes A, et al. Genome-wide signatures of complex introgression and adaptive 370 evolution in the big cats. Sci Adv. 2017:3(7):e1700299. https://doi.org/10.1126/SCIADV.1700299 371 Fouks B, Brand P, Nguyen HN, Herman J, Camara F, Ence D, Hagen DE, Hoff KJ, Nachweide S, 372 Romoth L, et al. The genomic basis of evolutionary differentiation among honey bees. Genome 373 Res. 2021:31(7):1203-1215. https://doi.org/10.1101/GR.272310.120/-/DC1 374 Gray NF. Ecology of nematophagous fungi: distribution and habitat. Annals of Applied Biology. 375 1983:102(3):501-509. https://doi.org/10.1111/J.1744-7348.1983.TB02721.X 376 Grundler MC, Rabosky DL. Rapid increase in snake dietary diversity and complexity following 377 the end-Cretaceous mass extinction. PLoS Biol. 2021:19(10):e3001414. 378 https://doi.org/10.1371/JOURNAL.PBIO.3001414 379 Hibbins MS, Gibson MJS, Hahn MW. Determining the probability of hemiplasy in the presence of 380 incomplete lineage sorting and introgression. elife. 2020:9:e63753. 381 https://doi.org/10.7554/ELIFE.63753 382 Hibbins MS, Hahn MW. Phylogenomic approaches to detecting and characterizing introgression. 383 Genetics. 2022:220(2):iyab173. https://doi.org/10.1093/genetics/iyab173 384 Hines HM, Rahman SR. Evolutionary genetics in insect phenotypic radiations: the value qof a 385 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 14 comparative genomic approach. Curr Opin Insect Sci. 2019:36:90-95. 386 https://doi.org/10.1016/J.COIS.2019.08.013 387 Jablonski D. Lessons from the past: evolutionary impacts of mass extinctions. Proc Natl Acad Sci 388 U S A. 2001:98(10):5393-5398. https://doi.org/10.1073/PNAS.101092598 389 Jarvis ED, Mirarab S, Aberer AJ, Li B, Houde P , Li C, Ho SYW, Faircloth BC, Nabholz B, 390 Howard JT, et al. Whole-genome analyses resolve early branches in the tree of life of modern 391 birds. Science. 2014:346(6215):1320-1331. https://doi.org/10.1126/science.1253451 392 Jiang XZ, Xiang MC, Liu XZ. Nematode-Trapping Fungi. Microbiol Spectr. 2017:5(1):10.1128. 393 https://doi.org/10.1128/microbiolspec.FUNK-0022-2016 394 Jin WT, Gernandt DS, Wehenkel C, Xia XM, Wei XX, Wang XQ. Phylogenomic and ecological 395 analyses reveal the spatiotemporal evolution of global pines. Proc Natl Acad Sci U S A. 396 2021:118(20):e2022302118. https://doi.org/10.1073/PNAS.2022302118/-/DCSUPPLEMENTAL 397 Katoh K, Standley DM. MAFFT Multiple Sequence Alignment Software Version 7: Improvements 398 in Performance and Usability. Mol Biol Evol. 2013:30(4):772-780. 399 https://doi.org/10.1093/MOLBEV/MST010 400 Korstian JM, Paulat NS, Platt RN, Stevens RD, Ray DA. SINE-Based Phylogenomics Reveal 401 Extensive Introgression and Incomplete Lineage Sorting in Myotis. Genes. 2022:13(3):399. 402 https://doi.org/10.3390/GENES13030399 403 L Rocha J, Silva P, Santos N, Nakamura M, Afonso S, Qninba A, Boratynski Z, Sudmant PH, 404 Brito JC, Nielsen R, et al. North African fox genomes show signatures of repeated introgression 405 and adaptation to life in deserts. Nat Ecol Evol. 2023:7(8):1267-1286. 406 https://doi.org/10.1038/S41559-023-02094-W 407 Letunic I, Bork P. 2021. Interactive Tree Of Life (iTOL) v5: an online tool for phylogenetic tree 408 display and annotation. Nucleic Acids Res. 2021:49(W1), W293-W296. 409 https://doi.org/10.1093/nar/gkab301 410 Li Y, Liu Z, Liu C, Shi Z, Pang L, Chen C, Chen Y, Pan R, Zhou W, Chen X, et al. HGT is 411 widespread in insects and contributes to male courtship in lepidopterans. Cell. 412 2022:185(16):2975-2987.e10. https://doi.org/10.1016/J.CELL.2022.06.014 413 Liu K, Zhang W, Lai Y , Xiang M, Wang X, Zhang X, Liu X. Drechslerella stenobrocha genome 414 illustrates the mechanism of constricting rings and the origin of nematode predation in fungi. BMC 415 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 15 Genomics. 2014:15(1):114. https://doi.org/10.1186/1471-2164-15-114 416 Lopes F, Oliveira LR, Kessler A, Beux Y, Crespo E, Cárdenas-Alayza S, Majluf P, Sepúlveda M, 417 Brownell RL, Franco-Trecu V, et al. Phylogenomic Discordance in the Eared Seals is best 418 explained by Incomplete Lineage Sorting following Explosive Radiation in the Southern 419 Hemisphere. Syst Biol. 2021:70(4):786-802. https://doi.org/10.1093/SYSBIO/SYAA099 420 Lowery CM, Fraass AJ. Morphospace expansion paces taxonomic diversification after end 421 Cretaceous mass extinction. Nat Ecol Evol. 2019:3(6):900-904. 422 https://doi.org/10.1038/S41559-019-0835-0 423 Marques DA, Meier JI, Seehausen O. A Combinatorial View on Speciation and Adaptive 424 Radiation. Trends Ecol Evol. 2019:34(6):531-544. https://doi.org/10.1016/J.TREE.2019.02.008 425 Meleshko O, Martin MD, Korneliussen TS, Schröck C, Lamkowski P, Schmutz J, Healey A, 426 Piatkowski BT, Shaw AJ, Weston DJ, et al. Extensive Genome-Wide Phylogenetic Discordance Is 427 Due to Incomplete Lineage Sorting and Not Ongoing Introgression in a Rapidly Radiated 428 Bryophyte Genus. Mol Biol Evol. 2021:38(7):2750-2766. 429 https://doi.org/10.1093/MOLBEV/MSAB063 430 Mendes FK, Hahn MW. Gene Tree Discordance Causes Apparent Substitution Rate Variation. Syst 431 Biol. 2016:65(4):711-721. https://doi.org/10.1093/SYSBIO/SYW018 432 Minh BQ, Schmidt HA, Chernomor O, Schrempf D, Woodhams MD, V on Haeseler A, Lanfear R, 433 Teeling E. IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the 434 Genomic Era. Mol Biol Evol. 2020:37(5):1530-1534. 435 https://doi.org/10.1093/MOLBEV/MSAA015 436 Mistry J, Chuguransky S, Williams L, Qureshi M, Salazar GA, Sonnhammer ELL, Tosatto SCE, 437 Paladin L, Raj S, Richardson LJ, et al. Pfam: The protein families database in 2021. Nucleic Acids 438 Res. 2021:49(D1):D412-D419. https://doi.org/10.1093/NAR/GKAA913 439 Nevado B, Wong EL Y , Osborne OG, Filatov DA. Adaptive Evolution Is Common in Rapid 440 Evolutionary Radiations. Curr Biol. 2019:29(18):3081-3086.e5. 441 https://doi.org/10.1016/J.CUB.2019.07.059 442 Parins-Fukuchi C, Stull GW, Smith SA. Phylogenomic conflict coincides with rapid 443 morphological innovation. Proc Natl Acad Sci U S A. 2021:118(19):e2023058118. 444 https://doi.org/10.1073/PNAS.2023058118 445 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 16 Pease JB, Haak DC, Hahn MW, Moyle LC. Phylogenomics Reveals Three Sources of Adaptive 446 Variation during a Rapid Radiation. PLoS Biol. 2016:14(2):e1002379-24. 447 https://doi.org/10.1371/JOURNAL.PBIO.1002379 448 R Core Team. R: A language and environment for statistical computing. 2023. 449 Rhodes JA, Baños H, Mitchell JD, Allman ES. MSCquartets 1.0: quartet methods for species trees 450 and networks under the multispecies coalescent model in R. Bioinformatics 451 2021:37(12):1766-1768. https://doi.org/10.1093/BIOINFORMATICS/BTAA868 452 Rivas-González I, Rousselle M, Li F, Zhou L, Dutheil JY , Munch K, Shao Y , Wu D, Schierup MH, 453 Zhang G. Pervasive incomplete lineage sorting illuminates speciation and selection in primates. 454 Science. 2023:380(6648):eabn4409. https://doi.org/10.1126/science.abn4409 455 Salichos L, Rokas A. Inferring ancient divergences requires genes with strong phylogenetic 456 signals. Nature. 2013:497(7449):327-331. https://doi.org/10.1038/NA TURE12130 457 Sepkoski JJ. Rates of speciation in the fossil record. Philos Trans R Soc Lond B Biol Sci. 458 1998:353(1366):315-326. https://doi.org/10.1098/RSTB.1998.0212 459 Shen XX, Steenwyk JL, Rokas A. Dissecting Incongruence between Concatenation- and 460 Quartet-Based Approaches in Phylogenomic Data. Syst Biol. 2021:70(5):997-1014. 461 https://doi.org/10.1093/SYSBIO/SYAB011 462 Silvestro D, Cascales-Miñana B, Bacon CD, Antonelli A. Revisiting the origin and diversification 463 of vascular plants through a comprehensive Bayesian analysis of the fossil record. New Phytol. 464 2015:207(2):425-436. https://doi.org/10.1111/NPH.13247 465 Sievers F, Higgins DG. Clustal Omega for making accurate alignments of many protein sequences. 466 Protein Sci. 2018:27(1):135-145. https://doi.org/10.1002/pro.3290 467 Song H, Wang Y, Shao H, Li Z, Hu P, Yap-Chiongco MK, Shi P, Zhang T, Li C, Wang Y, Ma P, 468 Vinther J, Wang H, Kocot KM. Scaphopoda is the sister taxon to Bivalvia: Evidence of ancient 469 incomplete lineage sorting. Proc Natl Acad Sci U S A. 2023:120(40):e2302361120. 470 https://doi.org/10.1073/pnas.2302361120 471 Steenwyk JL, Shen XX, Lind AL, Goldman GH, Rokas A. A Robust Phylogenomic Time Tree for 472 Biotechnologically and Medically Important Fungi in the Genera Aspergillus and Penicillium. 473 mBio. 2019:10(4):e00925-19. https://doi.org/10.1128/mBio.00925-19 474 Steenwyk JL, Lind AL, Ries LNA, Dos Reis TF, Silva LP, Almeida F, Bastos RW, Fraga da Silva 475 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 17 TFC, Bonato VLD, Pessoni AM, Rodrigues F, Raja HA, Knowles SL, Oberlies NH, Lagrou K, 476 Goldman GH,Rokas A. Pathogenic Allodiploid Hybrids of Aspergillus Fungi. Curr Biol. 477 2020a:30(13):2495-2507.e7. https://doi.org/10.1016/j.cub.2020.04.071 478 Steenwyk JL, Buida TJ 3rd, Li Y , Shen XX, Rokas A. ClipKIT: A multiple sequence alignment 479 trimming software for accurate phylogenomic inference. PLoS Biol. 2020b:18(12):e3001007. 480 https://doi.org/10.1371/journal.pbio.3001007 481 Steenwyk JL, Buida TJ, Labella AL, Li Y, Shen XX, Rokas A. PhyKIT: a broadly applicable 482 UNIX shell toolkit for processing and analyzing phylogenomic data. Bioinformatics. 483 2021;37(16):2325-2331. https://10.1093/bioinformatics/btab096 484 Steenwyk JL, Li Y , Zhou X, Shen XX, Rokas A. Incongruence in the phylogenomics era. Nature 485 Reviews Genetics. 2023:24(12):834-850. https://doi.org/10.1038/s41576-023-00620-x 486 Varga T, Krizsán K, Földi C, Dima B, Sánchez-García M, Sánchez-Ramírez S, Szöllő si GJ, 487 Szarkándi JG, Papp V , Albert L, et al. Megaphylogeny resolves global patterns of mushroom 488 evolution. Nat Ecol Evol. 2019:3:668-678. https://doi.org/10.1038/s41559-019-0834-1 489 Wang J, Street NR, Park EJ, Liu J, Ingvarsson PK. Evidence for widespread selection in shaping 490 the genomic landscape during speciation of Populus. Mol Ecol. 2020:29(6):1120-1136. 491 https://doi.org/10.1111/MEC.15388 492 Wu H, Wang Z, Zhang Y , Frantz L, Roos C, Irwin DM, Zhang C, Liu X, Wu D, Huang S, et al. 493 Hybrid origin of a primate, the gray snub-nosed monkey. Science. 2023:380(6648):eabl4997. 494 https://doi.org/10.1126/SCIENCE.ABL4997 495 Wu M, Kostyun JL, Hahn MW, Moyle LC. Dissecting the basis of novel trait evolution in a 496 radiation with widespread phylogenetic discordance. Mol Ecol. 2018:27(16):3301-3316. 497 https://doi.org/10.1111/MEC.14780 498 Yang E, Xu L, Yang Y , Zhang X, Xiang MC, Wang C, An Z, Liu XZ. Origin and evolution of 499 carnivorism in the Ascomycota (fungi). Proc Natl Acad Sci U S A 2012:109(27):10960-10965. 500 https://doi.org/10.1073/PNAS.1120915109 501 Yang Y , Yang E, An Z, Liu X. Evolution of nematode-trapping cells of predatory fungi of the 502 Orbiliaceae based on evidence from rRNA-encoding DNA and multiprotein sequences. Proc Natl 503 Acad Sci U S A. 2007:104(20):8379-8384. https://doi.org/10.1073/PNAS.0702770104 504 Yang Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol. 505 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 18 2007:24(8):1586-1591. https://doi.org/10.1093/MOLBEV/MSM088 506 507 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 19 Figure 1 Phylogeny of nematode-trapping fungi. Their phylogenic relationships were determined 508 using concatenated nucleotide sequences of the single-copy orthologous genes present in all 509 species. Dactylella cylindrospora, a non-NTF species, was used as the outgroup. Bootstrap values 510 were 100% on each node. Gene-concordance factors (gCF) values were calculated by IQ-TREE 511 and annotated on each node, with green indicating nodes greater than 60% and red indicating 512 nodes less than 60%. 513 514 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 20 Figure 2 Extensive conflict between the gene trees and the species tree. a. Densitree plot. Blue 515 represents the gene trees, and red represents the consensus tree inferred by the Densitree software, 516 which is consistent with the topology of the species tree. b. A plot resulting from 517 multi-dimensional scaling (MDS) analysis illustrates the topological differences between the gene 518 trees (denoted by blue dots) and the species tree (denoted by the red pentagram). 519 520 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 21 Figure 3 Origins of conflict between the gene trees and the species tree. a. ILS analysis based on 521 the Multispecies coalescent (MSC) analysis. Blue circles represent four-taxa scenarios in which 522 the topology can be explained solely by the ILS. Red triangles represent scenarios in which this 523 hypothesis is rejected, indicating that the topology is explained by other factors. The closer the 524 blue circles to the center of the triangle, the stronger the influence of ILS. b. Schematic 525 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 22 representation of D-statistic results. c. Reticulate phylogenetic tree inferred by Phylonet, with red 526 indicating gene introgression sites. When the number of hybridization events was set to 3, the tree 527 inferred by PhyloNet matched with the species tree, and the fit was optimal. d. Sankey diagram 528 depicting the suspected HGT events among NTF and the predicted sources of the genes. 529 530 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 23 531 Figure 4 Divergence nodes and cumulative branch lengths for the three NTF genera. a. 532 Topological structures of the three gene trees and their divergent branch lengths. b. Stacked bar 533 chart showing the proportions of the three types of gene tree topology inconsistent with the 534 species tree. c. Box plot of cumulative branch lengths for four types of gene trees. 535 536 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint 24 Figure 5 Functional enrichment analysis. Functional enrichment analysis of the ILS genes that are 537 linked to the divergence of three NTF lineages and display signs of significant positive selection. 538 The Gene Ontology (GO) terms enriched among those associated with the cell membrane system 539 and polarity division are shown. 540 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2024. ; https://doi.org/10.1101/2024.03.21.586083doi: bioRxiv preprint