Phylogenomic analyses of Alternaria section Alternaria: A high-resolution, genome-wide study of lineage sorting and gene tree discordance

ABSTRACT

The genus Alternaria contains a diversity of saprobic and pathogenic species that can be found in a wide range of environments. Alternaria is currently divided into 26 subgeneric sections, and the “small-spored” Alternaria section Alternaria includes many species that are economically important agricultural pathogens. Recognizing that a stable framework for systematics and species identification is essential for management and regulation purposes, this section has experienced much taxonomic debate and systematic revision in recent years. Molecular phylogenetic studies have challenged the reliability of using morphological characteristics to differentiate Alternaria species but have also suggested that commonly used molecular markers for fungal phylogenetics may not be sufficiently informative at this taxonomic level. To allow the assessment of molecular variation and evolutionary history at a genome-wide scale, we present an overview and analysis of phylogenomic resources for Alternaria section Alternaria. We review the currently available genomic resources and report five newly sequenced genomes. We then perform multiple comparative genomic analyses, including macrosynteny assessment and inference of phylogenetic relationships using a variety of data sets and analysis methods. Fine-scale, genome-wide phylogenetic reconstruction revealed incomplete lineage sorting and the genomic distribution of gene/species tree discordance. Based on these patterns, we propose a list of candidate genes that may be developed into informative markers that are diagnostic for the main lineages. This overview identifies gaps in knowledge and can guide future genome sequencing efforts for this important group of plant pathogenic fungi.

KEYWORDS:

• Alternaria alternata
• arborescens
• gaisen
• longipes
• whole-genome sequencing
• comparative genomics

INTRODUCTION

The genus Alternaria, characterized by chains of darkly pigmented, multicelled conidia with transverse and longitudinal septa (phaeodictyospores) and tapering apical beaks, was established by Nees in 1817 (Nees 1817). Over the next 150 years, other genera (e.g., Brachycladium, Macrosporium, Stemphylium, Ulocladium) were erected to accommodate variations in Alternaria-like phaeodictyospore morphologies (Lawrence et al. 2016). Decades of comprehensive revision of Alternaria and related genera based on morphological characteristics culminated in the recognition of a total of 275 Alternaria species. These species have been further divided into subgeneric groups based on broad patterns of morphological similarity (Simmons 1967, 2007).

Results of recent molecular phylogenetic studies have challenged the morphological basis of the accepted taxonomy of Alternaria and close relatives (Pryor and Gilbertson 2000; Peever et al. 2004; Hong et al. 2005; Andrew et al. 2009; Lawrence et al. 2013; Woudendberg et al. 2013, 2015; Armitage et al. 2015, 2020). In general, the morphologically distinct subgeneric species groups, now called sections, tend to be supported by molecular data. However, the phylogenetic relationships within species groups are often inconsistent with the subtle phenotypic characteristics used to delineate morphospecies. This pattern was also observed at the generic level: lack of phylogenetic differentiation led to the synonymization of 13 other genera into Alternaria (Woudenberg et al. 2013), thereby redefining and expanding the generic concept of Alternaria. Collectively, these studies suggest that morphological characteristics typically used to delineate species (e.g., conidium length, width, septation; chain structure; beak shape) may not reflect evolutionary relationships between taxa.

A subgeneric species group that has experienced much taxonomic debate and confusion is the “small-spored” Alternaria section Alternaria (Lawrence et al. 2013; Woudenberg et al. 2015). This section contains many species that cause disease of economically important crops, so consistent species delineation and identification are essential for management and quarantine purposes. Particular species are associated with specific host ranges and pathogenicity, due in part to the ability to produce host-specific toxins (Meena et al. 2017). Some species are included on regulated organism lists for international phytosanitary certification, further highlighting the need for stable taxonomy and nomenclature. Even in the hands of experienced mycologists, the morphological characteristics used to differentiate species in section Alternaria may overlap between species and be sensitive to variation in growth conditions. Different molecular studies have led to different systematics proposals for section Alternaria, potentially due to the use of different molecular markers for inferring phylogenetic relationships and the differentiation among taxa. When working with strains that are so closely related, the commonly used genes for fungal phylogenetics may not be informative enough at this taxonomic level. In such cases, adopting a phylogenomics approach and assessing variation at a genome-wide scale may provide the necessary discriminatory power.

Here, we explore a phylogenomics approach for Alternaria section Alternaria (also known as Alternaria alternata species group) to address questions of molecular systematics and evolutionary history. First, we provide a review of the genomic resources currently available for this section and release five newly sequenced genomes from isolates collected from cereal crops in Canada. We then perform a range of comparative genomic analyses, including macrosynteny assessment and inference of phylogenetic relationships using a variety of data sets and analysis methods. Fine-scale, genome-wide phylogenetic reconstruction allowed us to investigate incomplete lineage sorting and the genomic distribution of gene/species tree discordance. By viewing how patterns of concordance varied with genomic context, we determined the regions of the genome that supported the main lineages identified by whole-genome analyses. We then propose a list of candidate genes that can be developed into diagnostic markers for distinguishing the main lineages. This overview identifies gaps in knowledge and can guide future genome sequencing efforts for this important group of plant pathogenic fungi.

MATERIALS AND METHODS

Existing genomic resources.—

All genome assemblies for section Alternaria that were publically available on January 27, 2020, were downloaded. Assembly quality was assessed with QUAST (5.0.2; Mikheenko et al. 2018), and genomes with inconsistent genome sizes (>50 Mb) or low quality scores (N50 <10 kb) were omitted. For one additional isolate (EGS 90.391), raw short-read sequence data were obtained directly from the authors who published genome analyses previously (Woudenberg et al. 2015), resulting in 33 genomes (TABLE 1). Here, we present the original species names given to the strains at the time of genome submission but also indicate what their current names would be based on the revision of section Alternaria by Woudenberg et al. (2015).

Table 1. Summary of information and statistics for the analyzed section Alternaria genomes

