Genetic and Toxigenic Variability within Aspergillus flavus Population Isolated from Maize in Two Diverse Environments in Kenya

Aspergillus flavus is the main producer of carcinogenic aflatoxins in agricultural commodities such as maize. This fungus occurs naturally on crops, and produces aflatoxins when environmental conditions are favorable. The aim of this study is to analyse the genetic variability among 109 A. flavus isolates previously recovered from maize sampled from a known aflatoxin-hotspot (Eastern region, Kenya) and the major maize-growing area in the Rift Valley (Kenya), and to determine their toxigenic potential. DNA analyses of internal transcribed spacer (ITS) regions of ribosomal DNA, partial β-tubulin gene (benA) and calmodulin gene (CaM) sequences were used. The strains were further analyzed for the presence of four aflatoxin-biosynthesis genes in relation to their capability to produce aflatoxins and other metabolites, targeting the regulatory gene aflR and the structural genes aflP, aflD, and aflQ. In addition, the metabolic profile of the fungal strains was unraveled using state-of-the-art LC-MS/MS instrumentation. The three gene-sequence data grouped the isolates into two major clades, A. minisclerotigenes and A. flavus. A. minisclerotigenes was most prevalent in Eastern Kenya, while A. flavus was common in both regions. A. parasiticus was represented by a single isolate collected from Rift Valley. Diversity existed within the A. flavus population, which formed several subclades. An inconsistency in identification of some isolates using the three markers was observed. The calmodulin gene sequences showed wider variation of polymorphisms. The aflatoxin production pattern was not consistent with the presence of aflatoxigenic genes, suggesting an inability of the primers to always detect the genes or presence of genetic mutations. Significant variation was observed in toxin profiles of the isolates. This is the first time that a profound metabolic profiling of A. flavus isolates was done in Kenya. Positive associations were evident for some metabolites, while for others no associations were found and for a few metabolite-pairs negative associations were seen. Additionally, the growth medium influenced the mycotoxin metabolite production. These results confirm the wide variation that exists among the group A. flavus and the need for more insight in clustering the group.


