Widespread Occurrence of a CYP51A Pseudogene in Calonectria pseudonaviculata.

Calonectria pseudonaviculata and C. henricotiae are two closely related fungal species responsible for boxwood blight disease of ornamental shrubs (Buxus spp.) in the U.S. and Europe. A previous study has shown isolates of the latter species, which is restricted to Europe, to be less sensitive to tetraconazole, an azole fungicide. In this study, we have analyzed the CYP51 paralogs for polymorphism in 26 genomes, representing geographically disparate populations of C. pseudonaviculata (n = 19) and C. henricotiae (n = 7), from the U.S., Europe, Asia, and New Zealand. The presence of a CYP51A pseudogene and lack of a functional CYP51A paralog in all C. pseudonaviculata genomes examined is a novel discovery for fungi and could have implications for the evolution of resistance to antifungal chemicals.

Introduction

Calonectria pseudonaviculata and C. henricotiae are the causal agents of the invasive fungal disease boxwood blight in native and commercial Buxus species [1,2]. This disease manifests itself through foliar and stem lesions, potentially leading to plant death. C. pseudonaviculata causes the disease in Europe, the U.S., New Zealand, and Asia, whereas C. henricotiae has been documented only in Europe [3].

C. pseudonaviculata was originally identified as conspecific with C. henricotiae, though [3] showed that these are two distinct species based on a multi-gene phylogeny. Phenotypic differences were also reported, including higher thermotolerance in C. henricotiae [3].

Among the antifungal chemistries used to control boxwood blight, and a multitude of other fungal diseases in agriculture and medicine, is the azole chemical class, also known as demethylation inhibitors (DMIs) [4,5]. The azoles bind sterol 14α-demethylase, also known as CYP51, a fungal-specific enzyme which facilitates the synthesis of the fungal membrane sterol ergosterol, a critical component of fungal cell membranes [6,7]. By binding to the active site of this enzyme, the azoles are inhibitors of ergosterol biosynthesis, thereby weakening and disrupting the cell membrane.

Repeated exposure of fungal pathogens to azoles provides a significant selection pressure on populations of fungi [8-10]. Over time, individuals may develop reduced sensitivity and ultimately become differentially resistant to single azoles or to the entire class [10]. Fungi overcome inhibition of ergosterol biosynthesis by azoles via a variety of mechanisms involving CYP51. Many reports have shown that fungi incur nonsynonymous mutations in the CYP51 gene, which alter the conformation of the encoded active-site amino-acid residues. Such changes render the azoles unable to bind to the active site, thereby enabling continued ergosterol biosynthesis. Another commonly reported mechanism involves overexpression of the CYP51 gene in order to compensate for inhibition [10].

Some species of fungi have been shown to harbor CYP51 paralogs, or duplicated copies of the gene, which can accumulate mutations leading to gene loss or a change in protein structure [10,11]. These paralogs are denoted CYP51A, CYP51B, and CYP51C. Differential gene loss of the paralogs has resulted in only a single copy, either CYP51A or CYP51B, being left behind in diverse fungal lineages, with members of these lineages developing resistance and reduced sensitivity to azoles no matter which of the paralogs remain [10]. Recent studies suggest specific roles of the paralogs in the normal growth of Fusarium species, which are known to encode all three genes [11]. Individual gene knockouts of each of the Fusarium CYP51 paralogs resulted in increased sensitivity to different azole fungicides for CYP51A and CYP51B, but not CYP51C. Of the three paralogs, CYP51B was found to be most vital for normal growth.

Select genomes of C. pseudonaviculata and C. henricotiae species possess three CYP51 paralogs; however, in C. pseudonaviculata, CYP51A is a pseudogene due to a premature stop codon [12]. A widespread CYP51A pseudogene has been reported in another fungus, Rhynchosporium commune, where the presence or absence of a functional CYP51A paralog has a population-specific distribution [7,13]. To date, there are no reports of fungal species for which all isolates have only a CYP51A pseudogene and no functional CYP51A paralog. Based on in vitro assays on solid agar media, C. henricotiae displays significant reduced sensitivity to the azole tetraconazole when compared to C. pseudonaviculata, which could be explained by the presence of a functional CYP51A in the former and only a CYP51A pseudogene in the latter [3].

