The flow chart presented in Supplementary Figure S1 illustrates the main steps in our study. A total of 19 fungi were selected, including six rust fungi, three anther smuts, seven nonphytoparasitic fungi in Pucciniomycotina, and three nonbiotrophic fungi in Ustilaginomycotina. Among them, Austropuccinia psidii APSI1 (Aps), Cronartium quercuum f. sp. fusiforme G11 version 1.0 (Cqf) (Pendleton et al., 2014), Hemileia vastatrix XXXIII (Hvs), Melampsora larici-populina 98AG31 version 1.0 (Mlp) (Duplessis et al., 2011), Puccinia striiformis f. sp. tritici PST-130 (Pst) (Cantu et al., 2011), and Uromyces viciae-fabae I2 (Uvf) were selected as analytical data for representative rust fungi. Microbotryum intermedium 1389 BM 12 12 (Min) (Branco et al., 2017), Microbotryum lychnidis-dioicae p1A1 Lamole (Mld) (Perlin et al., 2015), and Microbotryum silenes-dioicae 1303 (Msd) were selected as phytoparasitic fungi for convergence analyses. The following nonphytoparasitic fungi were included for comparison: Cystobasidiopsis lactophilus JCM 7595 (Cla), Cystobasidium minutum MCA 4210 version 1.0 (Cmi) (Camiolo et al., 2019), Erythrobasidium yunnanense JCM 10687 (Eyu), Glaciozyma antarctica PI12 (Gan), Leucosporidiella creatinivora 62-1032 version 1.0 (Lcr), (Mondo et al., 2017), Rhodotorula graminis strain WP1 version 1.1 (Rgr) (Firrincieli et al., 2015), and Sporobolomyces roseus SR19 (Sro). Finally, Acaromyces ingoldii MCA 4198 (Ain) (Kijpornyongpan et al., 2018), Malassezia globosa CBS 7966 (Mgl) (Xu et al., 2007), and Pseudozyma antarctica T-34 (Pan) (Morita et al., 2013), which belong to three different classes in Ustilaginomycotina, were selected as the outgroup taxa.

The data of 11 selected fungi (Cqf, Mlp, Pst, Min, Mld, Cmi, Lcr, Rgr, Ain, Mgl, and Pan) are publicly available in the JGI MycoCosm portal (https://mycocosm.jgi.doe.gov/mycocosm/home). The protein and coding sequences (CDSs) of Msd were downloaded from the NCBI database. Genome assemblies without gene annotation for an additional seven fungi (Aps, Hvs, Uvf, Cla, Eyu, Gan, and Sro) publicly available in the NCBI database were downloaded. Then, protein-coding genes were predicted by BRAKER2 (Hoff et al., 2019). Detailed information for the selected fungi is available in Supplementary Table S1.

All amino acid sequences of the 19 fungi were analyzed to obtain orthologous groups by OrthoFinder with default parameters (Emms and Kelly, 2016). To determine the phylogenetic relationships between the selected fungi, we conducted a phylogenetic analysis using all the single-copy genes. Multiple sequence alignments of single-copy orthologous genes were performed using MAFFT (Katoh and Standley, 2013), and the gaps were removed using trimAl (Capella-Gutiérrez et al., 2009). Next, the genes were concatenated to generate a supergene sequence. The best-fitting model for phylogenetic analysis was estimated using ModelTest-NG (Darriba et al., 2020). With the input of these 19 supergene sequences, a phylogenetic tree was constructed using RAxML-NG software with 1,000 bootstrap replicates (Kozlov et al., 2019).

Only single-copy genes were utilized to identify convergent amino acid sites. First, ancestral states were reconstructed using CodeML. We set node “a” as the ancestral state of rust fungi, node “b” as the most recent common ancestor (MRCA) of rust fungi and other nonphytoparasitic fungi, node “c” as the ancestral state of anther smuts, and node “d” as the MRCA of anther smuts and other nonphytoparasitic fungi (Figure 1). We identified convergent amino acid sites between rust fungi and anther smuts with the following rules: 1) the amino acid residues of node a and node c were identical, 2) amino acid changes were inferred to have occurred between node a and node b, and 3) amino acid changes were inferred to have occurred between node c and node d. Noteworthily, detection of false and random convergence events may result from limited taxon sampling.

