Genomic and culture-based analysis of Cyclaneusma minus in New Zealand provides evidence for multiple morphotypes

Abstract

Cyclaneusma needle cast, caused by Cyclaneusma minus, affects Pinus species world wide. Previous studies suggested the presence of two distinct morphotypes in New Zealand, ‘verum’ and ‘simile’. Traditional mycological analyses revealed a third morphotype with clear differences in colony morphology and cardinal growth rates at varying temperatures. Genome sequencing of eight C. minus isolates provided further evidence of the existence of a third morphotype, named as ‘novus’ in this study. To further analyse these morphotypes, we predicted candidate effector proteins for all eight isolates, and also characterized a cell-death eliciting effector family, Ecp32, which is present in other pine phytopathogens. In concordance with their distinct classification into three different morphotypes, the number of Ecp32 family members differed, with patterns of pseudogenization in the ‘simile’ morphotype, and some members being found exclusively either in the ‘simile’ or ‘verum’ morphotypes. We also showed that the Ecp32 family proteins trigger cell death in non-host Nicotiana species, and, as previously demonstrated in other plant pathogens, the Ecp32 family proteins in C. minus adopt a β-trefoil fold. These analyses provide further evidence that the three morphotypes might be distinct species that need formal descriptions. Understanding the geographical range of different Cyclaneusma species and variations in virulence and pathogenicity will provide a better understanding of pine needle diseases and enable the development of more durable methods to control this disease.

Background

Cyclaneusma needle cast (CNC) is a disease of Pinus species in many parts of the world (Bednářová et al. 2013). CNC affects more than 30 pine species, including Pinus radiata in New Zealand and Australia, while in forests of Europe and North America, Scots pine (Pinus sylvestris) is prone to this disease (Merrill and Wenner 1996; Müller & Kurkela 2012). CNC results in premature needle loss, and therefore losses in growth and wood production. For New Zealand’s plantation forestry, the losses caused by CNC were estimated to be around $38 million per year (Bulman 2009). Fortunately, this impact has likely now been reduced through the planting of less susceptible genotypes. Regardless, the disease persists in some locations, and there are no alternative control measures (Bulman 1993).

Cyclaneusma needle cast is caused by the Leotiomycete fungus Cyclaneusma minus (Butin) DiCosmo, Peredo and Minter. C. minus was initially reported as Naemacyclus niveus (Gilmour 1959), which was later separated into two species, N. niveus and N. minor (Von Butin 1973). This separation was based on morphological characters on needles, the formation of the sexual and asexual stages in culture, and the relationship with different host species. Naemacyclus niveus and N. minor were later reassigned to the genus Cyclaneusma as Cyclaneusma niveum and C. minus, respectively (Dicosmo et al. 1983). The morphological features distinguishing C. niveum and C. minus were described in Von Butin (1973) and again in Dick et al. (2001), where C. niveum isolates had sickle-shaped pycnidiospores and much shorter apothecia, and C. minus had bacilliform-shaped pycnidiospores and longer apothecia. Further investigation of the C. minus isolates in New Zealand found that, based on apothecial length variation as well as fungal morphology in culture, C. minus should be separated into two morphotypes: tagged as C. minus ‘verum’ and C. minus ‘simile’ (Dick et al. 2001). A relationship was shown between apothecium dimensions from needles and these morphotypes. A multigene phylogenetic analysis of New Zealand and Australian isolates of Cyclaneusma demonstrated that the separation of isolates by morphotype was consistent with genetic analysis and suggested the need for species descriptions to formally re-classify, as well as pathogenicity studies to understand phenotypic differences between the morphotypes (Prihatini et al. 2014).

Hunter et al. (2016) developed molecular tools for differentiating the ‘verum’ and ‘simile’ morphotypes and looked at their distribution throughout New Zealand. The study found that ‘simile’ was the most common morphotype in New Zealand, and that molecular classification was much more accurate than morphological assessments. Amongst the isolates studied, one (NZFS3325) showed DNA sequence variability in the Internal Transcribed Spacer (ITS) region of its genome, when compared to the corresponding region in both the ‘simile’ and ‘verum’, which indicated that additional morphotypes might exist (Hunter et al. 2016). Rightfully identifying these different morphotypes would help to determine their current distribution as well as to predict future distribution under different climatic scenarios, in New Zealand and globally. Comparative studies between morphotypes can also determine whether a morphotype is more virulent than others and can have implications for disease diagnostics and forest management.

To further characterize the morphotypes of C. minus, we carried out genome sequencing for eight isolates: two belonging to morphotype ‘simile’ (Cms; NZFS759 and NZFS809), four belonging to morphotype ‘verum’ (Cmv; NZFS110, NZFS725, NZFS1800, and NZFS3276), and two isolates (Cmn; NZFS3305 and NZFS3325) that could not be differentiated into ‘verum’ or ‘simile’ using previously developed molecular tools. This genomic information will provide additional data to assist taxonomic classification. An additional benefit of genomic resources is that they can provide insights into the virulence on their hosts, through the prediction of candidate effector (CE) proteins. Effectors are usually small proteins secreted by pathogens into and around host cells to facilitate invasion, for example by suppressing host defence responses (Lo Presti et al. 2015; Rocafort et al. 2020). If the host has the ability to recognize such effectors, through its extracellular and intracellular resistance (R) proteins, immune responses are activated in order to stop pathogen infection (Jones and Dangl 2006; Cook et al. 2015). This recognition can lead to a hypersensitive response, which is generally characterized by a localized plant cell death (Win et al. 2012). Notably, pathogen effectors can be heterologously expressed in host and non-host plants, with those triggering a plant cell death response being candidates for recognition by a plant R protein.