Given species name at time of genome submission	Taxonomy notes^a	Isolate name	Sampling notes	Reference	Accession source	Data type	No. of contigs (>500 bp)	N50 (kb)	L50 (contigs)	Length (Mb)	GC%	Complete BUSCO genes (%)	Missing BUSCO genes (%)
A. alternata		ATCC 11680 (BMP 0238)	Allium, USA	Dang et al. 2015	mycocosm.jgi.doe.gov	Illumina	253	450.6	25	33.52	51.13	99.00	0.88
A. alternata		ATCC 34957	Sorghum, USA	Nguyen et al. 2016	LMXP01000000	PacBio	27	2828.0	5	33.50	50.95	98.94	0.94
A. alternata	A. alternata ex-type	ATCC 66891 (EGS 34.016, BMP 0269, CBS 916.96)	Peanut, India	Dang et al. 2015	mycocosm.jgi.doe.gov	Illumina	278	303.9	33	33.15	51.15	99.00	0.64
A. alternata		BMP 0270 (ATCC 66982, EGS 34.039)		Dang et al. 2015	mycocosm.jgi.doe.gov	Illumina	171	757.3	14	33.40	51.02	99.12	0.70
A. alternata		DET2001	Wheat, Canada	This study	JAAOQX010000000	Illumina	111	1206.5	11	33.47	51.03	99.36	0.59
A. alternata		DET2010	Wheat, Canada	This study	JAAOQV010000000	Illumina	216	1366.9	8	33.79	51.01	99.12	0.82
A. alternata		DET2019	Wheat, Canada	This study	JAAOQU010000000	Illumina	449	824.4	12	35.17	51.00	99.18	0.76
A. alternata	Apple pathotype	FERA 1177	Apple	Armitage et al. 2020	PDXD01000000	Illumina	735	536.6	17	35.63	51.11	99.18	0.82
A. alternata		JS-1623	Fir, Korea	Park et al. 2020	VZUP01000000	Illumina + PacBio	11	4961.7	3	33.67	50.97	92.85	6.21
A. alternata		KAS5516	Rye, Canada	This study	JAAOQT010000000	Illumina	174	1478.4	8	33.95	51.02	99.12	0.82
A. alternata		MOD1-FUNGI5	Grapes, USA	Gebru et al. 2020	SJDQ01000000	Illumina	1085	71.5	140	33.40	51.31	98.18	1.23
A. alternata	Peach pathotype	NAP07	Peach, Japan	Ushijima and Yamamoto 2019	BJEP01000000	Nanopore	19	3070.6	4	35.84	50.98	99.12	0.76
A. alternata		PPRI 21032	Onion, South Africa	Bihon et al. 2016	LSHC01000000	Illumina	530	244.4	41	33.01	51.30	98.01	1.76
A. alternata		SRC1lrK2f	Soil, USA	Zeiner et al. 2016	LXPP01000000	Illumina	91	1097.9	12	32.99	51.40	99.36	0.59
A. alternata	Tangerine pathotype	Z7	Tangerine	Wang et al. 2016	LPVP01000000	Illumina + PacBio	156	1182.7	9	34.36	50.99	99.36	0.59
A. arborescens		DET2008	Wheat, Canada	This study	JAAOQW010000000	Illumina	127	816.7	14	33.30	50.99	99.24	0.76
A. arborescens	Tomato pathotype, A. arborescens ex-type	EGS 39.128 (CBS 102605)	Tomato, USA	Hu et al. 2012	AIIC01000000	Illumina	493	314.3	30	33.79	50.93	99.30	0.64
A. arborescens		FERA 675	Asian Pear	Armitage et al. 2020	PDUP01000000	Illumina	325	862.6	14	33.94	51.06	99.24	0.70
A. arborescens		MOD1-FUNGI6	Apple, USA	Gebru et al. 2020	SJDP01000000	Illumina	459	302.5	30	33.75	51.10	99.12	0.70
A. arborescens		RGR 97.0013	Apple	Armitage et al. 2020	PDWY01000000	Illumina	287	684.7	17	33.80	50.89	99.36	0.64
A. arborescens		RGR 97.0016	Apple	Armitage et al. 2020	PEJP01000000	Illumina	339	528.3	23	33.77	51.08	99.30	0.70
A. citriarbusti	Moved to A. alternata	BMP 2343 (EGS 46.140)	Citrus, USA	Dang et al. 2015	mycocosm.jgi.doe.gov	Illumina	1269	48.0	212	33.85	50.98	96.95	1.88
A. fragariae	A. alternata f. sp. fragariae	BMP 3062	Strawberry	Dang et al. 2015	mycocosm.jgi.doe.gov	Illumina	779	78.2	126	33.13	50.99	97.95	1.29
A. gaisen	A. gaisen reference	EGS 90.391 (CBS 118488)	Asian Pear, Japan	Woudenberg et al. 2015	contact authors	Illumina	270	529.5	19	35.15	50.90	99.36	0.59
A. gaisen	Pear pathotype	FERA 650	Pear	Armitage et al. 2020	PDWZ02000000	Illumina + Nanopore	27	2110.0	6	34.35	51.14	99.36	0.59
A. limoniasperae	Rough lemon pathotype, moved to A. alternata	BMP 2335 (EGS 44.159)	Citrus	Dang et al. 2015	mycocosm.jgi.doe.gov	Illumina	1596	50.4	213	34.76	50.98	96.95	2.17
A. longipes	A. longipes reference	BMP 0313 (EGS 30.033, CBS 540.94)	Tobacco, USA	Dang et al. 2015	mycocosm.jgi.doe.gov	Illumina	1079	138.5	68	35.68	51.00	99.06	0.94
A. mali	Apple pathoype, moved to A. alternata	BMP 3063	Apple	Dang et al. 2015	mycocosm.jgi.doe.gov	Illumina	2697	21.2	511	33.61	50.78	91.03	5.80
A. mali	Apple pathoype, moved to A. alternata	BMP 3064	Apple	Dang et al. 2015	mycocosm.jgi.doe.gov	Illumina	1669	35.8	297	34.33	50.84	94.43	3.69
A. tangelonis	Tangerine pathotype, moved to A. gossypina	BMP 2327 (EGS 45.080)	Citrus	Dang et al. 2015	mycocosm.jgi.doe.gov	Illumina	1639	37.4	270	33.75	50.87	94.72	3.46
A. tenuissima	A. tenuissima reference, moved to A. alternata	BMP 0304 (ATCC 96828, EGS 34.015, CBS 918.96)	Dianthus, UK	Dang et al. 2015	mycocosm.jgi.doe.gov	Illumina	216	662.1	18	33.36	51.03	99.12	0.76
A. tenuissima	Moved to A. alternata	FERA 1082	Apple	Armitage et al. 2020	PDXA01000000	Illumina	345	547.4	21	33.94	51.14	99.12	0.76
A. tenuissima	Moved to A. alternata	FERA 1164	Apple	Armitage et al. 2020	PDXB01000000	Illumina	250	709.1	15	34.72	50.90	99.24	0.70
A. tenuissima	Apple pathotype, moved to A. alternata	FERA 1166	Apple	Armitage et al. 2020	PDXC01000000	Illumina + Nanopore	22	2583.9	6	35.70	51.08	99.12	0.76
A. tenuissima	Moved to A. alternata	FERA 24350	Asian pear	Armitage et al. 2020	PDXE01000000	Illumina	167	608.1	17	33.07	51.32	99.12	0.76
A. tenuissima	Apple pathotype, moved to A. alternata	FERA 635	Apple	Armitage et al. 2020	PDXF01000000	Illumina	912	318.6	25	36.06	50.92	99.12	0.70
A. tenuissima	Moved to A. alternata	FERA 648	Pear	Armitage et al. 2020	PDXG01000000	Illumina	124	1341.0	10	33.51	50.99	99.36	0.59
A. tenuissima	Apple pathotype, moved to A. alternata	FERA 743	Apple	Armitage et al. 2020	PDXH01000000	Illumina	788	622.1	15	35.92	51.05	99.12	0.70
Average							531.18	903.46	61.03	34.11	51.04	98.36	1.27

