The endophytic fungus NQ8GII4 was isolated by the Laboratory of Plant Pathology, School of Agriculture, Ningxia University, and deposited at China General Microbiological Culture Collection Center (CGMCC No. 19271). The strain was maintained on potato dextrose agar (PDA) medium and lyophilized filter paper for short and long term storage. For subsequent assays, this strain was growing on a PDA medium for 5 d at 25°C. The Agrobacterium tumefaciens strain GV3101 (hereafter GV3101) was used for agroinfiltration of plants, cultured on lysogeny broth medium at 28°C.

Gannon No. 3 is a variety of alfalfa, widely cultivated in Gansu, Ningxia, Qinghai, and Xinjiang Province in western China (Zhang, 2018), and was purchased from Ningxia Shanggu Agricultural and Animal Husbandry Development Co. (Ningxia, China). Nicotiana benthamiana seedlings were grown in a glasshouse at 14 h 22°C: 10 h 20°C, day: night, 72% relative humidity.

Three 5 mm NQ8GII4 mycelial disks were placed in Czapek medium and incubated on a rotary shaker at 175 rpm at 25°C for 4 days. The genomic DNA of NQ8GII4 was extracted from the mycelia using the Omage Fungal DNA Kit. Purified genomic DNA was quantified by TBS-380 fluorometer (Turner Biosystems Inc., Sunnyvale, CA). Short insert library (400 bp) was constructed and sequenced on an Illumina HiSeq X Ten Sequencing System (Shanghai Majorbio Bio-Pharm Technology Co., Ltd.), generating 150 bp paired-end reads. In order to make the subsequent assembly more accurate, low quality bases and sequencing adapter sequences were trimmed and filtered from the raw illumina reads using Trimmomatic v. 0.36 (Bolger et al., 2014). SOAPdenovo2 (Luo et al., 2012) was tested to assemble the NQ8GII4 genome. The gene set integrity was then evaluated with BUSCO v5.4.5 (Simão et al., 2015) and CEGMA v2.5 (Parra et al., 2007). The assembled genome sequence was deposited into the GenBank database under accession No. JAVIXY000000000.1.

The annotation of the genome of NQ8GII4 was performed with MAKER2 v2.31 (Holt and Yandell, 2011), Snap, Augustus v2.5.5,¹ and Evidence Modeler.² The protein sequences of the predicted genes were compared with NCBI Non-Redundant (NR), Kyoto Encyclopedia of Genes and Genomes (KEGG), Swiss-Prot protein, Pfam, Cluster of Orthologous Genes (COG), and Gene Ontology (GO) databases by BLASTp (BLAST+2.9.0 with E-value cut-off <1e−5). The repeats were identified by RepeatMasker v4.1.4.³ barrnap V0.9⁴ and tRNAscan-SE v2.0 (Lowe and Eddy, 1997) were used to predict the rRNA and tRNA contained in the genome.

OrthoVenn 3 (Sun et al., 2023) was used to analyze the orthologs pairs of 5 endophytes (NQ8GII4, E. festucae, S. indica, P. fici, Fo47), 5 pathogens (V. dahliae, P. oryzae, F. solani, Fol4287, C. higginsianum), and 2 saprophytes (A. nidulans, N. crassa) (Supplementary Table 1). The tools were executed with the E-value cutoff of 1e−5 for all-to-all protein similarity comparisons, and an inflation value of 1.5 for the generation of orthologous cluster using the Markov Cluster Algorithm. For evolutionary divergence data estimation, clustering, protein family selection and phylogenetic analysis were performed with OrthoVenn 3.

Three 5 mm NQ8GII4 mycelial disks were placed in Czapek medium and incubated on a rotary shaker at 175 rpm at 25°C for 4 days, the fermentation broth of containing 10⁷ CFU/mL of NQ8GII4 growth was obtained. To avoid the effect of microorganisms other than NQ8GII4 on gene expression in alfalfa. Seeds of healthy Gannon No.3 were cleaned in sterile water, rinsed with 75% ethanol for 30 s, and sterilized with 10% NaClO for 1 h, and washed four times with sterile water. Excess water was removed using sterile filter paper, and then the seeds were added to 70 mL of Hoagland’s nutrient solution in a 500 mL bottle. The seeds were germinated inside the bottles under the following conditions: 22°C with a light intensity of 4000 lx and a photoperiod of 16 h light and 8 h dark. After 25 days of incubation, the roots of the resulting seedlings were inoculated with 15 mL of fermentation broth containing 10⁷ CFU/mL of NQ8GII4 growth. All operations were carried out under sterile conditions, and the plants were grown in an incubator with 20 replicates.