In this study, genome sequencing was performed for eight C. minus isolates with the main objective of confirming if a third morphotype is present in New Zealand, and its relationship to other morphotypes. Moreover, the genome sequences of the different C. minus morphotypes were examined for genes encoding CE proteins. Here, we focused on homologues of the cell death-eliciting Ecp32 protein family, firstly identified in the tomato pathogen Fulvia fulva and subsequently in the pine pathogen Dothistroma septosporum (Mesarich et al. 2018; Tarallo et al. 2022). We compared their sequences, predicted tertiary structures, and screened these homologues for their ability to trigger cell death in Nicotiana benthamiana and N. tabacum using Agrobacterium tumefaciens-mediated transient transformation assays (ATTAs). By doing so, we aimed to provide novel insights into whether the C. minus morphotypes differ with respect to the conserved Ecp32 family, and to highlight potential differences in the molecular basis for how the different morphotypes interact with their hosts.

Results

C. minus isolates can be differentiated into three morphotypes based on molecular and morphological analyses

Eight C. minus isolates were collected from different parts of New Zealand and classified into two different morphotypes, ‘verum’ (Cmv) and ‘simile’ (Cms), according to molecular analysis (Table 1 and Additional file 1: Figure S1) (Hunter et al. 2016). However, two of these isolates, NZFS3305 and NZFS3325 (NZFS3325 was previously identified in Hunter et al. (2016)), could not be differentiated based on this molecular analysis and were, therefore, classified as a new C. minus morphotype, ‘novus’ (Cmn). The phylogeny suggested that morphotype Cmn is more closely related to Cms than to Cmv (Additional file 1: Figure S1).

Table 1 Isolates of Cyclaneusma minus used in this study

This classification was also supported by morphological analyses. Results for the temperature study are shown in Fig. 1. The widest temperature range for growth, as well as the fastest growth rate, was seen in the Cmv morphotype. It had an optimal growth at 25 °C, with a mean radial growth rate of 2.0 ± 0.7 mm/day at 25 °C. However, growth at our lowest (4°C) and highest (30°C) temperatures were also detected for this morphotype. The Cms and Cmn morphotypes shared similar optimal growth temperatures at 20°C and 22°C, respectively, while growth ceased at 4°C and 30°C for both. The Cms morphotype had a slightly faster mean radial growth rate of 0.96 ± 0.13 mm/day at its optimal growth temperature of 20°C, compared to that of the Cmn morphotype being 0.69 ± 0.15 mm/day at 22°C. Analysis of temperature study data using two-way ANOVA with Tukey’s post hoc showed that the radial growth per day of Cmv was statistically different from Cms and Cmn (p = 1.8 × 10⁶ and 2.9 × 10⁵ respectively), while there was no significant difference in the radial growth between Cms and Cmv (p = 0.65). Isolates and their corresponding morphotypes grown on 2% and 9% MEA at 22°C in the dark after 8 weeks are shown in Fig. 2. The colony morphology of each morphotype on 2% and 9% MEA at 22°C in the dark after 8 weeks was displayed in Additional file 2: Table S1.

Fig. 1

Effect of temperature on the growth of C. minus morphotypes. The ‘simile’ (blue), ‘verum’ (red), and ‘novus’ (green) morphotypes were assessed using mean radial growth, measured in millimetres per day with standard deviations

Fig. 2

Colony morphology of C. minus isolates. All isolates were grown at 22°C in the dark for 8 weeks

Differences in C. minus morphotype are supported by differences in genome size

The final genome assemblies for the eight C. minus isolates ranged in length from 31 to 43 Mb (Table 2 and Additional file 2: Table S2). The final set of gene models for each genome showed a completeness of > 96% with the retrieval of at least 728 complete BUSCO hits from a total of 758 fungal BUSCOs (Table 2 and Additional file 2: Table S2). The Cmv morphotype was predicted to have a larger genome size (based on assembly length) and higher number of gene models compared to the Cms and Cmn morphotypes (Table 2).

Table 2 Cyclaneusma minus genomes statistics

Prediction of candidate effector proteins from C. minus morphotypes

Effectors can have virulence functions during host infection by supressing plant defence responses or can be avirulence factors, in which case they are recognized by R proteins from the hosts. For this reason, following genome sequencing and assembly, secreted CE proteins were identified from the predicted protein set of all eight C. minus isolates to determine whether their CE protein repertoires were morphotype-specific. Strikingly, genomes from the same morphotype were predicted to have similar numbers of proteins and effectors, with Cmv having more predicted CEs than Cms or Cmn (Table 3). In all cases, a higher number of apoplastic effectors, when compared to cytoplasmic effectors, were predicted in all isolates analysed. Additional file 2: Table S3 shows all amino acid sequences of the CEs predicted from the isolates by using EffectorP.

Table 3 Effector protein prediction statistics across the eight C. minus isolates investigated in this study

The number of Ecp32 effector family members differs between the three morphotypes of C. minus, further supporting their differentiation

To investigate potential differences in the effector protein repertoires of C. minus morphotypes, we focused on the Ecp32 effector family that had been characterised in another fungal pine pathogen (Tarallo et al. 2022). An apoplastic protein from isolate NZFS809, Cms835, was firstly identified as an orthologue of DsEcp32-1 from D. septosporum and shown to be part of the Ecp32 effector family. With the Cms835 (NZFS809) amino acid sequence as reference, other Ecp32 family members in this and other isolates/morphotypes were identified using BLASTp and tBLASTn, respectively. Based on this analysis, it was determined that all isolates of all morphotypes had at least three members of the Ecp32 family of proteins (Fig. 3 and Additional file 2: Table S4).

Fig. 3