As a result of these insights, we sought to understand the distribution of the CYP51A pseudogene in C. pseudonaviculata isolates. Our specific objectives were to mine the CYP51 paralogs from genome sequences originating from a worldwide collection of C. pseudonaviculata isolates and compare these sequences with those of C. henricotiae in order to identify CYP51 pseudogenes and/or polymorphisms that could explain the differential sensitivity to tetraconazole reported by [3]. Our findings have implications for the evolution of the CYP51 pseudogenes and paralogs in fungal pathogens, which could inform future studies on the mechanisms underpinning resistance to azole antifungals.

Methods

Nineteen genomes of C. pseudonaviculata and seven genomes of C. henricotiae representing geographically disparate populations (Table 1) were generated at the Mycology and Nematology Genetic Diversity Laboratory at the U.S. Department of Agriculture. Genomic data were generated for the fungal isolates as described in [14] using the Illumina MiSeq platform. Sequence reads were further processed to remove adapter traces using cutadapt v1.14 [15]. Draft genome assemblies were generated from the processed paired-end reads using SPAdes v3.10.0 [16]. Gene models were predicted in the previously generated Illumina-PacBio hybrid assembly of isolate CT1 [3,17] using the AUGUSTUS web server with the Fusarium graminearum species training set [18].

Table 1.

Identification and origin of Calonectria pseudonaviculata and C. henricotiae isolates from which genomes were examined for this study.

Isolate Id	Alternate ID	Speciesa	Continent	Location	Year	Pseudogene
CT1	CpsCT1	C.ps.	United States	Connecticut	2011	Y
JAC13-14	CpsCT10	C.ps.	United States	Connecticut	2012	Y
JAC13-118	CB088	C.ps.	Europe	Belgium	2011	Y
JAC13-149	NC BB13	C.ps.	United States	North Carolina	2012	Y
JAC13-167	TU004	C.ps.	Asia	Turkey	2012	Y
JAC13-172	91.9.6 A	C.ps.	Asia	Iran	2013	Y
JAC13-182	STE-U 3399	C.ps.	New Zealand	New Zealand	1998	Y
JAC18-13	CDFAS582	C.ps.	United States	California	2017	Y
JAC18-14	CDFAS666	C.ps.	United States	California	2016	Y
JAC13-59	MDSH 7.7	C.ps.	United States	Maryland	2013	Y
JAC13-107	CB054	C.ps.	Europe	Belgium	2010	Y
JAC13-110	CB105	C.ps.	Europe	Belgium	2012	Y
JAC13-126	MM2013calTA7	C.ps.	Asia	Iran	2013	Y
JAC13-127	11-416-11	C.ps.	Europe	Slovenia	2011	Y
JAC14-13	CB002	C.ps.	Europe	Belgium	2008	Y
JAC14-15	RHS 21350	C.ps.	Europe	United Kingdom	2007	Y
JAC14-108	NY13-370.9a	C.ps.	United States	New York	2013	Y
JAC14-154	09-1762	C.ps.	Europe	France	2009	Y
JAC15-1	2015-7-36-0041	C.ps.	United States	Pennsylvania	2015	Y
JAC13-124	PD 011/04744201	C.h.	Europe	The Netherlands	2011	N
JAC13-131	P-10-5865	C.h.	Europe	Germany	2010	N
JAC13-147	DE017	C.h.	Europe	Germany	2011	N
JAC13-185	JKI 2100	C.h.	Europe	Germany	2007	N
JAC13-216	RHS 192076	C.h.	Europe	United Kingdom	2012	N
JAC14-101	NL018	C.h.	Europe	The Netherlands	2011	N
JAC14-47	NL019	C.h.	Europe	The Netherlands	2011	N

The presence of the CYP51A pseudogene is indicated for each isolate.

aC.ps. designates C. pseudonaviculata, while C.h. designates C. henricotiae.