Note. New genomes from Canadian isolates are shown in bold type.

^a Nomenclatural changes resulted from synonymization of original species names into A. alternata by Woudenberg et al. (2015).

Genome sequencing of new isolates.—

For sequencing of new genomes, we chose five section Alternaria isolates collected from rye or wheat in Canada between 2013 and 2018 (TABLE 1). Fungal isolates were grown on potato dextrose agar (potato extract, 4 g/L, dextrose, 20 g/L, agar 15 g/L) plates for 14 days at 25 C. Genomic DNA was extracted following Möller et al. (1992) with modifications (see SUPPLEMENTARY MATERIAL for details). The RNA polymerase second largest subunit (rpb2) gene was sequenced following polymerase chain reaction (PCR) amplification using primers and thermocycler conditions described in Woudenberg et al. (2013). Sequences of new isolates (data now shown) were compared with reference sequences for preliminary identification to species or species complex (under the taxonomic framework of Woudenberg et al. [2015]), later confirmed by genome sequencing described below.

Genomic DNA samples were normalized to 300 ng and sheared to approximately 350-bp fragment size using a Covaris M220 sonicator (Woburn, Massachusetts). The fragmented DNA was used as the template to construct PCR-free libraries with a NxSeq AmpFREE Low DNA Library Kit (Lucigen, Middleton, Wisconsin) according to the manufacturer’s instructions. Single-indexed libraries were sequenced with paired-end reads (2 × 150 bp) on a NextSeq 500 instrument (Illumina, San Diego, California) using a NextSeq High Output Reagent Kit (Illumina) according to the manufacturer’s recommendations. An average of 24.2 million reads were generated per sample (SUPPLEMENTARY TABLE 1).

Genome assembly.—

For new assemblies, raw reads were clipped of Illumina adapters and quality trimmed with Trimmomatic (0.38; Bolger et al. 2014). Assemblies were created with SPAdes (3.12; Bankevich et al. 2012), and contigs less than 500 bp in length were removed. See SUPPLEMENTARY MATERIAL for details. Contigs were also removed if BLASTn searches indicated that they were mitochondrial or bacterial in nature. The resulting assemblies had an average sequence coverage of 105× (SUPPLEMENTARY TABLE 1). The final versions of five new genome assemblies are available under National Center for Biotechnology Information (NCBI) accession numbers JAAOQT010000000 to JAAOQX010000000.

Gene prediction (and annotation).—