Schematic diagrams of the Ecp32 family from three different morphotypes of C. minus. One representative isolate from each morphotype was chosen: NZFS110 for morphotype ‘verum’ (Cmv), NZFS809 for morphotype ‘simile’ (Cms), and NZFS3305 for morphotype ‘novus’ (Cmn). Each rectangle represents a gene and orthologous genes between morphotypes are colour-coded in the same way as Fig. 4. Ψ refers to pseudogenes. The horizontal lines connecting the rectangles indicate the genes are on the same scaffold and adjacent in the genome. The percentages refer to full-length pairwise amino acid identity between orthologous proteins encoded by each gene from each representative morphotype (V = ‘verum’, S = ‘simile’, N = ‘novus’). The pink rectangle indicates the gene encoding the protein with a predicted intrinsically disordered region. An asterisk (*) next to the gene number indicates amino acid variants in the respective encoded proteins that occurred in one or more of the other Cmv isolates. NZFS3276 orthologues of 9109, 6323, and 2795 have variants I34T, V188I, and E40D, respectively. Orthologues of the other polymorphic protein showed that 5316 has a Q82H substitution in all three other Cmv isolates and a S215G substitution in isolates NZFS1800 and 3276

Although all isolates and morphotypes appeared to have proteins belonging to the Ecp32 family, the number of family members between each morphotype varied. A BLASTp search indicated that both the NZFS759 and NZFS809 genomes (both Cms) encode eight Ecp32 family proteins. Each of the Cmv isolates (NZFS110, NZFS725, NZFS1800, and NZFS3276) have six paralogues in their Ecp32 family, while Cmn isolates (NZFS3305 and NZFS3325) have only three. The Ecp32 genes from one representative isolate of each morphotype are illustrated in Fig. 3. Unlike what was observed for orthologues between D. septosporum and F. fulva, proteins from the Ecp32 family from each C. minus isolate did not group with the previously identified members of the Ds/FfEcp32 family (Additional file 1: Figure S2). As was also found for the D. septosporum and F. fulva Ecp32 proteins, all Ecp32 family members from C. minus have predicted signal peptides but no known functional domains.

Some Ecp32 family members are found exclusively in one morphotype. This is the case for a set of orthologous proteins (Cmv2795, Cmv1388, Cmv116, and Cmv3467) from the four Cmv isolates. These proteins have longer amino acid sequences than the rest of the proteins in this family. Normally, they range from 197 to 226 amino acids in length, however these four proteins each contain 268 amino acids. The additional amino acids form a predicted intrinsically disordered region (IDR) at the N-terminus, immediately following the putative signal peptide (Additional file 1: Figure S3).

The morphotypes also showed differences in pseudogenisation of Ecp32 genes. For Cms isolates NZFS809 and NZFS759, three of the eight Ecp32 genes were pseudogenes that were predicted to encode truncated proteins, with identical mutations in the homologues from these two isolates (NZFS809 is shown in Fig. 3). Two of the genes (Cms834 and Cms4060) had a premature stop codon at the position encoding amino acids 54 and 149, respectively, whilst Cms2573 had a cytosine insertion at nt position 213 in the coding sequence, causing a frameshift that also led to a premature stop codon. Interestingly, in both Cms isolates, Cms834 (and the NZFS759 orthologue) is on the same scaffold and adjacent to Cms835 (and the NZFS759 orthologue) (Fig. 3). Similarly, another gene cluster in this morphotype includes the Cms4060 pseudogene, alongside Cms4058 and Cms4059, which are predicted to encode functional proteins (Fig. 3). This is also true for the respective orthologues from NZFS759.

Given that the Cms and Cmn morphotypes are the most closely-related of the three morphotypes based on the multigene phylogeny generated in this study (Additional file 1: Figure S1), the difference in the number of Ecp32 gene family members between them was striking, although three of the eight Cms genes were pseudogenised (Fig. 3). Furthermore, the three genes that shared orthologues between Cms and Cmn showed sequence diversity, with predicted amino acid sequence identities of only 47, 65, and 74% between these morphotypes (Fig. 3). Together, these results support the identification of Cmn as a distinct morphotype from Cms and Cmv, and also suggest that several gene duplication events and possibly gene loss may have occurred in this effector family, either in C. minus or in an ancestral species.

The Ecp32 family from different C. minus morphotypes are sequence-divergent and some members are found exclusively in certain morphotypes

To further investigate differences between the sets of Ecp32 family members across the morphotypes and isolates of C. minus, we performed a protein sequence alignment of all full-length members from the Ecp32 family (Additional File 1: Figure S3). Family members that had identical amino acid sequences were excluded from this alignment, as were sequences encoded by pseudogenes. Figure 4 shows a phylogenetic tree constructed based on these alignments, including all family members identified, even those with identical amino acid sequences, but still excluding those encoded by pseudogenes.

Fig. 4

Phylogenetic tree of C. minus Ecp32 family members across eight isolates and three morphotyes of C. minus. Neighbour-joining tree obtained from the amino acid sequence alignment of Ecp32 proteins using Geneious Software v9.1.8. Bootstrap values are shown at the nodes as percentages. The scale bar represents 0.2 substitutions per site. The number following each morphotype ID (NZFS) refers to the protein number from the annotated genomes. The asterisk (*) indicates that the protein was incorrectly predicted, and the corrected protein sequence was used. The sequences encoded by predicted pseudogenes were excluded from the analysis. The black dot refers to the proteins that were expressed in N. benthamiana and N. tabacum using an Agrobacterium tumefaciens-mediated transient expression assay to assess their ability to elicit cell death. Cms: C. minus morphotype ‘simile’; Cmv: C. minus morphotype ‘verum’; Cmn: C. minus morphotype ‘novus’. Proteins were colour-coded according to their clustering in the phylogeny