The CYP51 paralog sequences were obtained by running local blastn searches of the above genome sequences using the CLC Genomics Workbench v11 (QIAGEN, Redwood City, CA, U.S.A.) [19]. Preliminary blastn searches were conducted using Fusarium graminearum CYP51 paralogs as query sequences [11]. Intron sequences were verified and trimmed from the assembled genome data guided by the gene prediction annotation of the CT1 genome. The sequences were aligned in MEGA v7 and translated into amino-acid sequences [20]. The CYP51 coding sequences were searched for early stop codons or nonsynonymous polymorphisms. The CYP51 paralog sequences for C. henricotiae were uploaded into GenBank with the following accessions: CYP51A, MN497050; CYP51B, MN497051; CYP51C, MN497052. The CYP51 paralog sequences for C. pseudonaviculata were uploaded into GenBank with the following accessions: CYP51A, MN497053; CYP51B, MN497054; CYP51C, MN497055.

A blastp search against the non-redundant protein sequence database at NCBI was performed in January and February of 2019 using the translated coding sequences for the CYP51 paralogs from isolates of both C. pseudonaviculata and C. henricotiae [19]. In both species, each paralog was found to have only one unique genotype and thus one unique amino acid sequence per species. Therefore, the blastp search was conducted using the coding sequences of CYP51 paralogs originating from a single genome of each species.

A phylogenetic analysis was carried out with amino acid sequences of the three paralogs, which were aligned and analyzed according to [9]. Briefly, amino acid sequences from this study were aligned to a subset of those from [6]. Sequences were aligned using the MAFFT (Multiple Alignment using Fast Fourier Transform) sequence alignment program [21]. A neighbor joining phylogenetic tree was generated with 1000 bootstrap replicates using MEGA v7 [20].

Results

The CYP51A gene had one predicted intron, which began at nucleotide 193 and ended at nucleotide 251. The CYP51B gene had two predicted introns, one beginning at position 247 and ending at position 355; the other beginning at position 556 and ending at position 630. The CYP51C gene had two predicted introns, the first beginning at position 424 and ending at position 478; the second beginning at position 1544 and ending at position 1591.

For the CYP51A amino acid sequence of C. pseudonaviculata, the blastp top hit was eburical 14-alpha-demethylase from Neonectria ditissima (KPM38499.1), with a query coverage of 99% and an E value of 0. The CYP51A sequence had 80% amino acid identity (404 of 506 amino acids) and 89% positive amino acid matches (453 of 506 amino acids) with respect to the N. ditissima protein. For CYP51B of C. pseudonaviculata, the top blastp hit was to a predicted protein from Nectria haematococca (XP_003054236.1), which had a query coverage of 100% and an E value of 0. The CYP51B sequence had 86% amino acid identity (455 of 527 amino acids) and 92% positive amino acid matches (488 of 527 amino acids) with respect to the N. haematococca protein. The top blastp hit for the CYP51C of C. pseudonaviculata was a hypothetical protein of Fusarium euwallaceae (RTE69955.1), which had a query coverage of 100% and an E value of 0. The CYP51C sequence had 83% amino acid identity (427 of 516 amino acids) and 91% positive amino acid matches (473 of 516 amino acids) with respect to the F. euwallaceae protein.

For the CYP51A amino acid sequence of C. henricotiae, the top blastp hit was to the same eburical 14-alpha-demethylase found in N. ditissima (KPM38499.1), which had a query coverage of 99% and an E value of 0. The CYP51A sequence also had 80% amino acid identity (404 of 506 amino acids) and 89% positive amino acid matches (453 of 506 amino acids) with respect to the N. ditissima protein. For CYP51B of C. henricotiae, the top blastp hit was to the same predicted protein in N. haematococca (XP_003054236.1), which had a query coverage of 99% and an E value of 0. The CYP51B sequence had 86% amino acid identity (455 of 527 amino acids) and 92% positive amino acid matches (488 of 527 amino acids) with respect to the N. haematococca protein. The top blastp hit for the CYP51C of C. henricotiae was also a hypothetical protein of F. euwallaceae (RTE69955.1), which had a query coverage of 99% and an E value of 0. The CYP51C sequence had 83% amino acid identity (428 of 516 amino acids) and 91% positive amino acid matches (473 of 516 amino acids) with respect to the F. euwallaceae protein.