The previously published genomes are not all in the same state of gene prediction and annotation. Different software/pipelines, or even different versions of the same pipeline, can introduce technical variation that may obscure direct comparison across genomes. To remove this potential source of bias, all assemblies were subjected to new gene prediction using funannotate (1.7.2; https://github.com/nextgenusfs/funannotate, see SUPPLEMENTARY MATERIAL for programs and parameters). In brief, gene prediction parameters were obtained from training with high-quality genomic and transcriptomic data that were available for the same isolate (SRC1lrK2f; Zeiner et al. 2016). Pretrained SRC1lrK2f parameters were used to guide gene prediction in the other genomes.

Genome completeness and macrosynteny.—

The BUSCO (benchmarking universal single-copy orthologs) metric was applied to assess genome completeness using BUSCO (4.0.6; Simão et al. 2015). Genome assemblies (nucleotide level) were compared with the ascomycota_odb10 ortholog data set, which contains 1706 conserved ascomycete orthologs. Whole-genome macrosynteny comparisons were performed with the Nucmer module of MUMmer (3.1; Kurtz et al. 2004; SUPPLEMENTARY MATERIAL). To restrict comparisons to only near-complete genomes, genomes with more than 30 scaffolds were excluded.

Phylogeny.—

Sequence alignment was performed with MUSCLE (3.8; Edgar 2004), and model selection and construction of maximum likelihood phylogenies were performed with IQ-TREE (2.0; Nguyen et al. 2015). Tree construction was performed under best-fit substitution models and partitions. Multispecies coalescent (MSC) analysis was implemented by ASTRAL-III (5.7.3; Zhang et al. 2018). See SUPPLEMENTARY MATERIAL for full details of phylogenetic analyses and tree comparison methods. Previous studies (Woudenberg et al. 2015; Armitage et al. 2020) have shown that the gaisen clade is the earliest diverging of the clades studied here, so it was used for rooting displayed trees.

RESULTS

Genome assembly and gene prediction.—

Summary statistics for the 38 genomes are shown in TABLE 1. Our newly generated assemblies for Canadian section Alternaria isolates were very similar to previously published genomes in regard to GC content (mean = 51.0%) and genome length (mean = 34.11 Mb). Overall, the quality of 38 analyzed genome assemblies was quite high, with a mean of 531 contigs and an N50 of 904 kb. As expected, genomes sequenced using long-read technologies such as PacBio (Menlo Park, California) or Oxford Nanopore (Oxford, UK) tended to have higher qualities than short-read Illumina-based genomes. For example, long-read assemblies had ~14× fewer contigs, and ~5× longer N50 values, than short-read assemblies.

To standardize our downstream comparisons, we subjected each assembly to ab initio gene prediction using the funannotate pipeline. An average of 12 788 gene models per genome were predicted (SUPPLEMENTARY TABLE 2), with long-read genomes having an average of 212 more gene models than short-read genomes.

Assessment of genome completeness.—

The BUSCO metrics were quite high (TABLE 1; SUPPLEMENTARY TABLE 2), indicating that most of these Alternaria assemblies were not missing significant genomic content. Based on the Ascomycota ortholog set (1706 genes), the average percentage of complete genes detected was 98.4%, with an average of only 1.3% of genes missing. Interestingly, the genome with the highest percent missing genes was JS-1623, which was sequenced using long reads and had the highest assembly quality scores. When this outlier was excluded, the long-read genomes had greater average BUSCO scores than the short-read genomes (99.2% versus 98.4% complete, respectively). Investigation of BUSCO genes that were missing from all 38 genomes (SUPPLEMENTARY TABLE 3) did not reveal any shared functional roles.

Genome structure and macrosynteny.—

Patterns of large-scale synteny along chromosomes were compared for the five genomes with near-complete assemblies. The 10 core chromosomes of the NAP07 strain have been matched to the 10 A. solani chromosomes, so the NAP07 genome was used as the reference genome for these analyses. Based on the four pairwise comparisons (FIG. 1), the overall structural arrangement of the genomes was fairly well conserved across strains. For example, although three NAP07 chromosomes (Chr01, Chr02, Chr04) were each matched to two scaffolds in ATCC 34957, there were no large-scale conflicts. This situation likely represents a lack of centromere coverage and failure to join the two chromosome arms in the ATCC 34957 assembly. The opposite situation appears in the JS-1623 comparison: NAP07 Chr01 and Chr08 equals JS-1623 Scaffold 1, and NAP07 Chr04 and Chr10 equals JS-1623 Scaffold 4. Although it is possible that JS-1623 has only eight core chromosomes, the more parsimonious explanation is that chromosomes were simply misjoined in the assembly process due to similarity of telomeric sequences. Despite some potential discrepancies observed for FERA 1166, the general conclusion is that large-scale duplications or rearrangements of core chromosomes were not common during the evolutionary history of the section Alternaria clade.

Figure 1. Circular plots showing macrosynteny between section Alternaria genomes sequenced with long-read technology. Four genomes were compared with the NAP07 genome. NAP07 chromosomes are color-coded and shown on the right side of each plot. Scaffolds from other genomes are in black and shown on the left side of each plot. Tick marks along chromosomes/scaffolds indicate 1 Mb distances

Species phylogeny (concatenated analyses).—

We used three data sets for inferring the species (organismal) phylogeny:

Data set 1—Orthologous proteins. The first approach was to extract all protein sequences and determine single-copy orthogroups with OrthoFinder (2.3.11; Emms and Kelly 2019; SUPPLEMENTARY MATERIAL). Consistent with macrosynteny results, OrthoFinder analyses of section Alternaria revealed very few gene duplication events: 95.9% of orthogroups experienced no duplications across the entire phylogeny. In total, protein sequences from 8750 single-copy orthologs shared by all taxa were used to construct a consensus phylogeny (SUPPLEMENTARY FIG. 1).

Data set 2—BUSCO genes. The second data set was composed of conserved, shared genes identified during BUSCO analyses. Nucleotide sequences for BUSCO genes present in all 38 genomes were extracted, and 983 genes were chosen for the final data set (SUPPLEMENTARY MATERIAL; SUPPLEMENTARY TABLE 4), with a cumulative alignment length of 1.97 Mb and 101.7 × 10³ (5.2%) variable sites. A maximum likelihood (ML) tree was constructed under the best-fit partition model (SUPPLEMENTARY FIG. 1).

Data set 3—Custom gene selection. The third data set was a custom selection of genes that were likely to be the most informative for phylogenetic reconstruction of our sample of strains from section Alternaria. Starting with a high-confidence gene set predicted with evidence from transcriptomic data, each of the 12 860 genes were searched for and extracted from the genomes. Only single-copy genes that were present in all 38 genomes, and contained sufficient phylogenetic signal (i.e., >3% phylogenetically informative sites; SUPPLEMENTARY MATERIAL), were retained. The final custom data set had 6209 genes, with a cumulative alignment length of 8.58 Mb and ~756 × 10³ (8.8%) variable sites. The resulting ML tree constructed under the best-fit partition model is presented in FIG. 2.

Figure 2. Maximum likelihood tree constructed from the custom data set of 6209 genes concatenated. Shading of squares at nodes represents the level of support each node received from bootstrapping, but also from other data sets or analyses methods (see inset). For shading, black indicates that the branch was present and significantly supported, gray indicates that the branch was present but not significantly supported, and white indicates that the branch was absent. Significant support was considered >95% in ML ultrafast bootstrapping for the BUSCO and custom data sets, >95% in posterior probability for the custom MSC data set, and presence in consenses tree for OrthoFinder data set. Numbers on branches indicate the gene concordance factor (GCF) for the 6209 genes trees compared with the displayed species tree. Current species name is indicated in brackets in taxon labels, if different from original species name. ^T indicates ex-type strain

Comparison of branch support values from the three data sets (FIG. 2) revealed that four main clades were identified by all analyses, hereafter referred to as the (i) “alternata/tenuissima” (alt/ten), (ii) “arborescens,” (iii) “longipes,” and (iv) “gaisen” clades. As predicted, the tree constructed from protein sequences (OrthoFinder) provided less resolution within clades, whereas the data set of 6209 custom genes provided the most resolution.

Species phylogeny (analyses of individual loci).—

Joint analyses of the 6209 concatenated genes included a very large amount of sequence data, so the resulting phylogeny was well resolved with nearly every branch receiving 100% bootstrap support. This approach averages over the signals from all genes but obscures the differences that may exist between the evolutionary histories of individual genes. When the sampled taxa are likely to span the species-population interface, one must accept that incomplete lineage sorting may cause discordance between gene trees (Pamilo and Nei 1988) and that gene trees may not match the species tree.

When gene tree discordance is likely to be high, we can apply multispecies coalescent (MSC) approaches that account for incomplete lineage sorting and estimate the species tree from the distribution of gene trees. To more comprehensively explore gene-tree space, and to account for gene-tree uncertainty, 10 nonparametric bootstrap ML trees were constructed for each of the 6209 genes under their respective best-fit models. MSC analysis of the resulting 62 090 gene trees with ASTRAL-III produced a species tree with topology similar to that from concatenated analyses (SUPPLEMENTARY FIG. 1).

Although all analyses support the existence of the four main clades (FIG. 2), the relationships among clades received mixed support. For example, the concatenated custom 6209-gene set suggests that the longipes clade is sister to the arborescens clade [longipes,arborescens], but most evidence supports a sister relationship between the longipes and alt/ten clades [longipes,alt/ten].

Gene concordance factor.—

To further explore gene/species tree concordance, we determined which branches in the species tree were present in each of the individual gene trees (gene concordance factor, GCF). As shown in FIG. 2, most of the main branches that received full support in the concatenated ML and MSC analyses were supported by a large proportion of the gene trees. The most notable exception was the branch for the alt/ten clade, which was supported by only 9.8% of individual gene trees. The placement of the longipes clade again received mixed support, with [longipes,alt/ten] and [longipes,arborescens] being supported by 12.4% and 11.1% of individual gene trees, respectively.

Genomic distribution of gene/species tree discordance.—

To visualize the relationship between gene/species tree concordance and genomic location, individual genes were scored for their support of the four main branches and the patterns were plotted against physical coordinates along chromosomes (FIG. 3). At the chromosome level, the proportion of analyzed genes that supported a particular main branch could vary quite drastically. For example, the gaisen clade was supported by 94.8% of genes from chromosome 6, but only 52.5% of genes from chromosome 8. When analyzed at a finer scale, gene/species tree discordance was not evenly distributed along individual chromosomes. Genomic clustering of discordance was most apparent for the gaisen clade where, for example, chromosomes 1 and 2 both had much of their gene/species discordance concentrated into regions of ~0.75 Mb in length (FIG. 3). Discordance for the arborescens and longipes clades was more evenly distributed, but they still had regions where continuous stretches of analyzed genes failed to support the clade branch. The distribution for the alt/ten branch was fundamentally different because gene/species discordance was the rule rather than the exception, but there were still peaks of genomically colocated concordance.

Figure 3. Plots of branch-specific concordance factor as a function of genomic location. Lines represent moving averages based on a sliding window of 20 genes. A concordance factor of 100 indicates that the branch was present in all included gene trees, whereas 0 indicates that the branch was absent in all included gene trees. The mean for the entire chromosome is indicated within each panel

We examined the gene trees located in some regions of high overall gene/species discordance to reveal the phylogenetic history of those regions (e.g., 34 genes in the ~154 kb region centered around the 5.9 Mb mark of the chromosome 1; 23 genes in the last ~185 kb at one end of chromosome 5). Although the topology of the trees varied quite drastically, a common observation was that the gaisen clade often intermingled with the alt/ten clade and the longipes clade often intermingled with the arborescens clade. ASTRAL MSC consensus trees produced from these discordant genomic regions are shown in SUPPLEMENTARY FIG. 2.

Putative diagnostic markers for detection and identification of main lineages.—

Now that the main lineages were circumscribed by genome-wide data, we searched for diagnostic molecular markers that were consistent with this framework. To support accurate identification, we chose molecular markers that contained sequence polymorphisms that could distinguish among all four of the main section Alternaria lineages (i.e., genes that had polymorphism arise, sort and fix between the divergence events giving rise to the lineages). We focused on regions of the genome that contained multiple consecutives genes that were fully concordant with the species phylogeny. TABLE 2 lists our selection of top 20 candidate diagnostic genes. We were able to find appropriate candidate genes from all chromosomes except chromosome 8, which is consistent with chromosome 8 having the lowest mean concordance factor for two of the four main branches (FIG. 3). Genes were chosen irrespective of their function, and information on predicted gene products provided no evidence for any obvious trends of functional enrichment.

Table 2. List of 20 candidate genes that are diagnostic for the four analyzed clades in section Alternaria

Locus	Related NCBI accession	Gene product	Chromosome	Gene midpoint (NAP07 coordinate)
ASA-01	XM_018527094	ClpP/crotonase (caseinolytic protease)	Chr01	804 003
ASA-02	XM_018532566	Glycoside hydrolase	Chr01	3 419 152
ASA-03	XM_018525396	Thermostable alpha-amylase	Chr01	4 422 548
ASA-04	XM_018526656	Hypothetical protein	Chr02	1 798 824
ASA-05	XM_018523574	Hypothetical protein	Chr02	2 280 898
ASA-06	XM_018525342	MFS general substrate transporter (major facilitator superfamily)	Chr02	5 311 887
ASA-07	XM_018527519	Hypothetical protein	Chr03	893 828
ASA-08	XM_018532806	RTA1-domain-containing protein	Chr03	1 251 119
ASA-09	XM_018524513	Hypothetical protein predicted sulfite efflux pump	Chr03	1 359 323
ASA-10	XM_018531969	Hypothetical protein	Chr04	819 652
ASA-11	XM_018533777	Class I glutamine amidotransferase-like protein	Chr04	2 294 759
ASA-12	XM_018527323	Homogentisate 1,2-dioxygenase	Chr05	650 876
ASA-13	XM_018534328	Hypothetical protein	Chr05	2 032 935
ASA-14	XM_018534844	Hypothetical protein	Chr06	773 854
ASA-15	XM_018524984	Non-ribosomal peptide synthetase C-terminal	Chr07	1 253 208
ASA-16	XM_018529476	Hypothetical protein	Chr07	2 178 554
ASA-17	XM_018527994	Acetyl-CoA synthetase-like protein	Chr09	1 007 905
ASA-18	XM_018524187	AhpC/TSA family protein-like protein (alkyl hydroperoxide reductase/thiol specific antioxidant)	Chr09	1 421 701
ASA-19	XM_018524743	Histone-fold-containing protein	Chr10	828 763
ASA-20	XM_018533177	zf-ZPR1-domain-containing protein (zinc-finger domain)	Chr10	1 038 833

For comparison, we also examined existing loci that previously have been applied to the molecular systematics of section Alternaria (Hong et al. 2005; Andrew et al. 2009; Woudenberg et al. 2015; Armitage et al. 2015). Gene trees built from extracted sequences for 15 of these loci (TABLE 3) were scored for the presence of main clades. The rDNA region is typically informative at deeper phylogenetic depths, so as expected, the 28S nuc rDNA, internal transcribed spacer (ITS), and 18S nuc rDNA markers provided very little resolution and did not support any of the four main clades. Only one of the 15 gene trees had the alt/ten clade present, consistent with our previous finding that less than 10% of loci support that branch. Loci that were designed specifically for Alternaria provided greater resolution, but only one locus (Tma22) supported all four main clades. Therefore, most existing loci may be able to place an unknown strain into a section, but section-specific molecular markers are required to distinguish lineages within section Alternaria.

Table 3. Existing markers and their ability to discriminate the four analyzed clades in section Alternaria.

				Gene tree supports clade?
Divergence level	Locus	Annotation	Chromosome	gaisen	alt/ten	arborescens	longipes
Deep	28S	28S nuclear rRNA (large subunit)	Chr07	N	N	N	N
Deep	ITS	Internal transcribed spacer (ITS1-5.8S-ITS2)	Chr07	N	N	N	N
Deep	18S	18S nuclear rRNA (small subunit)	Chr07	N	N	N	N
Medium	gapdh	Glyceraldehyde-3-phosphate dehydrogenase	Chr06	Y	N	N	Y
Medium	rpb2	RNA polymerase, second largest subunit	Chr05	Y	N	Y	Y
Medium	tef1	Translation elongation factor 1-alpha	Chr03	N	N	N	Y
Shallow	Alta1	Major allergen Alt a 1	Chr02	Y	N	Y	Y
Shallow	endoPG	Endopolygalacturonase	Chr01	Y	N	Y	Y
Shallow	KOG1058	Coatemer subunit beta	Chr10	Y	N	N	N
Shallow	KOG1077	AP-2 adaptor complex subunit alpha	Chr01	Y	N	Y	Y
Shallow	OPA10-2	Anonymous region (putative FAD-binding-like protein)	Chr01	Y	N	Y	Y
Shallow	OPA1-3	Anonymous region	Chr01	N	N	N	Y
Shallow	Pgs1	CDP-diacylglycerol-glycerol-3-phosphate 3-phosphatidyltransferase	Chr07	Y	N	N	N
Shallow	Rev3	Catalytic subunit of DNA polymerase zeta	Chr04	Y	N	Y	Y
Shallow	Tma22	Translation machinery associated protein	Chr05	Y	Y	Y	Y

Based on genome coordinates alone (SUPPLEMENTARY TABLE 5), many of these loci are not predicted to support the species tree because they are located in genomic regions with low gene/species tree concordance (FIG. 3). The one locus (Tma22) that supported the species tree was within a genomic region of high concordance and, coincidently, was directly adjacent to one of our chosen candidate diagnostic genes (ASA-13).

Variation in genome content.—

The previous analyses were based on genes that were present in the genomes of all sampled taxa. OrthoFinder analyses revealed that 6004 (38.8%) of 15 479 orthogroups were missing from at least one strain, and 99.3% (5964) of these were shared between multiple strains. Although macrosynteny results indicate that large-scale variation in core chromosome structure may be rare in section Alternaria, smaller-scale variation in genome content (i.e., one to a few genes) appears to be more common.

Although a detailed pan-genomic analysis is outside the scope of this study, we did search for genes or genomic regions that were potentially diagnostic for the major clades. Based on filtered OrthoFinder results, we identified 113 and 69 orthogroups whose presence or absence, respectively, was specific to one of the four major clades (i.e., present in all strains in one clade, and absent from all other strains; FIG. 4). Consistent with the relative branch lengths in the species phylogeny, the gaisen clade possessed the highest number clade-specific orthogroups (53), whereas the alt/ten clade had the lowest (4). Patterns of clade-specific orthogroup absence were very similar to that of orthogroup presence. The functional annotations of clade-specific orthogroups were investigated for the two clades with the largest samples sizes (alt/ten and arborescens). We identified putative functions or conserved protein domains for 55% of these genes (SUPPLEMENTARY TABLE 6), but no clear trends in function were observed. When considering orthogroups that were shared exclusively by all members of two clades (FIG. 4), the alt/ten and longipes clades shared the highest number (6), supporting a [longipes,alt/ten] sister relationship.

Figure 4. Numbers of orthogroups that were clade specific (present/absent). Numbers within circles indicate orthogroup patterns for single clades (e.g., present in all taxa in clade but absent from all other taxa/absent from all taxa in clade but present in all other taxa). Numbers between circles indicate orthogroup patterns shared between two clades

When the genomic locations of clade-specific genes were examined, little evidence of physical clustering was found. In total, 163 of the 182 (89.6%) genes were singletons that were not physically adjacent to other clade-specific genes. The mean and maximum numbers of physically consecutive clade-specific genes were 1.12 and 4, respectively, indicating that clade-specific genomic regions are very limited in size.

DISCUSSION

Review of genomic resources for Alternaria section Alternaria.—

Our first goal was to provide a brief review of the genomic resources that are currently available for Alternaria section Alternaria. Over the past 9 years, Alternaria researchers have generated and released genomes for strains sampled from a variety of hosts, with a focus on pathotypes infecting citrus, apple, and pear (see TABLE 1 for list of references). Although infection of cereals and contamination of harvested grains is a major problem caused by species in section Alternaria (Tralamazza et al. 2018; Tittlemier et al. 2019), there were no genomes available for strains sampled from major agricultural crops such as wheat, oat, barley, or rye. Here, we sequenced and assembled five new genomes from strains from wheat and rye, representing the first genomes produced for section Alternaria strains of Canadian origin.

The quality of the 38 analyzed genome assemblies was quite high, with some of the long-read-based genomes nearing the complete chromosome level, allowing for macrosynteny analyses and the detection of large-scale chromosomal rearrangements. Although we made only four pairwise comparisons, the structure of the genomes appeared to be fairly well conserved across section Alternaria. Note that we focused on core chromosomes and did not specifically analyze accessory chromosomes, also known as conditionally dispensable or supernumerary chromosomes. Although these accessory chromosomes are very interesting and often play important roles in conferring host-specific pathogenicity (Hu et al. 2012; Wang et al. 2019), a detailed analysis of their evolution and gene content will be performed in a future study. Recent work in other fungi have found that large-scale variation in genome content can be quite common and contribute to adaptive evolution (Möller and Stukenbrock 2017). A possible explanation for high macrosynteny observed here is that the depth of phylogenetic divergence within section Alternaria is relatively low, providing only a short time frame for major structural changes. When similar analyses are extended across additional sections of the Alternaria genus, more evidence of large-scale genome content variation may be found.

Phylogenomics of Alternaria section Alternaria.—

Another goal was to build a robust phylogeny for section Alternaria using a phylogenomics approach. Most previous phylogenetic work has relied on sequences from up to seven independent loci (Andrew et al. 2009; Lawrence et al. 2013; Armitage et al. 2015; Woudenberg et al. 2015), but a recent genome-based study included 500 BUSCO genes (Armitage et al. 2020). To fully explore the phylogenetic signals present in the genomic data, we compared multiple data sets and tree-building approaches.

An advantage of the OrthoFinder approach is that it allows sample-specific ortholog sets to be identified using established, user-friendly software packages. A disadvantage is that it relies heavily on reading frame and protein prediction, processes that may not be accurate or precise for the taxa under study. BUSCO genes are predesignated as single-copy orthologs across a large taxonomic group, in this case, all ascomycete fungi. Thus, the same set of genes can be used for molecular systematics of any ascomycete fungi, making comparisons across studies and taxa possible. However, BUSCO genes are by definition highly conserved, so their phylogenetic informativeness decreases as the taxonomic level decreases. Much like the OrthoFinder approach, BUSCO genes may not possess sufficient sequence variation to be useful at phylogenetic divergence levels found in section Alternaria (i.e., species-population interface). To maximize the phylogenetic informativeness of included loci, we also created a 6209-gene data set that was customized to section Alternaria. The advantages of this method are that it requires confident gene prediction for only one representative genome, which is often the case for many taxonomic groups, and it is well suited for large-scale studies that utilize mainly draft genomes without annotation.

One consistent pattern that was shared among all data sets and phylogenetic analyses was the well-supported existence of four main clades (FIG. 2): (i) alternata/tenuissima (alt/ten), (ii) arborescens, (iii) longipes, and (iv) gaisen clades. However, the relationships among the four distinct clades received mixed support. For example, the concatenated custom 6209-gene set suggests [longipes,arborescens], but most other analyses support [longipes,alt/ten]. The [longipes,alt/ten] grouping was recovered by Armitage et al. (2020) with a genome sample that overlaps significantly with ours, albeit using a different gene set and analysis method.

Gene/species tree discordance.—

Visualization of gene/species tree discordance (FIG. 3) revealed the dynamic nature of the evolutionary histories of specific regions along the length of chromosomes. Although the genomic regions displaying gene/species tree discordance could have arisen due to more complex evolutionary scenarios (e.g., horizontal gene transfer, genome duplication and loss), our findings of high macrosynteny and low rates of gene duplication provide strong evidence that our 6209 chosen loci are orthologous within our sample and therefore accurately reflect the evolutionary histories of the organisms. Since genes with low levels of informative polymorphism result in poorly resolved trees, we purposely excluded genes with low phylogenetic signal from our custom gene set. This screening criterion prevented low-signal genes from interfering with the high-signal genes during consensus and tree comparison analyses (i.e., a poorly supported polytomy should not detract from well-supported monophyly). Our approach is an extension of genealogical concordance concepts originally applied to species delineation using multilocus sequence data (Taylor et al. 2000; Dettman et al. 2003). Newly diverged species lose ancestral polymorphism via genetic drift, which may eventually lead to reciprocal monophyly. If effective population sizes are large, and time between divergence events is short, these random processes are expected to result in exclusive reciprocal monophyly for only a small proportion of genes (Pamilo and Nei 1988; Mehta and Rosenberg 2019).

The comparative levels of incomplete lineage sorting and resulting nonmonophyly can be estimated from the summed, genome-wide GCF of the section Alternaria clade-specific branches (FIG. 2). The alt/ten branch had the lowest GCF (9.8), and the internal branches within the alt/ten clade also had low GCF compared with other clades, suggesting that the alt/ten clade has had a larger historical effective population size. The alt/ten clade being large, widespread, and diverse is consistent with the frequency of occurrence of species in natural settings and surveys: the majority of small-spored Alternaria strains collected worldwide fall into the alt/ten clade (e.g., Zheng et al. 2015; da Cruz Cabral et al. 2017), although this may be due in part to difficulties with definitive taxonomic identification within section Alternaria. The other three clades are sampled less commonly and appear to have smaller effective population sizes, allowing for more rapid fixation of polymorphisms during divergence and higher probabilities of genealogical coalescence within clades. Although we included most publicly available genomes, the small sample size for some clades is a clear limitation of this study. The addition of new genomes may affect GCFs, but the magnitude of the effect, and which branches it would affect, would depend on the nature of the strains themselves. For example, adding a genome from a strain that is a first-generation hybrid between two clades would have a greater effect than a strain that falls cleanly into one clade. We look forward to expanding the genome sampling of underrepresented lineages and testing how well this framework holds up when challenged with more data.

Although no sexual state has been observed for species in section Alternaria, there is strong genetic evidence that recombination does occur in these taxa (Stewart et al. 2014; Meng et al. 2015). The high levels of gene tree discordance observed here are consistent with recombination within clades. In particular, the alt/ten clade displays high levels of within-clade diversity, and the within-clade branches typically have very low GCF values. Collectively, most evidence suggests that the alt/ten clade represents a single large, recombining species with wide geographic distribution. Further population-level studies are required to assess recombination in the other main clades.

Taxonomy and nomenclature of Alternaria section Alternaria.—

When erected by Lawrence et al. (2013), section Alternaria contained ~60 small-spored catenulate species, including well-studied species such as A. alternata, A. tenuissima, A. arborescens, and A. gaisen. Later, a comprehensive revision by Woudenberg et al. (2015) concluded that section Alternaria was composed of only 11 phylogenetic species and one species complex. Here, we describe the correspondence between named taxa and the four main clades identified in our phylogenomic analyses.

i) alternata/tenuissima (alt/ten) clade. Most of the 26 strains included in this clade were originally described as A. alternata or A. tenuissima (see TABLE 1). High plasticity and overlap among morphological characters makes reliable and consistent differentiation of these morphospecies very difficult. Our phylogenomic analyses failed to support a distinction between the two taxa, in agreement with most molecular phylogenetic studies involving these two morphospecies (Andrew et al. 2009; Armitage et al. 2015, 2020; Woudenberg et al. 2015; da Cruz Cabral et al. 2017). Highlighting the lack of phylogenetic divergence, Woudenberg et al. (2015) synonymized 35 morphospecies into A. alternata, including A. tenuissima. The name A. alternata was retained because it was the name of the type species for the genus and was well established in the research literature. Given the large body of supporting evidence, we also recommend that the name A. tenuissima be retired from use and be replaced by A. alternata instead.