The results affirm that proteins from the Ecp32 family are more expanded in some morphotypes than in others, and they cluster into seven separate groups (colour-coded to match Fig. 3). The clustered groups coloured in blue, yellow, and green are made up of eight protein members, with one paralogue belonging to each of the eight genomes analysed. The clustered group colour-coded in red has six members, with paralogues only belonging to Cms and Cmv, and no members belonging to Cmn. The pink group corresponds to Cmv proteins with a predicted disordered region of 73 additional amino acids, which are most closely related to green group proteins, despite their larger size.

A high level of sequence divergence was seen between the Ecp32 proteins from the three morphotypes. Between all orthologous Cms and Cmv Ecp32 protein family members, the pairwise amino acid identity varied between 25.1% and 71.4%. Likewise, between Cms and Cmn, the pairwise identity was between 32.1% and 68.9%, and between Cmv and Cmn it was 28.8% and 68.6%. Differences were also observed between isolates of the same morphotype, as the amino acid sequence identity among the Ecp32 proteins ranged from 25.3–50.6% (Cmv), 30.9–42.3% (Cms), and 34.4–46.8% (Cmn).

We next compared the locations of clustered Ecp32 genes with their positions in the phylogeny. All genes encoding Ecp32 proteins of the red and grey groups in the Cmv morphotypes are adjacent to each other in their respective genomes, as shown for Cmv5316 and Cmv5317 of isolate NZFS110 in Fig. 3; the corresponding proteins are more similar in sequence to each other than to any other Ecp32 protein from that isolate, with identities ranging from 26.8 to 50.6% (between Cmv5316 and Cmv5317) (Additional file 1: Figure S3). Together, these results suggest that these ‘red and grey’ family members may have arisen by an ancestral tandem duplication event. A similar linkage of these two family members was seen in the Cms morphotype (genes Cms4059 and Cms4060), although in that case the ‘grey’ Cms4060 gene was pseudogenised (Fig. 3). Although not shown in the tree, Cms4060 showed greatest similarity with the proteins colour-coded in grey (Additional file 1: Figure S3), which suggests a common ancestral origin for the red and grey-coloured gene pairs of the Cmv and Cms morphotypes. A third Ecp32 gene was also present in the red-grey gene cluster in the Cms morphotype: Cms4058 (orange) (Fig. 3). However, despite its close proximity to Cms4059 in the genome, Cms4058 was more similar in sequence to Cms835 (green) than to Cms4059 (red).

Together, these results highlight the expansion of the Ecp32 family in the three different morphotypes analysed, with possible gene duplication and the generation of new family members. In addition, sequence diversification was observed, not only between morphotypes, but also between family members from the same isolate.

Proteins from the Ecp32 family of different C. minus morphotypes trigger cell death in Nicotiana species

The next questions we addressed were whether the Ecp32 family proteins identified from the C. minus morphotypes could trigger cell death in Nicotiana species, similar to what was observed in D. septosporum and F. fulva (Tarallo et al. 2022), and whether there were differences in cell death-eliciting ability between morphotypes. The selected proteins tested by ATTAs are indicated with a black dot in Fig. 4. One protein from one representative isolate of each morphotype was selected from the clusters that had a higher number of family members.

In N. benthamiana, all proteins tested triggered certain degrees of cell death, except for Cmv937 (the Cmv representative from the blue cluster in Fig. 4) (Fig. 5). Among those tested proteins, proteins that consistently triggered cell death were those from the green cluster shown in Fig. 4 (Cmv6323, Cms835, and Cmn5237), and the yellow cluster (Cmv9109, Cms9849, and Cmn3003), with Cms9849 triggering a strong cell death response in all infiltration zones. Across all proteins tested, differences were seen between morphotypes (Fig. 5). The truncated protein Cms4060 (NZFS809) was also tested, which triggered a weak cell death response. Although this protein was not included in the phylogenetic tree shown in Fig. 4, it showed greatest similarity with the proteins colour-coded in grey, and an alignment between these proteins (Additional file 1: Figure S3) indicated conservation of the N-terminus (post signal peptide), which suggests that this region of the protein might be responsible for cell death elicitation, since Cmv5316 also triggered a cell death (Fig. 5).

Fig. 5

Ecp32 proteins from different morphotypes of C. minus triggered cell death in Nicotiana benthamiana. a Proteins were expressed in N. benthamiana using an Agrobacterium tumefaciens-mediated transient expression assay to assess their ability to elicit cell death. Representative images are shown (n = 12–30 infiltration zones) from at least three independent experiments. INF1, Phytophthora infestans elicitin positive cell death control; EV, empty vector, a negative no-cell death control; Cms: C. minus morphotype ‘simile’; Cmv: C. minus morphotype ‘verum’; Cmn: C. minus morphotype ‘novus’. Photos were taken 7 days after infiltration. b Graphs display the percentages of infiltration zones across the entire experiment that showed cell death in response to each different protein, divided into four categories: strong cell death, weak cell death, chlorosis, and no cell death. The proteins are shown as orthologous groups, with the different morphotypes represented for each. The coloured lines under the protein IDs refer to the different gene and protein groups in Figs. 3 and 4

In N. tabacum, yellow cluster proteins Cmv9109, Cms9849, and Cmn3003 (Fig. 4) again consistently triggered strong cell death, while all other proteins inconsistently triggered either weak cell death or chlorosis (Additional file 1: Figure S4). Interestingly, the truncated protein Cms4060 (NZFS809), as well as its paralogue Cmv5316, did not trigger cell death (Additional File 1: Figure S4).

Taken together, these results suggest that different Ecp32 family members from the same isolate differ in their ability to trigger a cell death response in non-host plants. There were also differences in the responses observed for orthologues from different morphotypes (e.g., proteins from the blue cluster).

Predicted tertiary structures of the C. minus Ecp32 family members suggest possible roles in virulence