Amino acid sequences encoded by the CYP51 gene paralogs were separated into three monophyletic clades along a phylogenetic tree with bootstrap values greater than 70% (Figure 1). The three clades consisted solely of CYP51A, CYP51B, or CYP51C sequences for various fungal lineages. The C. pseudonaviculata and C. henricotiae paralog sequences were paired together across each of the three clades with high bootstrap support, as expected for sister taxa [3].

Figure 1.

Neighbor-joining phylogenetic tree (1000 bootstrap replicates) of amino acid sequences of CYP51 proteins from diverse fungal lineages. Amino acid sequences from Calonectria pseudonaviculata and C. henricotiae were aligned with those of other fungal species from [6] excluding Fusarium graminearum accession AFN66169.1. Clades representing proteins encoded by CYP51 gene paralogs are indicated to the right of the tree. The scale bar indicates the substitution rate of amino acid residues per site.

All isolates of C. pseudonaviculata displayed a single CYP51A pseudogene, with premature stop codon in the CYP51A nucleotide sequence, thereby confirming the earlier report of this pseudogene (12; Table 2). The early TGA stop was the result of a nucleotide substitution (C688T). No functional CYP51A paralog was found from blast searches of the C. pseudonaviculata genomes. None of the C. henricotiae isolates displayed the premature stop codon in CYP51A.

Table 2.

Summary of single nucleotide polymorphisms (SNPs) from the CYP51 paralog coding sequences between Calonectria pseudonaviculata and C. henricotiae.

CYP51A		CYP51B		CYP51C
SNPa	Synonymous or nonsynonymousb	SNP	Synonymous or nonsynonymous	SNP	Synonymous or nonsynonymous
G109A	N, V37I	C384T	S, 128	G45A	S, 15
G207C	S, 69	A924G	S, 308	C54A	S, 18
G234T	N, K78N	C1218T	S, 406	T183C	S, 51
T496C	N, W166R			C204T	S, 68
T534G	N, H178Q			T619C	S, 207
A579G	S, 193			T753C	S, 251
G587A	N, R196K			G780A	S, 260
C652T	S, 218			A1020G	S, 340
T688C	N, *230R			A1430G	N, K477R
C704A	N, A235E
C758T	N, A253V
T847C	N, Y283H
G868A	N, V290M
A1049G	N, K350R
T1271C	N, V424A
T1366C	N, Y456H
G1418A	N, S473N
T1521C	S, 507
A1522C	N, K508Q

Amino acid substitutions resulting from nonsynonymous nucleotides are indicated.

aAny nucleotide or amino acid substitution is ordered as C. pseudonaviculata first and C. henricotiae second.

bS denotes synonymous polymorphism, while N denotes nonsynonymous polymorphism.

*Denotes the premature stop encoded by CYP51A from C. pseudonaviculata.

All 19 C. pseudonaviculata genomes shared a single genotype for CYP51B and CYP51C. Similarly, all seven C. henricotiae genomes shared a single genotype for these paralogs. With respect to isolates of C. henricotiae, C. pseudonaviculata exhibited polymorphism in both paralogs. In contrast, C. pseudonaviculata and C. henricotiae CYP51B gene paralogs encoded identical amino acid sequences (Table 2). In CYP51C, one nonsynonymous polymorphism was present.

Discussion

A CYP51A pseudogene was mined from the genomes of all 19 isolates of C. pseudonaviculata from four continents, confirming its widespread distribution. The premature stop codon occurred in the CYP51A pseudogene prior to the region encoding the heme-binding site residues of the CYP51A protein [22,23]. We hypothesize that this interrupts function of the CYP51A protein, thus restricting the potential total biosynthesis of ergosterol when compared to C. henricotiae, since we did not uncover evidence for a functional CYP51A paralog in C. pseudonaviculata. We speculate that this could partially explain the greater in vitro sensitivity of C. pseudonaviculata to tetraconazole when compared to C. henricotiae and perhaps be linked to differences in thermotolerance between these closely related species [3,7,10].