According to the infestation process of NQ8GII4 in the alfalfa root (Supplementary Figure 1), alfalfa seedlings were harvested at 0.5, 1, 6, and 14 days post inoculation (dpi), three biological replicates were collected at each treatment, and fifteen whole plants were pooled as a single replicate. The obtained plants were rapidly frozen with liquid nitrogen and maintained at −80°C until RNA was extracted. A1, A2, and A3 are alfalfa samples of NQ8GII4 inoculation at 12 hpi, B1, B2, and B3 are alfalfa samples of NQ8GII4 inoculation at 1 dpi, C1, C2, and C3 are alfalfa samples of NQ8GII4 inoculation at 6 dpi, D1, D2, and D3 are alfalfa samples of NQ8GII4 inoculation at 14 dpi, E1, E2, and E3 are mycelia of NQ8GII4.

Samples of A1 ∼ E3 total RNA were extracted using E.Z.N.A.^® Total RNA Isolation Kit (Omega, USA) according to the manufacturer’s recommendations. After integrity, quality, and purity checking, the purification of messenger RNA (mRNA), construction of complementary DNA (cDNA) libraries, cDNA end repair, adapter ligation, and cDNA amplification were performed by Majorbio Technologies Co., Ltd. (Shanghai, China). The final quantified libraries were sequenced on an Illumina Hiseq 2000 platform. RNA-Seq data generated in this study were mapped to the NQ8GII4 genome using HISAT2 (Kim et al., 2015), and transcriptome-based gene structures were obtained by String Tie (Pertea et al., 2015). RNA-Seq data with accession PRJNA1125526 were deposited in the NCBI databases.⁵

FPKM (fragments per kilobase of exon model per million mapped reads) (Trapnell et al., 2010) was used to standardize the expression levels. With the aid of the DESeq program (Love et al., 2014), differential expression analysis of A/B/C/D-vs-E was carried out. The cut-off value of P-value < 0.01 and at least 2 fold change (FC) was set as threshold for differential expression.

The identification and annotation of carbohydrate active enzymes (CAZymes) were performed on the dbCAN3 meta server (E-value < 1e−15, coverage > 0.35) (Zheng et al., 2023). Transporters, secreted peptidases, and G-protein-coupled receptors (GPCR) were annotated by a BLAST search of the Transporter Classification Database (TCDB) v2020091 (Saier et al., 2021), MEROPS (Rawlings et al., 2014), and GPCRdb (Pándy-Szekeres et al., 2023) using a BLASTp E-value cut-off of 1E−5. GPCR characterized by a seven transmembrane domain structure, and GPCR-like proteins were evaluated to verify the presence of seven transmembrane helices using TMHMM 2.0 (Krogh et al., 2001). GPCR-like proteins functional annotation was performed by blasting with the Pathogen Host Interaction (PHI) database (Urban et al., 2022). Cytochrome P450s were analyzed on the online tool at Majorbio Cloud Platform.⁶ SignalP V4.1 (Petersen et al., 2011), TMHMM 2.0, Big-PI Fungal Predictor (Eisenhaber et al., 2004), ProtComp V9.0 (Emanuelsson et al., 2000), and EffectorP v3.0 (Sperschneider and Dodds, 2022) were used to distinguish effector proteins in the NQ8GII4 genome. Gene clusters related to secondary metabolite biosynthesis were identified using the antibiotic and secondary metabolite analysis shell (antiSMASH V7.0) (Blin et al., 2021) online tool. The parameter setting was maintained as the default parameter values.