Considering that conserved sites permit a more accurate inference of the ancestral state to avoid false convergence (Xu et al., 2017), the convergent amino acid sites with high support were selected with the following rules: 1) the amino acid residues of all rust fungi and anther smuts were identical, 2) there was a significant difference in the usage of amino acids between phytoparasitic and nonphytoparasitic fungi (p values less than 0.0005, calculated by Fisher’s exact test), and 3) more than five nonphytoparasitic fungi in Pucciniomycotina and more than two nonphytoparasitic fungi in Ustilaginomycotina tended to use the same amino acid. In general, the ancestral states of node a and node c should be more accurate when the first rule is met. The second and third rules were used to improve the accuracy of the ancestral states of node b and node d.

To filter out noise resulting from chance amino acid substitutions, we then performed a Poisson test to verify whether the observed number of highly supported convergent sites of each gene was significantly greater than the number expected with random substitution under the JTT-fgene model, following the method of Zou and Zhang (2015) (conv_cal, https://github.com/ztzou/conv_cal). The Benjamini–Hochberg method was used to correct multiple testing (q < 0.05).

The aligned protein-coding sequences and the species tree obtained from the phylogenetic analysis were used to detect positive selection using the branch-site model in the CodeML program of PAML (Yang, 1997). For every single-copy gene, the branch-site model (model = 2 and NSsites = 2) of codon evolution was used. To identify genes under adaptive evolution related to phytoparasitism, the ancestral branch of rust fungi and the ancestral branch of anther smuts were designated as foreground branches, and the other branches were designated as background branches. We compared the model (which allowed sites to be under positive selection; fix_omega = 0) with the null model (in which sites may evolve neutrally or under purifying selection; fix_omega = 1 and omega = 1) with likelihood ratio tests (LRTs) performed by the CodeML program in PAML. The LRT statistics were evaluated by calculating the 2ΔLnL between the abovementioned models. The significance of LRT statistics was determined by a χ ² distribution, and p values were adjusted by following the Benjamini–Hochberg procedure. The alternative model was considered to be significantly better than the null model when q < 0.05.

MCMCTree analysis was applied to estimate divergence times among the 19 fungi (Yang, 1997). According to the TimeTree database (http://www.timetree.org/), ancestral Basidiomycota fungi diverged into Ustilaginomycotina and Pucciniomycotina 489 million years ago (MYA), while the MRCA of Puccinia and Melampsora diverged 114.5 MYA (Kumar et al., 2017). Using these two fossil calibration points, we estimated the divergence times of the 19 selected fungi based on all single-copy genes. After acquiring this ultrametric tree, orthologous gene group expansion and contraction were detected by CAFE v5.0 (Mendes et al., 2020). Only the multicopy gene group that contained at least one protein in no fewer than eight species was chosen for analysis. Then, the p-value was calculated for each orthologous group, and orthologous groups with p values less than 0.05 were considered to have accelerated rates of expansion and contraction. The groups that had specifically expanded and contracted in rust fungi and anther smuts were considered candidate gene groups. We also performed a manual check to identify genes specifically missing or present in selected rust fungi and anther smuts.

To study whether convergent evolution occurs in rust fungi and anther smuts, we selected 19 basidiomycetous fungi, including 16 fungi in Pucciniomycotina (six rust fungi, three anther smuts, and seven nonphytoparasitic fungi) and three nonphytoparasitic fungi in Ustilaginomycotina, for evolutionary analysis. We first identified all the orthologous gene groups in these fungi. In total, 19,479 orthologous gene groups were identified, including 530 single-copy orthologous gene groups and 1,276 multi-copy orthologous gene groups across all selected fungi (Supplementary Tables S2, S3). Consistent with the findings of previous studies (Cuomo et al., 2017), it is notable that rust fungi have very different gene contents from other selected fungi, harboring a larger number of unique paralogs (2,592, on average) and unclustered genes (7,004, on average) (Figure 2; Supplementary Table S2).

To identify the genomic signatures of convergent evolution between rust fungi and anther smuts, a maximum likelihood phylogeny was reconstructed on the basis of all the proteins encoded by the single-copy genes. The constructed genome-wide phylogenetic tree showed reasonable phylogenetic relationships. Within Pucciniomycotina, three anther smuts were clustered together with other Microbotryomycetes fungi. Compared with other selected Pucciniomycotina fungi, six rust fungi have occupied a basal position (Figure 1).

All 530 single-copy genes were used for genomic convergence analyses. Four hundred forty-four genes with 2,150 convergent amino acid substitutions between rust fungi and anther smuts were identified based on our initial criteria. Because genomic signatures of molecular convergence are generally not strong (Zou and Zhang, 2015), we surmised that random and false convergence events occurred frequently. Seventy-six genes with 84 highly supported sites that fit further criteria were chosen for the Poisson test. Overall, the Poisson test revealed that only two genes (OG0002989 and OG0003182) exhibited stronger signatures of convergent amino acid substitution than expected by chance (q < 0.05) (Supplementary Table S4).

These two convergent genes encode V-type H⁺-transporting ATPase subunit d (ATPeV0D) and ubiquitin-small subunit ribosomal protein S27Ae (RP-S27Ae). For ATPeV0D, two convergent substitutions occurred at site 17 (valine in nonphytoparasitic fungi to isoleucine in phytoparasitic fungi) and site 92 (lysine in nonphytoparasitic fungi to threonine in phytoparasitic fungi). These two sites are located in helix 1 and helix 5, respectively. Helix 1 is required for ATPase subunit assembly, and helix 5 contributes to the translocation of protons (Owegi et al., 2006). For RP-S27Ae, one convergent substitution occurred at site 57 (serine in nonphytoparasitic fungi to alanine in phytoparasitic fungi). This site is a phosphorylation site in the ubiquitin domain that regulates ubiquitin homeostasis (Figure 3) (Lee et al., 2017).

Across all 530 gene orthologues, we found convergent signals of positive selection in 51 genes during the evolution of rust fungi and anther smuts (Supplementary Table S5). To further characterize the functions of positively selected genes, we performed functional annotation with the Kyoto Encyclopaedia of Genes and Genomes (KEGG) database. In total, 43 of the 51 genes had KEGG annotations. Among them, 25 genes were associated with 27 pathways (Table 1). The genes were associated with very diverse pathways, including four main categories: metabolism, genetic information processing, environmental information processing, and cellular processes. There were eight genes in translation pathways (pathway ids: 00970, 03008, 03010, 03013, and 03015), which was the most abundant category. In addition, we found three positively selected genes associated with secretory pathways. They included a protein transport protein (Sec24, OG0003606), calnexin (CANX, OG0003436), and dolichyl-phosphate-mannose-protein mannosyltransferase (POMT, OG0002993). However, genes associated with other functions were relatively small in number.

To identify convergent gene gain or loss in rust fungi and anther smuts, an analysis was performed in CAFÉ. First, a clock-calibrated phylogenetic tree was constructed based on the concatenation of all single-copy genes (Supplementary Figure S2). A total of 4,559 orthologous gene groups were used for CAFÉ analysis by filtering (Supplementary Table S6). There were a few orthologous gene groups that showed a convergent pattern of expansion or contraction between rust fungi and anther smuts. Only two gene groups were specifically expanded in the selected rust fungi and anther smuts. One orthologous group (OG0000061) was annotated as oligopeptide transporters (OPTs), and the other group (OG0000075) was associated with vacuolar protease (Table 2). However, there were no orthologous gene groups that showed a convergent pattern of contraction. Furthermore, no genes were specifically absent or present in rust fungi and anther smuts.