Characterization of the transcriptional responses of Armillaria gallica 012m to GA3

Gastrodia elata needs to establish a symbiotic relationship with Armillaria strains to obtain nutrients and energy. However, the signaling cross talk between G. elata and Armillaria strains is still unclear. During our experiment, we found that the vegetative mycelium of Armillaria gallica 012m grew significantly better in the media containing gibberellic acid (GA3) than the blank control group (BK). To explore the response mechanism, we performed an RNA-sequencing experiment to profile the transcriptome changes of A. gallica 012m cultured in the medium with exogenous GA3. The transcriptome-guided differential expression genes (DEGs) analysis of GA3 and BK showed that a total of 1309 genes were differentially expressed, including 361 upregulated genes and 948 downregulated genes. Some of those DEGs correlated with the biological process, including positive regulation of chromosome segregation, mitotic metaphase/anaphase transition, attachment of mitotic spindle microtubules to kinetochore, mitotic cytokinesis, and nuclear division. These analyses explained that GA3 actively promoted the growth of A. gallica to some extent. Further analysis of protein domain features showed that the deduced polypeptide contained 41 candidate genes of GA receptor, and 27 of them were expressed in our samples. We speculate that GA receptors exist in A. gallica 012m. Comparative studies of proteins showed that the postulated GA receptor domains of A. gallica 012m have a higher homologous correlation with fungi than others based on cluster analysis.