As with many other plant pathogenic fungi, Alternaria strains that infect specific host plants or produce host-specific toxins have been described as formae speciales or pathotypes. In many cases, these strains were raised to the status of formally named species based on slight morphological differences (e.g., Simmons 1999a). Three examples of this are included in our alt/ten clade, all of which were named after their host associations: A. mali (apple), A. citriarbusti (citrus), and A. limoniasperae (lemon). Molecular studies of such named strains have shown that they are intermixed with other species and do not represent genetically distinct taxa (Peever et al. 2004; Rotondo et al. 2012). Accordingly, A. mali, A. citriarbusti, and A. limoniasperae were among the 35 morphospecies synonymized into A. alternata (Woudenberg et al. 2015). Our phylogenomic results support these large-scale synonymizations and the idea that apparent host-specific pathogenicity is a characteristic of individual strains rather than species.

For the sake of nomenclatural stability, we recommend that all strains in the alt/ten clade be referred to as A. alternata, with the recognition of host-specific pathotypes within the species.

ii) arborescens clade. A. arborescens is characterized morphologically by more frequent branching of the conidiophore and conidial chains, producing a three-dimensionally arborescent appearance (Simmons 1999b). This species was designated to contain strains causing stem canker of tomato (formerly known as the tomato pathotype, or A. alternata f. sp. lycopersici), but subsequent work has shown that host specificity is lacking and the species may be isolated from a wide range of hosts (Woudenberg et al. 2015). Phylogenomic analyses placed a strain named A. mali into this clade, likely representing an A. arborescens strain that was named A. mali simply because it was sampled from an apple host.