The tertiary structure of one C. minus Ecp32 family member from each of the clusters shown in Fig. 4 was predicted using AlphaFold2 (Additional file 1: Figure S5). The amino acid sequences of proteins from NZFS110 (Cmv) were used for these structural predictions since one Ecp32 family member from this isolate was present in all clusters, except for the orange one (Fig. 4). All of the structures had high pLDDT scores and TM-scores, indicative of highly confident predictions (Additional file 2: Table S5). As previously demonstrated for the Ecp32 proteins of D. septosporum and F. fulva, the Ecp32 family members from C. minus were predicted to be structurally similar to β-trefoil proteins, even though they had some differences in topology. The structures of Cmv937, Cmv6323, Cmv9109, Cmv5317, and Cmv2795 have two predicted disulphide bonds, while Cmv5316 has three.

Pairwise structural alignments between all of the predicted structures suggested that Cmv6323 (green group) and Cmv9109 (yellow group) are more similar to one another than any other structural pairwise alignment from the predicted structures for Cmv isolate NZFS110. Interestingly, Cmv6323 and Cmv9109 are part of the two groups that triggered the most consistent and strong cell death responses in N. benthamiana. This similarity, however, was not observed in the amino acid pairwise alignment of these proteins, for which the pair Cmv5316 and Cmv5317 have a higher percentage of identity than the Cmv6323 and Cmv9109 pair. The two predicted disulphide bonds shared across all predicted structures are conserved with the two predicted disulphide bonds present in all members of the Ecp32 family from D. septosporum and F. fulva (except Ds/FfEcp32-5, which does not have any cysteine residues), suggesting that this is a conserved feature in proteins belonging to this family (Additional file 1: Figure S5h) (Tarallo et al. 2022).

The structure of the truncated Ecp32 candidate Cms4060 was also predicted since the protein triggered weak cell death in the same manner as Cmv5316 (Fig. 5). Cms4060 has 77 amino acids less than the 226 expected for Cmv5316. The alignment between the predicted structures of Cms4060 and Cmv5316 (Additional file 1: Figure S5i) highlights the N-terminal conservation that was previously observed in the sequence alignment (Additional file 1: Figure S3). The predicted structure of Cms4060 (coloured green in the structural alignment) aligns with the region of Cmv5316 coloured in black, including one conserved disulphide bond (S-S1). The grey colouring represents the C-terminal region of Cmv5316 that is not present in Cms4060. Similarities in the structures of truncated protein Cms4060 and full-length Cmv5316, with conservation of cell death inducing activity, might provide insights into the regions of the proteins from the Ecp32 family that are responsible for plant receptor-mediated recognition.

Discussion

Cyclaneusma minus is an important pathogen of pines worldwide. However, very little is known about the biology of this pathogen. At first, differences in ascospore length and characteristics of the fungus grown on agar plates indicated that there were two morphological types of C. minus in New Zealand, termed C. minus ‘simile’ and C. minus ‘verum’ (Dick et al. 2001). However, with the development of molecular tools to differentiate these morphotypes, another morphotype was found to be present in the New Zealand C. minus population (Hunter et al. 2016). Of the 120 isolates tested in that study, only one isolate (NZFS3325) was found to differ from the simile and verum morphotypes. An additional isolate (NZFS3305) was subsequently identified (McDougal et al. 2016). Based on the symptoms observed on Pinus radiata needles, these morphotypes do not appear to differ in the symptoms that they cause, however this has not yet been studied in detail. Host genotype is known to strongly influence symptom development of CNC (Suontama et al. 2019; Ismael et al. 2020). With climate change and seasonal variation shifting, it becomes more important to understand the biology of these morphotypes, in order to enable mapping of likely distributions and elucidating their virulence and pathogenicity under future climatic conditions.

Through morphological and molecular analysis, we have confirmed the presence of a third C. minus morphotype in New Zealand, here classified as ‘novus’. Growth temperature data suggested that the optimum is different for ‘verum’ compared to ‘simile’ and ‘novus’. This in turn indicated that sporulation and germination may occur at different temperatures, which could influence the geographical distribution of these morphotypes, highlighting the importance of comparative studies between them. It has been proposed previously that ‘verum’ and ‘simile’ are likely different species, and that these should be formally described (Prihatini et al. 2014). The discovery of this new morphotype novus, that is different again to the other morphotypes, provides further impetus for the need of formal descriptions.

To help shed light on how genetic differences in morphotypes might affect C. minus virulence or pathogenicity, the genomes of isolates from three morphotypes of C. minus were sequenced. These are the first genomic resources generated for Cyclaneusma and have enabled the prediction and comparison of CEs. The effector prediction pipeline generated a first list of possible effector proteins for C. minus. Because effectors are often described as SSCPs (Lo Presti et al. 2015), these features were used to identify C. minus CEs.

Some pathogens are known to secrete effectors that occur in families (Lo Presti et al. 2015; Denton-Giles et al. 2020; Rocafort et al. 2022). Following the observed similarity between DsEcp32-1 (part of the Ecp32 family from D. septosporum) and Cms835 (from isolate NZFS809), it was assessed to what extent the Ecp32 protein family was also present in this and the other seven isolates of C. minus. The results indicated that different morphotypes have different numbers of family members, with the morphotypes also showing differences in pseudogenisation of Ecp32 genes. This was also observed in the Ecp32 family from D. septosporum and D. pini (both pine pathogens) and F. fulva, as they have four, three, and five family members respectively, with the presence of pseudogenes in some isolate genomes (Tarallo et al. 2022).