Introduction
Gastrodia elata Blume is commonly called Tian ma in Chinese and is primarily distributed in the mountainous areas of East Asia, Southeast Asia and Oceania (Ahn et al. 2007;Chen, Ji and Zhu, 1999;Liu et al. 2015;Tang and Eisenbrand 1992;Zhan et al. 2016). Modern pharmacological experiments have demonstrated that extracts of G. elata or its active compounds possess wide-reaching biological activities, including antitumor, antivirus, memory-improving, antioxidation and anti-aging actions (Heo et al. 2007;Hu and Shen, 2014;Huang and Wang 1985). As a perennial, non-photosynthetic orchid, G. elata needs to establish a symbiotic relationship with Armillaria strains at late stages to obtain nutrients and energy (Chen et al. 2020;Lan and Li, 1994;Wang et al. 2007;Zou, Bai, Lv, Xu and Sun, 2010). The species of genus Armillaria were well known as phytopathogens that caused serious root rot disease on natural forests, tree plantations for timber production, orchards, vineyards and gardens (Heinzelmann et al. 2019;Koch,  Wilson, Sense, Henkel and Aime, 2017;Labbé et al. 2017;Onsando et al. 1997). However, several species of Armillaria have been identified to be able to form symbiosis with G. elata (Cha and Igarashi 1995;Cha & Igarashi 1996;Guo, Wang, Xue, Zhao and Yang, 2016;Sekizaki et al. 2008) and affect both the yields and the quality of the tubers of Tian ma (Guo et al. 2016).
It has been known for many years that phytohormones and their signaling network play a central role in the regulation of plant growth and development during plant life cycle, and recent research has revealed that these signaling pathways are probably present in non-plant species (Bari and Jones 2009;Grant and Jones 2009). Various microbes can produce or affect the phytohormone levels of their hosts; furthermore, auxins, cytokinins (CKs), GAs, abscisic acid (ABA) and ethylene (ET) had been isolated from microbial fermentation (Chanclud and Morel 2016;R. Eichmann et al. 2021;Ludwig-Müller 2015). On the other side, many works have disclosed that host hormones play an important role in plant-microbe interaction. IAA had been proved to be able to induce the production of invasive filaments in Saccharomyces cerevisiae (Prusty et al. 2004) and increase the antimicrobial activity of Streptomycetaceae strains isolated from Arabidopsis roots (van der Meij et al. 2018). CKs improve the progression of fungal biotrophic pathogens and stimulate beneficial AMF interactions in pea (Chanclud and Morel 2016). Exogenous ABA induces vast gene expression changes in the endophyte Aspergillus nidulans . GAs inhibit nodulation and AMF formation, while BRs seem to support it (Foo et al. 2016(Foo et al. , 2013. Historically, GA was primarily known to be a prominent type of phytohormone that regulates many aspects of plant development and physiology (de Lucas et al. 2008;Olszewski et al. 2002). Although GA biosynthesis by fungi has been documented and reviewed, few studies have focused on the signaling pathways of GAs in fungi (Hedden 2020). Through transcriptome analysis and phytohormone quantification, Takeda et al. revealed the accumulation of GA in Japanese lotus root lines infected with AM fungi, indicating the importance of de novo synthesis of GA in arbuscular mycorrhizal development (Takeda et al. 2015). It was reported that exogenous GA3 significantly increased hyphopodium formation at the epidermis, thus leading to the promotion of fungal colonization and arbuscule formation in the root cortex (Tominaga et al. 2020). Furthermore, GAs could promote the germination of conidia and the growth of the mycelium of Neurospora crassa (Barker and Tagu 2000). However, no specific studies have been carried out at the transcriptome level to clarify the responses of Armillaria starin to GA3 and how they relate to the symbiotic relationships.
Some research works have investigated the mechanisms of its symbiotic relationship with G. elata (Li, Guo and Lee, 2020;Tsai et al. 2016). Yuan et al. observed that strigolactone, as a phytohormone, could promote the branching and development of Armillaria mellea hyphae (Yuan et al. 2018). However, there is still limited knowledge on the symbiotic molecular mechanism between the Armillaria strain and G. elata. In previous works, we sequenced the genome of A. gallica 012m isolated from the tuber of G. elata collected from a plantation field, and investigated the symbiotic relationship with G. elata (Zhan et al. 2020). During the experiment, we found that the vegetative mycelium of A. gallica 012m grew significantly better in media containing GA3 than the control. The experimental phenomenon implied that GA3 might be an important signal molecule for A. gallica 012m. Here, we show the genome-guided analysis of differential expression between additional GA3 and BK. These studies provide a foundation to characterize GA3 signaling in A. gallica 012m.

Sampling of A. gallica 012m
In a previous study, we isolated A. gallica 012m from G. elata tubers in 2017 (Zhan et al. 2020). In this experiment, we activated the mycelium of A. gallica 012m with a modified Czapek-Dox medium with an addition of 2 mg/L GAs or not and incubated at 25 ℃ for 15 d in the dark (Zhan et al. 2020). The fresh mycelium was placed in RNA later, stored at room temperature for 24 h, then moved to −20 ℃ for storage until RNA isolation.
Three pellets from liquid cultivation were inoculated in liquid media of GAs and BK in each bottle. Each medium had four parallels and was cultured in a 160 r/min shaker at 25 °C for 16 d. After the cultivation was completed, the pellets were collected and weighed by filtration through a Buchner funnel. Using IBM SPSS v23 (Tanner-Smith and Tipton 2014) software, a paired-sample T test was performed on the blank group, and the GA groups were used to detect the difference in hormone treatment on the growth of A. gallica 012m.

RNA extraction, library preparation and sequencing.
Total RNA was isolated from fresh mycelium using the RNeasy mini kit (Qiagen, Hiden, Germany) following the manufacturer's instructions. The cDNA library was sequenced using the Illumina HiSeq (Caporaso et al. 2012) with a double-end sequencing strategy in Novogene Bioinformatics Technology, Beijing, China. The original data was stored in the National Center for Biotechnology Information (NCBI) database with SRA research registration number PRJNA759758.

Reads assembly and transcriptomes annotation
The high-quality clean data of samples were obtained by removing the sequence of adapters or filtering low-quality reads with Trimmomatic v.0.36 (Bolger et al. 2014). Then, the quantification of clean reads was estimated using FastQC v0.11.9 (Brown et al. 2017). Next, Hisat2 v2.1.0 (Hical Indexing for Spliced Alignment of Transcripts), a sequence alignment software, was used to mapping the chain-specific transcriptome data sequence with A. gallica genome to obtain the SAM data. Next, the sam file was converted into bam file by using Samtools software as the input file of StringTie v2.1.0 (Pertea et al. 2015). FPKM values (Fragments Per Kilobase of exon model per Million mapped fragments) were calculated. The gene expression matrix was initially filtered to filter out genes with low expression levels, and finally genes that could be used for subsequent analysis were obtained.

Gene ontology (GO) functional classification
The most detailed information on gene function was provided by enrichment analyses which facilitated us to screen for unique insights into GA3 responses. The read count matrix was imported into R v3.6.3. The edgeR of R package was used for differential gene analysis, with an FDR < 0.05 and |log2FC|> 1. Then, for functional annotation and classification, all transcripts and their corresponding genes were compared by eggNoG-mapper (Huerta-Cepas et al. 2017). The enriched GO terms were ordered based on three categories, including biological process (BP), molecular function (MF) and cellular component (CC).

Functional analysis of growth-related genes of DEGs
Gene sequences with biological process function were extracted. Then the protein domains were annotated using Pfam platform (http:// pfam. xfam. org/) to analyze the function of the protein domain.

Prediction and expression of GA3 receptors in A. gallica 012m
To obtain entire GA receptor members in the genome of A. gallica 012m, we annotated the genome of A. gallica 012m using InterProScan (Paysan-Lafosse et al. 2022) and obtained candidate GA receptor genes. The expressed sequences were selected with a significance threshold of E-value ≤ le-5 and their conserved domains were detected by SMART (http:// smart. embl-heide lberg. de/).

Construction of receptor sequence evolutionary tree
We constructed an evolutionary tree including the screened sequences and the downloaded receptor protein sequences. All sequences were imported into MEGA7 (Kumar et al. 2016), and the ClustalW2 was used for multiple amino acid sequence alignment. The phylogenetic tree was constructed by NJ (neighbor-joining) method, and the number of calculations was 1000.

Morphological characteristics of A. gallica grown under GA
It was found that GA3 (2 mg/L) could promote the growth of A. gallica 012m in solid medium, as shown in Fig. 1.

Effects of plant growth substances on the biomass of A. gallica 012m
Compared to blank control group, the biomass of A. gallica 012m was significantly increased in those mediasupplemented with different concentrations of GA3. Among them, 2 mg/L GA made largest biomass of A. gallica 012m (Fig. 2, Table S1). In the three hormone-liquid medium, 2 mg/L GA made Armillaria's largest biomass of 0.7741 ± 0.0470 g/L. According to the biological expression data of different concentrations and blank control group, we found that GA3 could increase the biomass of A. gallica 012m by 56 ± 9%. Using IBM SPSS v23 software, the result showed that the increase of its biomass was significant (Table S1).

Evaluation of transcriptome sequencing data.
We obtained the above 10.0 Gb raw data by sequencing on the Illumina HiSeq platform, which could be used for further analysis after quality control and comparative quantification. 3 repeat transcriptome samples of 2 mg/L GA3 and BK were sequenced by using Illumina HiSeq platform to obtain 27,631,384; 30,479,924; 25,566,599; 29,533,772; 30,365,200 and 26,125,609 pairs of PE reads respectively (Table 1).
After eliminating the reads with low-quality adapter sequences, an average of 21.916.601 to 28.961.759 clean reading pairs was chosen from GA3 and BK, respectively. Read mapped percentage of clean data from all samples was higher than 91.74%, and read mapped percentages of BK3 308 Page 4 of 13 were the highest. Additionally, 88.01%-91.48% of pure readings were successfully mapped to the reference genome of A. gallica 012m.

PCA analysis
Dimensionality reduction was performed by principal component analysis, and the distribution of sample points displayed on the two-dimensional or three-dimensional plane. According to the position of the points, it can also be seen whether the samples belonging to the same group are together and whether the samples between different groups are separated. The results are shown in Fig. 3. It can be seen from the figure that the samples of different conditions are clearly distinguished, and the distance between biological replicates is close, indicating that the consistency of biological replicates and the difference between different groups are better (Fig. 3).

Enriched GO terms of upregulated and downregulated DEGs
In the GA3 vs. BK comparison, 1309 genes were differentially expressed, including 361 upregulated genes (Supplementary Table. S1) and 948 downregulated genes ( Fig. 4; Supplementary Table. S2). To better visualize the functions and pathways of DEGs obtained from RNA-Seq, we  compared sequences of A.gallica 012m with the GO database. The GO database's annotation result is shown in Fig. 5.

Genes involved in the growth of A. gallica 012m
A total of 1309 differential genes were annotated in the comparison of GA3 vs BK, including 361 upregulated genes and 948 downregulated genes. To further explore the mechanism of stimulating the growth of A. gallica 012m, we compared the biological growth process of DEGs. The result showed that 68 differential genes were expressed in the upregulated genes, and 50 differential genes were expressed in the downregulated genes. In the upregulated genes, the genes controlling the growth of A. gallica 012m mainly involved cell cycle, signal transaction and cellular protein metabolism, such as positive regulation of chromosome segregation (GO:0051984), positive regulation of mitotic metaphase / anaphase transition (GO:0045842), positive regulation of attachment of mitotic spindle microtubules to the kinetochore (GO: 1,902,425) and so on. In the downregulated genes, some negative regulation genes were involved, such as negative regulation of cell cycle (GO:0045786) related to the cellular process, negative regulation of signal transaction (GO:0009968) related to signaling, negative regulation of cell aging (GO:0090344) related to development process. Those differential expressions explain that GA3 could promote the growth of A. gallica 012m to some extent.

Growth-related genes of DEGs analysis.
To further study the molecular mechanism of hormones promoting the growth of Armillaria, we analyzed the domain of genes involved in biological process regulation. The result showed that some genes are involved in regulating cell growth and metabolism (Table. 2). In the upregulated DEGs, Armga012mGene25682 contains HMG (high mobility group) domains, which are involved in the regulation of DNA-dependent processes, and Armga012m-Gene10992 containing S_TKc domains (serine/threonine protein kinases, catalytic domain), which play a role in a multitude of cellular processes, including division, proliferation, apoptosis and differentiation. Some DEGs have BTB domains (Broad-complex, Tramtrack and Bric a brac), GAL4 (GAL4-like Zn(II)2Cys6 (or C6 zinc) binuclear cluster DNA-binding domain), Methyltransf_11 domain and AAP2 (amino acid permease 2), which play an active role in the growth and development of cells. In the downregulated DEGs, Armga012mGene03731 and Armga012mGene12792 have protein domains of the WD40 domains (WD or betatransducin repeats), which can regulate the process of apoptosis. In addition, Armga012mGene06892 has BRCT protein domains (breast cancer carboxy-terminal domain), directly occluding this peptide-binding groove or disrupting essential conserved BRCT core-folding determinants. In short, the DEGs analysis could explain that GA enables the promotion of the growth of A. gallica 012m to some extent.

Prediction and expression of GA3 receptors in A. gallica 012m
GAs are natural complex biomolecules initially identified as secondary metabolites in the fungus Gibberella fujikuroi with strong implications in plant physiology, and GAs receptors are a class of proteins that play a vital role in the growth system. Forty-one candidate genes containing GA receptor domains were identified, and 27 speculated GAs receptor domains of A. gallica 012m were predicted using SMART (Table. 3) and expressed. The predicted GA receptor proteins were annotated, and 27 speculated GA receptor domains had the Abhydrolase_3 domain. Those predicted GA3 proteins have similar domains to the determined GA3 receptor proteins, and the expression levels of candidate GA3 receptor protein sequences were different (Fig. 6) in A. gallica 012m. Thus, we speculated that exogenous GA3 affected the growth of A. gallica 012m through GA receptors.

Genes containing GA receptor domain in A. gallica 012m
We downloaded the sequence information of GA receptors in fungi, peanut, soybean and rice from GenBank (Supplementary Table. S3) with a total of 52 sequences. Finally, by comparing the 52 genes with GA receptor domain and A. gallica 012m, an analysis of the phylogenetic relationship is shown in Fig. 7. The result showed that A. gallica 012m (Armgao012mGene20317) not only had homology with Valsa mali var. pyri (GenBank: KUI53266.1), but also had a high homology correlation with Colletotrichum gloeosporioides (GenBank: KAF3800814.1). The phylogenetic tree analysis idicated the results that A. gallica 012m might possess GA receptor domains.

Discussion
Plant-microbe interaction is essential to the balance in an ecosystem compared to the other microbial interactions (Sharma et al. 2020). Plants and microorganisms have coexisted and coevolved for millions of years, during which time both plants and microorganisms have been exposed to many signaling molecules produced and released by each other (Eichmann, Richards and Schafer, 2021). During our work on A. gallica culture, we attempted to investigate the responses of A. gallica to GA3. The results show that rhizomorph production of A. gallica 012m was simulated significantly by GA3. In another coincidental research, Jacobs et al. reported that the rhizomorph production of A. gallica was seriously inhibited by paclobutrazol (PBZ), known as a GA inhibitor. Their experiments showed that only 4% of colonies growing in PBZ plates formed rhizomorphs compared to 100% of colonies in control plates (Jacobs and Berg 2000). Those results suggest that GA3 might influence the rhizomorph production of A. gallica as a signaling molecule. The DEGs analysis of GA3 vs. BK showed that there were 361 upregulated genes and 948 downregulated genes. In addition, enrichment analyses of GO terms were performed for DEGs using EggNoG-mapper. One interesting finding was that the genes regulating growth showed upregulated expression in the GO analysis of differentially expressed genes, such as "regulation of biological process", "positive regulation of biological process" and "biological process". These analyses explained that GA3 actively promoted the growth of A. gallica to some extent.
Increased analysis of the signaling pathway revealed that the GA receptor members in A. gallica 012m were predicted in its genome and RNA-Seq based on BLASTp search with the full-length GA receptor protein sequence from UniProt. The results showed that candidate genes containing GA receptor domain of A. gallica 012m has better homologous correlation with those of fungi than plants.
Interestingly, more and more phytohormones receptors, including bacteria, fungi, algae, animal, etc. (Bodrato et al. 2009;Frebortova and Frebort 2021;Galindo-Trigo, Gray, and Smith, 2016;Herivaux et al. 2017;Jiang and Guo 2010;Ju et al. 2015;Lenasi et al. 2002;Samanovic et al. 2018;Zwanenburg and Pospisil 2013), had been predicted based on specific responses and gene sequences in non-plant species. Wang et al. had identified PcrK, a membrane-bound histidine kinase from the phytopathogenic bacterium Xanthomonas Campestris, as a bacterial receptor that specifically detects the plant cytokinin 2-isopentenyladenine (2iP) to promote adaptation to oxidative stress (Wang et al. 2017). In another case, SynEtr1 from Synechocystis PCC 6803 had been characterized as a functional ethylene receptor (Hernandez-Prieto et al. 2012;Klahn et al. 2015;Lacey and Binder 2016;Narikawa et al. 2011;Ramakrishnan and Tabor 2016;Song et al. 2011). Moreover, some CKs and ET receptor homologs were identified in zooxanthellae, diatoms, chromerids, free-living pseudofungus, filamentous marine protists and brown algae (Bidon et al. 2020;Kabbara et al. 2019;Papon and Binder 2019). From this point of view, phytohormones as inter-kingdom signaling molecules could influence several non-plant species and the receptors could exist in those species. Taken together, we speculate that GA could be a type of signaling molecule and there may be some receptors in A. gallica 012m.

Conclusions
In summary, GA3 improved the mycelium growth of A. gallica 012m significantly in both solid and liquid media. Several DEGs identified through RNA-seq analysis were related to biological process, including "regulation of biological process", "positive regulation of biological process" and "biological process". These analyses explained that GA3 actively promoted the growth of A. gallica to some extent. In this study, analysis of protein domain features showed that deduced polypeptide contained 41 candidate genes of the GA receptor domain. Comparative studies of proteins showed that the postulated GA receptor domains of A. gallica 012m have a higher homologous correlation with fungi than others based on cluster analysis. We speculate that GA receptors exist in A. gallica 012m.

Author contributions All authors reviewed the manuscript.
Funding This work was financially supported by the National Natural Science Foundation of China (Nos. 81860624). The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript. The following grant information was disclosed by the authors: National Natural Science Foundation of China: 81860624.

Data availability
The following information was supplied regarding data availability: data is available at SRA: PRJNA759758.

Declarations
Competing interests The authors declare that there are no competing interests, and the manuscript was approved by all authors for publication.
Ethical approval We hereby certify that our research did not involve human participants and/or animals.