Quantitative, real-time RT-PCR (qRT-PCR) assays were conducted with the same total RNA samples that were used for RNA-seq studies. We chose 10 DEGs at random for qRT-PCR investigation to confirm the accuracy of the RNA-Seq data. TUB was selected as an internal standard to calculate the relative expression levels. This analysis was conducted using the qTOWER 2.2 system (Jena, Germany) with TB Green Premix Ex Taq II FAST qPCR mix (Takara, Japan). qRT-PCR was carried out with 12.5 μL of 2 × TB Green Premix Ex Taq II FAST qPCR mix, 1.0 μL of sense and antisense primers (2.5 μM), 1.0 μL of cDNA in a final volume of 25 μL. The cycling program was set as follows: 94°C for 10 min, followed by 40 cycles of 94°C for 5 s, and 60°C for 10 s. Each plate was repeated thrice in independent runs for all reference and selected genes. Data were analyzed using qPCRsoft 4.1 quantitative PCR software. All primers are listed in the Supplementary Table 2. Primers were synthesized by Sangon Biotech Co., Ltd. (Shanghai, China). The 2^–ΔΔCT method was used for relative quantification of gene expression and transformed into log₂FC.

Coding regions of 25 genes encoding putative effectors with signal peptides (SPs) were amplified from a NQ8GII4 complementary DNA (cDNA) library with gene-specific primers (Supplementary Table 2) using Hieff Canace^® Plus High-Fidelity DNA Polymerase (Yeasen, China). The amplicons were subsequently ligated with potato virus X (PVX) vectors (ClaI and SalI for PVX vector). All constructs were validated by sequencing in Sanger (Sanger Biotech, Shanghai, China).

The constructs were transformed into GV3101 using hot shock method. After selection with selective antibiotics, individual colonies validated by PCR were cultured in lysogeny broth medium at 28°C in a shaking incubator at 220 rpm for 48 h. The bacteria were then pelleted by centrifugation and suspended in MES buffer (10 mM magnesium chloride (MgCl₂), 10 mM MES, 200 μM acetosyringone, pH 5.7) in the dark for 3 h at room temperature (RT) before infiltration. For infiltration, suspended GV3101 cells were adjusted to a final OD600 of 0.6. GV3101 cell suspension was infiltrated into plant leaves using a syringe without a needle. To determine cell-death-inducing activity of the proteins, PVX constructs were agroinfiltrated into plants. Symptom development was monitored visually 5 d post agroinfiltration (dpa). GV3101 transformed with BAX1, and INF1 were used as a positive control. GV3101 transformed with PVX was used as negative control. In addition, the bacterial suspension were infiltrated into N. benthamiana 12 h before infiltrating GV3101 bearing the INF1 or BAX1-expressing construct.

A Maximum Likelihood tree was generated with amino acid sequences of A, KS, and effector with the program MEGA v.11 using the Jones Taylor Thornton model (JTT), 100 replicates for bootstrap analysis. A volcano plot was drawn by Hiplot.⁷ Heat map, stacked bar chart, circular bar plot and phylogenomic tree were drawn by Chiplot⁸ (Xie et al., 2023; Li J. et al., 2023; Ji et al., 2022). The venn diagram was drawn by TBtools (Chen et al., 2023).

The genome of NQ8GII4 was assembled from the data generated by the Illumina sequencing platform. After quality control of raw data, 4,997.1 Mb data were generated. The assembled NQ8GII4 genome consisted of 2,064 scaffolds with N50 length of 144.83 Kb and a combined size of 50,827,403 bp (48.47 Mb) (Supplementary Table 3). The average GC ratio was 53.93%. A total of 15,149 protein-coding genes were predicted with a length of 28,820,034 bp. The average gene length was 1902.44 bp, and the gene density was 307 genes per 1 Mb of the genome (Figure 1). These predicted genes account for 56.70% of the genome. About 14,279 (94.25%), 10,276 (67.83%), 9,605 (63.40%), 9,428 (62.23%), 5,665 (37.40%), and 4,404 (29.07%) of the predicted genes presented homologies with known functions in the NR, Pfam, Swiss-Model, GO, COG, and KEGG databases, respectively. In addition, 233 tRNA and 50 rRNA were predicted in the assembly genome.

The repeat elements was 5.36% (2,727,425 bp) in the NQ8GII4 genome, including 4.72% (2,400,760 bp) as interspersed repeats and 0.64% (326,665 bp) as tandem repeats. The most abundant repetitive elements was the DNA transposons (1.42%, 720,961 bp), followed by long terminal repeat (LTR) (0.27%, 138,531 bp), non-LTR retrotransposons LINEs (long interspersed nuclear elements) (0.39%, 199,653 bp), rolling circle (RC) (0.72%, 367,651 bp). Among the tandem repeats, a total 7,526 (321,471 bp) microsatellite were identified, accounting for 0.63% of the genome. In addition, 41 (0.01%, 5,194 bp) satellite were identified. 4,414 disrupted genes were identified in NQ8GII4 genome, suggesting that the repeat elements have had an important influence on the dynamics of the evolution of the NQ8GII4 genome (Supplementary Table 4).