The diversification of gene families can occur due to tandem and/or segmental duplications (Leister 2004), which both seem to be present in the expansion of this particular effector family from C. minus. For example, the genes encoding proteins Cms4058, Cms4059, and Cms4060 (truncated) from NZFS809 are on the same scaffold and adjacent to one another. This is different from what was observed in D. septosporum and F. fulva, where the locations of the Ecp32 genes are conserved in matching chromosomes between the two species but none of these family members are localized on the same chromosome (Mesarich et al. 2023). The differences in the number of family members and pseudogenes in each morphotype might reflect pathogen virulence and host range, as the selection pressure imposed by host R genes on plant pathogenic fungi can lead to an expansion and sequence divergence of an effector gene family (Pendleton et al. 2014; Franceschetti et al. 2017; Denton-Giles et al. 2020).

Also consistent with the predicted structures of Ecp32 proteins from D. septosporum and F. fulva, the predicted structures generated for representative members of the C. minus Ecp32 family were shown to adopt a β-trefoil fold. This is a fold present throughout different kingdoms that is involved in many different processes and interactions between organisms (Žurga et al. 2015). The most common types of proteins that share this fold are trypsin inhibitors and lectins, many being virulence factors in plant-pathogenic fungi (Schubert et al. 2012; Varrot et al. 2013; Juillot et al. 2016). It remains to be determined whether these β-trefoil proteins from C. minus have any virulence functions during infection. It would be interesting to determine the structures of other CEs from C. minus and identify to what extent this fold is present in the structures of any of those proteins. For example, the Ecp20 family, another cell death-eliciting family present in both D. septosporum and F. fulva (Tarallo et al. 2022), could not be detected in C. minus through primary sequence similarity searches. It remains possible, however, that searches in the predicted structures of all C. minus CEs might identify proteins that adopt the characteristic of Ecp20 family proteins, thereby expanding the repertoire of common effector types between these foliar fungal pathogens.

Proteins from the Ecp32 family in C. minus were found to have members that differed in their ability to trigger cell death in the non-host species N. benthamiana and N. tabacum. This difference was also observed in members of the Ecp32 family from D. septosporum and F. fulva (Tarallo et al. 2022). Whether the identified sequence variation contributes to the virulence functions of these proteins, as well as cell death elicitation, remains to be determined. Of note, other experiments, such as reactive oxygen species measurements and the analysis of defence marker gene expression in Nicotiana species, could be performed in the future to determine if these cell death responses are actually due to activation of the plant immune system. Similarly, the CE proteins should be transiently expressed in BAK1/SOBIR1-silenced or knock-out plants of Nicotiana species, to determine whether BAK1 and SOBIR1 (extracellular co-receptors involved in transducing defence response signals following apoplastic effector recognition) (Wang et al. 2018) are required for the elicitation of cell death by Ecp32 family members. The differences in cell death-triggering activity observed between different family members might be due to differences in the composition and positioning of amino acids on the surface on the proteins, for example as a result of insertions, deletions and/or substitutions of residues, which occur more often in loops than in α-helices and β-sheets (Fiser et al. 2000). In some cases, loops in tertiary structures of proteins are known to be involved in protein-protein interactions, recognition sites, signalling cascades, and ligand binding (Fetrow 1995), all of which could explain why some proteins from the C. minus Ecp32 family trigger more consistent cell death than others. Also, sequence and structural diversification observed for the proteins in this family could have a role in evasion of host immunity, another possible explanation for why this family has been expanded in C. minus.

Cms4060 (truncated) and Cmv5316 (full-length) triggered the same degree of cell death in N. benthamiana, but no cell death was observed in N. tabacum. This suggests that the N-terminal part of the protein (that is also conserved between Cms4060 and the other Ecp32 proteins) might contain the region, or conserved epitope, that is either recognized by a plant extracellular immune receptor in N. benthamiana or it interacts with plasma membrane proteins to trigger cell death. Previous studies have demonstrated the possibility of swapping loop regions from non-functional proteins for other loops from functional proteins (Wolfson et al. 1991; Pardon et al. 1995). Considering the consistent cell death triggered by Cmv9109, Cms9849, and Cmn3003, it would be informative to try to determine differences in the loop regions of these proteins that might be responsible for cell death, in comparison with other proteins that trigger inconsistent or no cell death. It would then be possible to determine if there are differences between the C. minus morphotypes in how these proteins function.

Conclusions

The generation of new genome sequences of eight C. minus isolates supported taxonomic classification of C. minus ‘novus’ as a third morphotype in the New Zealand population of this pathogen. It is important to understand the geographical range and variations in virulence and/or pathogenicity of these morphotypes, as well as the relationships between host genotype and fungal morphotypes, all of which have implications for diagnostics and biosecurity for the forest industry. The formal description of these morphotypes as separated species will be necessary to support forest biosecurity and disease management. It was also possible to identify CE proteins that might be the candidates for future investigation in the C. minus-pine interaction. We also showed that, just as in D. septosporum (another pine pathogen) and F. fulva, the Ecp32 family is present in C. minus, and the evidence suggests that some members from this family might be recognized by plant immune receptors. This, combined with the fact that these proteins are conserved at the primary sequence and tertiary structure levels, suggest that they have important and/or conserved roles in fungal virulence. Amino acid sequence searches revealed that the Ecp32 protein family is widespread across other fungal species (Tarallo et al. 2022) and, therefore, might be an important core effector of plant-pathogenic fungi. Furthermore, effector identification and analysis from the Cyclaneusma genomes has provided novel means to compare these morphotypes and demonstrated genomic-based differences that corroborate the hypothesis that these are in fact separate species. Disease resistance based on core effectors that are vital for a pathogen’s ability to cause disease is more likely to be durable. Knowledge of core effectors from different C. minus morphotypes can enable the development of molecular tools that can both assist in better diagnostics and understanding of the pathogen population as well as complement breeding programmes.

Methods

Microorganisms and plants

