Comparative transcriptome analysis of juniper branches infected by Gymnosporangium spp. highlights their different infection strategies associated with cytokinins
BMC Genomics volume 24, Article number: 173 (2023)
Gymnosporangium asiaticum and G. yamadae can share Juniperus chinensis as the telial host, but the symptoms are completely different. The infection of G. yamadae causes the enlargement of the phloem and cortex of young branches as a gall, but not for G. asiaticum, suggesting that different molecular interaction mechanisms exist the two Gymnosporangium species with junipers.
Comparative transcriptome analysis was performed to investigate genes regulation of juniper in responses to the infections of G. asiaticum and G. yamadae at different stages. Functional enrichment analysis showed that genes related to transport, catabolism and transcription pathways were up-regulated, while genes related to energy metabolism and photosynthesis were down-regulated in juniper branch tissues after infection with G. asiaticum and G. yamadae. The transcript profiling of G. yamadae-induced gall tissues revealed that more genes involved in photosynthesis, sugar metabolism, plant hormones and defense-related pathways were up-regulated in the vigorous development stage of gall compared to the initial stage, and were eventually repressed overall. Furthermore, the concentration of cytokinins (CKs) in the galls tissue and the telia of G. yamadae was significantly higher than in healthy branch tissues of juniper. As well, tRNA-isopentenyltransferase (tRNA-IPT) was identified in G. yamadae with highly expression levels during the gall development stages.
In general, our study provided new insights into the host-specific mechanisms by which G. asiaticum and G. yamadae differentially utilize CKs and specific adaptations on juniper during their co-evolution.
Gymnosporangium asiaticum and G. yamadae are obligate biotrophic pathogens that cause pear-rust disease and apple-rust disease, respectively, hindering the development of orchard industry [1,2,3]. In addition, their infections also cause swelling of leaf axil and stems of the important medicinal plant juniper (Juniperus chinensis), threatening the cultivation of juniper [4,5,6]. Both G. asiaticum and G. yamadae are demicyclic and heteroecious rust fungi, because they parasite on gymnosperms as telial host, while selecting dicotyledons as aecial hosts [4, 7]. Although they can parasitize on one branch of the juniper at the telial stage, the symptoms are quite different (Fig. 1). Usually, G. asiaticum produces tongue-shapes or wedge-shaped telia directly at the base of needle leaves with slightly swelling (Fig. 1a), while G. yamadae infection results in galls forming in twigs (Fig. 1b) with enlargement of the phloem and cortex (Fig. 1c, d). After spring rains, the teliospores of G. asiaticum germinate preferentially and about a month earlier than G. yamadae, then G. asiaticum airborne basidiospores infect Pyrus, Chaenomeles and Cydonia species while G. yamadae airborne basidiospores infect only Malus species without cross-infection [7,8,9]. Similarly, the mature telia of G. yamadae emerge from the galls and germinate successively after several rains (Fig. 1b). The life cycle characteristics of the two Gymnosporangium species on junipers suggest a potentially special host selection mechanism, unlike other rust fungi whose telial hosts are angiosperms.
In the last few years, various studies have focused on the host selection mechanisms of heteroecious rust fungi. Related studies on the phylogeny of Pucciniales have demonstrated that the aecial stage is under more selective pressure and stronger evidence of co-diversification with hosts than the telial stage, and the results reflect that the telial hosts are less constrained for Pucciniales species [10,11,12]. However, phylogenetic analyses did not prove clear host specificity of Gymnosporangium species. The overlapping host range of telial stages within Gymnosporangium might be driven by multiple speciation mechanisms [8, 13]. Further comparative transcriptomic studies have revealed that stage-specific proteins (e.g., effector proteins and carbohydrate-active enzymes) secreted by heteroecious rust fungi show high correlation with the hosts that probably determine the host-specific selection [14, 15]. Until now, several transcriptome analyses have been performed to investigate the molecular mechanisms underlying the host specificity of Gymnosporangium-juniper interaction, from the perspective of rust pathogen. As expected, Gymnosporangium species parasitizing on junipers exhibited a conserved genetic expression program corresponding to energy, translation and signal transduction processes . In spite of the composition of CAZymes is similar across different categories, indicating a series of common genes between G. asiaticum and G. yamadae that ensure the ability to infect juniper , the clearly different symptoms on infected junipers between the two Gymnosporangium species imply their potentially different infection strategies.
Evidence showed that pathogens and insect species induce abnormal plant tissues, such as swollen gall, club-shaped root and leaf tumors, providing shelter and nutritional resource for these parasites and symbionts [16,17,18,19]. The formation mechanisms of abnormal host plant tissues have been investigated in recent years. Phytohormones, especially indole-3-acetic acid (IAA) and cytokinins (CKs), play key roles in regulating abnormal plant tissue development . Plant parasitic fungi, such as Taphrina deformans , Exobasidium gracile , Moniliophthora perniciosa [23, 24] and Ustilago maydis [19, 25], produce IAA or CKs participating in hyperplasia of host leaf, which contribute to their infections. Certain insect species synthesize IAA and CKs that reprogram and modify plant tissue to a nutrient-rich shelter for them [26,27,28,29]. In addition to the effects of parasites-derived phytohormones, other effector molecules including MicroRNAs , secreted proteins [28, 31] and small RNAs  are deemed to be required to regulate the abnormal growth of host plant. The bottom line is that these organisms that provoke plant tissue proliferation are all dependent on living plants, which may reflect the common survival strategies of these parasites or symbiotic organisms. Rust fungi can induce malformation, swelling and even galls on plant hosts, the symptom is relatively more typical in Gymnosporangium species [7, 13, 33, 34]. In G. juniper-virginianae-cedar pathosystem, the rust fungi synthesized CKs that might facilitate the gall formation , which enlighten us postulate that CKs derived by G. yamadae might be an important factor that causes galls and results in the different symptoms on juniper between G. asiaticum-juniper and G. yamadae-juniper pathosystems.
The mechanism of cytokinins biosynthesis is conserved in plants . On the one hand, most isopentenyladenine (iP)- and trans-zeatin (tZ)-type cytokinins synthesis are synthesized de novo by the activity of adenylate isopentenyltransferases (IPTs); on the other hand, a large fraction of the cis-zeatin (cZ)-type cytokinins is derived from tRNA degradation, which depends on the catalysis of tRNA-isopentenyltransferase (tRNA-IPTs) [35, 36]. However, all CKs presented in fungi are consistent with origins from within the tRNA degradation pathway . In recent years, investigations of microorganism emphasized the importance of tRNA-IPTs in tRNA degradation pathway of CKs production. The first characterization of tRNA-IPT gene in filamentous fungi was carried out in Claviceps purpurea and was proved essential for cZ biosynthesis and contributing to the formation of iP and tZ in the fungus . Research on Magnaporthe oryzae  and U. maydis [25, 37] also demonstrated that tRNA-IPTs mainly contribute to cZ-type CKs production and involve in virulence of these plant pathogens. Several previous studies have also shown that rust fungi, including Puccinia recondita f. sp. tritici  and G. juniper-virginianae  are able to produce CKs, while the molecular mechanisms of rust-derived CKs as well as their roles in rust-plant interactions remain unclear.
In this study, we executed RNA-seq analysis of juniper branch tissues after the two Gymnosporangium spp. infections at different stages and revealed similar and different genetic responses of juniper in the interaction between G. asiaticum and G. yamadae. The temporal transcriptional analysis of the gall showed an extensive and dynamic gene regulation pattern that is G. yamadae dominated. Furthermore, we obtained evidence that G. yamadae-derived CKs depended on the catalysis of tRNA-isopentenyltransferase (tRNA-IPT) identified in G. yamadae are most likely contribute to galls formation and development, which suggested the distinct infection strategies between G. asiaticum and G. yamadae. Taken together, the study provides new insights into the host-specificity mechanism of G. asiaticum and G. yamadae.
Experimental design and RNA-seq results
Total RNA for each sample from three independent stages was extracted and subsequently sequenced as a dedicated analysis pipeline (Additional file 1: Fig. S1). Average of 560, 549, 411, and 459 million clean reads per sample were obtained from infected juniper tissues, and 460 and 409 million clean reads per sample from healthy branches for which collected accordingly at the late stage of G. asiaticum and G. yamadae, respectively (Table 1). The proportion of reads mapped ranged from 87.21% for average in healthy branches sample to 56.65% for average in the late stage of infection branches sample, which reflects the dynamic development of plant tissue and rusts (Table 1).
We finally obtained 235,589 unigenes for juniper and 57,513 unigenes for the two Gymnosporangium species after separated assembled unigenes between juniper and the two Gymnosporangium species (Additional file 2: Table S1). The proportional distribution of homology of Juniperus chinensis and the two Gymnosporangium species unigenes in NR were shown in Additional file 1: Fig. S2. We assessed the overall reproducibility of the data by combining principal component analysis (PCA) and Pearson correlation analysis, and eliminated III_GY_J2 because of its unclear separation among independent conditions based on gene expression (Additional file 1: Fig. S3). Overall, the transcriptome generated by RNA-seq approach is suitable for subsequent analysis between different samples.
In total, 40,581 and 27,173 CDSs were determined using in juniper and the two Gymnosporangium species transcriptomes, respectively. Ultimately, a total of 32,200 (79.35%) unigenes of juniper and 25,623 (94.30%) unigenes of rust were annotated in at least one database, and 9,065 (22.34%) and 6,975 (25.67%) unigenes, respectively, were annotated in all of the above databases (Additional file 2: Table S2). A total of 4,193 differentially expressed genes (DEGs) (3,216 up-regulated and 977 down-regulated) were identified in III_GA_J, likewise 4,843 (2,813 up-regulated and 2,030 down-regulated), 9,522 (7,023 up-regulated and 2,499 down-regulated) and 6,311 DEGs (3,116 up-regulated and 3,195 down-regulated) were identified in I_GY_J, II_GY_J and III_GY_J, respectivly (Fig. 2, Additional file 2: Table S3). These results indicated that the widespread variation in gene expression were induced in the infected plant tissues.
Transcriptional responses of juniper tissues infected by G. asiaticum and G. yamadae
The Venn diagram (Additional file 1: Fig. S4) showed that 1,391 DEGs were commonly up-regulated and 580 DEGs were commonly down-regulated in III_GA_J sample and III_GY_J sample. KEGG pathway enrichment analysis of these DEGs (Fig. 3) showed that the commonly up-regulated genes were mostly enriched in ‘transport and catabolism’, ‘transcription’, ‘spliceosome’ pathways. Conversely, the commonly down-regulated genes were predominantly enriched in ‘energy metabolism’, ‘photosynthesis’, and ‘ribosome’ pathways (Fig. 3a). Furthermore, the DEGs enriched in the 'carbohydrate metabolism', 'enrgy metabolism' and 'biosynthesis of other secondary metabolites' pathways were up-regulated in III_GY_J sample but down-regulated in III_GA_J sample (Fig. 3b).
We noted that there 27 DEGs (Type_a) were downregulated in III_GA_J sample while upregulated in III_GY_J sample, 226 DEGs (Type_b) were upregulated in III_GA_J sample while downregulated in III_GY_J sample (Additional file 1: Fig. S4). Based on Swissport annotations, half of DEGs in Type_a relate to plant immunity, such as berberine bridge enzyme-like 18 and probable aspartic proteinase GIP2; some genes relate to plant hormone synthesis and cell growth and development, such as 1-aminocyclopropane-1-carboxylate oxidase, 12-oxophytodienoate reductase 3, protein SOSEKI 1 and protein ASMMETRIC LEAVES 2 (Additional file 2: Table S4). In Type_b DEGs, most genes were related to plant immune defense, such as leucine-rich repeat receptor-like serine/threonine-protein kinase, cysteine-rich receptor-like protein kinase 10 and aspartic proteinase nepenthesin-2 (Additional file 2: Table S4). However, plant defense-related genes were differently regulated in juniper after infection by G. asiaticum and G. yamadae.
MapMan analysis (Fig. 3c) showed that a large number of genes were significantly regulated in signaling, proteolysis and secondary metabolism categories during G. asiaticum and G. yamadae treatments. However, G. yamadae infection inhibited transcription of the cascades associated with the plant defense mechanism after recognition by the corresponding receptors, while fewer DEGs identified during G. asiaticum infection (Fig. 3c, Additional file 2: Table S5). More importantly, most G. yamadae DEGs were downregulated, especially genes involved in signaling (262 DEGs) and PR-proteins (64 DEGs), compared to G. asiaticum DEGs (66 and 5 DEGs in signaling and PR-proteins) (Fig. 3c, Additional file 2: Table S5). Hormone signaling-related genes, including ethylene, auxins, ABA and JA were downregulated in III_GY_J sample while highly expressed III_GA_J sample (Fig. 3c).
G. yamadae induces widespread transcriptomic changes in galls
The functional enrichment analyses in GO databases showed that mainly up-regulated genes were mostly classified into cellular component and biological process categories. The common GO terms, including ‘nucleus’ and ‘kinase activity’ between I_GY_J and II_GY_J; ‘cytoplasm’ and ‘intracellular membrane-bounded organelle’ between II_GY_J and III_GY_J; ‘cellular process’ and ‘ion binding’ between I_GY_J and III_GY_J (Fig. 4a). KEGG pathways analysis indicated that DEGs were commonly significantly enriched in the photosynthesis-related pathways, with most of them being down-regulated across all three stages. The rest of the enriched pathways, however, varied among the three stages. In I_GY_J sample, pathways such as ‘glycine, serine and threonine metabolism’, ‘RNA polymerase’, ‘steroid biosynthesis’ were significantly enriched; in II_GY_J sample, pathways representing ‘folding, sorting and degradation’, ‘transcription’ and ‘metabolism of cofactors and vitamins’ were significantly enriched; in III_GY_J sample, ‘energy metabolism’, ‘plant hormone signal transduction’ and ‘alpha-linolenic acid metabolism’ pathways were predominantly enriched, which mainly involved with up-regulated genes (Fig. 4b).
The heatmaps showed that during the middle stage (II_GY_J) of gall development, genes associated with photosynthesis, sugar metabolism, plant hormones and defense-related pathways were predominantly up-regulated at (Fig. 5, Additional file 2: Table S6). The most significant up-regulated genes in II_GY_J sample involved in photosynthesis are the light-harvesting chlorophyll protein complexes LHCB4, psaG in photosystem I, cytochrome b6f complex petC, psbW in photosystem II and photosynthetic electron transport petF (Fig. 5a, Additional file 2: Table S6). Additionally, there is another notable cluster of up-regulated genes in sugar metabolism pathway, including beta-glucosidase (which demonstrated a 13-fold increase), and glycogen synthase and trehalose 6-phosphate synthase (both of which had an eightfold increasein the II_GY_J sample) (Fig. 5b, Additional file 2: Table S6). Genes involved in jasmonic acid (JA), salicylic acid (SA) and ethylene (Eth) signal transduction pathways showed a common expression pattern, which were up-regulated gradually from the early stage to the middle stage and repressed at the late stage of the gall development (Fig. 5c). In contrast, a set of genes in cytokinins (CKs), abscisic acid (ABA) and auxin signal transduction pathways were particularly up regulated II_GY_J sample, such as auxin influx carrier, cytokinins receptor and ABA receptor (Fig. 5c, Additional file 2: Table S6). Plant defense-related genes were extensively up-regulated in II_GY_J sample, especially some typical plant immune signal recognition receptors such as calcium-dependent protein kinase and calmodulin, as well as some bacterial induced defense genes (Fig. 5d, Additional file 2: Table S6).
Identification of CKs biosynthesis-related gene tRNA-IPT in G. yamadae
Gall tissue and the telia of G. yamadae contained substantially higher concentrations of CKs than healthy young branches, while concentrations of CKs in G. asiaticum telia were slightly higher than healthy young branches (Fig. 6a). To further test the ability of CKs production in G. yamadae, we examined the changes of CKs content during the teliospores germination. The result showed (Fig. 6b) that the CKs content was higher at 8 h than at 4 h during the teliospores germination.
The features of CKs distribution in the gall of G. yamadae facilitate our study of the possibility that G. yamadae-derived CKs. We searched the KOG annotation results in our previously sequenced transcriptomes of G. asiaticum and G. yamadae teliospores and identified c25847_g1 and c18939_g3 annotated as “tRNA delta (2)-isopentenyl pyrophosphate transferase”. However, tRNA-IPT identified in this study was obtained from only one candidate, TRINITY_DN3125_c0_g1_i1, in juniper transcriptome database (Additional file 2: Table S7). We confirmed c18939_g3 and TRINITY_DN3125_c0_g1_i1 were one gene by protein sequence alignment, and identified TRINITY_DN3125_c0_g1_i1 as G. yamadae tRNA-IPT (GytRNA-IPT), as well as G. asiaticum candidate tRNA-IPT (c25847_g1) and two juniper candidate tRNA-IPTs (TRINITY_DN46394_c0_g1_i2 and TRINITY_DN16719_c0_g1_i1) were identified in this study to examine the relatedness to other tRNA-IPTs from a phylogenetically distinct group of organisms (Additional file 2: Table S7). An unrooted phylogenetic tree was constructed after trimming and aligning the tRNA-IPT sequences (Fig. 7). The phylogenetic tree was divided into six clusters, which were eukaryotic tRNA-IPTs, plant adenylate IPTs, eukaryotic origin plant tRNA-IPTs, prokaryotic origin plant tRNA-IPTs, fungal adenylate IPT and bacterial adenylate IPTs. In the phylogenetic tree, GytRNA-IPT was grouped with c25847_g1 in the cluster of eukaryotic tRNA-IPTs of plant pathogenic fungi. Alignment of GytRNA-IPT and c25847_g1 with other biotrophic fungal tRNA-IPTs, revealed conservation of ATG/GTP binding site, DMAPP binding site and a typical eukaryotic tRNA-IPTs zinc-finger domain in C-terminal (Additional file 1: Fig. S5). We noted that a chloroplast target sequence (18–54 amino acids) was identified within the N-terminal extension, and two nuclear localization signals (KKLWNEHVQSKRHRS, KRHRSASRPRKSYQHHKK) were predicted for the C-terminus of the protein.
The expression level of GytRNA-IPT was up-regulated in the middle stage (II_GY_J) relative to in the early stage (I_GY_J) and the late stage (III_GY_J) of gall development (Fig. 8). The RT-qPCR results of GytRNA-IPT were consistent with the FPKM values in the transcriptome data. These results indicated that the synthesis of CKs is stronger in gall tissues than in the branch tissues where galls are formed.
The two Gymnosporangium species induced specific transcriptional responses during the biotrophic interactions with juniper
Generally, genes associated with photosynthesis and disease resistance in host plants are inhibited by rust pathogens during their infections, which are well known to be an indispensable infection strategy of obligate biotrophic pathogens [41,42,43]. In the present study, juniper genes related to photosynthesis and metabolism of terpenoids and polyketides, linolenic acid and α-linolenic acid in both G. asiaticum and G. yamadae infected tissues, which suggesting general responses of plant after pathogen infection [16, 23, 44]. However, up-regulated genes in G. yamadae induced gall tissues enriched in various pathways related to substance synthesis and metabolism, suggesting a more active metabolic activity in G. yamadae induced gall tissues than in G. asiaticum infected juniper tissues.
The plant responses to biotic stress including recognition of the initial signal from pathogen, followed by triggering a series transcriptional cascade and finally inducing the production of defense molecules [45, 46]. At the very late infection stage of G. asiaticum on juniper, downstream signal transduction and defense mechanism were extensively activated whereas bypassed initial signal recognition. However, major genes, particularly related to signal recognition, transduction and PR-proteins encoding were inhibited in G. yamadae induced gall tissues. It reflected that the differences of the two Gymnosporangium species infection strategies originally occurred at the molecular level of signal recognition in their interactions with juniper.
Temporal transcriptional analysis revealed a G. yamadae-manipulated gene regulation pattern associated with gall development
Insects-induced leaf galls as resource sinks and shelters for these living organisms, for which exhibited wide range of gene expression changes, referring to the inhibition of genes related to photosynthesis and plant immune defense response and activation of genes related to nutrient metabolism [16, 17, 29, 47]. In the present study, genes associated with photosynthesis, chlorophyll synthesis and porphyrin and chlorophyll metabolism pathways were suppressed at the early and the late stages of the gall development . Unexpectedly, a specific set of genes was significantly up-regulated at the middle stage (II_GY_J) and involved in light-harvesting and excitation transfer category (psaG, psbW and LHCB4), electron transporting (petF) and redox (petC), which may essentially constructed a cooperative framework of the two photosystems [48,49,50]. Our study revealed that photosynthesis-related genes were not continuously suppressed during gall development, and the enhanced photosynthetic activity may be responsible to produce more nutrients during the maturity of G. yamadae teliospores .
Genes involved in ‘glycine, serine and threonine metabolism’ pathway at the early stage (I_GY_J) of the gall development may suggest the nutrient reserves in gall tissues. More importantly, the increase of soluble sugar content might contribute to the elevated osmotic, thus increasing the water content in gall tissues . Subsequently, along with G. yamadae teliospores development and maturity, vesicular transport-related genes and additional starch and sucrose synthesis-related genes were up-regulated at the middle stage (II_GY_J) and numbers of energy metabolism-related genes were up-regulated in the late stage (III_GY_J), which suggested a general function of galls as a carbohydrate source for the growth of G. yamadae teliospores. Correspondingly, rapid changes in the water and nutrient content in the gall lead to an increase in osmotic pressure, which induces the expression of ABA synthesis-related genes and an increase in ABA content, especially during the formation of large numbers of teliospores (II_GY_J) [33, 52]. Besides, genes involved in auxin signaling transduction pathway were also specifically up-regulated in II_GY_J, which possibly contributed to teliospores germination [20, 53, 54]. These results collectively reflected an autotrophic characteristic in G. yamadae induced gall tissues and indicated that the gall probably provide nutrients and water for teliospores of G. yamadae during the long-time development .
Many studies demonstrated that insects may have to suppress plant resistance systems during gall formation and development for better survival [16, 29]. However, we noted that a large number of plant defense response-related genes were activated during the middle stage (II_GY_J) of gall development. In addition, JA and SA signal transduction pathways were progressively activated at the middle stage (II_GY_J) and triggered defense-related response in juniper, which may be a result of the stronger interaction between G. yamadae and juniper along with the gall development and the teliospores maturation. Ultimately, the delayed activation of juniper defense responses at the late stage (III_GY_J) reflected a common infection strategy of biotrophic and hemibiotrophic plant pathogens [42, 55,56,57].
The main distinguishing features between other parasite-gall systems are the observation of gene expression profiles and the capture of dynamical changes in multiple developmental stages of gall, consistent with the phenotypic development of gall and the interaction between G. yamadae and juniper [29, 57, 58]. Furthermore, differences in the genetic regulation of host plant physiological responses, referring to secondary metabolism and stress defense mechanisms, maybe caused by phenotypic disparities [59,60,61]. Gymnosperms, especially conifers, are quite different from angiosperms in terms of leaf physical morphology and substance metabolism [61,62,63,64]. Consequently, the molecular mechanisms of fungi gall initiation and development remains an open question due to the lack of study on fungi-induced galls development and the difficulties in studying non-model plants .
G. yamadae-derived cytokinins synthesized related with tRNA-IPT involved in the regulation of gall development
CKs plays an important role in the information and development of galls induced by insects and in the colonization and virulence of plant pathogenic fungi [20, 25, 38, 65]. While it is unclear whether there is rust-derived cis-zeatin. A recent study confirmed that cis-zeatin was detected in G. juniper-virginianae induced gall tissues as well as the wet telia, but not in control cedar branchlets . To verify the effect of G. yamadae derived CKs on the induction of galls, we firstly examined the CKs content in juniper and G. yamadae, and found that galls and the telia contained significantly higher concentrations of CKs than the control sample (Fig. 6), consistent with the result in G. juniper-virginianae-galls system . Considering the close relationship between G. yamadae and G. juniper-virginianae from a phylogenetic point of view , we proposed that cis-zeatin synthesized by the tRNA degradation pathway is responsible for a large number of CKs in G. yamadae consistent with demonstrating that the expression pattern of GytRNA-IPT gene was consistent with the developmental characteristics in G. yamadae induced galls. Although we were unable to verify direct regulation of the tRNA-IPT gene in the synthesis of CKs in G. yamadae due to the lack of genetic transformation in obligate biotrophic rust, the above results suggest that G. yamadae-derived CKs were essential in the gall formation and development. More importantly, G. yamadae-derived CKs in the galls may also play a role in nutrient mobilization and delaying senescence in juniper young branches, allowing the biotrophic G. yamadae to establish and grow in the host [24, 65].
We identified the orthologues in G. yamadae (GytRNA-IPT) from our published data which contains all the typical characteristics of eukaryotic tRNA-IPTs [66, 67] and shows a close phylogenetic relationship with C. purpurea tRNA-IPT and U. maydis tRNA-IPTs whose functions have been well understood, which suggested their similar functions [25, 38]. The derivation of two juniper candidate tRNA-IPTs (TRINITY_DN46394_c0_g1_i2 and TRINITY_DN16719_c0_g1_i1) from both eukaryotic origin plant tRNA-IPTs and prokaryotic origin plant tRNA-IPTs indicates a functional differentiation of juniper tRNA-IPTs. The phylogenetic tree also showed that tRNA-IPTs from obligated biotrophic fungi, such as Austropuccinia psidii, Puccinia striiformis f. sp. tritici and Melampsora larici-populina were clustered in eukaryotic tRNA-IPTs clade (Fig. 7). Further studies focusing on the functions of those tRNA-IPT genes would be instructive in understanding the role of fungal-derived CKs in plant-fungi interactions.
A surprising finding was that GytRNA-IPT contains a predicted chloroplast transit peptide within the N-terminal extension. So far, tRNA-IPTs with chloroplast transit peptide have only been found in Arabidopsis . We confirmed that tRNA-IPT of M. oryzae also have a chloroplast transit peptide for N-terminus (Additional file 2: Table S7) whereas have no description in the previous study . One possibility is that tRNA-IPT in fungi can be transferred to function in host plants, similar to effectors.
Differentially utilizing of CKs providing a new insight into host selection mechanisms of rust fungi
Studies on the molecular mechanisms of rust-host plant specific interactions have shown that rust secretes effector proteins that largely determine host selection [15, 69, 70]. Nonetheless, the molecular mechanisms of rust-host plant interactions are mostly constrained from the perspective of effector proteins due to the lack of a genetic transformation system for rust. In contrast, in more plant-pathogen systems, secondary metabolites show a strong association with host specificity . As for Alternaria alternata, a conditionally dispensable chromosome controls toxins production that has a contribution to host-specific pathogenicity . At the genome level of Suillus, both secondary metabolites and pathways play a key role in enhancing host specificity by deactivating of reactive oxygen species . Based on the current experimental data, we propose that G. yamadae induced galls determined by cZ-type CKs synthesized via GytRNA-IPT catalysis cannot be omitted in G. yamadae-juniper specific interactions. Interestingly, CKs are also accumulated in alternate host lobes without gall formation during the G. yamadae-apple interaction . Combined with previous studies [28, 32], it is strongly suggested that there are other factors (e.g., effectors and small RNAs) in G. yamadae-juniper system that determine gall formation and development along with cZ-type CKs. The different infection characteristics of G. yamadae in juniper and apple leaves reflected the complicated mechanisms underlying the host-specific interactions of heteroecious rust.
Gymnosporangium species have a special host alternation pattern with its telial stage occurring on gymnosperms and aecial stage on angiosperms, which is contrary to most rusts parasitizing on gymnosperms during the aecial stage, such as Cronartium ribicola and M. larici-populina [8, 74, 75]. The majority of Gymnosporangium species have a narrow range of telial host selection restricted to only one or two genera, which usually results in different Gymnosporangium species parasitizing on the same tree, such as G. asiaticum and G. yamadae [4, 8, 13]. The two co-evolved Gymnosporangium species may exhibit strategies to avoid interspecific competition. Galls induced by G. yamadae on juniper branches ensure adequate nutrient and water during the long-time germination of the teliospores, thus providing more opportunities to transfer to alternative host. On the other hand, the adaptive trait of G. yamadae induced galls may be beneficial to avoid the competition of nutrients and water with G. asiaticum teliospores that also parasites on juniper. We hypothesized that the two Gymnosporangium species have evolved different infection strategies, i.e., different utilization of phytohormones, during their long evolutionary period, for the harmonious and common utilization of limited host resources. Although exactly how G. yamadae obtained the ability to induce the gall during the evolution remains unknown. Future work should emphasize the regulation of juniper cytokinin levels through pharmacological approaches, potentially disrupting teliospore formation as a strategy to manage rust disease.
In this study, we explore gene regulatory changes in juniper following G. asiaticum and G. yamadae infection through comparative transcriptome analysis to further understand the development and host specificity of Gymnosporangium species symptoms. Temporal transcriptional analysis revealed the extensive reprogramming of gene expression in gall tissues induced by G. yamadae. More importantly, we proposed that G. yamadae derived CKs depended on GytRNA-IPT regulation is an essential factor in the gall formation, which reflects a conservative characteristic of biotrophic pathogenic fungi to use their own phytohormones for infection. The different infection strategies of G. asiaticum and G. yamadae remind us that the factors determining the host-specific interaction of rust fungus are not single and may relate to secondary metabolites such as phytohormones.
Juniper branch morphology and anatomy
The upper surface and longitudinal sections of G. asiaticum telial on juniper leaves and G. yamadae telial on juniper branches were photographed by a digital camera (Leica, M205FA). To illustrate the morphology and anatomy of the gall tissues, infected branches of G. yamadae at the middle stage of gall development, with full-grown gall and obvious teliospores, were utilized. Healthy juniper branches and the gall tissues were collected simultaneously from a single juniper tree. The branches or galls were cut into longitudinal sections by hand or into 0.2 mm thick sections using a freezing microtome (Leica CM1950) for cross-sections. Longitudinal sections of branch were visualized under a stereomicroscope (Leica, M205FA); cross-sections of healthy branch and gall tissues were photographed by light microscopes (Leica, DM3000 and Leica, DM2500).
Experimental setup and sample collection
Juniper leaves and young branches infected by G. asiaticum and G. yamadae, including the early stage of gall development without the teliospores formation (I_GY_J), the middle stage with the teliospores break through the gall (II_GY_J) and the late stage of gall development with the matured teliospores (III_GA_J; III_GY_J), and healthy leaves and young branches (0_GY_J and 0_GA_J) as control were collected from Northwest A&F University, Yangling, Shaanxi and Haidian park, Beijing, China, between April and May 2021. All plant material were identified by Prof. Yingmei Liang. All samples were deposited in the Mycological Herbarium, Museum of Beijing Forestry University Mycologicum, Academiae (BJFC). Deposition number were: BJFC-R01700 (0_GA_J), BJFC-R02162 (III_GA_J), BJFC-R01692 (0_GY_J), BJFC-R02376 (I_GY_J)), BJFC-R02382 (II_GY_J), BJFC-R03656 (III_GY_J). Figure 1 and Additional file 1: Fig. S1 illustrate the experimental design. Sections of juniper tissues colonized by rust were collected and removed apparent fungal materials. Three biological replicates (i.e., J1, J2, J3) were harvested from each juniper, frozen separately in liquid nitrogen and stored at − 80 ℃ prior to RNA extraction.
RNA extraction, RNA-sequencing and assembly
Total RNA was isolated from 100 mg of infected leaf and branch tissues using the RNA Easy Fast Plant Tissue Kit (Tiangen, Beijing) according to the manufacturer’s instructions. RNA quality and quantity were checked with a Qubit® RNA Assay Kit with a Qubit® 2.0 Fluorometer (Life Technologies, CA, USA). As shown in Additional file 1: Fig. S1, we constructed the cDNA libraries according to the standard protocols developed by Illumina (Illumina, San Diego, CA, USA) following the manufacturer’s recommendations. In each library, 150 bp paired-end reads were generated from Illunima NovaSeq 6000 platform (Illumina, San Diego, CA, U.S.A.). The raw sequenced reads were processed by FASTQ . Next, the clean reads were used for de novo assembly using the Trinity software . The kallisto software  was used to estimate the expression abundance of the assembled transcripts with default parameters. Further cluster analysis by using CD-HIT software package  retained only one of the transcripts whose similarity exceeded the threshold. Such sequences removed from redundancy based on similarity are defined as unigenes.
In order to assess the quality and homology of the transcriptome data, we mapped the unigenes of J. chinensis to the genomes of Pinus lambertiana (GCA_001447015.2_Sugar_pine_JHU_assembly_genomic.fna), Pinus taeda (GCA_000404065.3_Ptaeda2.0_genomic.fna), Picea abies (GCA_900067695.1_Pabies01_genomic.fna), Picea glauca (GCA_000411955.6_PG29_v5_genomic.fna) and Gnetum montanum (Gnetum.final.fa) through BLASTn and BLASTp  with default parameters in NCBI. The unigenes of G. asiaticum and G. yamadae were annotated as described in the previous study . The unigenes that without comparison results were further annotated by a subset of plants and fungi from the NR database (NCBI non-redundant protein sequences). All unigenes that can be aligned were used to generate the transcriptomes of J. chinensis, G. asiaticum and G. yamadae (Additional file 1: Fig. S1).
TransDecoder (https://github.com/TransDecoder/TransDecoder/wiki) was used to predict the coding sequences (CDS) for all unigenes. Gene function was annotated in KO (KEGG Ortholog database), GO (Gene Ontology), Pfam (Protein family), and NR databases by EggNOG-emapper  and in Swiss-Prot (A manually annotated and reviewed protein sequence database), KOG/COG (Clusters of Orthologous Groups of proteins) databases by diamond , respectively.
Gene expression analysis
The FPKM expression levels and counts for all unigenes were estimated in each replicate by RSEM . Principal component analysis (PCA) was performed from read counts via the plotPCA function in R package . The similarity between samples at the expression levels was visualized by calculating the Pearson correlation coefficient between samples. Different expression genes (DEGs) was performed by DEseq2  and unigenes with a significant p-value (< 0.05) and more than a twofold change in transcript level were deemed as differentially expressed. Volcano plots of DEGs and heatmaps of gene expression profiles were generated using TBtools software . To derive expression patterns of genes in the different gall development stages, unigenes expressed at all stages were selected across all the stages with p-value < 0.05.
Functional analysis of DEGs
The annotated DEGs were used in functional enrichment analysis. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment analysis  and GO category annotation of significant DEGs were performed by clusterProfiler  with adjusted p-value < 0.05, and the TBtools software  was used to test the statistical enrichment of differentially expressed genes. In addition, significant DEGs were classified to functional categories using the pathway analysis program MapMan  with default parameters.
CKs extraction and quantitative analyses
Total cytokinins were extracted from 100 mg of healthy branches, G. yamadae infected young branches (galls) and the telial of G. asiaticum and G. yamadae. Measured teliospores were artificially germinated by immersing them in distilled water . Then the supernatant was took at the middle (4 h) and late (8 h) time stages to measure the CKs content. Briefly, all samples were ground in liquid nitrogen, and then lysed by PBS and RIPA buffer, respectively. The obtained 1 mg/mL homogenates were used for CKs content determination. Afterward, 10µL supernatant were used for quantitative analyses by using Plant Cytokinin (CTK) Enzyme-Linked Immunoassay Kit (WuHan Saipei Biotechnology, WuHan) to detect plant samples and using Cytokinin (CTK) Enzyme-Linked Immunoassay Kit (WuHan Saipei Biotechnology, WuHan) to detect teliospores samples performed in Multiskan Go (Thermo scientific, USA), according to the manufacturer’s instructions.
Identification of potential cytokinins biosynthetic genes in G. asiaticum and G. yamadae and phylogeny analyses
Based on the previous published data , potential tRNA-isopentenyl transferase genes were identified in telial transcriptomes of G. asiaticum and G. yamadae. TBtools Reciprocal BLAST was then conducted to search for homologous genes in the transcriptome of this study. Combined the results from previous studies [25, 66], we added 1 tRNA-IPT from Cholletotrichum gloeosporioides (hemitropic fungi) and 3 tRNA-IPTs from rust (biotrophic fungi) through a BLASTp search  in NCBI. There were 39 tRNA-IPT protein sequences listed in Additional file 2: Table S7 and the identified tRNA-IPT orthologues of G. asiaticum and G. yamadae were aligned using MAGA 6.0 using Neighbor-joining algorithm. Bootstrap values were determined using 1000 replications. Protein sequences were aligned by multiple alignment with SnapGene 6.1 and visualized using GeneDoc . Localization signals of G. yamadae tRNA-IPT protein were predicted by using LOCALIZER .
RT-qPCR analysis was conducted to validate genes expression profiles between infected leaf samples and healthy leaf samples. Primers were designed by Primer 3  and synthesized by RuiBioTech (Beijing, China). Primer information is presented in Additional file 7: Table S7. The RT-qPCR was performed in CFX96™ SYBR® Real-TimePCR System (Bio-Rad Laboratories, Hercules, CA, USA), as described previously . Transcript expression levels were normalized with ubiquitin-conjugated enzyme E3 . Data logging and the determination of Cq values were done using Bio-Rad CFX Manager 2.1 (Bio-Rad Laboratories, Hercules, CA, USA).
Availability of data and materials
The data that support the findings of this study are available from the corresponding author upon reasonable request. RNA-Seq raw data were obtained from the Sequence Read Archive (https://www.ncbi.nlm.nih.gov/bioproject/) under the accession numbers SRR23018171-SRR23018188.
Differentially expressed genes
Fragment per kilobase per million mapped reads
- 0_GA_J and 0_GY_J:
Juniper healthy branch tissues
Juniper branch tissues infected by G. yamadae at the early stage of gall development without the teliospores formation
Juniper branch tissues infected by G. yamadae the middle stage with the teliospores break through the gall
Juniper branch tissues infected by G. yamadae the late stage of gall development with the matured teliospores
Juniper branch tissues infected by G. asiaticum at the stage with the matured teliospores
Kyoto Encyclopaedia of Genes and Genomes
Principal component analysis
Hiratsuka Y, Kern FD. A revised taxonomic account of Gymnosporangium. Mycologia. 1973;66(1):206.
Cummins GB, Hiratsuka Y. llustrated genera of rust fungi. Minneapolis: Burgess Pub. Co.; 2003.
Li BH, Wang CX, Dong XL, Li HY, Sun LJ, Zhang HH. Research progress in apple diseases and problems in the disease management in China. Plant Prot. 2013;39(2):46–64.
Tao SQ, Cao B, Tian CM, Liang YM. Comparative transcriptome analysis and identification of candidate effectors in two related rust species (Gymnosporangium yamadae and Gymnosporangium asiaticum). BMC Genomics. 2017;18(1):651.
Park SA, Jegal J, Chung KW, Jung HJ, Noh SG, Chung HY, Ahn J, Kim J, Yang MH. Isolation of tyrosinase and melanogenesis inhibitory flavonoids from Juniperus chinensis fruits. Biosci Biotechnol Biochem. 2018;82(12):2041–8.
Khamis AS, Chai LC. Chemical and antimicrobial analyses of Juniperus Chinensis and Juniperus Seravschanica essential oils and comparison with their methanolic crude extracts. Int J Anal Chem. 2021;2021:9937522.
Tao SQ, Cao B, Kakishima M, Liang YM. Species diversity, taxonomy, and phylogeny of Gymnosporangium in China. Mycologia. 2020;112(5):941–73.
Zhao P, Liu F, Li YM, Cai L. Inferring phylogeny and speciation of Gymnosporangium species, and their coevolution with host plants. Sci Rep. 2016;6:29339.
Tao SQ, Cao B, Morin E, Liang YM, Duplessis S. Comparative transcriptomics of Gymnosporangium spp. teliospores reveals a conserved genetic program at this specific stage of the rust fungal life cycle. BMC Genomics. 2019;20(1):723.
Aime MC, Bell CD, Wilson AW. Deconstructing the evolutionary complexity between rust fungi (Pucciniales) and their plant hosts. Stud Mycol. 2018;89:143–52.
Leppik EE. Some viewpoints on the phylogeny of rust fungi. V. Evolution of biological specialization. Mycologia. 1965;57(1):6–22.
Leppik EE. Some viewpoints on the phylogeny of rust fungi VI Biogenic radiation. Mycologia. 1967;59(4):568–79.
Zhao P, Qi XH, Crous PW, Duan WJ, Cai L. Gymnosporangium species on Malus: species delineation, diversity and host alternation. Persoonia. 2020;45:68–100.
Petre B, Lorrain C, Stukenbrock EH, Duplessis S. Host-specialized transcriptome of plant-associated organisms. Curr Opin Plant Biol. 2020;56:81–8.
Duplessis S, Lorrain C, Petre B, Figueroa M, Dodds PN, Aime MC. Host adaptation and virulence in heteroecious rust fungi. Annu Rev Phytopathol. 2021;59(1):403–22.
Takeda S, Yoza M, Amano T, Ohshima I, Hirano T, Sato MH, Sakamoto T, Kimura S. Comparative transcriptome analysis of galls from four different host plants suggests the molecular mechanism of gall development. PLoS ONE. 2019;14(10): e0223686.
Nabity PD, Haus MJ, Berenbaum MR, DeLucia EH. Leaf-galling phylloxera on grapes reprograms host metabolism and morphology. Proc Natl Acad Sci U S A. 2013;110(41):16663–8.
Rolfe SA, Strelkov SE, Links MG, Clarke WE, Robinson SJ, Djavaheri M, Malinowski R, Haddadi P, Kagale S, Parkin IA, et al. The compact genome of the plant pathogen Plasmodiophora brassicae is adapted to intracellular interactions with host Brassica spp. BMC Genomics. 2016;17:272.
Bruce SA, Saville BJ, Neil Emery RJ. Ustilago maydis produces cytokinins and abscisic acid for potential regulation of tumor formation in maize. J Plant Growth Regul. 2011;30(1):51–63.
Chanclud E, Morel JB. Plant hormones: a fungal point of view. Mol Plant Pathol. 2016;17(8):1289–97.
Cisse OH, Almeida JM, Fonseca A, Kumar AA, Salojarvi J, Overmyer K, Hauser PM, Pagni M. Genome sequencing of the plant pathogen Taphrina deformans, the causal agent of peach leaf curl. mBio. 2013;4(3):e00055-00013.
Jahn L, Hofmann U, Ludwig-Muller J. Indole-3-acetic acid is synthesized by the endophyte Cyanodermella asteris via a tryptophan-dependent and -independent way and mediates the interaction with a non-host plant. Int J Mol Sci. 2021;22(5):2651.
Teixeira PJ, Thomazella DP, Reis O, do Prado PF, do Rio MC, Fiorin GL, Jose J, Costa GG, Negri VA, Mondego JM, et al. High-resolution transcript profiling of the atypical biotrophic interaction between Theobroma cacao and the fungal pathogen Moniliophthora perniciosa. Plant Cell. 2014;26(11):4245–69.
Costa JL, Paschoal D, Dasilva EM, Silva JS, Docarmo RM, Carrera E, Lopez-Diaz I, Rossi ML, Freschi L, Mieczkowski P, et al. Moniliophthora perniciosa, the causal agent of witches’ broom disease of cacao, interferes with cytokinin metabolism during infection of Micro-Tom tomato and promotes symptom development. New Phytol. 2021;231(1):365–81.
Morrison EN, Emery RJN, Saville BJ. Fungal derived cytokinins are necessary for normal Ustilago maydis infection of maize. Plant Pathol. 2017;66(5):726–42.
Yamaguchi H, Tanaka H, Hasegawa M, Tokuda M, Asami T, Suzuki Y. Phytohormones and willow gall induction by a gall-inducing sawfly. New Phytol. 2012;196(2):586–95.
Murakami R, Ushima R, Sugimoto R, Tamaoki D, Karahara I, Hanba Y, Wakasugi T, Tsuchida T. A new galling insect model enhances photosynthetic activity in an obligate holoparasitic plant. Sci Rep. 2021;11(1):13013.
Korgaonkar A, Han C, Lemire AL, Siwanowicz I, Bennouna D, Kopec RE, Andolfatto P, Shigenobu S, Stern DL. A novel family of secreted insect proteins linked to plant gall development. Curr Biol. 2021;31(9):1836-1849 e1812.
Takeda S, Hirano T, Ohshima I, Sato MH. Recent progress regarding the molecular aspects of insect gall formation. Int J Mol Sci. 2021;22(17):9424.
Lu S, Sun YH, Amerson H, Chiang VL. MicroRNAs in loblolly pine (Pinus taeda L.) and their association with fusiform rust gall development. Plant J. 2007;51(6):1077–98.
Cambier S, Ginis O, Moreau SJM, Gayral P, Hearn J, Stone GN, Giron D, Huguet E, Drezen JM. Gall wasp transcriptomes unravel potential effectors involved in molecular dialogues with oak and rose. Front Physiol. 2019;10:926.
Ren B, Wang X, Duan J, Ma J. Rhizobial tRNA-derived small RNAs are signal molecules regulating plant nodulation. Science. 2019;365(6456):919–22.
Moffatt K, Flores CM, Andreas P, Kisiala A, Emery RJN. Phytohormone profiling reveals fungal signatures and strong manipulation of infection cycle in the Gymnosporangium juniper-virginianae dual-host plant system. Botany. 2018;96(1):57–65.
Yun HY, Hong SG, Rossman AY, Lee SK, Lee KJ, Bae KS. The rust fungus Gymnosporangium in Korea including two new species G monticola and G unicorne. Mycologia. 2009;101(6):790–809.
Sakakibara H. Cytokinins: activity, biosynthesis, and translocation. Annu Rev Plant Biol. 2006;57:431–49.
Miyawaki K, Tarkowski P, Matsumoto-Kitano M, Kato T, Sato S, Tarkowska D, Tabata S, Sandberg G, Kakimoto T. Roles of Arabidopsis ATP/ADP isopentenyltransferases and tRNA isopentenyltransferases in cytokinin biosynthesis. Proc Natl Acad Sci U S A. 2006;103(44):16598–603.
Morrison EN, Knowles S, Hayward A, Thorn RG, Saville BJ, Emery RJ. Detection of phytohormones in temperate forest fungi predicts consistent abscisic acid production and a common pathway for cytokinin biosynthesis. Mycologia. 2015;107(2):245–57.
Hinsch J, Galuszka P, Tudzynski P. Functional characterization of the first filamentous fungal tRNA-isopentenyltransferase and its role in the virulence of Claviceps purpurea. New Phytol. 2016;211(3):980–92.
Chanclud E, Kisiala A, Emery NR, Chalvon V, Ducasse A, Romiti-Michel C, Gravot A, Kroj T, Morel JB. Cytokinin production by the rice blast fungus is a pivotal requirement for full virulence. PLoS Pathog. 2016;12(2): e1005457.
Hu GG, Rijkenberg FHJ. Ultrastructural localization of cytokinins in Puccinia recondita f.sp. tritici-infected wheat leaves. Physiol Mol Plant Pathol. 1998;52:79–94.
Santos SA, Vidigal PMP, Guimaraes LMS, Mafia RG, Templeton MD, Alfenas AC. Transcriptome analysis of Eucalyptus grandis genotypes reveals constitutive overexpression of genes related to rust (Austropuccinia psidii) resistance. Plant Mol Biol. 2020;104(4–5):339–57.
Duplessis S, Major I, Martin F, Seguin A. Poplar and pathogen interactions: Insights from populus genome-wide analyses of resistance and defense gene families and gene expression profiling. Crit Rev Plant Sci. 2009;28(5):309–34.
Dorostkar S, Dadkhodaie A, Ebrahimie E, Heidari B, Ahmadi-Kordshooli M. Comparative transcriptome analysis of two contrasting resistant and susceptible Aegilops tauschii accessions to wheat leaf rust (Puccinia triticina) using RNA-sequencing. Sci Rep. 2022;12(1):821.
Tang J, Dunshea FR, Suleria HAR. LC-ESI-QTOF/MS characterization of phenolic compounds from medicinal plants (hops and juniper berries) and their antioxidant activity. Foods. 2019;9(1):7.
Saur IML, Huckelhoven R. Recognition and defence of plant-infecting fungal pathogens. J Plant Physiol. 2021;256: 153324.
Zhou JM, Zhang Y. Plant immunity: danger perception and signaling. Cell. 2020;181(5):978–89.
Shih TH, Lin SH, Huang MY, Sun CW, Yang CM. Transcriptome profile of cup-shaped galls in Litsea acuminata leaves. PLoS ONE. 2018;13(10): e0205265.
Nelson N, Yocum CF. Structure and function of photosystems I and II. Annu Rev Plant Biol. 2006;57:521–65.
Garcia-Cerdan JG, Kovacs L, Toth T, Kereiche S, Aseeva E, Boekema EJ, Mamedov F, Funk C, Schroder WP. The PsbW protein stabilizes the supramolecular organization of photosystem II in higher plants. Plant J. 2011;65(3):368–81.
Ilikova I, Ilik P, Opatikova M, Arshad R, Nosek L, Karlicky V, Kucerova Z, Roudnicky P, Pospisil P, Lazar D, et al. Towards spruce-type photosystem II: consequences of the loss of light-harvesting proteins LHCB3 and LHCB6 in Arabidopsis. Plant Physiol. 2021;187(4):2691–715.
Chen YP, Chen T, Zhang XM, Zhang YF, An LZ. The relation of seasonal changes in water and organic osmotica to freezing tolerance in the leaves of Sabina. Bull Botan Res. 2008;28(3):336–41.
Wilkinson S, Davies WJ. ABA-based chemical signalling: the co-ordination of responses to stress in plants. Plant Cell Environ. 2002;25(2):195–210.
Teruko N, Kaori T, Yukiko K, Tadako M. Effect of auxin and gibberellin on conidial germination in Neurospora crassa ii. “Conidial density effect” and auxin. Plant Cell Physiol. 1982;8:1363–9.
Teruko N, Yukiko K, Eiko T, Noriko T, Tadako M. Effects of auxin and gibberellin on conidial germination in Neurospora crassa. Plant Cell Physiol. 1978;4:705–9.
Doehlemann G, Wahl R, Horst RJ, Voll LM, Usadel B, Poree F, Stitt M, Pons-Kuhnemann J, Sonnewald U, Kahmann R, et al. Reprogramming a maize plant: transcriptional and metabolic changes induced by the fungal biotroph Ustilago maydis. Plant J. 2008;56(2):181–95.
Adhikari TB, Balaji B, Breeden J, Goodwin SB. Resistance of wheat to Mycosphaerella graminicola involves early and late peaks of gene expression. Physiol Mol Plant Pathol. 2007;71(1–3):55–68.
Martinson EO, Werren JH, Egan SP. Tissue-specific gene expression shows a cynipid wasp repurposes oak host gene networks to create a complex and novel parasite-specific organ. Mol Ecol. 2022;31(11):3228–40.
Schultz JC, Stone GN. A tale of two tissues: Probing gene expression in a complex insect-induced gall. Mol Ecol. 2022;31(11):3031–4.
Liesche J, Vincent C, Han X, Zwieniecki M, Schulz A, Gao C, Bravard R, Marker S, Bohr T. The mechanism of sugar export from long conifer needles. New Phytol. 2021;230(5):1911–24.
Stull GW, Qu XJ, Parins-Fukuchi C, Yang YY, Yang JB, Yang ZY, Hu Y, Ma H, Soltis PS, Soltis DE, et al. Gene duplications and phylogenomic conflict underlie major pulses of phenotypic evolution in gymnosperms. Nat Plants. 2021;7(8):1015–25.
De La Torre AR, Piot A, Liu B, Wilhite B, Weiss M, Porth I. Functional and morphological evolution in gymnosperms: a portrait of implicated gene families. Evol Appl. 2020;13(1):210–27.
Han X, Turgeon R, Schulz A, Liesche J. Environmental conditions, not sugar export efficiency, limit the length of conifer leaves. Tree Physiol. 2019;39(2):312–9.
Hu XG, Liu H, Zhang JQ, Sun YQ, Jin YQ, Zhao W, El-Kassaby YA, Wang XR, Mao JF. Global transcriptome analysis of Sabina chinensis (Cupressaceae), a valuable reforestation conifer. Mol Breeding. 2016;36(7):99.
Canovas FM, Avila C, Canton FR, Canas RA, de la Torre F. Ammonium assimilation and amino acid metabolism in conifers. J Exp Bot. 2007;58(9):2307–18.
Walters DR, McRoberts N. Plants and biotrophs: a pivotal role for cytokinins? Trends Plant Sci. 2006;11(12):581–6.
Kakimoto T. Biosynthesis of cytokinins. J Plant Res. 2003;116(3):233–9.
Smaldino PJ, Read DF, Pratt-Hyatt M, Hopper AK, Engelke DR. The cytoplasmic and nuclear populations of the eukaryote tRNA-isopentenyl transferase have distinct functions with implications in human cancer. Gene. 2015;556(1):13–8.
Kasahara H, Takei K, Ueda N, Hishiyama S, Yamaya T, Kamiya Y, Yamaguchi S, Sakakibara H. Distinct isoprenoid origins of cis- and trans-zeatin biosyntheses in Arabidopsis. J Biol Chem. 2004;279(14):14049–54.
Beckerson WC, Rodriguez de la Vega RC, Hartmann FE, Duhamel M, Giraud T, Perlin MH. Cause and effectors: Whole-genome comparisons reveal shared but rapidly evolving effector sets among host-specific plant-castrating fungi. mBio. 2019;10(6):2391–2319.
Zhao J, Duan W, Xu Y, Zhang C, Wang L, Wang J, Tian S, Pei G, Zhan G, Zhuang H, et al. Distinct transcriptomic reprogramming in the wheat stripe rust fungus during the initial infection of wheat and barberry. Mol Plant Microbe Interact. 2021;34(2):198–209.
Borah N, Albarouki E, Schirawski J. Comparative methods for molecular determination of host-specificity factors in plant-pathogenic fungi. Int J Mol Sci. 2018;19(3):863.
Hatta R, Ito K, Hosaki Y, Tanaka T, Tanaka A, Yamamoto M, Akimitsu K, Tsuge T. A conditionally dispensable chromosome controls host-specific pathogenicity in the fungal plant pathogen Alternaria alternata. Genetics. 2002;161(1):59–70.
Lofgren LA, Nguyen NH, Vilgalys R, Ruytinx J, Liao HL, Branco S, Kuo A, LaButti K, Lipzen A, Andreopoulos W, et al. Comparative genomics reveals dynamic genome evolution in host specialist ectomycorrhizal fungi. New Phytol. 2021;230(2):774–92.
Liu JJ, Sturrock RN, Sniezko RA, Williams H, Benton R, Zamany A. Transcriptome analysis of the white pine blister rust pathogen Cronartium ribicola: de novo assembly, expression profiling, and identification of candidate effectors. BMC Genomics. 2015;16:678.
Hacquard S, Petre B, Frey P, Hecker A, Rouhier N, Duplessis S. The poplar-poplar rust interaction: insights from genomics and transcriptomics. J Pathog. 2011;2011: 716041.
Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–90.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52.
Freedman AH, Clamp M, Sackton TB. Error, noise and bias in de novo transcriptome assemblies. Mol Ecol Resour. 2021;21(1):18–29.
Fu L, Niu B, Zhu Z, Wu S, Li W. CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics. 2012;28(23):3150–2.
Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215(3):403–10.
Huerta-Cepas J, Forslund K, Coelho LP, Szklarczyk D, Jensen LJ, von Mering C, Bork P. Fast genome-wide functional annotation through orthology assignment by eggNOG-Mapper. Mol Biol Evol. 2017;34(8):2115–22.
Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2015;12(1):59–60.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
Chen C, Chen H, Zhang Y, Thomas HR, Frank MH, He Y, Xia R. Tbtools: an integrative toolkit developed for interactive analyses of big biological data. Mol Plant. 2020;13(8):1194–202.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Schwacke R, Ponce-Soto GY, Krause K, Bolger AM, Arsova B, Hallab A, Gruden K, Stitt M, Bolger ME, Usadel B. MapMan4: A refined protein classification and annotation framework applicable to multi-omics data analysis. Mol Plant. 2019;12(6):879–92.
Tao SQ, Auer L, Morin E, Liang YM, Duplessis S. Transcriptome analysis of apple leaves infected by the rust fungus Gymnosporangium yamadae at two sporulation stages. Mol Plant Microbe Interact. 2020;33(3):444–61.
Nicholas K, Nicholas H. GeneDoc: a tool for editing and annotating multiple sequence alignments. 1997. http://www.pscedu/biomed/genedoc.
Sperschneider J, Catanzariti AM, DeBoer K, Petre B, Gardiner DM, Singh KB, Dodds PN, Taylor JM. LOCALIZER: subcellular localization prediction of both plant and effector proteins in the plant cell. Sci Rep. 2017;7:44598.
Untergasser A, Cutcutache I, Koressaar T, Ye J, Faircloth BC, Remm M, Rozen SG. Primer3–new capabilities and interfaces. Nucleic Acids Res. 2012;40(15): e115.
Link TI. Testing reference genes for transcript profiling in Uromyces appendiculatus during urediospore infection of common bean. PLoS ONE. 2020;15(8): e0237273.
Shao C, Lao W, Liang Y. Reference genes selection of Gymnosporangium yamadae during the interaction with apple leaves. J Fungi (Basel). 2022;8(8):830.
The author thanks the lab members for assistance and anonymous reviewers for their thorough and constructive comments on the manuscript.
This work was supported by the National Natural Science Foundation of China (31870628). The funding bodies did not contribute intellectually to conceiving of the idea, planning of experiments, interpreting the data or writing the manuscript.
Ethics approval and consent to participate
All the samples of junipers (Juniperus chinensis L.) with rust pathogens G. asiaticum were obtained from Northwest A&F University, Yangling, Shaanxi, China. All the collecting procedures were carried out under the permission and assistance of Northwest A&F University. All the samples of junipers (Juniperus chinensis L.) with rust pathogens G. yamadae were obtained from Haidian park, Beijing, China with the permission of Haidian park. All methods were performed in accordance with the relevant guidelines/regulations/legislation.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure S1. The pipeline of RNA-seq and bioinformatic analysis used in the study. The plant unigenes were mapped to five publicly available databases of Pinus lambertiana, P.s taeda, Picea abies, P. gluca, Gnetum montanum in National Center for Biotechnology Information. The rust fungi unigenes were annotated as described in the methods. Figure S2. The proportional distribution by species of Juniperus chinensis (a) and Gymnosporangium species (b) unigenes with homology in NR database. Figure S3. Assessment of RNA-seq data reproducibility. (a), Principal component analysis based on gene expression level of healthy leaves and branches (0_GA_J and 0_GY_J), Gymnosporangium asiaticum infected leaves with the the matured teliospores (III_GA_J), G. yamadae infected branch tissues at the early (I_GY_J), middle (II_GY_J) and late stage (III_GY_J) of the gall development, showing the clear separation of the six tested samples and the proximity of biological replicates. (b), Sample correlation matrix of replicates of healthy leaves and branches (0_GA_J and 0_GY_J), G. asiaticum infected leaves with the the matured teliospores (III_GA_J), G. yamadae infected branch tissues at the early (I_GY_J), middle (II_GY_J) and late stage (III_GY_J) of the gall development based on gene expression. Figure S4. Statistics of annotated different expression genes. Venn diagram showing the number of up- and down-regulated genes (|log2Fold-Change| > 1, p-value < 0.05) in III_GA _J and III_GY_J samples compared with the control, respectively. Type_a and Type_b refer to juniper unigenes regulated reversely in III_GA _J and III_GY_J samples. Figure S5. Amino acid sequence alignment of tRNA-IPT homologs. The tRNA-IPTs from Saccharomyces cerevisiae (NP_014917.3), Claviceps purpurea (CCE29200.1), Sporisorium reilianum (CBQ67583.1), Ustilago maydis (XP_011386632.1), Gymnosporangium asiaticum (c25847_g1) and Gymnosporangium yamadae (TRINITY_DN3125 c0_gI_i1). The dark to light colors indicate amino acids conservation from high to low. Black lines indicate ATG/GTP binding site, DMAPP binding site and zinc-finger motif; the two putative nuclear localization signals of GytRNA-IPT (TRINITY_DN3125 c0_gI_i1) and GatRNA-IPT (c25847_g1) are highlighted by green and red boxes; black box indicate chloroplast localization signal of GytRNA-IPT.
Table S1. The statistics of unigene classification. Unigenes assigned with “juniper” species were determined as juniper unigenes, whie assigned with “fungi” species were determined as the two Gymnosporangium spp. unigenes. Table S2. Functional annotation of Juniperus chinensis and the Gymnosporangium spp. unigenes in selected six public databases (Nr, KEGG, GO, Swissprot, KOG, PFAM). Table S3. Differentially expressed genes in III_GA, I_GY, II_GY and III_GY. Table S4. Swissport annotations of Type_a and Type_b DEGs in the Venn diagram. Table S5. MapMan annotations of different expression genes (|log2Fold-Change| ≥ 5, p <0.01) in III_GA and III_GY samples. Table S6. FPKM and KO number of photosynthesis, sugar metabolism, plant hormone and defense-related genes in I_GY, II_GY and III_GY samples. Table S7. The tRNA-isopentenyltransferase proteins used for sequence alignment and phylogenetic tree creation.
About this article
Cite this article
Shao, C., Tao, S. & Liang, Y. Comparative transcriptome analysis of juniper branches infected by Gymnosporangium spp. highlights their different infection strategies associated with cytokinins. BMC Genomics 24, 173 (2023). https://doi.org/10.1186/s12864-023-09276-7
- Gymnosporangium asiaticum
- Gymnosporangium yamadae