This is the first report to our knowledge of a ubiquitous CYP51 pseudogene without evidence for a co-occurring functional CYP51A paralog. Brunner et al. [13] examined isolates of Rhynchosporium commune from barley fields across the globe and found the presence of a CYP51A pseudogene to be absolute; however, the presence of a functional CYP51A paralog was found in many cases to co-occur with the pseudogene, especially in countries where prior introduction of azole fungicides has occurred. Further, these populations of R. commune were shown to consist of both members with only the CYP51A pseudogene and those with a functional CYP51A paralog in addition to the pseudogene [13]. Historical sampling by [7] reported that during the twentieth century R. commune populations at an experimental site in the UK displayed the CYP51A pseudogene only at high frequencies and that a functional CYP51A paralog began to increase in frequency during the mid-1980’s, a phenomenon termed “paralog re-emergence”. This occurred during a period which coincided with introduction of new azole fungicides to the United Kingdom mirroring the findings of [13] which noted that the functional CYP51A paralog appeared to be undergoing significant positive selection, likely due to intragenic recombination, in countries with prior azole exposure (2016). In both studies, the presence of the functional CYP51A paralog was definitively associated with a dramatic reduction in sensitivity to azoles.

In our analyses of genomes from nineteen isolates of C. pseudonaviculata from across four continents, we did not uncover a functional CYP51A gene paralog, though all isolates of C. henricotiae possessed a CYP51A paralog and lacked the CYP51A pseudogene. Perhaps with deeper sampling of populations from worldwide collections of C. pseudonaviculata, CYP51A paralog reemergence will be documented. LeBlanc et al. [24] characterized a global collection of C. pseudonaviculata, a subset of which were applied to this study, with simple sequence repeat markers and concluded that the widespread clonality in this pathogen could be explained by multiple introductions of a single clonal lineage. Given the lack of sexual recombination reported for this fungus, the likelihood of paralog reemergence of CYP51A in C. pseudonaviculata is unclear, though results presented here should facilitate future efforts to document this phenomenon [24,25].

Both Brunner et al. [13] and Hawkins et al. [7] concluded that the presence of a functional CYP51A paralog alone was associated with a dramatic decrease in sensitivity to azoles, which was thus not attributed to polymorphism in the CYP51B paralog. The absence of nonsynonymous substitutions in CYP51B between the two species is further support that the lack of a functional CYP51A in C. pseudonaviculata may at least partially explain the previously reported discrepancy in sensitivity to tetraconazole [3]. A single nonsynonymous polymorphism was identified in CYP51C between the two species, but evidence for a role of this gene in sensitivity to azoles has not been reported for fungi. Future work will examine expression levels of the CYP51 paralogs in response to tetraconazole exposure, in addition to promoter analysis, in order to screen for single nucleotide polymorphisms or repetitive DNA elements, which have been linked to overexpression of CYP51A in plant pathogenic fungi [26]. Additionally, RNA-Seq transcriptomic analyses could be applied to identify detoxification genes from the fungal xenome that share a coordinated role in the reduced sensitivity to an azole fungicide reported for C. henricotiae [9,27].

In summary, our results contribute to the knowledge base on the evolution of CYP51 gene paralogs in the ascomycete fungi, which may have implications for treating boxwood blight and other fungal diseases. Consequently, the possibility of paralog reemergence of CYP51A may exist for C. pseudonaviculata, which could reduce the efficacy of azoles in treating future boxwood blight outbreaks. In turn, the reported clonal nature of the pathogen, the ubiquity of the CYP51A pseudogene and apparent lack of a functional CYP51A paralog in C. pseudonaviculata may place constraints on CYP51B mutation under selection of azoles applied in the field. Finally, research on this pathosystem presents a unique future opportunity to examine additional genetic and genomic factors underpinning differential antifungal sensitivity in two recently diverged species.