Cyclaneusma minus isolates NZFS110, NZFS725, NZFS759, NZFS809, NZFS1800, NZFS3305, NZFS3276, and NZFS3325 (Table 1) were used in this study. Each was grown on a sterile cellophane sheet overlaying 2% (w/v) Difco™ Malt Extract Agar (Becton Dickinson & Company, New Jersey, USA) and incubated at 22°C in the dark. Once the mycelium filled half of the plate, it was scraped into a sterile 15 mL Falcon tube and stored at − 80°C in preparation for freeze-drying and genomic DNA (gDNA) extraction.

Escherichia coli DH5α was used for CE gene cloning, and A. tumefaciens GV3101 was used for ATTAs. Nicotiana tabacum Wisconsin 38 and N. benthamiana were used as model non-host plants for ATTAs. Plants were grown in a temperature-controlled room at 22ºC and 12/12 light/dark photoperiod with 80–85 µmol/m²/s.

Colony morphology analysis

For colony morphology assessments, three 5 mm diameter plugs from the edge of an actively growing culture of each C. minus isolate were placed in the middle of 2% and 9% MEA plates. Plates were incubated as above and assessed for colony morphology at eight weeks post-inoculation using a standard colour chart used to describe colours in colonies from Rayner (1970). For the temperature study, a 5 mm diameter plug from each isolate was placed onto an independent 2% MEA plate as before, each with three replicates, and incubated at 4, 10, 17, 20, 22, 25, 28, or 30 °C in the dark. Two measurements at right angles for each replicate were taken at 14 days post-inoculation. The mean and standard deviation of the radial growth rate per day for each isolate were then calculated for each incubation temperature. Differences between radial growth per day at different temperatures were analysed using a two way ANOVA, containing temperature (as a character) and morphotype. Significant terms were followed up by applying a multiple-comparison procedure based on Tukey-adjusted contrast. The cardinal temperatures for each morphotype were determined from these data.

Genome sequencing, assembly, and annotation

Total gDNA for genome sequencing was extracted from freeze-dried mycelia of C. minus isolates using a DNeasy Plant Maxi Kit (Qiagen, USA), according to the manufacturer’s instructions. For each isolate, 200 mg of tissue was added to Lysing Matrix C tubes containing 1.0 mm silica spheres (MP Biomedical, California, USA) and homogenised to a fine powder using a FastPrep-24™ Instrument (MP Biomedical) with two cycles of 20 s at 5 m/sec. The extracted gDNA was concentrated by ethanol precipitation and resuspended in TE buffer (10 mM Tris and 1 mM EDTA, pH 8.0). The gDNA concentrations were quantified by spectrophotometry using both the Denovix DS-11 FX (DeNovix, Delaware, USA) and the DeNovix dsDNA Broad Range Assay (DeNovix, Delaware, USA).

The gDNA from eight isolates of C. minus (Table 1) was used for whole-genome sequencing using Illumina short-read sequencing technology (Illumina, California, USA). Paired-end libraries with insert sizes of 327–353 bp were sequenced on a NovaSeq 6000 system (Illumina, California, USA).

De novo genome assemblies of filtered reads were generated using SPAdes v3.14.0 (Nurk et al. 2013) with default settings. Contigs were assembled into scaffolds with redundans v0.13c (Pryszcz and Gabaldón 2016). Genome assembly quality was measured in QUAST (Gurevich et al. 2013). Assembly completeness was estimated with Benchmarking Universal Single-Copy Orthologs (BUSCO) comparisons to the fungal dataset (Simão et al. 2015).