The phylogeny was calibrated using three calibration points. A. nidulans and N. crassa diverged approximately 300 ∼ 400 MYA, P. oryzae and C. higginsianum diverged approximately 207 ∼ 339 MYA, V. dahliae and C. higginsianum diverged approximately 100 ∼ 150 MYA (Hacquard et al., 2016; O’Connell et al., 2012). NQ8GII4 and F. oxysporum diverged approximately 60 ∼ 70 MYA, NQ8GII4 and F. solani diverged relatively recently (45 ∼50 MYA), suggested that NQ8GII4 is closely related to F. solani (Geiser et al., 2021). It also suggested that each lifestyle has evolved independently of the others (Figure 2A).

An orthovenn 3.0 analysis identified 17,618 ortholog groups (127,204 genes) that were clustered between NQ8GII4 and the other 11 fungal genomes [comprising phytopathogens (V. dahliae, P. oryzae, F. solani, Fol4287, and C. higginsianum), endophytes (E. festucae, S. indica, P. fici, and Fo47), and saprophytes (A. nidulans, and N. crassa)]. On average, each ortholog group contained approximately 7.22 genes. No clusters were identified that were shared by all symbiotic fungi, but were absent from the genomes of pathogens and saprophytes (Figures 2B, C).

To understand transcriptional reprogramming during the mutualistic symbiosis of NQ8GII4 and alfalfa, gene profiles of NQ8GII4 were studied during infection at 0.5, 1, 6, and 14 dpi. In total, 2,091 genes were identified as DEGs. Among these genes, 100 (7 upregulated and 93 downregulated), 143 (75 upregulated and 68 downregulated), 597 (56 upregulated and 541 downregulated), and 1,714 (661 upregulated and 1,053 downregulated) were identified in 0.5, 1, 6, and 14 dpi. Groups of 8, 57, 282, and 1,413 were identified as unique to 0.5, 1, 6, and 14 dpi. A total of 52 DEGs were commonly expressed at 0.5, 1, 6, and 14 dpi (Supplementary Figure 4A).

To validate the reliability of the RNA-Seq results, we randomly selected 10 DEGs for qRT-PCR analysis. The trends in expression levels of the selected genes were similar to those found in the RNA-Seq, indicating the reliability of the DEGs analyzed by RNA-Seq (Supplementary Figure 4B).

The gene profile of NQ8GII4 at 0.5 dpi was similar to 1 dpi. Downregulated DEGs were mainly enriched in GTs, FCWDEs, S08A, steroid-14α-demethylase (CYP51), and the sugar porter family. One DHA1 and one DHA2 were upregulated at 1 dpi. Gene09380, annotated as GH12 protein, a putative cellulase (Wang et al., 2018), was upregulated at 1 dpi, suggesting that NQ8GII4 begins infecting by degrading the alfalfa cell wall. At 6 dpi, upregulated DEGs were enriched in plant cell wall-degrading enzymes (PCWDEs), such as enzymes degrading cellulose, hemicellulose, and pectin. At 14 dpi, upregulated DEGs were enriched in the group of effectors, and CYPs BM-3 (CYP102). Two different PTH11-like GPCRs were significantly upregulated at 6 and 14 dpi, respectively. These findings suggest that NQ8GII4 respond to different environment signals at 6 and 14 dpi. The core biosynthetic gene of secondary metabolite cluster 9 (NRPS) were downregulated at 0.5, 1, 6, and 14 dpi, and the core biosynthetic gene of secondary metabolite cluster 4 (T1PKS) was significantly downregulated at 14 dpi. 3, 12, and 31 putative effector genes were upregulated at 1, 6, and 14 dpi, respectively. Of these putative effector genes, 1, 5, and 11 genes upregulated at 1, 6 and 14 dpi, respectively, were also annotated as CAZyme-encoding genes. Futhermore, autophygy related genes, such as gene08729 (homolog of ATG1), gene00631 (homolog of ATG2), gene07294 (homolog of ATG11), and others, were downregulated at 6 and 14 dpi (Figure 11).

