Functional dynamics reveal the response of the crabapple (Malus sp.) phyllosphere microbiome to Gymnosporangium yamadae infection

9 The plant microbiome is crucial for maintaining plant health under disease stress. While 10 considerable research has focused on belowground systems, the phyllosphere microbiome remains 11 underexplored, particularly regarding its functional responses to pathogen infection. In this study, 12 we investigated the phyllosphere microbiome dynamics of crabapple ( Malus ‘Kelsey’) infected by 13 Gymnosporangium yamadae using metatranscriptomic sequencing across six stages of rust disease 14 progression. Our analysis revealed a general increase in the diversity of fungi, bacteria and viruses 15 in infected leaves. Notably, fungal diversity negatively correlated with bacterial diversity, 16 reflecting competitive interactions during disease development. Microbial taxa in diseased leaves 17 exhibited heightened expression activity compared to healthy leaves, with fungi progressively 18 dominating the microbial community. Functional co-occurrence networks of the phyllosphere 19 microbiome in infected leaves were more complex than in healthy leaves, suggesting adaptive 20 reorganization in response to pathogen invasion. Differentially expressed genes at each stage were 21 significantly enriched in carbohydrate metabolism pathways and gene regulation-related functions, 22 enabling functional adaptations to rust diseases. Random forest modeling identified key microbial 23 transcripts associated with pathogen abundance, including beneficial microbes like Saitozyma 24 podzolica, which secretes glucan-degrading enzymes, and Mortierella elongata, involved in sterol 25 biosynthesis and plant resistance. Conversely, Alternaria alternata emerged as a major 26 pathobiome contributor, secreting enzymes that degrade plant cell wall components (e.g., pectin, 27 cellulose, lignin, and xylan), and engaging in MAPK signaling pathways critical for pathogenesis. 28 Our findings underscore the vital role of the phyllosphere microbiome in mediating 29 plant-pathogen interactions and shaping disease progression, providing a foundation for 30 microbiome-based strategies to enhance plant resilience. 31 Importance 32 Our findings reveal the dynamic shifts in the expression patterns and adaptive functional strategies 33 of crabapple phyllosphere microbiome in response to the Gymnosporangium yamadae infection. 34 We identified several key microbes that may play vital roles in pathogenesis and speculated on 35 their roles and functions in plant-pathogen interactions. In conclusion, our study highlights the 36 potential of the phyllosphere microbiome in regulating plant health. 37