The strong phylogenomic support observed here for the distinct arborescens clade is consistent with previous molecular phylogenetic studies (e.g., the “arborescens clade” identified by Armitage et al. [2020]). Woudenberg et al. (2015) found that multilocus sequence data could not resolve the relationships between A. arborescens and ex-type strains from three other closely related morphospecies (A. cerealis, A. geophila, and A. senecionicola) and collectively considered them as the A. arborescens species complex. Support for an A. arborescens species complex is found when comparing the GCF values for the within-clade branches: values for the arborescens clade are greater than that for alt/ten clade, suggesting that the former clade has lower rates of recombination and higher levels of phylogenetic substructure. We note that uneven sampling across the main clades may introduce bias, so further sampling of the arborescens clade is required before meaningful comparisons of within-clade diversity can be made. Our newly sequenced strain DET2008 appears to be early diverging and may represent a previously undescribed species, pending increased sampling and analyses. Furthermore, better representation of the other species within the complex are needed to determine whether they represent multiple distinct species or one larger species. Until these issues are resolved, we follow Woudenberg et al. (2015) and recommend using the name A. arborescens species complex.

iii) longipes clade. The two strains comprising this clade are from closely related, sister species. One strain was from A. longipes, representing the tobacco pathotype of small-spored Alternaria. The second strain was originally identified as A. tangelonis (tangerine pathotype), however, A. tangelonis was synonymized into A. gossypina based on comparative analyses of ex-type strains (Woudenberg et al. 2015). The placement of a tangerine pathotype in the longipes clade highlights the fact that citrus and apple pathotypes are typically polyphyletic and cannot be assumed to represent distinct species (Armitage et al. 2020).