We carried out transient expression of 25 full-length effectors in N. benthamiana. The positive control using vector containing the BAX or INF gene, and PVX vector was used as the negative control. It was shown that, FnEG1 (Gene00104), FnEG3 (gene00314), FnEG7 (gene00635), FnEG8 (gene00639), FnEG168 (gene10147), FnEG214 (gene11857), and FnEG218 (gene11959) could trigger cell death in N. benthamiana. By contrast, PVX vector did not induce plant cell death, suggesting the cell death was significantly activated by these effectors. BAX and INF strongly induced plant cell death at 5 d post inoculation (dpi), whereas FnEG100 (gene07554) could fully suppress INF or BAX-triggered plant cell death (Figure 12). FnEG1, FnEG7, FnEG8, FnEG100, FnEG168, FnEG214, and FnEG218 contained AA9, CFEM, Cerato-platanin, SUN, AltA, Pectata lyase, and GH16 domain, respectively (Supplementary Figure 5). In RNA-seq analysis of NQ8GII4-infected alfalfa roots, FnEG168 was upregulated at 1 and 14 dpi, while FnEG1, FnEG8, and FnEG214 were upregulated at 14 dpi.

Molecular phylogeny, whole-genome alignment, and divergence date estimates indicate that 5 endophytes examined evolved independently from each other. We then compared gene families — such as CAZymes, transporters, CYPs, GPCRs, secreted peptidases, effectors, and secondary metabolite gene clusters—among the 12 fungi studied, we found different endophytes possessed different structures of gene family. The gene family profiles of E. festucae and S. indica were similar to those of saprophytes. Epicholoë could act as saprophytes to enhance nutrient cycling in grassland ecosystems (Song et al., 2021). S. indica, displaying both biotrophic and saprotrophic characteristics, is particularly useful for studying the evolution from saprotrophism to biotrophism (Lahrmann et al., 2015). Endophytes like NQ8GII4, Fo47, and P. fici share similarities with pathogens. P. fici was first identified as a pathogen of Ficus carica. Mycoviruses and the expansion of pectinase encoding genes might have participated in the transition of lifestyles (Wang et al., 2015; Wang et al., 2023). F. oxysporum causes vascular wilt diseases in over 100 different crops (Constantin et al., 2020). The nonpathogenic strain Fo47 is the most studied biocontrol strain of F. oxysporum (Aimé et al., 2013). Redkar et al. (2021) demonstrated that pathogenic and nonpathogenic F. oxysporum strains share a set of conserved pathogenicity factors, underscoring that the nonpathogenic F. oxysporum strains are closely related to pathogenic strains. Notably, NQ8GII4 and F. solani are closely related taxa that diverged 45 to 50 MYA. We speculated NQ8GII4 evolved from a Fusarium phytopathogen ancestor.

GTs are enzymes that build complex carbohydrates from activated sugar donors and are involved in essential activities of fungal cells, such as cell wall synthesis, the glycogen cycle, and the trehalose cycle (O’Connell et al., 2012). The fungal cell wall is primarily composed of chitin, glucans, mannans, and glycoproteins. Fungal cell wall-degrading enzymes (FCWDEs), which are responsible for the breakdown or modification of the fungal cell wall, are crucial for fungal development (Lyu et al., 2015). CYP51, known as sterol-14α-demethylase, is a critical enzyme in ergosterol biosynthesis, influencing cell membrane fluidity and hyphal morphology (Tsybruk et al., 2023). The classes of GTs, CYP51, and FCWDEs enriched abundant downregulated DEGs at 0.5, 1, and 6 dpi. Given the roles of GTs, CYP51 and FCWDEs in fungal biology, this down regulation could play a role in maintaining the balance of the mutualistic relationship by preventing overgrowth of the fungus.