metatranscriptome, phyllosphere microbiome, rust disease, Malus sp. 38 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint 1 Introduction 39 The plant microbiome encompasses an entire community of microorganisms and their 40 associated molecular spectrum within plant environments, forming a diverse and complex 41 ecological network that significantly impacts plant health, resilience and defense mechanisms [1–5]. 42 Current research underscores that plant fitness and immunity are not shaped solely by individual 43 pathogens or beneficial species but by the collective interactions within the microbiome [6]. 44 Beneficial microbial communities, especially within the rhizosphere and phyllosphere, play 45 critical roles by modulating immune res ponses, competing with pathogens, and producing 46 metabolites that enhance defense [3,7–9]. For instance, Pseudomonas and Bacillus in the rhizosphere 47 produce antibiotics that directly suppress pathogens [10,11]. Furthermore, plants can actively engage 48 in recruiting beneficial microbes, particularly under pathogen pressure, to alleviate stress and 49 restructure microbial communities assemblies to enhance resilience [12,13]. For example, the 50 rhizosphere microbiome of the Ralstonia solanacearum-resistant tomato variety Hawaii 7996 has 51 a higher abundance of beneficial microbes, such as Flavobacteriaceae , Sphingomonadaceae and 52 Pseudomonadaceae compared to the susceptible variety Monkeymaker, providing a natural 53 microbial buffer against pathogens [14]. Similarly, Arabidopsis thaliana recruits protective bacteria, 54 such as Microbacterium, Stenotrophomonas, and Xanthomonas, which collectively induce 55 systemic immunity against downy mildew [15]. Likewise, in wheat, exposure to high pathogen 56 pressure from Fusarium pseudograminearum leads to an enrichment of Stenotrophomonas 57 rhizophila, functioning as an early defense alert in the root endosphere [16]. This observed increase 58 in protective microbial abundance in various plant-pathogen systems aligns with the studies 59 showing that stressed plants under pathogen pressure can engage in a “cry for help”, attracting 60 beneficial microbes that bolster plant defenses [15,17,18]. 61 As research advances, scientists increasingly recognize the potential of leveraging 62 microbiome functions for sustainable crop protection and yield improvement [3]. Functional 63 studies have linked specific gene expression profiles to pathogenic responses, further deepening 64 our understanding of microbiome interactions [19]. For example, during Fusarium wilt in peppers, 65 genes associated with detoxification, biofilm formation, and chemotaxis were significantly 66 enriched, suggesting targeted adaptations in the root endosphere, a coordinated microbial response 67 to bolster plant defenses [20]. Similarly, Rhizoctonia solani infection in sugar beets led to an 68 increase in Chitinophagaceae and Flavobacteriaceae, which produced antifungal enzymes and 69 metabolites like phenazines, polyketides, and siderophores, providing an add itional line of defense 70 against pathogens [21]. 71 However, not all microbial members benefit plant health, some can be detrimental by forming 72 harmful partnerships with pathogens, disrupting plant resilience, and facilitating disease 73 progression [6,22,23]. For example, Verticillium dahlia infection was shown to stabilize certain 74 microbial networks within the bulk soil and the rhizosphere, with viruses becoming as central 75 players in the ‘pathobiome’, potentially aiding pathogen invasion and colonization [24]. Similarly, 76 studies on the diseases affecting tomato, tobacco and other plants have shown that certain 77 microbes associated with nematode pathogens release enzymes that degrade plant cell walls, 78 which aids pathogen entry into roots and exacerbates disease progression [25,26]. This multifaceted 79 relationship underscores the role of the plant microbiome in pathogenesis cannot be viewed as a 80 straightforward dichotomy of pathogenic or protective effects, but must be analyzed in context [27]. 81 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint Despite the insights gained from studies on the rhizosphere, most research on 82 plant-microbiome-pathogen interactions has focused on the belowground microbiome, particularly 83 the rhizosphere. In contrast, the aboveground microbiome, or phyllosphere-constituting the above 84 ground plant microbiome-remains relatively underexplored despite it holds significant potential as a 85 critical line of defense against airborne pathogens [28–30]. Under abiotic stress, the Arabidopsis 86 mutant mfec exhibited significantly reduced phyllosphere microbiome diversity, which was linked 87 to the loss of immune and MIN7 pathways, resulting in dysbiosis and plant symptoms. This 88 underscores the critical role of phyllosphere microbiome homeostasis in maintaining plant health 89 [29]. Increasing evidence also suggests the assembly strategies of the phyllosphere microbiome are 90 dynamic, with strategies varying in response to different pathogen pressures [31]. For instance, the 91 phyllosphere microbiome of citrus plants infected with Diaporthe citri exhibited lower community 92 evenness and a more complex co-occurrence network, with the emergence of new microbial taxa 93 [30]. Similarly, temporal analysis of the phyllosphere and rhizosphere microbiomes in 94 Phytophthora sojae -infected and Septoria glycines -infected soybean plants revealed that the 95 structure of the phyllosphere microbial communities was more responsive to pathogen infection 96 and disease progression, with a noticeable increase in the abundance of saprophytic fungi [32]. 97 During disease progression, the phyllosphere microbial communities in apple, wheat and tobacco 98 were observed to exhibit higher diversity and more complex co-occurrence networks [8,33–35]. 99 However, while many of the studies to date has focused on identifying the microbial components of 100 the disease-induced phyllosphere microbiome, the functional responses of these communities and 101 their direct links to plant health have been largely overlooked [36]. 102 Crabapple trees (Malus spp.), cherished in landscaping for their attractive shape and leaf color, 103 face a significant threat from the rust fungus Gymnosporangium yamadae , which severely 104 diminishes their ornamental value [37]. G. yamadae, a heteroecious fungal pathogen, requires two 105 different hosts (Malus spp. and Juniperus chinensis) to complete its infection cycle and produces 106 four morphologically distinct spores [38,39]. In early spring, brownish telia break through the 107 epidermis of the teilal host (J. chinensis ), forming a bright yellow gelatinous mass and releasing 108 haploid basidiospores that infect Malus species. Initial infection manifests as chlorotic spots and the 109 development of yellowish droplets (spermogonia) on the upper surface of the infected leaves [38]. 110 The biotrophic and unculturable nature of G. yamadae complicates traditional pathogen-host 111 studies, prompting researchers to pivot towards investigating the response patterns of microbiome 112 [13]. Recent studies have highlighted the essential role of microbiome in modulating plant 113 responses to biotrophic pathogens. For example, studies of the phyllosphere microbiome in 114 cucumber infected with powdery mildew cucumber revealed significant differences in bacterial 115 alpha diversity across varying disease severity, with more severely diseased leaves exhibiting 116 higher bacterial diversity and an increased abundance of certain beneficial microbes [40]. Similarly, 117 in Arabidopsis thaliana infected with the downy mildew pathogen Hyaloperonospora 118 arabidopsidis, three bacterial species were specifically enriched in the rhizosphere, collectively 119 inducing systemic resistance and promoting plant growth [15]. In the studies of apple and G. 120 yamadae biotrophic interactions, a meta-transcriptome analysis was conducted to compare the 121 phyllosphere fungal communities in infected and uninfected apple leaves ( M. domestica cv. Fuji) 122 at 10- and 30-days post-inoculation (dpi). This analysis revealed a significant shift in the 123 community composition occurs during the later stages of infection, with a notable increase in the 124 abundance of Alternaria and Fonsecaea species at 30 dpi [34]. Subsequent research expanded on 125 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint these findings with a comprehensive time-course analysis of fungal and bacterial community 126 dynamics, as well as key leaf metabolites, in two crabapple varieties, ‘Flame’ and ‘Kelsey’, over 127 six distinct stages of rust disease progression [8]. The integrative study highlighted that the leaves 128 regulate disease progression by secreting specific metabolites, which help mediate the enrichment 129 of potential beneficial microbes, supporting the “cry for help” strategy, where plants under stress 130 recruit microbial allies for defense [8]. A complementary investigation explored the diversity and 131 structural dynamics of endophytic microbial communities in apple ( Malus domestica cv. Gala) 132 leaves across various stages of rust infection. This study utilized amplicon sequencing to offer 133 foundational insights into the predicted functional profiles of endophytic fungi and bacteria, 134 shedding light on the microbial shifts associated with disease progression [41]. While these studies 135 have advanced our understanding of the diversity and structural changes in the phyllosphere 136 microbiome of Malus species under G. yamadae infection, a comprehensive understanding of the 137 dynamic functional profiles of these microbial communities throughout the disease course is still 138 lacking. Such insights are critical for uncovering the adaptive responses of Malus phyllosphere 139 microbiome to G. yamadae infection and understanding their impact on host plant health. 140 This study aims to systematically investigate the functional response of the phyllosphere 141 microbiome in crabapple leaves infected by G. yamadae at various stages of lesion expansion using 142 meta-transcriptomic technology. Our main objectives were to (i) assess the diversity and 143 composition of phyllosphere microbial transcriptomes at various stages of lesion expansion; (ii) 144 pinpoint differentially expressed microbial genes at each stage of infection and examine their 145 dynamic expression profiles over time; (iii) elucidate the functional roles and activities of 146 phyllosphere microbiota and evaluate their potential interactions with plant disease processes and 147 plant health. 148 2 Results 149 2.1 The taxa transcriptomes diversity and structure patterns in the phyllosphere 150 with infection 151 The phyllosphere of crabapple leaves harbors a complex microbial community whose 152 diversity and structure undergo significant shifts in response to G. yamadae infection. Based on 153 the analysis of transcript expression levels and taxonomic annotations, the microbial community 154 in the phyllosphere was categorized into bacterial, fungal, viral, and archaeal transcriptomes, with 155 species-level expression quantified as FPKM values (Table S1). Archaea were excluded from 156 further transcriptome diversity and structural analyses due to their low diversity. 157 The alpha diversity analysis of the phyllosphere microbiome during G. yamadae infection 158 revealed distinct trends in bacterial, fungal, and viral transcriptomes. Overall, fungal 159 transcriptomes consistently exhibited greater alpha diversity than that of bacterial and viral 160 transcriptomes (Figure 1a, S1a). For bacterial transcriptomes, diseased leaves generally exhibited 161 higher species richness compared to healthy leaves. The Shannon index in healthy leaves peaked at 162 the third stage and declined to its lowest value by the sixth stage. In early to mid-disease stages 163 (Stages 1 to 4), healthy leaves displayed greater bacterial diversity than diseased leaves. However, 164 as the lesions expanded, the Shannon index and Pielou’s evenness index of bacterial 165 transcriptomes in diseased leaves incrementally increased, peaking at the sixth stage (Figure 1a, 166 S1b). For fungal transcriptomes, healthy leaves did not display any distinct patterns in diversity 167 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint indices. However, G. yamadae infection initially resulted in a decline in species richness, which 168 then rose to peak at the final stage. Interestingly, the Shannon index, as well as the Pielou’s 169 evenness index, declined steadily in diseased leaves, an opposite trend compared to bacterial 170 diversity (Figure 1, S1b; Table S2a). The correlation analysis between the Shannon indices of 171 bacterial and fungal communities through the simple linear model confirmed the inverse 172 relationship between the Shannon indices of bacterial and fungal communities, with bacterial 173 diversity increasing the fungal diversity decreasing in diseased leaves (Figure 1b). The viral alpha 174 diversity did not exhibit any significant patterns across stages or between leaf conditions (Figure 175 S1a; Table S2a). 176 177 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint Figure 1 Phyllosphere transcriptomes alpha diversity and structural changes of bacteria and fungi 178 under G. yamadae invasion. (a) The alpha diversity of bacterial and fungal transcriptomes in 179 diseased leaves and healthy leaves across six developmental stages of crabapple rust disease. Box 180 plots show the range of estimated values between the 25th and 75th percentiles, with the median, 181 minimum, and maximum observed values within each dataset. Different letters indicate 182 statistically significant differences was determined using one-way ANOV A with Tukey-HSD post 183 hoc test ( P < 0.05) or Kruskal-Wallis with Wilcoxon’s test ( P < 0.05). (b) The correlation of 184 Shannon index between bacterial and fungal transcriptomes. Solid lines represent the results of the 185 simple linear regression models, while the surrounding grey band represents the 95% confidence 186 interval. (c) PCoA of bacterial and fungal transcriptomes based on the Bray-Curtis distance matrix, 187 R² and P were calculated using PERMANOVA test. Asterisks show the P significance level, * P < 188 0.05, **P < 0.01, ***P < 0.001 and ****P < 0.0001 and NS denotes no statistical significance. 189 We found that phyllosphere bacterial, fungal and viral transcriptome compositions were 190 significantly influenced by three main factors: leaf condition (healthy or diseased), sampling stage, 191 and their interaction. Among these, leaf condition emerged as the most significant determinant of 192 transcriptome composition across all microbial groups (bacterial, fungal, and viral; Figure 1c, S1c; 193 Table 1, S2b). To better understand the compositional differences, we employed Principal 194 Coordinate Analysis (PCoA) based on Bray-Curtis distance. This analysis highlighted the distinct 195 clustering of microbial transcriptomes across conditions and stages. For fungal and viral 196 transcriptomes, the primary axis (PCoA 1) accounted for 46% and 48% of the total variance, 197 respectively, delineating the overall differences in transcriptome composition. Separation based on 198 leaf condition (healthy vs. diseased) was observed along the secondary axis (PCoA 2), which 199 explained an additional 33% of the variation. Further analysis showed that bacterial 200 transcriptomes in diseased leaves exhibited the highest Bray−Curtis dissimilarity during the early 201 stage of infection (Stage 1). In contrast, viral communities displayed peak dissimilarity during the 202 final stage of disease progression (Figure S1d). 203 Table 1 The effects of leaf condition and sampling stage on the structure of bacterial and fungal 204 transcriptomes based on PERMANOVA with 999 permutations. 205 Bacteria Fungi R² F P R² F P Leaf condition (LC) 0.33476 27.4413 0.001 0.44416 50.7015 0.001 Sampling stage (SS) 0.15774 2.5861 0.004 0.17891 4.0845 0.001 Interaction (LC & SS) 0.21473 3.5204 0.001 0.16668 3.8054 0.001 2.2 Overview of the composition and differential transcripts of phyllosphere 206 transcriptomes 207 We characterized the composition and differential taxonomic profiles of phyllosphere 208 transcriptomes in response to G. yamadae infection (Figure 2). The infection substantially altered 209 the transcriptomic landscape of the crabapple phyllosphere. In healthy leaves, bacterial 210 transcriptomes consistently dominated, typically accounting for 50% or more of the total 211 transcriptomes across all stages (Figure 2a). Except for the first stage and the sixth stage, the 212 relative expression abundance of microbial groups followed this order: bacteria > fungi > viruses. 213 Conversely, in diseased leaves, fungal transcript abundance progressively increased, peaking at 214 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint 96.2% in the final stage. Concurrently, the relative expression abundance of bacterial and viral 215 transcriptomes decreased, although bacterial transcripts remained more abundant than viral ones 216 throughout the progression. Additionally, archaeal transcripts were exclusively detected in the 217 final diseased stage. 218 219 Figure 2 Overview of the composition and differential transcripts of phyllosphere transcriptomes. 220 (a) The composition of phyllosphere transcriptomes in different stages, showing the relative 221 abundance of bacterial, fungal, and viral transcriptomes. (b) Taxonomic compositions of the 222 bacterial transcriptome at the order level and the fungal transcriptome at the family level, 223 highlighting the predominant groups in both healthy and diseased leaves at different stages of 224 infection. (c) Significantly differentially expressed transcripts in crabapple leaves under different 225 conditions. Differential expression analysis was performed using a generalized linear model 226 (GLM) to identity transcripts showing significant differences at each developmental stage ( P < 227 0.05, FDR corrected). 228 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint Taxonomical classification of the bacterial transcriptomes at the order level and fungal 229 transcriptomes at the family level reveals key shifts in microbial composition based on 230 transcriptomic profiles, highlighting the top 20 taxa transcripts by relative expression abundance 231 (Figure 2b). In diseased leaves, the dominant bacterial taxa included Enterobacterales (63.4%), 232 Rickettsiales (12.1%), Pseudomonadales (6.6%), Thiotrichales (5.1%) and Sphingobacteriales 233 (3.7%). In the healthy leaves, the most abundant bacterial orders were observed for 234 Enterobacterales (51.5%), Sphingobacteriales (27.9%), Micrococcales (4.3%), Cytophagales 235 (3.9%) and Pseudomonadales (3.9%). Several bacterial orders, including Thiotrichales, Nostocales, 236 Bacteroidales, Geodermatophilales, Rickettsiales, Xanthomonadales and Burkholderiales, showed 237 significant upregulation in diseased samples compared to healthy leaves (Figure 2c, S2). 238 Conversely, Cytophagales, Vibrionales, Rhodospir illales, Rhizobiales and Bacillales, were 239 significantly downregulated in diseased samples (P < 0.05). 240 For fungal transcriptomes, the family Pleosporaceae dominated both diseased and healthy 241 leaves, accounting for 48.2% and 64.6%, respectively (Figure 2b). Other families with higher 242 relative abundance in diseased samples included Ophiocordycipitaceae and Didymellaceae. 243 Fungal families represented by Dermateaceae were notably upregulated in diseased leaves 244 compared to healthy leaves and clustered together in the hierarchical analysis (Figure 2c). Unlike 245 the bacterial communities, fungal transcripts generally showed increased expression during rust 246 infection, with relatively fewer taxa downregulated. Interestingly, some bacterial taxa exhibited 247 stage-specific differential expression patterns. Sphingobacteriales and Micrococcales were 248 significantly upregulated in the third diseased stage and downregulated in the early and late stages. 249 However, no similar stage-dependent trends were observed in fungal communities. 250 2.3 Changes in phyllosphere gene expression upon G. yamadae infection 251 To determine the gene expression response during each stage of G. yamadae infection, we 252 quantified and visualized the number of expressed genes (Figure 3a). Overall, the total number of 253 expressed genes, as well as stage-specific gene expression, was lowest at the initial stage of rust 254 infection and peaked at the final stage. Pairwise comparisons of phyllosphere transcriptomes 255 between diseased and healthy crabapple leaves were c onducted at each stage to identify 256 differentially expressed genes (DEGs). Transcripts with significant differential expression 257 (log2Fold-Change ≥ 1 or ≤ -1, adjusted P -value < 0.05) were categorized DEGs (Figure 3b; Table 258 2, S3a). As the rust spots expanded, there was a marked increase in gene expression, with the 259 number of upregulated genes consistently surpassing the number of downregulated genes at each 260 stage. 261 Table 2 The results of differential gene expression analysis in crabapple phyllosphere at each 262 stage. 263 Stage Upregulated a Downregulated b Total 1 194 0 194 2 339 40 379 3 298 78 376 4 382 67 449 5 226 129 355 6 306 168 474 a Genes that were significantly upregulated in diseased leaves (log 2Fold-Change ≥ 1, adjusted 264 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint P-value < 0.05) compared to healthy leaves were classified as upregulated genes. 265 b Genes that were significantly downregulated in diseased leaves (log 2Fold-Change ≤ -1, adjusted 266 P-value < 0.05) compared to healthy leaves were classified as downregulated genes. 267 268 Figure 3 Expression profiling of phyllosphere target genes infected with G. yamadae . (a) The 269 overlap of expressed genes across different stages of rust disease, showing top 22 sets. (b) 270 Differentially expressed genes (DEGs) between diseased and healthy leaves at each stage of 271 infection. Red dots indicate upregulated genes in diseased samples compared to healthy samples 272 (log2Fold-Change ≥ 1, adjusted P-value < 0.05), green dots indicate downregulated genes in 273 diseased samples compared to healthy samples (log 2Fold-Change ≤ -1, adjusted P-value < 0.05), 274 and the yellow dots denote genes with no significant difference between diseased and healthy 275 samples. (c) Expression profile of DEGs clusters through c-means clustering. The color scale 276 indicates the degree of membership of each gene to the respective clusters based on expression 277 patterns, with darker colors representing higher membership to a given cluser. 278 To further investigate the characteristics of the DEGs that actively respond to G. yamadae 279 infection, we performed clustering analysis using the c-means clustering method and visualized 280 the expression profiles of each cluster (Table S3b; Figure 3c). The six resulting clusters exhibited 281 distinct expression patterns and taxonomic compositions. In particularly, Cluster 5 exhibited high 282 expression in healthy leaf tissues. The Klebsiella genus was the predominant transcript in both 283 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint diseased and healthy samples, accounting for 61% and 22.9% respectively, with all identified as 284 Klebsiella pneumoniae. Additionally, Sphingobacteriaceae bacterium was a major contributor in 285 healthy leaves, comprising 51% of the total transcripts. Clusters 1, 3 and 4 showed similar 286 expression trends but differed significantly in their taxonomic composition. Cluster 1 contained 287 genes with highest relative expression abundance in diseased leaves, primarily from viruses and 288 bacteria, such as Foveavirus (33.3%), Acinetobacter (12.9%), and Klebsiella (11.3%). In contrast, 289 Clusters 3 and 4 were predominantly composed of fungal transcripts. In cluster 3, Monilinia 290 accounting for 10.2% of the relative expression abundance in diseased samples, whereas Cluster 4 291 was mainly composed of Alternaria, which represented 40.5% of the total transcripts. 292 2.4 Dynamic functional profiles of phyllosphere transcriptomes across different 293 stages 294 To investigate the impact of rust infection and disease progression on the functional attributes 295 and activities of the crabapple phyllosphere microbiota, we annotated the transcripts for functional 296 analysis (Figure 4, S3, S4). Of the target unigenes retrieved from the sequence and bioinformatics 297 data, 29.5%, 52.7% and 2.8% were assigned to the KEGG (Kyoto Encyclopedia of Genes and 298 Genomes), GO (Gene Ontology) and CAZy (Carbohydrate-Active enZYmes) databases, 299 respectively. 300 301 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint Figure 4 Functional annotation of phyllosphere transcripts in crabapple leaves based on the Kyoto 302 Encyclopedia of Genes and Genomes (KEGG), Gene Ontology (GO) and Carbohydrate-Active 303 enZYmes (CAZy) databases. (a) and (b), The KEGG pathway and GO functional module 304 enrichment analysis of all DEGs based on the hypergeometric test, respectively. The X-axis 305 represents the enrichment significance of the pathway or module ( P < 0.05, FDR corrected, log2 306 transformed). The circle and triangle sizes represent the number of enriched genes. (c) The relative 307 abundance of significantly differentially expressed carbohydrate-active enzymes (CAZy) across 308 rust disease stages. Differential expression analysis based on generalized linear model (GLM) was 309 used to identity CAZy families showing significant differences at each diseased stage (DESeq2; P310 ＜ 0.05, FDR corrected). 311 After filtering out annotations unrelated to phyllosphere microorganisms, we performed 312 KEGG and GO enrichment analysis on DEGs at each stage and identified significantly enriched 313 pathways (Figure 4a). Overall, the active pathways in phyllosphere microorganisms of crabapple 314 under infection by G. yamadae were primarily involved in primary metabolic processes, with 315 almost all KEGG classifications increasing due to rust disease (Figure 4a, S3). Although the 316 metabolism-related classifications did not exhibit overall significant temporal trends in diseased 317 leaves, some specific metabolic pathways showed significant patterns compared to healthy 318 samples. For example, the pentose and glucuronate interconversions pathway was enriched almost 319 throughout the entire disease progression (Figure 4a). Conversely, transcripts associated with RNA 320 methyltransferase activity, RNA helicase activity, dioxygenase activity and viral process were 321 significantly downregulated during the complete course of rust disease (Figure 4b). Specifically, at 322 the onset of the disease, certain transcripts related to carbohydrate metabolism were notably 323 enriched, such as those involved in ascorbate and aldarate metabolism category (Figure 4a). 324 Simultaneously, pathways involved in the degradation of aromatic compounds were active during 325 the second stage of infection (Figure 4a). Transcripts associated with cellular components 326 affecting gene expression and regulation, such as nucleosome, as well as those related to protein 327 kinase and dimerization activity, were significantly upregulated in the early stages of rust disease 328 (Figure 4b). As lesions expanded, pathways including steroid biosynthesis, galactose metabolism 329 and MAPK signaling were significantly enriched, alongside sustained activity of nucleosome 330 (Figure 4a, 4b). Interestingly, a few categories showed significant upregulation in the late stages of 331 rust disease, such as the aminoacyl-tRNA biosynthesis pathway related to translation, and genes 332 associated with protein kinase activity also exhibited significant enrichment in the latest stage of 333 the disease progression. In addition to the pathways and molecular functions mentioned above, 334 which were mostly downregulated across all diseased stages, certain functions were significantly 335 downregulated only in the late stages of rust disease, including carbohydrate metabolism, lipid 336 metabolism, amino acid metabolism, and apoplastic activity (Figure 4a). Notably, some pathways 337 exhibited temporal patterns, with the ascorbate and aldarate metabolism pathway enriched in the 338 first stage of G. yamadae infection but suppressed in the fifth stage. 339 To gain a more detailed resolution of specific functions associated with metabolism-related 340 pathways, we searched for carbohydrate-active enzymes (CAZy) database. Similar to the KEGG 341 analysis results, the diversity of KEGG orthology (KO) and CAZy families showed similar 342 temporal patterns, and G. yamadae infection also led to an increase expression of all CAZy 343 families (Figure S4a, S4b). After identifying for differentially expressed CAZymes across various 344 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint disease stages, we focused on 11 CAZy families actively involved in regulation (Figure 4c, S4c). 345 Remarkably, certain glycoside hydrolases, potentially related to fungal cell wall component 346 degradation, exhibited increased activity in diseased leaves. These include enzymes involved in 347 the degradation of chitin (GH5), glucans (GH5, GH16 and GH13) and mannans (GH5 and GH31). 348 The relative expression abundance of these glycoside hydrolases (GH5, GH16, and GH13) 349 increased as the lesions expanded and became more enriched in the late stages of rust disease 350 (Figure 4c). Additionally, CAZy families that may interact with the host plant were also enriched 351 in the infected leaves at different stages, including enzymes that participated in the degradation of 352 plant cell wall components such as pectin (PL3 and CE8), cellulose (AA3), lignin (AA3), and 353 xylan (CE5) (Figure 4c, S4c). Among these, the AA3 and CE5 families were upregulated in early 354 and middle stages of rust disease, similar to the patterns of glycoside hydrolases. In contrast, PL3 355 and CE8 families were upregulated in the late stages (Figure S4c). 356 2.5 Dynamic changes in the complexity of phyllosphere microbiome functional 357 co-occurrence networks as lesions expanded 358 359 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint Figure 5 The co-occurrence networks of phyllosphere microbiome functional genes in diseased 360 and healthy samples at different stages of rust disease. (a) Co-occurrence networks of functional 361 genes in diseased and healthy samples at early, middle, and late stages of infection. The nodes in 362 the networks are colored according to functional annotations derived from the KEGG database. 363 Positive correlations between genes are represented by red edges, while negative correlations are 364 indicated by blue edges. (b) Topological properties of the functional gene co-occurrence networks 365 for both diseased and healthy samples across different stages of disease progression. 366 To assess how lesion expansion impacts interactions among functional genes in the 367 phyllosphere, we conducted co-occurrence network analysis and calculated topological properties 368 for both healthy and diseased leaves at various stages (Figure 5). In all groups, functional genes 369 related to microbial metabolism were the most predominant category (Figure 5a). The infection of 370 G. yamadae altered the complexity of the phyllosphere co-occurrence networks, with network 371 complexity indices showing opposite trends between healthy and diseased samples. Specifically, 372 diseased leaves exhibited higher values in the number of nodes, number of edges, and average 373 degree indices compared to healthy leaves, with these indices gradually increased as lesions 374 expanded (Figure 5b). The modularization index, initially higher in diseased during the early 375 stages of infection than in healthy samples, decreased progressively throughout the rust disease 376 progression, eventually falling below the levels seen in healthy leaves by the middle and late 377 stages of infection. In contrast, network density and clustering coefficient indices were higher in 378 healthy leaves at early stages, but elevated in diseased samples as the infection advanced to the 379 middle and late stages. 380 2.6 Contribution of phyllosphere functional genes to the pathogenesis of 381 crabapple rust disease 382 To identify transcripts actively involved in the pathogenesis of crabapple rust disease, we 383 conducted a random forest analysis, focusing on transcripts related to G. yamadae (order 384 Pucciniales) and their contribution to disease progression (Figure 6; Table S4). Out of the 34 385 transcripts identified as significant important for predicting Pucciniales abundance ( P < 0.05), 386 Alternaria alternata emerged as the most predominant species. These key transcripts were 387 primarily involved in several critical biological processes and pathways: the MAPK signaling 388 pathway, aminoacyl-tRNA and steroid biosynthesis, as well as the production of CAZy families 389 that target plant cell wall components. Specifically, enzymes that degrade pectin (PL3), xylan 390 (CE5), and cellulose (AA3) were linked to this functional activity. Additionally, Pseudovirgaria 391 hyperparasitica was found to contribute significantly to the production of pectin-degrading 392 enzymes (CE8), which are involved in breaking down the plant cell wall. Transcripts involved in 393 steroid biosynthesis process were primarily attributed to Mortierella elongata , and these 394 transcripts were identified as the most important contributor to changes in Pucciniales abundance. 395 Other metabolism-related transcripts also showed strong significance for Pucciniales dynamics. 396 For instance, Microbotryum intermedium played a role in the degradation of aromatic compounds 397 and ascorbate and aldarate metabolism, Jaapia argillacea was involved in galactose metabolism, 398 and Termitomyces sp. J132 played a role in pentose and glucuronate interconversions. Interestingly, 399 several transcripts associated with fungal cell wall degradation also showed significant importance 400 in the model, including genes involved in the degradation of glucans (e.g., from Saitozyma 401 podzolica and Alternaria tenuissima ) and mannans (e.g., from Alternaria tenuissima and Jaapia 402 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint argillacea). These findings underscore the significance of microbial interactions and enzymatic 403 activities in contributing to the pathogenesis and progression of crabapple rust disease. 404 405 Figure 6 Importance ranking of phyllosphere microbial functional genes in the severity of 406 crabapple rust disease determined by random forest analysis. The bar plot on the left shows the 407 importance ranking of functional genes, with the X-axis representing the percentage increase of 408 mean square error (MSE), and the y-axis represents functional gene identifiers. The higher the 409 percentage increase in MSE, the more important the gene is in predicting pathogen severity. The 410 dot plot on the right displays the specific functions associated with these genes, where the x-axis 411 represents the functional categories and the y-axis lists the gene identifiers. All functional genes 412 shown in this plot were statistically significant in the random forest model ( P < 0.05). 413 3 Discussion 414 3.1 Rust disease shapes phyllosphere microbiome assembly in crabapple 415 The assembly and composition of the plant microbiome are intricately linked to the plant 416 health, with pathogen infections often resulting in dramatic shifts in microbial composition and 417 functional strategies across different plant-pathogen systems, potentially influencing the host 418 resistance or susceptibility to disease [3]. In our study, we observed a notable increase in species 419 richness across bacteria, fungi, and viruses, which occurred following G. yamadae infection in 420 crabapple leaves, coupled with more active microbial expression patterns in diseased leaves 421 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint (Figure 1a, S1a, S1b). These shifts could reflect a complex interplay, where the plant may recruit 422 antagonistic microbes to counter the pathogen or, alternatively, pathogen-associated microbes 423 could be facilitating disease progression by occupying ecological niches vacated by stressed or 424 damaged host tissues [8]. 425 An intriguing observation in our study was the negative correlation between bacterial and 426 fungal diversity in response to pathogen infection (Figure 1a,1b), suggesting competitive 427 interactions between these two major microbial groups, potentially driven by pathogen-induced 428 shifts in the phyllosphere microbiome. As G. yamadae colonizes the leaf surface and expands 429 lesions, this could create selective pressures that favor the proliferation of certain fungal species, 430 potentially outcompeting the bacterial populations. Such a competitive dynamic would be 431 consistent with the ecological imbalance induced by the rust pathogen, where opportunistic fungi 432 associated with the pathogen gain a competitive edge. Further support for this competitive 433 dynamic is provided by changes in community evenness: as lesions expanded, bacterial species 434 displayed more uniform expression abundance (increased evenness), while the fungal species 435 showed greater divergence in expression abundance (decreased evenness) (Figure S1b). We 436 postulate that this divergence in fungal community evenness may be induced by pathogen 437 invasion, which disrupts the host's defense mechanisms. This allows certain epiphytic fungi, 438 particularly those associated with the pathogen, to rapidly colonize and infect the compromised 439 host tissues, thereby putting bacterial communities at a competitive disadvantage, especially in 440 terms of nutrient acquisition [8]. This hypothesis is reinforced by the observed patterns in the 441 relative expression abundance across the microbial groups in the phyllosphere: bacteria 442 consistently dominated across all stages in healthy leaves, whereas fungi progressively emerged as 443 the predominant microbial group in diseased leaves, with their relative abundance peaking as the 444 lesions expanded (Figure 2a). These findings underscore the increasing importance of fungi—both 445 opportunistic pathogens and symbiotic microbes—in the pathogenesis process. As fungal 446 populations gain prominence in diseased tissues, they may not only exacerbate the disease but also 447 shift the overall microbial landscape, influencing the host's response to infection and potentially 448 altering disease outcomes. 449 3.2 Shifts in microbial community composition and functional groups 450 PCoA results revealed that microbial community composition in the phyllosphere was 451 strongly influenced by leaf condition, with distinct microbial signatures found in diseased leaves 452 versus healthy leaves (Figure 1c, S1c, S1d). Interestingly, during the early and middle stages of 453 infection, bacterial diversity was lower than in healthy leaves at corresponding stages, however, in 454 the late stages of infection, bacterial diversity in the diseased leaves surpassed that in healthy 455 leaves (Figure 1a). This suggests an adaptive restructuring of the microbial community over the 456 course of the infection, where specific taxa may have been displaced or outcompeted initially, but 457 later stages of infection provide new ecological niches, allowing for the gradual re-establishment 458 of bacterial communities. This restructuring may reflect the shifting availability of nutrients or 459 other ecological factors that favor certain bacterial taxa over time [18]. 460 The infection of G. yamadae significantly altered the composition of the phyllosphere 461 microbial community, and the profiles of relative expression abundances revealed notable shifts 462 (Figure 2). While the dominant bacterial taxa at the order level were relatively similar in both 463 healthy and diseased samples, certain transcripts with low cumulative relative abundance 464 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint exhibited significant upregulation or downregulation in diseased leaves (Figure 2b, 2c, S2). This 465 suggests that pathogen infection triggers a distinct microbial response in terms of gene expression, 466 likely reflecting shifts in microbial activity associated with disease progression. In the fungal 467 community, the family Pleosporaceae was the dominant group in both disease and healthy leaves, 468 with most fungal transcripts generally upregulated under pathogen infection, particularly when 469 compared to bacterial community. This suggests a greater involvement of fungi in the host 470 response to pathogen infection, possibly through interactions that ex acerbate or modulate 471 pathogen development [42]. 472 We identified several microbes in both the bacterial and fungal communities with potential 473 antagonistic characteristics (Figure 2, 3). For instance, certain members of the order Nostocales, 474 such as Nostoc calcicola and Nostoc linckia , have been documented as antagonists capable of 475 inhibiting Fusarium oxysporum -induced wilt in tomatoes [43]. Additionally, strains of Klebsiella 476 from the order Enterobacterales have been shown to suppress rust lesions in coffee leaves [44]. The 477 genus Methylobacterium, enriched in the healthy phyllosphere, has been associated with improved 478 plant growth and reduced disease incidence [45]. These antagonistic microbes, by limiting the 479 growth or activity of pathogens, might help mediate plant defense responses and support overall 480 plant health during pathogen exposure. Moreover, we observed groups with potential antagonistic 481 effects within the fungal community. For example, members of Mortierellaceae family can 482 promote plant growth and seed production of Arabidopsis thaliana by mediating the upregulation 483 of plant’s hormone production and activating its defense responses against pathogens [46]. In the 484 present study, the transcripts of these microbes were significantly enriched in diseased samples, 485 suggesting their potential role in enhancing the host's ability to resist G. yamadae infection, where 486 beneficial fungi or bacteria suppress pathogen proliferation and bolster host defenses [8]. 487 In addition to these beneficial microbes, we also observed that fungal communities that 488 gradually dominated the phyllosphere with lesion enlargement consisted primarily of well-known 489 opportunistic pathogens. The upregulated taxa in this category included members of the families 490 Diaporthaceae, Clavicipitaceae, and Mycosphaerellaceae [47–49]. This suggests that these fungal 491 groups may be opportunistically colonizing the highly vulnerable tissues of diseased leaves, 492 exploiting the conditions created by G. yamadae infection to further promote disease spread. The 493 opportunistic nature of th ese pathogens underscores the complexity of plant-microbe interactions 494 in the context of disease, where both beneficial and harmful microbes dynamically interact and 495 influence disease outcomes [8]. 496 3.3 Functional attributes and microbial network complexity 497 Functional analysis of the phyllosphere microbiome revealed G. yamadae infection 498 significantly influenced the functional attributes and activities of microbial communities, with a 499 noticeable shift in the complexity of functional co-occurrence networks compared to those of 500 healthy leaves. Specifically, the functional co-occurrence networks in G. yamadae-infected leaves 501 exhibited more complex patterns, aligning with the increased network complexity observed in our 502 prior co-occurrence network analyses at the taxonomic level (Figure 5) [8]. This shift suggests an 503 adaptive response by the microbial community to pathogen invasion, characterized by the 504 reorganization of microbial interactions. Previous studies have highlighted that pathogen infection 505 often leads to the formation of more intricate microbial networks a part of a defense strategy 506 against the invading pathogen [50,51]. For example, more complex microbial networks have been 507 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint reported in the phyllosphere of Diaporthe citri -infected citrus plants and in both the aboveground 508 and belowground compartments of chili pepper affected by Fusarium wilt disease, compared to 509 their healthy counterparts [20,30]. These findings support the notion that increased network 510 complexity may be a general response mechanism of the phyllosphere microbiome to external 511 stress, aimed at maintaining community stability and enhancing functional redundancy under 512 pathogen pressure. In our study, the increased functional network complexity in G. 513 yamadae-infected samples likely reflects a coordinated response by the phyllosphere microbiome 514 to counteract the pathogen’s effects. We hypothesize that this complexity reflects enhanced 515 interactions among microbial functions that could be protective or antagonistic against the pathogen. 516 Key topological metrics, such as modularity—a measure that often reflects the degree of 517 specialized, compartmentalized interactions within a microbial community-further illustrate this 518 dynamic shift. During the early stages of G. yamadae infection, we observed higher modularity in 519 the functional networks of infected leaves compared to healthy ones (Figure 5b), which suggests 520 that the microbial gene functions in diseased leaves were more stably correlated, potentially 521 reflecting positive co-regulation with distinct microbial subgroups. These compartmentalized 522 responses may be part of a plant defense, where specialized groups of microbes act in concert to 523 reinforce the plant’s defense mechanisms against pathogen invasion. However, as the disease 524 progressed, this stable co-regulation is disrupted, becoming less pronounced compared to healthy 525 leaves, which evolved more stable functional networks over time [51]. 526 3.4 Metabolic adaptations and enzymatic activities 527 The functional enrichment analysis provided critical insights into the mechanisms underlying 528 the observed shifts in the phyllosphere functional co-occurrence networks during G. yamadae 529 infection. For instance, DEGs were substantially enriched in carbohydrate metabolism pathways, 530 including pentose and glucuronate interconversions pathway, ascorbate and aldarate metabolism, 531 and galactose metabolism (Figure 4a). These findings suggest an upregulation of energy 532 metabolism in the phyllosphere microbiome in response to the stress induced by G. yamadae. This 533 increased activity in carbohydrate and energy-related pathways may be essential for microbial 534 survival and proliferation during pathogen attack, potentially supporting the microbiome's overall 535 resilience under stressful conditions. Additionally, the pathogen also modulated microbial gene 536 expression patterns, with nucleosome activity and dimerization activity being prominent in the 537 early diseased stages, followed by the upregulation of the aminoacyl-tRNA biosynthesis pathway 538 in later stages (Figure 4a, 4b). This dynamic regulation highlights the microbiome functional 539 adaptation throughout the rust disease progression, likely aimed at enhancing survival and 540 mitigating the adverse impact of pathogen [35]. 541 Another key finding was the upregulation of degradation pathways involved in breaking 542 down plant-produced antimicrobial substances. Specifically, we observed significant enrichment 543 of genes associated with the degradation of aromatic compounds, such as those in the ascorbate 544 and aldarate metabolism pathways (Figure 4a). This is consistent with our previous study on 545 metabolic profiles of the host plants used in the same experiment setup, further emphasizing the 546 complex interactions within the plant-pathogen-microbiome model [8]. Such interactions are 547 integral to the plant's response to pathogen attack, as plants often recruit beneficial 548 microorganisms with antagonistic functions against pathogens. Rhizoctonia solani -infected beet 549 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint roots are enriched with Chitinophagaceae and Flavobacteriaceae, which secrete fungal cell 550 wall-degrading enzymes to inhibit pathogen infection [21]. Similarly, in our study, the secretion of 551 fungal cell wall degrading enzymes (GH5, GH16, GH13, and GH31) that target cell wall 552 components (chitin, glucan, and mannan) during the early and middle stages of rust disease 553 supports this viewpoint (Figure 4c, S4c). These enzymes target key fungal cell wall components 554 such as chitin, glucan, and mannan, providing a defensive mechanism against pathogen by directly 555 undermining the structural integrity of the invading fungi. These findings underscore the critical 556 role of microbial enzymatic activities in shaping plant-microbe interactions during G. 557 yamadae infection. The secretion of fungal cell wall-degrading enzymes and the overall metabolic 558 reprogramming of the microbiome may be vital strategies by which the phyllosphere microbiota 559 responds to infection, either through direct antagonistic interactions with the pathogen or by 560 supporting the host's immune responses. 561 3.5 Role of beneficial microbes and pathobiome in disease dynamics 562 Random forest model further revealed important microbial taxa and their associated 563 functional activities that significantly contribute to the regulation of microbial functions and 564 pathogen dynamics during G. yamadae infection. Notably, we identified key microbes like Jaapia 565 argillacea and Saitozyma podzolica , whose transcripts were essential in shaping microbial 566 responses and modulating G. yamadae abundance (Figure 6). S. podzolica, a yeast known for its 567 plant growth-promoting (PGP) abilities, has been shown to combat pathogens like Fusarium 568 oxysporum f. sp. Melongenae through the secretion of fungal cell wall-degrading enzymes and 569 antifungal metabolites [52]. Additionally, we observed the enrichment of transcripts from beneficial 570 microorganisms like Mortierella elongata , known for its involvement in steroid biosynthesis and 571 its influence on G. yamadae abundance (Figure 6). M. elongata has been frequently identified as a 572 PGP microorganism in soil and the rhizosphere, enhancing hormone production (e.g., IAA) and 573 bolstering plant resistance to pathogens [53,54]. These recruited beneficial microbes likely play dual 574 roles, both in promoting plant health and in directly inhibiting pathogen proliferation. 575 In contrast to the beneficial microbes, we also detected the presence of certain pathogenic 576 partners, or ‘pathobiomes’, that likely assist G. yamadae in facilitating the progression of disease, 577 highlighting the complexity of the plant-pathogen-microbiome interplay [23]. These 578 microorganisms, typically part of the resident microbiota in healthy leaves, can be exploited by 579 pathogens during infection. This phenomenon highlights how pathogens can "hijack" commensal 580 or symbiotic microorganisms, converting them into allies that aid in overcoming host defenses and 581 establishing infections. One key mechanism involves symbiotic microorganisms capable of 582 producing plant cell wall-degrading enzymes, which can enhance pathogen colonization by 583 breaking down structural components of the host plant tissues. Such activities have been observed 584 in the infections of tomato and tobacco plants, where these enzymes facilitate pathogen entry and 585 spread [25,26]. Our analysis of CAZy indicated a significant enrichment of enzymes like pectate 586 lyase (PL3), carbohydrate esterase (CE8), and acetyl xylan esterase (CE5) (Figure 4c). These 587 enzymes are involved in the degradation of pectin and xylan, key polysaccharides that form part of 588 the plant cell wall. Their upregulation throughout various stages of G. yamadae infection suggests 589 that the pathogen and its associated microbiome are actively degrading plant cell walls, weakening 590 host tissues and facilitating deeper colonization. 591 Despite functional redundancy observed across different microbial taxa, Alternaria alternata 592 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint was identified as a primary contributor to the functional enrichment of microbial activities [55,56]. 593 Remarkably, A. alternata also participated in the MAPK signaling pathway, which is crucial for 594 regulating spore formation, resistance to oxidative and osmotic stress, as well as for its 595 pathogenesis in citrus [57]. The results of GO functional module enrichment analysis further 596 revealed active dynamics of phyllosphere microorganisms during the late stages of disease in the 597 plant apoplast (Figure 4b). The plant apoplast, a key physical barrier against pathogens, play a 598 significant role in triggering defense responses [28]. Thus, we speculate that the dynamic shifts 599 observed in microbial activities during late-stage infection may reflect increased opportunism by 600 certain microbes, exploiting the deteriorating phyllosphere environment. However, the roles of 601 these microorganisms in contributing to or mitigating disease progression require further 602 experimental validation to elucidate their precise functions and interactions within this complex 603 microbial ecosystem. 604 In summary, our study underscores the critical role of the phyllosphere microbiome in 605 mediating plant-pathogen interactions and shaping disease outcomes. These insights provide a 606 foundation for developing microbiome-based strategies to enhance plant health and resilience 607 against pathogens. 608 4 Materials and Methods 609 4.1 Sample collection 610 Crabapple ( Malus ‘Kelsey’) leaves were collected from trees located at the south gate of 611 Olympic Park, Beijing (40°N latitude, 116.38 ° E longtitude) between June and September 2021. 612 The sample process targeted both healthy (non-infected) and diseased ( Gymnosporangium 613 yamadae-infected) leaves from the same trees across six distinct stages of rust disease progression, 614 ranging from the formation of spermogonia to the maturation of aecia. To ensure uniformity and 615 control for genetic variation, three biological replicates were collected for each condition and 616 stage. The leaf sample were immediately transported to the laboratory on dry ice immediately and 617 stored at -80/i3 until RNA extraction. 618 4.2 RNA extraction and metatranscriptomic sequencing 619 To prepare for RNA extraction, the leaf tissues were ground into a fine powder using 620 RNase-free mortars and pestles with liquid nitrogen. Total RNA was extracted from the 621 phyllosphere using the Fecal RNA Extraction Kit (Majorbio, Shanghai, China), following the 622 manufacturer's procedure. The integrity and concentration of the extracted RNA were measured 623 with NanoDrop 2000 spectrophotometer (Thermo Scientific, MA, USA) and an Agilent 5300 624 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). Ribosomal RNA (rRNA) was depleted 625 using the RiboCop rRNA Depletion Kit for Mixed Bacterial Samples (Lexogen, USA), focusing on 626 preserving mRNA for library preparation. For sequencing, 200 ng for each sample’s RNA was 627 used for library preparation with the Illumina® Stranded mRNA Prep, Ligation (Illumina, San 628 Diego, CA, USA). Paired-end sequencing was carried out on the Illumina Novaseq 6000 at 629 Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China). 630 4.3 Bioinformatics and statistical analysis 631 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint 4.3.1 Data preprocessing and quality control 632 Raw sequencing data were processed to ensure high quality and minimize contamination 633 before further analysis. The sequencing data were processed using fastp v0.19.6 [58] for quality 634 control, trimming adapter sequences from both 3’ and 5’ ends of reads. Readers shorter than 50 bp 635 or those with an average base quality score below 20 were discarded to ensure high-quality data, 636 and reads containing ambiguous bases (denoted as ‘N’) were also removed. Potential contaminant 637 reads aligned to the host plant genomes were removed, including Malus domestica 638 (GCF_002114115.1), Malus baccata (GCA_006547085.1), Malus sylvestris (GCF_916048215.2), 639 and Malus sieversii (GCA_020795835.1) using BWA [59]. To further eliminate non-target 640 sequences, rRNA contamination was filtered using SortMeRNA [60], a tool designed to remove 641 ribosomal RNA sequences. After cleaning the data, de novo transcript assembly was performed 642 using Trinity v2.2.0 [61], a widely used tool for assembling RNA-seq data into longer transcript 643 sequences without requiring a reference genome. Transcripts with lengths ≥ 300base pairs were 644 retained for further analysis. The assembled transcripts were then de-replicated using CD-HIT 645 v4.6.1 [62], setting an identity threshold of 0.95 and a minimum coverage of 0.9, which reduced 646 redundancy and ensured that the final set of transcripts represented unique functional genes. The 647 longest sequence was used as representative unigenes for downstream analysis. 648 4.3.2 Taxonomic annotation 649 The assembled unigenes were aligned against the NCBI NR (non-redundant) database using 650 DIAMOND v0.8.35 [63] with an e-value cutoff of 1e-5 to ensure reliable matches. Sequences 651 identified as belong to Viridiplantae (plants), Metazoan (animals) and the pathogen (Pucciniales) 652 were excluded from the analysis, allowing the focus to remain on the targeted microbial 653 communities (Table S5). Transcript abundance was estimated using the RSEM v1.3.2 [64], 654 generating FPKM (fragments per kilobase of transcript per million mapped reads) values for each 655 sample, ensuring accurate comparisons of microbial activity in the phyllosphere microbiome. 656 4.3.3 Diversity analysis 657 Alpha diversity indices, including Richness, Shannon and Pielou’s evenness, as well as beta 658 diversity metrics such as Bray-Curtis dissimilarity, were calculated using the vegan package [65] in 659 R. For statistical comparisons of alpha diversity indices (bacterial, fungal, and viral transcriptomes) 660 between diseased and healthy leaves, one-way ANOVA with Tukey-HSD post hoc test or 661 Kruskal-Wallis with Wilcoxon’s test was employed using the multcomp package [66]. The Shannon 662 index values for bacterial and fungal transcriptomes were correlated using linear regression 663 models (Generalized Linear Models, GLM), implemented with the ggpmisc package [67]. Beta 664 diversity was assessed using Bray-Curtis distance matrices for different stages and leaf infection 665 conditions, which were visualized through principal coordinates analysis (PCoA). Permutational 666 multivariate analysis of variance (PERMANOVA) was applied to test the significant differences 667 in community composition between sample groups using the vegan package [65]. 668 4.3.4 Differential expression and functional analysis 669 The expression matrix of the unigenes was filtered to retain genes expressed in all the 670 consistently expressed across three biological replicates for each treatment. Genes that were 671 expressed at each stage of rust disease was counted, and the intersection of these genes were 672 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint visualized using the UpSetR package [68]. Differentially expressed genes (DEGs), transcripts and 673 carbohydrate-active enzymes (CAZymes) were identified using the DESeq2 package [69]. Pairwise 674 comparisons were made between adjacent sampling stages, applying a generalized linear model 675 (GLM) approach with a significance threshold of P < 0.05 after false discovery rate [FDR] 676 correction. Expression data, transformed to Z-scores and log-transformed FPKM values, were 677 calculated using the c-means clustering algorithm in the TCseq package [70]. For functional 678 annotation, DIAMOND [63] was used to annotate with the KEGG and GO databases, and hmmscan 679 v3.3.2 [71] was used for the CAZy database, with the alignment parameter e-value set to 1e-5 for 680 both. The dynamic abundance patterns of KEGG classifications and CAZy families were presents 681 based on normalized functional geng FPKM values and analyzed using one-way ANOVA test. 682 Functional enrichment analysis of DEGs was performed using the clusterProfiler package [72]. The 683 unigenes obtained from KEGG and GO annotations were used as the background gene set to 684 identify significantly upregulated and downregulated functions associated with DEGs at each time 685 point, with FDR correction applied using the Benjamini-Hochberg method. The Spearman 686 correlation of functional genes was calculated using the Hmisc package [73], and topological 687 properties and the functional co-occurrence network were analyzed and visualized using Gephi [74] 688 as disease lesions expanded. To identify functional genes most predictive of pathogen abundance, 689 a random forest model we constructed used the rfPermute package [75]. The relative abundance of 690 Pucciniales was set as the response variable, with the expression levels of other functional gens 691 serving as predictor variables. This model was employed to pinpoint microbial functional 692 activities potentially linked to the pathogen’s colonization success. 693 Acknowledgments 694 This study was supported by the National Natural Science Foundation of China (grant 32101527). 695 Author Contributions 696 Q.X. conducted the experiments, analyzed the data, and edited the manuscript. Y.Z. conducted the 697 experiments and edited the manuscript. S.T. designed the experiments, provided the experimental 698 conditions, and contributed to the editing and review of the manuscript. All authors have read and 699 approved the final version of the manuscript for publication. 700 Data Availability Statement 701 The sequencing raw data have been uploaded to the National Genomics Data Center. 702 Conflict of Interest Statement 703 The authors declare no conflict of interest. 704 705 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint

706 1. Vandenkoornhuyse P, Quaiser A, Duhamel M, Le Van A, Dufresne A. The importance of the 707 microbiome of the plant holobiont. New Phytologist 2015;206(4):1196–206. 708 2. Simon JC, Marchesi JR, Mougel C, Selosse MA. Host-microbiota interactions: from holobiont 709 theory to analysis. Microbiome 2019;7(1):5. 710 3. Trivedi P, Leach JE, Tringe SG, Sa T, Singh BK. Plant–microbiome interactions: from 711 community assembly to plant health. Nat Rev Microbiol 2020;18(11):607–21. 712 4. Berg G, Rybakova D, Fischer D, Cernava T, Vergès MCC, Charles T, et al. Microbiome definition 713 re-visited: old concepts and new challenges. Microbiome 2020;8(1):103. 714 5. Chialva M, Lanfranco L, Bonfante P. The plant microbiota: composition, functions, and 715 engineering. Current Opinion in Biotechnology 2022;73:135–42. 716 6. Vayssier-Taussat M, Albina E, Citti C, Cosson JF, Jacques MA, Lebrun MH, et al. Shifting the 717 paradigm from pathogens to pathobiome: new concepts in the light of meta-omics. Front Cell 718 Infect Microbiol 2014;4. 719 7. Gu S, Wei Z, Shao Z, Friman VP, Cao K, Yang T, et al. Competition for iron drives phytopathogen 720 control by natural rhizosphere microbiomes. Nat Microbiol 2020;5(8):1002–10. 721 8. Zhang Y, Cao B, Pan Y , Tao S, Zhang N. Metabolite-Mediated Responses of Phyllosphere 722 Microbiota to Rust Infection in Two Malus Species. Microbiology Spectrum 723 2023;11(2):e03831-22. 724 9. Zhou X, Zhang J, Khashi u Rahman M, Gao D, Wei Z, Wu F, et al. Interspecific plant interaction 725 via root exudates structures the disease suppressiveness of rhizosphere microbiomes. Molecular 726 Plant 2023;16(5):849–64. 727 10. Pieterse CMJ, Zamioudis C, Berendsen RL, Weller DM, Van Wees SCM, Bakker PAHM. 728 Induced systemic resistance by beneficial microbes. Annu Rev Phytopathol 2014;52:347–75. 729 11. Hacquard S, Spaepen S, Garrido-Oter R, Schulze-Lefert P. Interplay Between Innate Immunity 730 and the Plant Microbiota. Annu Rev Phytopathol 2017;55:565–89. 731 12. Noman M, Ahmed T, Ijaz U, Shahid M, Azizullah, Li D, et al. Plant–Microbiome Crosstalk: 732 Dawning from Composition and Assembly of Microbial Community to Improvement of Disease 733 Resilience in Plants. International Journal of Molecular Sciences 2 021;22(13):6852. 734 13. Pereira LB, Thomazella DPT, Teixeira PJPL. Plant-microbiome crosstalk and disease 735 development. Current Opinion in Plant Biology 2023;72:102351. 736 14. Kwak MJ, Kong HG, Choi K, Kwon SK, Song JY , Lee J, et al. Rhizosphere microbiome 737 structure alters to enable wilt resistance in tomato. Nat Biotechnol 2018;36(11):1100–9. 738 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint 15. Berendsen RL, Vismans G, Yu K, Song Y, De Jonge R, Burgman WP, et al. Disease-induced 739 assemblage of a plant-beneficial bacterial consortium. The ISME Journal 2018;12(6):1496–507. 740 16. Liu H, Li J, Carvalhais LC, Percy CD, Prakash Verma J, Schenk PM, et al. Evidence for the 741 plant recruitment of beneficial microbes to suppr ess soil-borne pathogens. New Phytologist 742 2021;229(5):2873–85. 743 17. Yuan J, Zhao J, Wen T, Zhao M, Li R, Goossens P, et al. Root exudates drive the soil-borne 744 legacy of aboveground pathogen infection. Microbiome 2018;6(1):156. 745 18. Rolfe SA, Griffiths J, Ton J. Crying out for help with root exudates: adaptive mechanisms by 746 which stressed plants assemble health-promoting soil microbiomes. Current Opinion in 747 Microbiology 2019;49:73–82. 748 19. Mannaa M, Seo YS. Plants under the Attack of Allies: Moving towards the Plant Pathobiome 749 Paradigm. Plants 2021;10(1):125. 750 20. Gao M, Xiong C, Gao C, Tsui CKM, Wang MM, Zhou X, et al. Disease-induced changes in 751 plant microbiome assembly and functional adaptation. Microbiome 2021;9(1):187. 752 21. Carrión VJ, Perez-Jaramillo J, Cordovez V, Tracanna V, de Hollander M, Ruiz-Buck D, et al. 753 Pathogen-induced activation of disease-suppressive functions in the endophytic root microbiome. 754 Science 2019;366(6465):606–12. 755 22. Bass D, Stentiford GD, Wang HC, Koskella B, Tyler CR. The Pathobiome in Animal and Plant 756 Diseases. Trends in Ecology & Evolution 2019;34(11):996–1008. 757 2 3 . L v T , Z h a n C , P a n Q , X u H , F a n g H , W a n g M , e t a l . P l a n t p a t h o g e n e s i s : T o w a r d 758 multidimensional understanding of the microbiome. iMeta 2023;2(3):e129. 759 24. Yan ZZ, Wang J, Liang J, Batista BD, Liu H, Xiong C, et al. Evidence of distinct response of soil 760 viral community to a plant infection and the disease pathobiome. Journal of Sustainable 761 Agriculture and Environment 2023;2(4):382–7. 762 25. Lu P, Shi H, Tao J, Jin J, Wang S, Zheng Q, et al. Metagenomic insights into the changes in the 763 rhizosphere microbial community caused by the root-knot nematode Meloidogyne incognita in 764 tobacco. Environmental Research 2023;216:114848. 765 26. Tian BY, Cao Y, Zhang KQ. Metagenomic insights into communities, functions of endophytes 766 and their associates with infection by root-knot nematode, Meloidogyne incognita , in tomato 767 roots. Sci Rep 2015;5(1):17087. 768 27. Malacrinò A, Abdelfattah A, Berg G, Benitez MS, Bennett AE, Böttner L, et al. Exploring 769 microbiomes for plant disease management. Biological Control 2022;169:104890. 770 28. V orholt JA. Microbial life in the phyllosphere. Nat Rev Microbiol 2012;10(12):828–40. 771 29. Chen T, Nomura K, Wang X, Sohrabi R, Xu J, Yao L, et al. A plant genetic network for 772 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint preventing dysbiosis in the phyllosphere. Nature 2020;580(7805):6 53–7. 773 30. Li PD, Zhu ZR, Zhang Y , Xu J, Wang H, Wang Z, et al. The phyllosphere microbiome shifts 774 toward combating melanose pathogen. Microbiome 2022;10(1):56. 775 31. Liu H, Brettell LE, Singh B. Linking the Phyllosphere Microbiome to Plant Health. Trends in 776 Plant Science 2020;25(9):841–4. 777 32. Díaz-Cruz GA, Cassone BJ. Changes in the phyllosphere and rhizosphere microbial 778 communities of soybean in the presence of pathogens. FEMS Microbiology Ecology 779 2022;98(3):fiac022. 780 33. Seybold H, Demetrowitsch TJ, Hassani MA, Szymczak S, Reim E, Haueisen J, et al. A fungal 781 pathogen induces systemic susceptibility and systemic shifts in wheat metabolome and 782 microbiome composition. Nat Commun 2020;11(1):1910. 783 34. Tao SQ, Auer L, Morin E, Liang YM, Duplessis S. Transcriptome Analysis of Apple Leaves 784 Infected by the Rust Fungus Gymnosporangium yamadae at Two Sporulation Stages. Molecular 785 Plant-Microbe Interactions 2020;33:444–461. 786 35. Peng D, Wang Z, Tian J, Wang W, Guo S, Dai X, et al. Phyllosphere bacterial community 787 dynamics in response to bacterial wildfire disease: succession and interaction patterns. Front 788 Plant Sci 2024;15. 789 36. Lemanceau P, Blouin M, Muller D, Moënne-Loccoz Y . Let the Core Microbiota Be Functional. 790 Trends in Plant Science 2017;22(7):583–95. 791 37. Guo L, Cao Y , Quan J, Liu BY. Crabapple in China: past, present and future. Acta Horticulturae 792 2019;1263:55–60. 793 38. Kern FD. A revised taxonomic account of Gymnosporangium. Pennsylvania State University 794 Press, State College. PA 1973. 795 39. Cummins G, Hiratsuka Y. Illustrated genera of rust fungi. APS Press, St Paul, MN 2003. 796 40. Yue C, Du C, Wang X, Tan Y, Liu X, Fan H. Powdery mildew-induced changes in phyllosphere 797 microbial community dynamics of cucumber. FEMS Microbiology Ecology 798 2024;100(5):fiae050. 799 41. Li Y , Tao S, Liang Y. Time-Course Responses of Apple Leaf En dophytes to the Infection of 800 Gymnosporangium yamadae. Journal of Fungi 2024;10(2):128. 801 42. Collinge DB, Jensen B, Jørgensen HJ. Fungal en dophytes in plants and their relationship to 802 plant disease. Current Opinion in Microbiology 2022;69:102177. 803 43. El-Sheekh MM, Deyab MA, Hasan RSA, Abu Ahmed SE, Elsadany AY. Biological control of 804 Fusarium tomato-wilt disease by cyanobacteria Nostoc spp. Arch Microbiol 2022;204(1):116. 805 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint 44. Shiomi HF, Silva HSA, Melo IS de, Nunes FV, Bettiol W. Bioprospecting endophytic bacteria 806 for biological control of coffee leaf rust. Sci agric (Piracicaba, Braz) 2006;63:32–9. 807 45. Murugaiyan S, Ramasamy K. Isolation and Characterization of Tomato Leaf Phyllosphere 808 Methylobacterium and Their Effect on Plant Growth. International Journal of Current 809 Microbiology and Applied Sciences 2017;6:2121–36. 810 46. Vandepol N, Liber J, Yocca A, Matlock J, Edger P, Bonito G. Linnemannia elongata 811 (Mortierellaceae) stimulates Arabidopsis thaliana aerial growth and responses to auxin, ethylene, 812 and reactive oxygen species. PLOS ONE 2022;17(4):e0261908. 813 47. Chen SF, Morgan DP, Michailides TJ. Botryosphaeriaceae and Diaporthaceae associated with 814 panicle and shoot blight of pistachio in California, USA. Fungal Diversity 2014;67(1):157–79. 815 48. White J, Sullivan R, Moy M, Patel R, Duncan R. An overview of problems in the classification 816 of Plant-parasitic Clavicipataceae. Studies in Mycology 2000;45:95–105. 817 49. Thomma BPHJ, Van Esse HP, Crous PW, De Wit PJGM. Cladosporium fulvum (syn. Passalora 818 fulva), a highly specialized plant pathogen as a model for functional studies on plant pathogenic 819 Mycosphaerellaceae. Molecular Plant Pathology 2005;6(4):379–93. 820 50. Faust K, Raes J. Microbial interactions: from networks to models. Nat Rev Microbiol 821 2012;10(8):538–50. 822 51. Wei Z, Gu Y , Friman VP, Kowalchuk GA, Xu Y, Shen Q, et al. Initial soil microbiome 823 composition and functioning predetermine future plant health. Science Advances 824 2019;5(9):eaaw0759. 825 52. Das S, Rabha J, Narzary D. Assessment of soil yeasts Papiliotrema laurentii S-08 and 826 Saitozyma podzolica S-77 for plant growth promotion and biocontrol of Fusarium wilt of brinjal. 827 Journal of Applied Microbiology 2023;134(11):lxad252. 828 53. Li F, Chen L, Redmile-Gordon M, Zhang J, Zhang C, Ning Q, et al. Mortierella Elongata’s 829 Roles in Organic Agriculture and Crop Growth Promotion in a Mineral Soil. Land Degradation & 830 Development 2018;29. 831 54. Ozimek E, Hanaka A. Mortierella Species as the Plant Growth-Promoting Fungi Present in the 832 Agricultural Soils. Agriculture 2021;11(1):7. 833 55. Louca S, Polz MF, Mazel F, Albright MBN, Huber JA, O’Connor MI, et al. Function and 834 functional redundancy in microbial systems. Nat Ecol Evol 2018;2(6):93 6–43. 835 56. Paasch BC, He SY . Toward understanding microbiota homeostasis in the plant kingdom. PLOS 836 Pathogens 2021;17(4):e1009472. 837 57. Chung KR. Mitogen-activated protein kinase signaling pathways of the tangerine pathotype of 838 Alternaria alternata. MAP Kinase 2013;2:e4. 839 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint 58. Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. 840 Bioinformatics 2018;34(17):i884–90. 841 59. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. 842 Bioinformatics 2009;25(14):1754–60. 843 60. Kopylova E, Noé L, Touzet H. SortMeRNA: fast and accurate filtering of ribosomal RNAs in 844 metatranscriptomic data. Bioinformatics 2012;28(24):3211–7. 845 61. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo 846 transcript sequence reconstruction from RNA-seq using the Trinity platform for reference 847 generation and analysis. Nat Protoc 2013;8(8):1494–512. 848 62. Fu L, Niu B, Zhu Z, Wu S, Li W. CD-HIT: accelerated for clustering the next-generation 849 sequencing data. Bioinformatics 2012;28(23):3150–2. 850 63. Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat 851