iv) gaisen clade. A. gaisen encompasses strains that are specifically pathogenic on Japanese pear. This clade included a genome that was accessioned as A. fragariae, named after the fact the strain was isolated from strawberry. However, this is not a formally named species and should actually be A. alternata f. sp. fragariae (Simmons 1999b). Recent work has found that some strawberry and pear pathotypes are morphologically and phylogenetically indistinguishable and should both be included within A. gaisen (Nishikawa et al. 2019). Once again, the apparent misplacement of taxa can be explained by erroneous identifications based on host association.

Our study highlights how nomenclature at the point of submission to public sequence databases typically does not get updated when taxonomic revisions occur. Of the 32 genomes obtained from public databases, 14 (44%) of them are accessioned under species names that have been synonymized or changed (sensu Woudenberg et al. 2015). However, this is not just a problem for retroactive correction of names in existing database accessions. Revisions in nomenclature or taxonomy can take time to be incorporated into the research community. For example, even though evidence shows that A. tenuissima is not a distinct species, and was synonymized into A. alternata in 2015, the tenuissima epithet continues to be frequently used in publications even today.

Future directions.—

Our list of candidate diagnostic molecular markers was developed for distinguishing among the four main section Alternaria clades identified in our phylogenomic analyses. Appropriate marker choice is critical for accurate and reliable methods of detection and identification, and these tools are critical for dealing with quarantine and restricted species that impact international trade. Future work is required to test how well these newly proposed markers perform when challenged by a wider range of section Alternaria strains, including strains from other species not studied here. Ideally, we hope to identify markers that are consistent within clades but also diagnostic across the entire section Alternaria.