Microorganisms have the ability to produce a vast array of bioactive secondary metabolites, most of which are derived from a polyketide and/or non-ribosomal peptide parent compound(s). Formation of a polyketide parent compound is catalyzed by a PKS, while formation of a non-ribosomal parent compound is catalyzed by an NRPS (Hansen et al., 2015). Luo et al. (2020) discovered the deletion of the V. dahliae NRPS gene VdNPS reduced virulence of the fungus on cotton and tobacco. Lu et al. (2022) demonstrated the Curvularia lunata PKS gene Clpks18 was required for production of the toxin methyl 5-(hydroxymethyl) furan 2-carboxylate as well as wild-type levels of virulence of the fungus on maize. In the current study, the core biosynthetic gene of secondary metabolite cluster 9 (NRPS) were downregulated at 0.5, 1, 6, and 14 dpi, and the core biosynthetic gene of secondary metabolite cluster 4 (T1PKS) was significantly downregulated at 14 dpi. We speculate that this may be related to the synthesis of some pathogenic toxins. However, additional experiments are required to determine whether suppression of toxin production contributes to the endophytic lifestyle of NQ8GII4.

Plants have developed two primary defense systems to counteract microbial attacks. The first, known as MAMP-triggered immunity (MTI), involves the recognition of microbe-associated molecular patterns (MAMPs) by pattern recognition receptors (PRRs) (Rebaque et al., 2021). Successful pathogens can disrupt MTI by delivering effectors, leading to effector-triggered susceptibility (ETS) (Ngou et al., 2022). In response, plants deploy resistance proteins (R) to detect these effectors, initiating effector-triggered immunity (ETI). This dynamic interaction creates a balance between plants and microorganisms (Dodds and Rathjen, 2010).

The putative NQ8GII4 effectors FnEG1, FnEG8, FnEG168, and FnEG214 triggered cell death in N. benthamiana in a Agrobacterium-mediated transient expression assay. The corresponding genes were upregulated at 14 dpi in infected alfalfa roots. FnEG1, FnEG8, FnEG168, and FnEG214 contained AA9, Cerato-platanin, AltA, and Pectata lyase domain, respectively. The function of these proteins in triggering plant immunity has been widely reported. MoHrip1, a member of the AltA family, is recognized as a MAMP by the plant immune system. When tobacco plants were treated with MoHrip1, they exhibited necrosis, increased reactive oxygen species, and upregulated defense-related genes (PR1a, PR5, HSR203J, and HIN1), along with enhanced resistance to Botrytis cinerea and Pseudomonas syringae pv. Tabaci (Zhang et al., 2017). The recombinant proteins of Epl-1 (Cerato-platanin) from Trichoderma asperellum ACCC30536 could alter signal transduction pathways of SA, JA, and IAA in Populus davidiana × P. alba var. pyramidalis (PdPap), promoting seedling growth and pathogen resistance against Alternaria alternata (Yu et al., 2018). Purified VdPEL1 (pectate lyase) triggered defense responses and conferred resistance to B. cinerea and V. dahliae in tobacco and cotton plants (Yang et al., 2018). Zarattini et al. (2021) treated Arabidopsis with Thermothielavioides terrestris purified AA9 protein. The phenomenon, such as an increased level of ethylene, jasmonic and salicylic acid hormones, along with deposition of callose in the cell wall, was observed. Therefore, we speculate that the establishment of the mutualistic relationship between NQ8GII4 and alfalfa involves production, by the fungus, of effectors that trigger plant immunity, which in turn slows growth of the fungus.

Autophagy is a programmed cell degradation mechanism that is ubiquitous in eukaryotic cells. Organisms can remove and degrade excess biological macromolecules and self-damaged cells that are produced in biological processes through autophagy and use the degradation products to provide energy. In fungi, autophagy plays an important role in regulating growth, development, and pathogenicity (Jiao et al., 2021). Wang et al. (2022) demonstrated that FvAtg4 and FvAtg8 contribute to F. verticillioides pathogenicity to maize by regulating the autophagic pathway to control lipid turnover, fumonisin biosynthesis, and pigmentation during its infectious cycle. Jiao et al. (2022) concluded that the SsAtg1 gene of Sclerotia sclerotiorum is involved in response to nutritional stresses that govern mycelial growth and is essential for sclerotia formation, contributing to virulence and the development of compound appressoria. In our study, a large number of autophagy-related genes, such as ATG1, ATG2, ATG11, and others. were downregulated at 6 and 14 dpi. This finding suggests that the autophagy pathway was inhibited during the process of NQ8GII4 symbiosis. Therefore, we hypothesized the endophytic characteristics of NQ8GII4 might be related to the inhibition of autophagy pathway.