INTRODUCTION
The filamentous fungus Aspergillus flavus is a cosmopolitan soil-born saprophytic organism with opportunistic parasitic behaviors to plants, animals and humans. A. flavus is the most common causal agent of aflatoxin (AF) contamination in food and feed, and the main producer of carcinogenic aflatoxin B1 (AFB1), (International Agency for Research on Cancer (IARC), 2002). A. flavus is known as the second cause of human invasive aspergillosis after A. fumigatus (Denning, 1988;Latge, 1999;Hedayati et al., 2007). Variability in the A. flavus phenotype exists, for example, in sclerotia formation, culture characteristics, and ability to produce aflatoxins. Reports illustrate that other metabolites can be produced with an unraveled toxicological profile (Pildain et al., 2008;Nicholson et al., 2009;Ehrlich and Brian, 2014). Culture and molecular methods have been used to group isolates into aflatoxigenic and non-aflatoxigenic strains (Zarrin and Erfaninejad, 2016), and to investigate the variability in fungal species and subspecies (Geiser et al., 1998a(Geiser et al., ,b, 2000Frisvad et al., 2007;Jurjević et al., 2015). Frisvad et al. (2005) raised A. flavus var. parvisclerotigenus to another species level, based on its atypical S-type sclerotia, profile of excreted metabolites and difference in rDNA-sequence. Geiser et al. (2000) subdivided A. flavus into Group 1 and Group 11 on the basis of sclerotia type and aflatoxin profile. Several authors have given evidence that A. flavus sensu lato may consist of several species (Geiser et al., 1998a(Geiser et al., ,b, 2000Pildain et al., 2008). The variation within the A. flavus isolates suggests the need for an extended taxonomic delineation for the benefit of the development of effective and safe biocontrol application technologies.
Among the mycotoxins produced by A. flavus, the production of aflatoxins have been studied more with the objective of understanding the diversity of the fungus, management of aflatoxin production by applying biocontrol technique and flow of genes involved in aflatoxin synthesis within a population. The aflatoxin biosynthetic pathway in A. flavus and A. parasiticus, involves approximately 30 genes clustered together in a 70-kb DNA region of the fungal genome on chromosome 111 roughly 80-kb away from telomere (Yu, 2012). The order of genes within the cluster is reported to be highly conserved wihin Aspergillus section Flavi. Among the 30 genes identified the functions of 19 have been assigned (Yu and Ehrlich, 2011). PCR detection of aflatoxin biosynthetic gene presence or expression has been used as diagnostic tool for aflatoxigenic fungi in selected commodities (Geisen, 2007). Sequence variability and deletions in various genes/regions of the aflatoxin biosynthetic cluster have also been used to determine the polyphyletic assemblages of A. flavus genotypes (Chang et al., 2005(Chang et al., , 2006Mohankumar et al., 2010;Gonçalves et al., 2012;Al-Wadai et al., 2013). Gene deletions and presence of single nucleotide polymorphisms (SNPs) in aflatoxin biosynthetic genes have been associated with A. flavus inability to produce aflatoxins. (Solorzano et al., 2014) observed genetic diversity among a population of A. flavus in a field in Mississipi over years and reported mutations among toxigenic isolates resulting into loss of ability to produce aflatoxins. Horizontal transfer of this gene cluster has also been documented (Susca et al., 2017). In this study the presence of the structural genes aflD (nor-1) which encodes the ketoreductase needed for the conversion of the 1'-keto group in norsolorinic acid (the first stable aflatoxin precursor) to the 1 ′ -hydroxyl group of averantin; aflP (omtA) involved in the conversion of sterigmatocystin (ST) to O-methylsterigmatocystin (OMST) and demethylsterigmatocystin (DMST) to dihydro-O-methylsterigmatocystin (DHOMST); aflQ (ordA) involved in the conversion of OMST to AFB 1 and AFG 1 and of DHOMST to AFB 2 and AFG 2 ; and the regulatory gene (aflR) are determined in A. flavus strains collected from Eastern and Rift valley regions of Kenya.
Kenya has experienced several fatal acute aflatoxicosis outbreaks following consumption of contaminated food and feed (Okoth, 2016). Most of these outbreaks have occurred in the Eastern part of the country, after consumption of home-grown maize. This region is reported to harbor a large population of the S-strain of A. flavus (Okoth et al., 2012). This study sought to establish the variability patterns of A. flavus strains previously isolated from the Eastern and Rift Valley regions of Kenya using sequence analysis of the ITS region, and parts of the β-tubulin and calmodulin genes and metabolite profiles (Okoth et al., 2012). The Rift Valley is the major commercial maize-growing area of Kenya, and no aflatoxicosis outbreaks have been reported there. Such information is important given that aflatoxin production is a polygenic variable trait that may experience outcrossing between two strains during reproduction. Reproduction in A. flavus is mostly asexual, however sexual reproduction is also known to occur together with non-sexual recombination among isolates that are highly similar. The mechanisms are not clear, and can only be deduced by population studies to unravel the diversity of the fungus. This is the first comprehensive analysis of A. flavus isolates in Kenya using molecular analysis and metabolic profiling. The possible variations within A. flavus isolates is explained here using the sequence analysis of the conserved ITS region of the rDNA gene cluster and parts of the calmodium and β-tubulin genes, presence of four genes in the aflatoxin biosynthetic pathway and metabolite profile.

Source of Aspergillus Flavus Isolates
The isolates used in this study were previously recovered from maize collected from farmers in the Rift valley (n = 255) and the Eastern regions of Kenya (n = 258) in 2010, from a larger study carried out by Okoth et al. (2012) and Nyongesa et al. (2015). The isolates were identified as A. flavus using morphological and cultural characteristics. The isolates were preserved on silica gel at the University of Nairobi. From this collection, 109 isolates [n = 53 (Eastern Kenya), n = 56 (Rift Valley)] were randomly selected, and grown at 28 • C on potato dextrose agar (PDA).

Reagents and Chemicals
The individual mycotoxin solid calibration standards (1 mg) of aflatoxin B1(AFB1), aflatoxin B2 (AFB2), aflatoxin G1 (AFG1), aflatoxin G2 (AFG2), and zearalanone (ZAN) (internal standard) were obtained from Sigma Aldrich (Bornem, Belgium) Zearalanone was used as internal standard to quantify the AFs, but also to feature as identification criterium for the other metabolites. One of the identification criterium is the relative retention time. The retention time of the internal standard is compared to the retention time of the metabolites, and should remain constant.
Frontiers in Microbiology | www.frontiersin.org profile. A loop of mycelia and spores of the isolates scraped from 9 mm PDA plugs were inoculated on 9 cm diameter Petri plates, and incubated at 28 • C for 7 days in the dark. Using a 9 mm diameter sterile cork borer, 9 plugs were harvested uniformly from the plates into 50 mL propylene tubes. The agar was submitted to acetonitrile maceration (5 mL) overnight (15 h). Then, metabolites were extracted from the matrix with 10 mL of water and 20 mL of dichloromethane while shaking for 30 min. The mixture was centrifuged (3,000 g, 15 min), 10 mL from the bottom face was evaporated and resuspended with 1 mL of methanol/acetonitrile/water (30/30/40, v/v/v).

LC-MS/MS Methodology
An LC-MS/MS method was developed and validated for the analysis of 32 analytes ( Table 1). A Waters Acquity UPLC system coupled to a Xevo TQS mass spectrometer (Waters, Milford, MA, USA) was used to analyse the samples, equipped with MassLynx R (version 4.1) and QuanLynx R (version 4.1) software (Waters, Manchester, UK) for data acquisition and processing. A ZORBAX Eclipse XDB C18-column (1.8 µm, 100 × 2.1 mm) was applied (Agilent Technologies, Diegem, Belgium). The mobile phase consisted of water/methanol (95/5, v/v (A)) and methanol/water (95/5, v/v (B)), both adjusted with 5 mM ammonium acetate and 0.1% formic acid, at a flow rate of 0.4 mL/min. The gradient elution programme started at 100% mobile phase A for 0.5 min. Then, the mobile phase B increased with a linear increase to 99% in the 20 min. The mobile phase B was kept at 99% for 5 min. The mobile phase linearly decreased till 0% for 1 min. Mobile phase A (100%) was running for 4 min. The duration of each LC run was 30 min, including re-equilibration. The capillary voltage was 3 kV, and nitrogen was applied as spray gas. Source and desolvation temperatures were set at 120 • and 400 • C, respectively. The argon collision gas pressure was 9 × 10 −6 bar, the cone gas flow 35 L/h and the desolvation gas flow 800 L/h. For increased sensitivity and selectivity, the instrument was operated in the selected reaction monitoring (SRM) mode, and two SRM transitions were monitored for each analyte. The SRM-transitions for every analyte are described in Table 1. Matrix-matched calibration plots were constructed for the determination of the AFB1, AFB2, AFG1, and AFG2. ZAN was used as internal standard in the multi-mycotoxin analysis. The limit of detection (LOD) was calculated as three times the standard error of the intercept, divided by the slope of the standard curve; the limit of quantification (LOQ) was computed in a similar way, except for the standard error, which was by a factor of six. The calculated LOD and LOQ were verified by the signal-to-noise ratio (s/n), which was more than 3 and 10, respectively, in accordance with the IUPAC guidelines (IUPAC, 1995).

DNA Extraction
Conidia of the A. flavus isolates grown on PDA for 5 days were inoculated into 250 mL Erlenmeyer flasks containing 100 mL of Wickerham medium (40 g/L glucose, 5 g/L peptone, 3 g/L yeast extract, 3 g/L malt extract). Incubation was carried out at 25 • C under shaking conditions (150 rpm) for 2 days. Then, mycelia were harvested by filtration and washed with 0.5 M ethylenediamine tetra-acetic acid (EDTA) & sterile distilled water (dH 2 O), and freeze-dried at -70 • C for DNA extraction. The mycelia were ground into a fine powder using a pestle and mortar. The powder (∼100 mg) was transferred into a 2 mL sterile eppendorf tube containing 2 mm diameter glass beads. Then, 650 µL of lysis CTAB (cetyltrimethylammonium bromide) buffer was added and the mixture was ground in a tissue miller at a frequency of 30/s for 3 min after which the samples were FIGURE 2 | Co-occurrence of the observed metabolites (A CYA and B YESA).
A positive association means that in case one metabolite is observed it is very likely that also the other metabolite will be present. A negative association means that in case one metabolite is observed it is very unlikely that also the other metabolite will occur. A random association means that nothing can be said concerning the co-occurrence of both metabolites.
incubated at 65 • C for 1 h. Proteins were precipitated by adding 600 µL phenol, inverted gently and centrifuged at 14000 rpm for 20 min. The top aqueous layer was pipetted into a new tube. Afterwards, 600 µL phenol/chloroform (25/24, v/v) was added, inverted gently and centrifuged at 14000 rpm for 20 min. This procedure was repeated twice with 600 µL of chloroform. The clean supernatant was transferred to a new 1.5 mL eppendorf tube to which 60 µL 3M sodium acetate (pH 8) & 800 µL isopropanol were added, and incubated overnight at 4 • C for precipitation. The samples were centrifuged at 14000 rpm for 10 min at 4 • C. The supernatant was discarded, after which the DNA pellet was washed twice with 1 mL 70% ethanol, and centrifuged again at 14000 rpm for 5 min at 4 • C. The DNA pellet was dried in an oven with the eppendorf tube lid open at 55 • C for 30 min. The pellet was re-suspended in 80 µL low salt TE buffer, and 5 µL of RNase (10 mg/mL) was added to remove any RNA contamination. Depending on the yield 1 µL on the 1% agarose gel was transferred, and stored at −20 • C. 50-100 ng of DNA was checked for intactness on 1% agarose. The ratio 260/280 = 1.6-1.8 and 260/230 = 1.8-2.2 obtained showed that the DNA was of food quality, and this was stored at −20 • C.

PCR Amplification and Sequencing
PCR analysis of the internal transcribed spacers (ITS1-5.8S-ITS2 cluster) regions of the ribosomal DNA gene cluster, and parts of the B-tubulin and calmodulin genes was done for each isolate using the primer pairs shown in Table 2. PCR amplifications were performed on 25 µL of a reaction mixture containing MgCl 2free reaction buffer, 50 mM MgCl 2 , 10 mM dNTP mix, 10 µM of each primer, 5 U/µL Taq DNA polymerase and 5 ng/µL of template DNA. The AccuPower R Taq PCR PreMix (Bioneer, Korea) was preferred as the Taq in the premix is not a high fidelity enzyme. PCR was carried out as follows: (1) one step at 94 • C for 3 min; (2) 30 cycles of the following three steps: 1 min at 94 • C, 1 min at 57 • C, 1 min at 72 • C, and (3) one final 10-min step at 72 • C. The PCR products were separated by 1.2% agarose gel electrophoresis in a Tris-base, acetic acid and EDTA buffer, and stained with ethidium bromide. PCR product purification was done using QIAquick PCR Purification Kit (Qiagen) and sequencing performed at Inqaba Ltd (Cape Town, Republic of South Africa). The obtained DNA sequences were assembled and edited using CLC Main Workbench 7.7.3 (https://www.qiagenbioinformatics. com/), and then compared with the sequences that are deposited in the GenBank (http://www.ncbi.nlm.nih.gov/). The sequences of each of the three genes were aligned separately with references in the NCBI database, and the SNPs/Indels were observed. Phylogenetic trees were constructed using maximum likelihood with 1000 bootstraps replicates using MEGA v 6. Sequences were deposited at GenBank under accession numbers indicated on Supplementary Tables 1-3.

Detection of Aflatoxin Genes
All A. flavus isolates were tested for the presence of structural genes aflD, aflQ, aflP, and regulatory genes aflR, using some of the primers described by Gallo et al. (2012), listed in Table 2, and optimized in the laboratory conditions. The aflD gene encodes an enzyme that catalyzes the conversion of the first stable aflatoxin biosynthesis intermediate, norsolorinic acid to averantin in both A. flavus and A. parasiticus (Bennett, 1981;Papa, 1982 (Cleveland, 1989) and in A. flavus (Yu et al., 1998).

Data Analysis
For statistical evaluation, the R software package (R core Team, 2014) version 2.15.3 was used. To cluster the isolates based on their toxin profile and to cluster the toxins, a hierarchical clustering, based on the simple matching distance, which is a measure for the dissimilarity between binary sample sets, was used. The data are represented as heatmaps, where to color codes indicate the presence (green) or absence of a certain metabolite or toxin. Furthermore, based on the expected and observed co-occurrence of the toxin, pairwise associations were classified as negative, positive or random. To verify whether the clustering of the isolates based on the toxin profile on CYA and YES was similar, correlation analysis was performed applying the Mantel test.

Multi-Analysis of Aspergillus Metabolites
To cluster the isolates based on their toxin profile and to cluster the toxins, a hierarchical clustering, based on the simple matching distance, which is a measure for the dissimilarity between binary sample sets, was used. The data are represented as heatmaps, where color codes indicate the presence (green) or absence of a certain metabolite or toxin. Furthermore, based on the expected and observed co-occurrence of the toxin, pairwise associations were classified as negative, positive or random. In Figure 1 the row dendrogram shows the clustering of the different isolates grown on either CYA or YESA medium. Isolates with a similar toxin profile were grouped together in the dendrogam. The dendrogram shows the associations between metabolites, metabolites with a similar occurrence pattern for the different isolates cluster together. To verify whether there was a correlation between the distances between isolates based on their toxin profile produced on CYA and the distance between isolates based on their toxin profile on YESA, a Mantel test with 1000 permutations was used (Figure 2). This test checks for significant correlations between distance matrices by a permutation procedure. The mantel correlation was 0.31, which was significant at α = 0.05. This means that the clustering of the isolates based on their toxin profile on CYA is correlated with the clustering of the species based on their toxin profile on YESA.

Identification and Analysis of Polymorphism in A. flavus Calmodulin Gene
During the analysis of calmodulin gene sequences, 563 characters were used, 17 of which were informative. Of the three genes used, calmodulin gene had the most SNPs and indels (Table 4). Among the isolates that amplified for this gene, 1, 30 and 37, all from Eastern Kenya, were identified as A. minisclerotigenis in agreement with the results in section Identification and analysis of polymorphism in A.flavus b-tubulin gene above. The remaining isolates were 100% identical to A. flavus sequences in NBCI with exception of isolate 90, which was 100% identical to A. parasiticus sequences in the genebank. Phylogenetic tree-based on calmodulin gene data indicates 2 major clades, namely A. parasiticus-forming one clade while A. flavus and A. minisclerotigenis forming a second clade. Several sub-clades for A. flavus are observed indicating a wide variation among A. flavus species (Figure 4).

Identification and Analysis of Polymorphism in A. flavus ITS Region
Four hundred and sixty characters were used in the analysis of ITS gene sequences. The ITS region was conserved in most regions for these isolates. Five SNPs were observed ( Table 5).
All isolates were identified as A. flavus with identical similarity (100%), excepted strains 1 and 28 which possessed a heterozygous loci, and matched (99%) to A. flavus. Isolate 90 was 100% identical to A. transmontanensis and A. parasiticus sequences in the NCBI genebank. Phylogenetic analysis of the sequence data of the ITS region shows isolates 1 and 28 formed a clade with a closer relationship with A. parasiticus (Figure 5). Variation exists among A. flavus isolates; those with SNPs (1, 28, 90) are closer in relationship.

DISCUSSION
Use of molecular markers in this study separates a population initially named as A. flavus using morphological and cultural (AFPA) characteristics into two species, namely A. flavus, The mycotoxin metabolic profile of the Aspergillus strains was highly variable suggesting possible high levels of genetic recombination among members of this species (Taylor et al., 1999). For example, isolate 90, which was identified as A. parasiticus produced CPA, which is not a known characteristic for this species. The same applied to some isolates of A. minisclerotigenes producing only AFB's, while some isolates of A. flavus produced both AFB's and AFG's. A. minisclerotigenes is known to produce both AFB's and AFG's while some isoltes of A. flavus have been classified as producing only AFBs (Pildain et al., 2008) while others produce both (Okoth et al., 2012). Similarly, Olarte et al. (2012) demonstrated high AF heritability, genetic variability and recombination in the AF gene cluster, analogous to spontaneous recombination during sexual reproduction in natural populations.
Knowledge of relative stability of characters among isolates of this species is important in the selection and use of non-toxigenic strains for the biological control of toxigenic ones. Further, A. flavus is an important taxa in food safety and as human pathogen. However, information of the relative phylogenetic distances at which these characters are acquired and lost, is yet to be established.
It was not always possible to distinguish aflatoxigenic species from non-aflatoxigenic species using PCR detection of AF biosynthesis genes, as some of the results were conflicting with the mycotoxin metabolite analysis. Both Gallo et al. (2012) and Levin (2012) reported similar observation and attributed this to inter-and intra-specific genetic mutations within the primers' targeted binding site. Fakruddin et al. (2015) tested 15 isolates of A. flavus for aflatoxigenic genes (aflD, aflM, aflP, aflQ, aflS, aflO, and aflR), and only one isolate possessed all of the 7 genes tested. The remaining isolates varied widely in terms of which genes among the 7 they possessed. The authors explained this to the inability of the primers used to amplify the genes. A few studies have reported a correlation between AF gene expression and AF production using RT-PCR and real time PCR targeting aflD, aflQ, aflO, aflP, aflR, and aflS genes (Sweeney et al., 2000;Jamali et al., 2013;Davari et al., 2015;Fakruddin et al., 2015;Baquiao et al., 2016). Mahmoud (2015) reported a correlation between PCR amplification of four aflatoxin biosynthetic pathway genes (aflD, aflM, aflP, and aflQ) with the AF production capability of A. flavus isolates from stored peanuts in Egypt. The same authors, however, noted that the expression of aflD and aflQ was a good marker for differentiating toxigenic from atoxigenic isolates. Accinelli et al. (2008) and Lourenco et al. (2007) did not find any correlation between profiles of aflatoxigenic A. flavus isolates and AFB1 concentrations in the soil. Jamali et al. (2013) recorded a correlation with AF produced. This study confirms the lack of correlation between the presence of the aflD, aflQ, aflP, aflR genes & AFB1 and AFB2 production, and promotes the usefulness of using both molecular and mycotoxin metabolite profiles in identification of the isolates. Though microbial and immunological methods are considered costly, time-consuming and may produce false-positives and/or false-negative results, a combination of these techniques together with molecular methods yields reliable results. Molecular biology methods combined with morphology and physiology reveal the diversity within the A. flavus species, and the need for further studies to enable identification of varieties. To date, the only clear morphological variation is the size of sclerotia. Small sclerotial A. flavus and A. minisclerotigenes are known to produce large amounts of aflatoxins, and these are most frequent in aflatoxinhot spot Eastern region of Kenya.

AUTHOR CONTRIBUTIONS
SO: Conceived and designed the experiments; SO, MD, AV, JD, SL, SD, MK, and JH: Generated data; JN, MD, and SO: Analyzed data; MD, SO, and SD: Wrote the manuscript.