The five newly sequenced genomes were all from isolates sampled from cereal crops of Canadian origin, but they did not form a monophyletic clade that corresponded to host or geography. Given the observed diversity among these five isolates, it is clear that our sampling of section Alternaria biodiversity in Canada is far from complete, and we will continue to characterize the species of Alternaria infecting a wide range of agricultural hosts.

Our analyses included whole genome sequences for as many section Alternaria species as were available; however, section Alternaria contains at least seven other species for which genome sequences do not currently exist. To characterize the biodiversity of this important group of fungi, genome sampling should be expanded to encompass the full phylogenetic diversity of small-spored Alternaria. Increased sampling of the arborescens clade will help to determine whether it represents a species complex, and proper population samples will allow more comprehensive tests for recombination.

ACKNOWLEDGMENTS

We thank Allen Xue and Yuanhong Chen for providing access to grain samples for fungal isolation; Keith Seifert for providing a fungal strain; Joyce Woudenberg and Ewald Groenewald for sharing sequence data; Kasia Dadej and staff at the Agriculture and Agri-Food Canada (AAFC) Molecular Technologies Laboratory for sequencing support; and Hai Nguyen for reviewing an early draft of the manuscript.

supplementary material

Supplemental data for this article can be accessed on the publisher’s Web site.

Additional information

Funding

This study was supported by Agriculture and Agri-Food Canada under grant J-002272.