2015;12(1):59–60. 852 64. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without 853 a reference genome. BMC Bioinformatics 2011;12:323. 854 65. Oksanen J, Kindt R, Legendre P, O’Hara B, Stevens MHH, Oksanen MJ, Suggests M. The 855 vegan package. Community Ecology Package 2007;10: 631–637. 856 66. Hothorn T, Bretz F, Westfall P. Simultaneous Inference in General Parametric Models. 857 Biometrical Journal 2008;50(3):346–63. 858 67. Aphalo PJ, Slowikowski K, Mouksassi S. ggpmisc: Miscellaneous Extensions to ‘ggplot2’ 859 2024. 860 68. Gehlenborg N. UpSetR: A More Scalable Alternative to Venn and Euler Diagrams for 861 Visualizing Intersecting Sets. 2015. 862 69. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq 863 data with DESeq2. Genome Biol 2014;15(12):550. 864 70. Gu L. TCseq: time course sequencing data analysis. 2023. 865 71. Potter SC, Luciani A, Eddy SR, Park Y, Lopez R, Finn RD. HMMER web server: 2018 update. 866 Nucleic Acids Research 2018;46(W1):W200–4. 867 72. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool 868 for interpreting omics data. Innovation (Camb) 2021;2(3):100141. 869 73. Harrell Jr. Hmisc: Harrell Miscellaneous version 5.1-3 from CRAN. 2024 870 74. Bastian M, Heymann S, Jacomy M. Gephi: An Open Source Software for Exploring and 871 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint Manipulating Networks. Proceedings of the International AAAI Conference on Web and Social 872 Media 2009;3(1):361–2. 873 75. Archer E. rfPermute: Estimate Permutation p-Values for Random Forest Importance Metrics. 874 2011. 875 876 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 11, 2025. ; https://doi.org/10.1101/2025.02.11.637274doi: bioRxiv preprint