Gene annotations were predicted in each masked genome assembly using the Funannotate v1.8.9 pipeline (Palmer and Stajich 2020). Gene prediction was carried out using predicted protein sequences of Elytroderma deformans CBS183.68 (available at: https://mycocosm.jgi.doe.gov/Elyde1/Elyde1.home.html) as a guide. At first, protein sequences of E. deformans were mapped to the C. minus assembly file and the orthologs from Cyclaneusma were extracted. These orthologous sequences were used to train Augustus v3.3.3 (Stanke et al. 2004), followed by three rounds of iterative gene prediction. Genes were predicted by Augustus v3.3.3, GeneMark ES v4.68 (Ter-Hovhannisyan et al. 2008), snap 2006-07-28 (Korf 2004), and GlimmerHMM v3.0.4 (Majoros et al. 2004), while transfer RNAs were predicted by tRNAscan-SE v2.0.9 (Chan et al. 2021). EVidenceModeler v1.1.1 (Haas et al. 2008) was used to generate consensus gene models. Putative functions of gene products were determined by searches to the MEROPS v12.0 (Rawlings et al. 2018), UniProt DB 2021_03 (The UniProt 2021), dbCAN v10.0 (Yin et al. 2012), PFAM v34.0 (Mistry et al. 2021), and BUSCO v2.0 (Simão et al. 2015) databases.

Phylogenetic trees

The phylogenetic positions of the C. minus strains were assessed by constructing phylogenetic trees using concatenated nucleotide sequences from the ITS of ribosomal DNA, nuclear LSU (nLSU), and translation elongation factor 1α (tef-1). An initial alignment of these barcode genes from Lophodermium conigenum NZFS790 to the assembled genomes was done using LASTZ (Harris 2007) and Geneious R10 (Kearse et al. 2012) to find the C. minus orthologs. Orthologs were then extracted and aligned to DNA sequences from other Cyclaneusma species, as well as Lophodermium and Strasseria species (Additional file 2: Table S6). The alignment for the tree was generated using a cost-matrix of 51% with pairwise alignment or multiple alignment gap opening penalty of 12, plus a gap extension penalty of 0.3. The DNA sequences were aligned in Geneious R10. Nucleotide sequences not aligned at the 5’ and 3’ extremities of alignments were trimmed, and internal gaps maintained. The multi-gene phylogenetic trees were created by a maximum likelihood algorithm with 1000 replicates and distances obtained using the General Time Reversible model using Mega X (Kumar et al. 2018).

Effector prediction

Following genome assembly and annotation of the eight C. minus isolates, a previously described pipeline (Hunziker et al. 2021) was used to identify CE proteins from all genomes, based on their full sets of predicted proteins. Briefly, signal peptides were predicted using SignalP v3.0 (Bendtsen et al. 2004). To exclude membrane-bound proteins, TMHMM v2.0 (Krogh et al. 2001) and BIG-PI (Eisenhaber et al. 2004) were used. Predicted proteins of ≤ 300 amino acids in length with ≥ 4 cysteine residues were then selected and classified as small, secreted cysteine-rich proteins (SSCPs). EffectorP v3.0 (Sperschneider and Dodds 2021) was then used to select proteins classified as apoplastic or cytoplasmic effectors, with candidates that were classified as non-effectors being discarded.

The amino acid sequence of DsEcp32-1, a member of the Ecp32 effector family from D. septosporum (Tarallo et al. 2022), was firstly used to identify the likely orthologue in C. minus isolate NZFS809. This C. minus sequence was then used to identify possible Ecp32 family homologues in the other seven isolates of C. minus, through reciprocal BLASTp and tBLASTn searches in conjunction with an E-value cut-off of < 10⁻⁵.

Protein sequence alignments were performed using Clustal Ω (Sievers et al. 2011) in Jalview (Waterhouse et al. 2009). Phylogenetic trees were constructed with the neighbour-joining method based on 1000 bootstrap replicates using Geneious Software v9.1.8 (Kearse et al. 2012) in conjunction with CE protein sequence alignments. IDRs were predicted by querying protein sequences in the Predictor of Natural Disordered Regions (PONDR® VLXT) server (Romero et al. 2001). A prediction score ≥ 0.8 was considered a significant IDR prediction. Protein tertiary structure predictions of Ecp32 family members were performed using AlphaFold2 in conjunction with ColabFold (Jumper et al. 2021; Mirdita et al. 2022). Two values were used to assess the confidence of the predictions: the predicted local-distance difference test (pLDDT) score, which ranges from 0 to 100, with values closer to 100 indicating a high confidence in prediction, and the global superposition metric template modelling score (TM-score), ranging from 0 to 1, with a value of ≥ 0.5 indicating a similar fold between structures. The Dali server (Holm 2020) was used to identify proteins with structural similarity in the Research Collaboratory for Structural Bioinformatics Protein Data Bank. Here, a Dali Z-score of ≥ 2 was used to infer structural similarity. Protein tertiary structures were visualized and rendered in PyMOL v2.5, with alignments carried out using the CEalign tool (DeLano 2002).

Agrobacterium tumefaciens-mediated transient expression assays (ATTAs)

Cyclaneusma minus ATTA expression vectors were generated using the method described by Guo et al. (2020). The expression vector used was pICH86988, which contains a N. tabacum N-terminal PR1α signal sequence for apoplast secretion followed by an N-terminal 3xFLAG tag for western blotting detection (Weber et al. 2011). CE gene sequences of C. minus, excluding the nucleotide sequence encoding the predicted signal peptide, were synthesized by Twist Bioscience (California, USA), and included BsaI recognition sites and 4 bp overhangs specific for Golden Gate modular cloning (Engler et al. 2009). CE genes were also cloned into a version of the vector that lacked the sequence encoding the PR1α signal peptide, but that contained only a start codon. Constructs were transformed into competent cells of E. coli by electroporation and plasmids were extracted using an E.Z.N.A.® Plasmid DNA Mini Kit I (Omega Bio-tek, Georgia, USA). Correct assemblies were confirmed by sequencing in association with the Massey Genome Service (Palmerston North, New Zealand). Expression vectors were then transformed into competent cells of A. tumefaciens by electroporation (Guo et al. 2020).

To perform the ATTAs, transformed cells of A. tumefaciens containing the CE genes were inoculated into lysogeny broth supplemented with 50 µg/mL kanamycin, 10 µg/mL rifampicin, and 30 µg/mL gentamycin and grown overnight at 28ºC in the dark until an OD₆₀₀ of 0.4–2.0 was reached. Cells were collected by centrifugation at 2500×g for 5 min and resuspended in 1 mL of infiltration buffer (10 mM MgCl₂:6H₂O, 10 mM MES-KOH, 100 µM acetosyringone (Sigma-Aldrich, Missouri, USA). The OD₆₀₀ was measured again, and sterile water added to reach a final OD₆₀₀ of 0.5. These cultures were then incubated for 3 h at room temperature before infiltration. An A. tumefaciens transformant carrying an expression vector for the extracellular elicitin INF1, from Phytophthora infestans (Kamoun et al. 1997), was used as a positive control and the empty pICH86988 vector was used as a negative control. The abaxial surface of N. benthamiana and N. tabacum leaves was infiltrated with the appropriate A. tumefaciens solution using 1 mL needleless syringes. The screening was done with at least three biological replicates. Each biological replicate consisted of two plants, each with two leaves infiltrated, resulting in at least 12 infiltration zones. Plants were monitored for up to seven days, at which time photographs were taken with a Nikon D7000 camera. To assess the plant response caused by each expressed protein, infiltration zones were divided into a strong or weak cell death response, chlorosis, or no cell death. A strong cell death response was considered when it was indistinguishable from the response elicited by INF1 and covered the entire infiltration zone, while a weak response was considered when the protein elicited a response across approximately 50% of the infiltration zone when compared with INF1. Chlorosis was considered when leaf discoloration occurred, but no apparent cell death was observed. A no-cell death response was indicated when the infiltration zone had the same lack of a response as the negative empty pICH86988 vector control.