Mutation of the LRG1 Rho-GAP gene is responsible for the hyper branching C-variant phenotype in the quorn mycoprotein fungus Fusarium venenatum A3/5

Abstract

Background

Quorn mycoprotein, a protein-rich meat alternative, is produced through large-scale fermentation of the fungus Fusarium venenatum. However, a major challenge during F. venenatum fermentation is the consistent appearance of mutants called colonial variants (C-variants). These C-variants have a highly branched morphology, which ultimately lead to a less desirable final product and early termination of the fermentation process. This study aimed to identify the genetic mutations responsible for C-variant morphology.

Results

We first isolated both C-variant and wild-type strains from commercial fermentation samples and characterised radial growth rates on solid media. Whole genome sequencing facilitated the identification of mutations in a gene called jg4843 in 11 out of 12 C-variant isolates, which were not observed in the wild-type isolates. The jg4843 gene was identified as the ortholog of LRG1, a Rho-GTPase activating protein that regulates the Rho1 signalling pathway affecting fungal growth. Notably, the mutations in jg4843 were primarily located in the RhoGAP domain responsible for LRG1 activity. To confirm the role of these mutations, we used CRISPR/Cas9-mediated homology-directed recombination to introduce the C-variant mutations into the wild-type isolate, which successfully recapitulated the characteristic C-variant morphology.

Conclusions

This study identified mutations in the LRG1 ortholog jg4843 as the genetic cause of C-variant morphology in commercial fermentation F. venenatum isolates. Understanding this genetic basis paves the way for developing strategies to prevent C-variants arising, potentially leading to more efficient and sustainable production of Quorn mycoprotein.

Background

The successful large-scale cultivation of Quorn mycoprotein has revolutionised the alternative protein market, offering an environmentally friendly and relatively more sustainable alternative to traditional meat production. The production of mycoprotein, pioneered by Quorn, has now expanded, with numerous emerging companies entering the market and further exploring mycoprotein as a sustainable protein source. Mycoprotein is derived from the mycelial biomass of the filamentous fungus Fusarium venenatum, cultivated through a controlled continuous flow fermentation process [1, 2]. The combination of the sparsely branched F. venenatum A3/5 strain and the production process results in the formation of biomass with a desirable meat-like texture and composition [1, 2]. However, the commercial production of mycoprotein is hampered by the consistent appearance of highly branched mutant strains during each fermentation campaign [1]. Originally named colonial variants (C-variants) due to their compact colony morphology, these C-variants are characterised by excessive and irregular hyphal branching patterns, resulting in dense mycelial masses that negatively affect the organoleptic quality of the final product [3, 4]. If left to continue, the C-variants would eventually out compete the wild-type A3/5 population. This limits the efficiency and capacity of production, and each fermentation campaign can only run for a limited time (usually 4–6 weeks) before having to be stopped [3]. By delaying or preventing the appearance of these C-variants through strain development, each fermentation campaign could be prolonged, leading to a more efficient and sustainable Quorn production process [5]. However, despite previous preliminary work identifying several candidate genes, including the central gene of this study [6], the causal mutation(s) responsible for the appearance of these highly branched strains remain unknown, therefore, understanding the genetic basis underlying the occurrence of C-variants is crucial for optimising mycoprotein production.

Commercial mycoprotein fermentations are run with all nutrients (including glucose) in excess [1, 2]. It has been understood for many years that continuous culture under constant conditions confers a selection pressure for higher growth rates [7, 8]. The mutations occurring in C-variants are likely to result in a selective advantage over the parental strain in this environment due to a shorter doubling time. Moreover, Simpson et al. [9] demonstrated that C-variants isolated from commercial mycoprotein fermentation samples displayed growth rate advantages over the wild-type isolate in both conditions of nutrient excess and various nutrient limitations. In contrast, C-variants isolated from lab-based glucose-limited, continuous flow cultures had growth rate advantages over the wild-type isolate in glucose-limited conditions, but not when glucose was in excess [10].

Previous studies in Neurospora crassa have identified a large number of genes associated with alterations in hyphal morphology [11]. It is equally likely in F. venenatum that numerous genes determine morphology and that mutations in any one gene may have a profound effect on the fitness of the organism, especially in suspension culture. While the causal mutation(s) that lead to C-variant formation in F. venenatum are not currently known, classical genetics work, particularly through heterokaryon formation of different C-variant isolates, has been used to determine the number of loci involved [4, 9]. These analyses have identified multiple different ‘complementation groups’, dependent on the origin of the C-variant isolates. Strains isolated from glucose-limited continuous cultures were assigned to three complementation groups [4], whereas C-variants isolated from commercial fermentation samples were assigned to a single complementation group [9]. Irrespective of the isolate source, all strains were found to be caused by recessive mutations [4, 9]. Furthermore, preliminary work identified several C-variant candidate genes in isolates from commercial fermentations [6]. Although C-variants are grouped by a common phenotype, these genetic complementation analyses, together with the differences observed in growth rate advantages in different nutrient conditions, suggests that there are separable groups of strains indicating that either different loci or different gene variants are responsible for the C-variant phenotype dependent on isolate origin.

Several interventions to prevent or delay the appearance of C-variants during mycoprotein production have been previously identified. Many of these strategies involved altering the selection pressure. For example, running the fermenter at low dilution rates [12], reducing the pH from 5.8 to 4.5 [13], changing the nitrogen source [14], periodically altering the limiting nutrient source [3, 13], or isolating more stable strains of F. venenatum [3, 15]. Hyphal branching in F. venenatum is also altered by a range of chemical interventions, including choline [16], cAMP and cGMP [17]. However, since all these strategies negatively affect other aspects of F. venenatum growth, none of these interventions have been found to be commercially viable and, as such, are not currently implemented in mycoprotein production. To mitigate the impacts of C-variant appearance it is, therefore, crucial to better understand the molecular mechanisms underlying the appearance of these highly branched mutants.

To date, only four genome assemblies are publicly available for F. venenatum, representing approximately 13,000 coding sequences within a 38 Mb genome, the most contiguous of these consists of 6 contigs [18]. However, genomic data is not currently available for C-variant isolates. In this study, we have investigated the genetic cause of C-variant formation in F. venenatum through an integrative omics approach, using whole genome sequencing, transcriptomic analysis, and genome editing. We aimed to unravel the genetic mutations driving the hyper branching phenotype by performing genetic variant analysis on 12 C-variant and 12 post-fermentation wild-type strains isolated from distinct fermentation campaigns, separating mutations that also may be selected for in fermentation from those specifically responsible for the C-variant phenotype. We identified genetic variants present only in C-variant genomes and, by incorporating transcriptomic data, we identified a single gene that was exclusively mutated and/or downregulated in all C-variant isolates. Finally, CRISPR/Cas9-mediated homology-directed recombination (CRISPR-HDR) was employed to validate the causal effects of C-variant SNPs in a wild-type background. These findings highlight that a single gene is likely to be responsible for the appearance of C-variants in commercial mycoprotein production.

Methods

Fungal strains and culture conditions

F. venenatum (ATCC 20334, A3/5) was acquired from Marlow Foods Ltd, and hereafter, will be referred to as wild-type (WT). It and its derived strains were conserved as a mixture of mycelial fragments in 25% glycerol at − 80 °C, and routinely grown in/on potato dextrose agar/broth (PDA/PDB) at 28 °C. Escherichia coli DH5α was grown on/in LB agar/broth (with or without appropriate antibiotic selection) at 37 °C for routine gene cloning and vector construction.

Strain isolation from mixed fermentation samples

C-variant (Cv) and post-fermentation wild-type (pfWT) strains were isolated from mixed cultures acquired from independent fermentation campaigns conducted at Marlow Foods Ltd. (Billingham, UK). Mixed fermentation samples were plated onto PDA and incubated for 7 days at 28 ^oC. After this initial growth period, three 10 mm fungal plugs were used to inoculate 200 ml carboxymethyl cellulose sodium salt (CMC) media (1.5% carboxymethyl cellulose sodium salt, 0.1% NH₄NO₃, 0.1% KH₂PO₄, 0.05% MgSO₄.7H₂O, 0.1% yeast extract), which were grown on a rotary shaker at 25 ^oC for 7 days. The culture was filtered through 2 layers of sterile Miracloth (Millipore), and the filtrate was centrifuged at 4000 ×g for 10 min. After discarding the supernatant, conidia were resuspended in sterile water and counted. Approximately 100 spores were plated onto PDA and grown for 12 h at 28 ^oC. Following incubation, several single spores were transferred to fresh PDA and grown for 7 days at 28 ^oC. Cv colonies exhibit a more compact structure and have slower radial growth rates as compared to the WT isolate [4]. Therefore, these morphological characteristics were used to classify colonies as either Cv or pfWT.

Colony radial growth assay

F. venenatum isolates were grown on PDA for 10 days at 28 °C. Following this 10-day growth, 5 mm plugs were transferred to fresh PDA. Following a 2-day growth period, colony diameters were measured in two 90^o planes using digital callipers at the same time daily for the next five days. Colony radial growth rates (K_r) were calculated as the average growth rate per day over the five-day period. Each radial growth assay was performed three times. A one-way ANOVA was performed to investigate the relationship between radial growth rate and isolate in R v4.2.3. Additionally, a Dunnett’s post hoc test was conducted using the DescTools library to compare each isolate to the WT. Data were plotted using the python library seaborn v0.12.2 [19]. Fungal colonies were imaged using a Canon EOS 4000D DSLR camera. Hyphal morphology was assessed using a Leitz Dialux 20EB microscope with a Leica EC3 camera.

DNA extraction and sequencing

Three 3 mm plugs were used to inoculate 50 ml PDB cultures, which were grown for 7 days on a rotary shaker at 28 ^oC. Mycelial biomass was harvested on filter paper under vacuum, followed by freeze-drying. DNA extraction was performed using a Macherey-Nagel NucleoSpin Plant II kit (Thermo Fisher Scientific, 11912262), following manufacturer’s recommended protocol for fungal DNA extraction. DNA quantification was performed using a Qubit (Thermo Fisher Scientific), and DNA purity was measured using a Nanodrop 1000 spectrophotometer (Thermo Fisher Scientific). For long-read sequencing, approximately 1 µg of high molecular weight genomic DNA from the WT was used with the SQK-LSK108 Ligation Sequencing Kit (Oxford Nanopore Technologies), following the manufacturer’s protocol. The long-read libraries were sequenced on R9.4.1 Spot-On Flow cells (FLO-MIN106) using the GridION X5 platform (Oxford Nanopore Technologies) set to high accuracy base calling. For Illumina sequencing, the WT and first batch of isolates (Cv 1–6, pfWT 1–6) were sequenced in-house on an Illumina MiSeq. Libraries were prepared as follows: DNA was sheared using the Covaris M220 with microTUBE-50 (Covaris, 520166) and size selected using the Blue Pippin (Sage Science). Illumina libraries were constructed with a PCR-free method using NEBNext End Repair (NEB, E6050S), NEBNext dA-tailing (NEB, E6053S) and Blunt T/A ligase (NEB, M0367S) modules. Libraries were multiplexed and sequenced using Illumina Miseq v3 2 × 300 bp PE (Illumina, MS-102-3003). The second batch of isolates (Cv 7–12, pfWT 7–12) were sent to Novogene (Cambridge, UK) for library preparation and sequencing. DNA sequencing data were deposited at the NCBI under the Bioproject ID PRJNA1151243.

Genome assembly

Long read sequencing data were quality controlled using NanoPlot v1.30.1 [20] (Figure S1) and adapters were trimmed using Porechop v0.2.4 [21], with default parameters. Long reads were assembled using NECAT v0.0.1_update20200803 [22] using a genome size of 38 Mb and all other parameters as default. Following assembly, error correction with long read data was performed using one iteration of Racon v1.4.20 [23], followed by one iteration of Medaka v1.5.0 [24] using the r941_min_high_g360 model. Illumina paired-end reads were quality controlled using FastQC v0.11.9 [25] (Figure S2), and adapters and low-quality regions were trimmed using Trimmomatic v0.39 [26]. Error correction with Illumina data was performed by aligning the short reads to the corrected long-read assembly using Bowtie2 v2.4.4 [27] to inform three iterations of polishing using Pilon v1.24 [28]. Assembly statistics were generated using a custom python script, and single copy ortholog benchmarking was performed using BUSCO v5.2.2 [29] with the hypocreales_odb10 database. The final genome assembly is deposited at the NCBI (MXQG00000000).

Variant calling pipeline

DNA sequencing reads were quality controlled using FastQC v0.11.9 [25], and adapters and low-quality regions were trimmed using Trimmomatic v0.39 [26] using a sliding window of 4 and minimum Phred score of 20. The first 10 nucleotides were trimmed and reads less than 100 nucleotides and unpaired reads were discarded. The WT genome was indexed, and reads were aligned using BWA mem v0.7.17 [30] using default settings. Alignments were then converted to binary files, sorted, and indexed using SAMtools v1.17 [31]. VCF files were generated using BCFtools v1.17 [31] mpileup command with the -m 5 option, and variants were called using BCFtools call with the --ploidy 1 and -c options. Variants were filtered using the BCFtools view command, excluding those with a quality score less than 40, an average per sample coverage less than 10, and variants present in the WT reads. Variants were then annotated using SnpEff v5.0e [32]. CMPlot v4.5.1 [33] was used to visualise the distribution and density of variants across the genome.

RNA extraction and sequencing

Three 3 mm plugs were used to inoculate 50 ml PDB cultures, with four biological replicates per isolate. The cultures were grown at 28 ^oC on a rotary shaker and were harvested at 40 h post inoculation during the exponential growth phase, as previously determined [6]. Mycelial biomass was harvested and washed on filter paper under vacuum, blot dried and flash frozen in liquid nitrogen prior to extraction. All samples were ground using a sterile pestle and mortar, and total RNA was extracted using a modified Macherey-Nagel NucleoSpin RNA kit (Thermo Fisher Scientific), whereby prior to use of the NucleoSpin columns, samples were ground in TRIzol reagent and ethanol precipitated. Due to issues with sample degradation, at least three biological replicates per isolate were submitted for library preparation and paired-end RNA sequencing, which was performed by Novogene (Cambridge, UK) on an Illumina HiSeq 4000 platform. RNA sequencing data were deposited at the NCBI under the Bioproject ID PRJNA1151243.

RNA-seq analysis pipeline

RNA sequencing reads were quality controlled using FastQC v0.11.9 [25] (Figure S3). Adapters and low-quality regions were trimmed using Trimmomatic v0.39 [26] with a sliding window of 4 and minimum Phred score of 20. The first 10 nucleotides were trimmed and reads less than 100 nucleotides and unpaired reads were discarded. The WT transcriptome was indexed, and reads were quantified using Salmon v1.6.0 [34]. Differential expression analysis was performed using the R package DESeq2 v1.40.2 [35], and using contrasts between the WT and each Cv and pfWT isolate. By default, DESeq2 uses the Benjamini-Hochburg correction to determine the false discovery rate (FDR) [35]. An FDR below 0.05 was used to identify differentially expressed genes (DEGs). Variance stabilizing transformation (VST) was used to normalise the raw read counts and a principal component analysis (PCA) was performed. A heatmap of DEG log2 transformed fold change (log2FC) data was produced using the python library seaborn v0.12.2 [19]. Shared DEGs across all isolates were plotted using the UpsetR v1.4.0 package in R [36]. Two lists representing all DEGs in the pfWT isolates and the Cv isolates were generated using custom bash commands, and were used to determine gene ontology (GO) enrichment of biological processes using the topGO v2.50.0 package in R [37].

Genome annotation

A repeat library specific to the WT was constructed using RepeatModeler v2.0.4 [38], which was then used with RepeatMasker v4.1.5 [39] to mask repeat regions throughout the genome. The softmasked assembly was indexed and RNA-seq reads were aligned using HISAT2 v2.2.1 [40] with default settings. Gene prediction was performed using BRAKER2 v2.1.6 [41], using the -etpmode setting with the RNA-seq data and the fungi OrthoDB [42] as evidence. Completeness of the annotation was assessed using BUSCO v5.2.2 [29], using default parameters and the gene families set defined for the hypocreales_odb10 database.

Functional annotation

Functional annotation of predicted proteins was performed using InterProScan v5.65 [43]. OrthoFinder v2.5.4 [44] was used to identify orthologous genes in WT assembly with proteomes from model fungal species, Aspergillus nidulans FGSC A4 (GCF_000011425.1), Candida albicans SC5314 (GCF_000182965.3), N. crassa OR74A (GCF_000182925.2), and Saccharomyces cerevisiae S288C (GCF_000146045.2). Protein domain prediction was performed using the ScanProsite tool [45], and lollipop plots were plotted using the R package trackViewer [46].

CRISPR/Cas9-mediated homology-directed recombination (CRISPR-HDR)

SgRNA cloning

Cloning of the sgRNA into the pFC332 CRISPR/Cas9 expression vector was performed as previously described [47]. Briefly, sgRNA PAM sites were identified using Geneious Prime v2023.0.4 (https://www.geneious.com) and the 5 S rRNA promoter-sgRNA scaffold (Table S1) was synthesised as a ‘gBlock Gene Fragment’ (IDT). This scaffold was cloned into the USER site of pFC332, which also contains an Aspergillus niger codon-optimized Cas9, hygromycin resistance, and AMA1 components. Two 190 bp homology-directed recombination (HDR) donor oligonucleotides were synthesised as ‘Alt-R HDR Donor Oligos’ (IDT), one oligo contained the LRG1 SNP identified in Cv6 and the other, the LRG1 SNP identified in Cv11. Both oligos contained a synonymous mutation in the sgRNA PAM site to prevent recognition and cleavage of the integrated DNA. The sequences of all primers, sgRNA scaffold, and HDR donor oligonucleotides used in this study are listed in Table S1.

Protoplast isolation and transformation

Protoplast isolation and transformation of the WT was performed as previously described [47]. Briefly, the protoplasts were prepared using 25 mg ml^− 1 Driselase (Sigma, D9515), 50 µg ml^− 1 Chitinase (Sigma, C6137) and 1 mg ml^− 1 Lyticase (Sigma, L4025) at 28 °C, 80 rpm for 2 h. Following digestion, the protoplasts were washed twice in STC (1.2 M Sorbitol, 10 mM Tris-HCl pH 8.0, 50 mM CaCl₂) and resuspended to 1 × 10⁸ protoplasts ml^− 1 in STC with 8% PEG 4000. For transformation, 100 µl of protoplasts was mixed with 5 µg plasmid DNA (with 5 µg HDR oligo, where required) and incubated at room temperature for 20 min. Then, 1 ml of STC with 40% PEG 4000 was added and incubated at room temperature for a further 20 min. Following these incubations, the transformation mix was added to 5 ml TB3 medium (0.3% yeast extract, 0.3% acid-hydrolysed casein, 1% D-glucose) supplemented with 100 µg ml^− 1 hygromycin B and incubated at 21 ^oC overnight. Each culture was combined with 10 ml Top Agar (TB3 with 18% Sorbitol, 1.5% agar), cooled to 50 ^oC and supplemented with 100 µg ml^− 1 hygromycin B, and then poured into 90 mm petri dishes. Following an 8 hour incubation at room temperature, a further 10 ml Top Agar supplemented with 100 µg ml^− 1 hygromycin B was added to each plate, and plates were incubated at room temperature. Colonies from transformation plates were transferred to PDA supplemented with 75 µg ml^− 1 hygromycin B after 4 days, and then incubated for 2 weeks at 28 °C.

CRISPR-HDR mutant confirmation

Genomic DNA for PCR was extracted from transformation colonies by macerating a small piece of mycelium in 20 µl of 0.5 M NaOH and incubating at room temperature for 20 min. This was then diluted with 400 µl 0.1 M Tris-HCl pH 8.0, and 5 µl was used in a 50 µl PCR reaction using Phusion HF DNA polymerase (NEB) and the primers LRG1_HDR_Chk_F and LRG1_HDR_Chk_R (Table S1), which bound 139 bp upstream and 143 bp downstream of the sgRNA target site, respectively. Amplicons were purified using a Monarch PCR & DNA Cleanup kit (NEB) and sent to Genewiz (Azenta Life Sciences) for Sanger sequencing using primer LRG1_HDR_Chk_F.

Results

Isolation and characterisation of fermenter strains

To understand the genetic basis of the hyper branching C-variant morphology, six C-variant strains (Cv 1–6) were isolated from samples taken at the end of independent commercial fermentation campaigns (Fig. 1A). Since the DNA of each nucleus within every fermentation campaign will accumulate mutations over time, to facilitate the identification of mutations occurring only in C-variants, six post fermentation wild type-like strains (pfWT 1–6) were also isolated. Concurrent with previous work, all isolates were selected based on their colony appearance whilst growing on an agar medium, the Cv isolates developing densely formed colonies compared to the ‘fluffier’ appearance of the WT and pfWT colonies (Fig. 1A). Following visual assessment, Cv isolates were found to exhibit increased levels of hyphal branching, whereas the hyphal morphology of the pfWT isolates was similar to that of the WT (Fig. 1B). The growth of these twelve isolates was further characterised by determining colony radial growth rates (K_r) over a seven-day period. Most of the pfWT isolates showed similar radial growth rates to the WT, with K_r of approximately 11 mm day^− 1 (Fig. 1C). However, pfWT1 had a significantly increased radial growth rate (p = 0.0179), with an average K_r of 12.7 mm day^− 1 (Fig. 1C). All of the Cv isolates had significantly reduced radial growth rates (p < 0.001), with K_r ranging between 3 and 4 mm day^− 1 (Fig. 1C), characteristic of highly branched F. venenatum [13].

Fig. 1

Isolation and characterisation of F. venenatum isolates from commercial fermentation samples. (A) Representative colony morphology of WT, C-variant (Cv) and post-fermentation WT (pfWT) isolates after growth on PDA for seven days at 28 °C. (B) Representative hyphal morphology of all batch 1 isolates after growth in PDB for 40 h at 28 °C. Scale bar represents 100 μm. (C) Boxplot of radial growth rate for isolates averaged over a five-day period (n = 3). *** p < 0.001, * p < 0.05

To facilitate genomic investigations, a chromosome-level genome assembly of the production F. venenatum A3/5 strain was produced. Long read Oxford Nanopore data (37x) was used to generate an initial assembly, which was subsequently polished using high quality Illumina data (66x). The resulting gapless assembly consisted of four chromosomes, with lengths ranging between 8 Mb and 12 Mb, and a 92.5 kb mitochondrial genome. This chromosome-level assembly had a similar overall size to the other publicly available F. venenatum assemblies and was comparable to the complete F. graminearum genome (Table S2). Furthermore, compared to the next best F. venenatum A3/5 assembly (GCA_900007375.1, Table S2), this new assembly had no unplaced contigs, contained no unknown sequences and slightly improved the single copy ortholog completeness based on the BUSCO Hypocreales database (Table S2). Gene prediction resulted in 14,211 genes, with an average gene coding length of 1,411 bp and an average of 2.9 exons per gene, comparable to other available annotations (Table S2).

Genetic and transcriptomic characterisation of commercial isolates

Illumina whole genome DNA sequencing data was also generated for the six Cv and six pfWT isolates. Genetic variants were called against the NIAB v2 WT genome assembly, and a total of 24 SNPs at 20 positions and eight INDELs were identified across all twelve isolates. Of the eight identified INDELs, seven occurred in repetitive or intergenic regions of the genome with low coverage so were not investigated further. However, one INDEL caused a disruptive inframe deletion in the gene jg4630 in pfWT6 (Table S3A). Of the 24 SNPs, six were identified across the pfWT isolates and 18 were identified across the Cv isolates (Table 1, Table S3A). No SNPs identified in the pfWT isolates were found in the Cv isolates, nor vice versa. Of the six SNPs identified in the pfWT isolates, three were upstream of genes, two were synonymous mutations, and one was a missense mutation; whilst in the Cv isolates, 11 SNPs were missense mutations, five occurred in upstream regions, one was a synonymous mutation, and one SNP occurred in an intronic region (Table 1, Table S3A). To visualise the distribution of these mutations across the genome, the variants across all isolates were aggregated in 100 kb windows across the genome and the data plotted for each chromosome (Fig. 2A). This identified two hotspots that were frequently mutated in the Cv isolates, one on chromosome 1 in the region 4,100,000–4,199,999 and the other on chromosome 4 in the region 7,200,000–7,299,999 (Fig. 2A). Interestingly, Cv6 only contained a single SNP, which was located in the hotspot on chromosome 1 (Table S3A).

Table 1 Quantification of types of SNPs identified in Cv and pfWT sisolates

Fig. 2

Heatmap of the distribution and density of all identified SNPs throughout the F. venenatum genome. The number of SNPs per 100 kb region across each chromosome were plotted for all of the isolates in batch 1 (A) and batch 2 (B). The blue triangles indicate SNPs that were only identified in pfWT isolates, all other SNPs were found only in Cv isolates

To determine the effect of these mutations on gene expression, we performed RNA-seq and differential gene expression analysis on all twelve fermenter isolates and the WT, each grown to mid-exponential phase, as determined by preliminary experiments [6]. The total number of differentially expressed genes (DEGs) varied across all isolates and did not correlate with morphological classification (Table S4). However, there was minimal overlap in DEGs between pfWT and Cv isolates (Fig. 3A). Most DEGs were either unique to or shared within each morphological group, with only a small subset common to both (Fig. 3A). Gene Ontology (GO) analysis indicated that DEGs in pfWT isolates were enriched for RNA binding and processing, whereas Cv isolates showed enrichment in GTP binding and GTPase activity (Fig. 3B). An enrichment in translation initiation was observed in both groups (Fig. 3B).

Fig. 3

(A) Upset plot showing the numbers of DEGs shared across all batch 1 isolates (B) Dotplot highlighting significantly enriched biological processes following gene ontology (GO) analysis. (C) Heatmaps of log2FC and FDR of mutated genes. Heatmap of log2 fold change (log2FC) of mutated genes in batch 1 isolates. Letters denote the presence of a SNP and indicate the type of mutation (M = missense, U = upstream, S = synonomous, I = intron). Black boxes highlight significant differentially expressed genes (FDR < 0.05)

When correlated with the genetic mutations, none of the five mutated genes in the pfWT isolates had significant changes in expression compared to the WT (Fig. 3C, Table S5). However, correlations between SNPs and expression were observed in five out of the ten genes mutated in the Cv isolates (Fig. 3C, Table S5). Mutations upstream of genes jg4215 (Cv2), jg6073 (Cv3) and jg6271 (Cv1), an intronic mutation in jg6573 (Cv2), and missense mutations in jg10978 (Cv4 and Cv5) had no significant effects on the expression of these genes (Fig. 3C, Table S5). Conversely, an upstream mutation in Cv2 and a missense mutation in Cv1 corresponded to a significant decrease in the expression of jg4843, a putative Rho-type GTPase-activating protein (Table S7), which was also downregulated in Cv3 (Fig. 3C, Table S5). The expression of jg64, a cellulose binding domain-containing protein (Table S7), was significantly upregulated in all Cv isolates except for Cv3, which had a mutation upstream of this gene that corresponded with WT levels of expression (Fig. 3C, Table S5). The missense mutation in the putative serine/threonine protein kinase gene jg12728 in Cv3 (Table S7) correlated with a significant downregulation of expression, yet this gene was also significantly downregulated in Cv1, which did not have any mutations in jg12728 (Fig. 3C, Table S5). Missense mutations in jg13054, a gal11 coactivator domain-containing protein (Table S7), in Cv1 and Cv3 were associated with a decrease in expression compared to the WT, yet missense mutations in this gene in Cv2 and Cv5 were not (Fig. 3C, Table S5). Finally, the SNP upstream of jg13105, a putative oxalate decarboxylase/oxidase (Table S7), in Cv2 correlated with a significant increase in expression in this gene, while Cv4 and Cv6 had significantly downregulated levels of expression compared to the WT (Fig. 3C, Table S5). These initial genomic analyses identified two mutational hotspots unique to the Cv isolates, one of which contained mutations in a gal11 coactivator domain-containing protein (jg13054) and the other, a putative Rho-type GTPase-activating protein (jg4843), found in four and five Cv isolates, respectively. Notably, Cv6 was the only Cv isolate containing a single SNP, which was located in jg4843. Transcriptomic analysis revealed significant downregulation of both genes in certain Cv isolates. Additionally, GO analysis indicated an enrichment in GTP binding and GTPase activity. Together, these findings suggest a potential role for jg4843 and/or jg13054 in C-variant morphology.

Confirmation of C-variant mutations in a second batch of commercial isolates

To gain further genetic evidence for the mutations responsible for C-variant morphology, a second batch of isolates were recovered from additional independent commercial fermentation campaigns (Fig. 4A). Six pfWT isolates (pfWT 7–12) and six Cv isolates (Cv 7–12) were identified based on the observed colony morphology on agar (Fig. 4A) and hyphal morphology was visually confirmed (Fig. 4B). The radial growth rates for these isolates were also characterised (Fig. 4C). All the pfWT isolates had similar radial growth rates to the WT, with K_r of between 9 and 11 mm day^− 1 (Fig. 4C). The Cv isolates in the second batch were more variable than the first, with K_r ranging between 2 and 6 mm day^− 1 but all were found to have significantly slower radial growth rates than the WT (p < 0.001, Fig. 4C).

Fig. 4

Isolation and characterisation of F. venenatum isolates from additional commercial fermentation samples. (A) Representative colony morphology of WT, C-variant (Cv) and post-fermentation WT (pfWT) isolates after growth on PDA for seven days at 28 °C. (B) Representative hyphal morphology of all batch 2 isolates after growth in PDB for 40 h at 28 °C. Scale bar represents 100 μm. (C) Boxplot of radial growth rate for additional isolates averaged over a five-day period (n = 3). *** p < 0.001

Whole genome sequencing data was generated for this second batch of isolates and variant analysis was performed as above. A total of 29 SNPs at 28 positions and 13 INDELs were identified across all twelve isolates. Eleven of the INDELs were identified in low coverage, repetitive regions so were not investigated further. The INDEL identified in the batch 1 isolate pfWT6, causing a disruptive inframe deletion in jg4630, was also found in the batch 2 isolates pfWT8 and pfWT9 (Table S3B). Another INDEL was identified in Cv10 that resulted in a frameshift mutation in jg5568 (Table S3B). Of the 29 SNPs, eight were identified across all the pfWT isolates and 21 were identified across the Cv isolates (Table 1, Table S3B). None of the SNPs were shared between the pfWT and Cv groups. In the pfWT isolates, five of the eight SNPs occurred upstream of genes, two were synonymous mutations, and one was a missense mutation (Table 1, Table S3B). Of the 29 SNPs identified in the Cv isolates, 14 were missense mutations, five were found in upstream regions, one was a synonymous mutation, and one SNP resulted in a premature stop codon (Table 1, Table S3B). Again, all of the variants across all batch 2 isolates were aggregated and plotted in 100 kb windows across the genome (Fig. 2B). The previously identified mutational hotspot at the region 4,100,000–4,199,999 on chromosome 1 was also enriched for this second batch of Cv isolates, whilst the hotspot at the region 7,200,000–7,299,999 on chromosome 4 was not (Fig. 2B). Following further investigation, only one gene in this region on chromosome 1 (jg4843) was mutated, encoding a putative Rho-type GTPase-activating protein (Table S7), with a single SNP in five out of the six batch 1 Cv isolates and all six of the batch 2 isolates. Most of these mutations were concentrated in the 3’ end of exon 3 of the gene, with a single SNP 22 bp upstream of the start codon, which was identified in two Cv isolates (Cv2 and Cv8), one from each batch (Fig. 5A). Furthermore, the presence of this upstream SNP in Cv2 correlated with a significant downregulation of this gene (Fig. 3C). Interestingly, jg4843 was also significantly downregulated in Cv3, the only Cv isolate without mutations in this gene (Fig. 3C).

Fig. 5

Distribution of SNPs throughout the LRG1 gene (A) and protein (B). Blue circles indicate SNPs identified in batch 1 isolates and red circles indicate SNPs identified in batch 2 isolates

An orthology search against model fungal species identified jg4843 as the F. venenatum ortholog of a LIM and RhoGAP domain-containing protein LRG1 (Table S7). LRG1 is a Rho1 GTPase activating protein (GAP) and acts as a negative regulator for the Rho1 GTPase protein, functioning through the RhoGAP domain [48]. Following a protein domain search, three LIM domains (corresponding to amino acids 95–157, 159–220, and 465–529) and a RhoGAP domain (amino acids 795–1000) were identified in jg4843 (Fig. 5B), which we have called FvLRG1. Most of the missense mutations identified in FvLRG1 occurred in the RhoGAP domain (Fig. 5B). Eight isolates had a missense mutation within the RhoGAP domain, representing seven different amino acid substitutions throughout the 205 amino acid domain, with the S952F mutation occurring in Cv6 and Cv7 (Fig. 5B). One missense mutation, a leucine to proline substitution (L730P) occurred upstream of the RhoGAP domain in Cv12 (Fig. 5B). Notably, all batch 2 Cv isolates, as well as five out of six batch 1 Cv isolates, carried mutations in FvLRG1, most of which were concentrated in the 3’ end of exon 3, in the key functional RhoGAP domain. Together, these data strongly suggest that disruption of FvLRG1 is the causal driver of the hyper branching C-variant phenotype in commercial fermentation isolates.

Introduction of a single SNP confirms effects of LRG1 mutation

To confirm the role of FvLRG1 in C-variant morphology, we utilised a CRISPR/Cas9-mediated homology-directed recombination (CRISPR-HDR) approach to separately introduce two different SNPs identified in the Cv isolates into the FvLRG1 gene in a WT background (Fig. 6A). This allowed us to assess the role of these mutations independent of the other mutations occurring in the Cv isolates. The missense mutation S952F shared by two different Cv isolates, Cv6 and Cv7, and the nearby N958T identified in Cv11 were chosen for further analysis. Both mutations occurred towards the C-terminal region of the RhoGAP domain (Fig. 5B), and their proximity to each other allowed the use of a single guide RNA.

Fig. 6

Morphology and radial growth rates of CRISPR-HDR colonies. (A) Schematic diagram depicting the CRISPR-HDR process. The pFC332 vector, containing the Cas9 (blue), sgRNA (purple), hygromycin resistance (pink) and AMA1 (brown) components, is transformed into the fungus alongside the single strand HDR oligo. The Cas9-sgRNA complex cuts the LRG1 gene (green) at the target sequence, and the DNA is repaired using the HDR oligo to incorporate the SNP (red) into the genome. (B) Representative colony and hyphal morphology of Empty Vector, LRG1 sgRNA (gRNA), LRG1 sgRNA with Cv6 HDR oligo (gRNA + Cv6 HDR), and LRG1 sgRNA with Cv11 HDR oligo (gRNA + Cv11 HDR) transformants. Colonies were imaged after growth on PDA for seven days at 28 °C and hyphae were imaged after growth in PDB for 40 h at 28 °C. Number of colonies displaying representative morphology indicated below images. Scale bar represents 100 μm. (C) Boxplot of radial growth rate for WT, Cv6, Cv11 and selected transformant colonies averaged over a five-day period (n = 3). *** p < 0.001

Following transformation, four empty vector and 10 sgRNA, sgRNA with Cv6 (S952F) HDR oligo, and sgRNA with Cv11 (N958T) HDR oligo colonies were picked and transferred to PDA for morphological characterisation. When the empty vector was transformed into the WT background, all four of the isolated transformant colonies grew like the non-transformed WT, with similar K_r of approximately 9 mm day^− 1 (Fig. 6B and C). Additionally, with the sgRNA construct alone, there was no change in colony morphology, hyphal branching or radial growth rate compared to the WT (Fig. 6B and C). However, when the S952F and N958T were introduced into the FvLRG1 gene in the WT background, a dramatic change in colony morphology, hyphal branching and radial growth rate was observed (Fig. 6B and C). In eight out of the ten selected colonies for the Cv6 (S952F) HDR oligo and nine out of the ten Cv11 (N958T) HDR oligo colonies, the colony and hyphal morphology mirrored that of the original Cv isolates, with K_r between 2 and 3 mm day^− 1 (Fig. 6B and C). The presence of the expected mutations in the transformant colonies used in the radial growth rate assays was confirmed using Sanger sequencing (Figure S4). These data confirmed that mutations in FvLRG1 are sufficient to cause the characteristic hyper branching C-variant phenotype.

Discussion

Quorn mycoprotein, made from the biomass of the filamentous fungus F. venenatum, is a key source of sustainable protein [5]. Whilst the commercial fermentations produce several tons of mycoprotein per hour, these campaigns can only run continuously for a few weeks before highly branched C-variants appear in the population and eventually out compete the parent strain [1, 2]. If continuous cultivation could be prolonged by delaying the appearance of the C-variant phenotype, the production process could become more efficient and sustainable. In this study, we have taken a step closer to this goal by identifying the previously unknown, causal mutations responsible for C-variant morphology in commercial mycoprotein fermentations.

C-variant colonies were isolated from multiple independent commercial fermentation campaign samples, which all displayed characteristic colony and hyphal growth (Figs. 1 and 4) as has been previously reported [4, 9, 10, 13]. A total of 12 post-fermentation isolates with WT-like morphology were also isolated to control for natural mutation rates that occur under continuous culture conditions [7, 8]. Although Wiebe et al. [4] and Simpson et al. [9]. suggested that between one to three genes may be responsible for the C-variant morphology, to date, there is limited knowledge of the causal genetic mutations. Following whole genome sequencing of the Cv and pfWT isolates, mutations were detected in gene jg4843 within the genomes of 11 of the 12 Cv isolates. Remarkably, one Cv isolate (Cv6) differed from the WT isolate at a single base pair, which resulted in a missense mutation in jg4843. Mutations in this gene did not occur in the pfWT isolates, and therefore, were only found in the Cv isolates. Notably, only Simpson et al. [9]. used strains isolated from the Quorn production plant, from which a single complementation group was identified. It is therefore possible that gene jg4843 is responsible for this single complementation group in this study, and that additional mutations are responsible for C-variants from other environments, such as the glucose-limited conditions used by Wiebe et al. [4]. Although the isolates identified in this study represent the diversity observed across commercial campaigns, it is still unknown whether multiple different C-variant isolates exist in each campaign or whether a single C-variant predominates.

The jg4843 gene encodes the F. venenatum ortholog of LRG1, a Rho-type GTPase-activating protein that acts as a form of molecular switch, cycling between the active GTP-bound and inactive GDP-bound forms of Rho1 [49]. The small GTPase Rho1 is essential for β-glucan synthesis and the cell wall integrity pathway [50], and is critical for cell wall remodelling and hyphal branching in several fungal species, including A. nidulans [51], Fusarium oxysporum [52], N. crassa [53], and Ustilago maydis [54]. Furthermore, RhoA, the ortholog of Rho1 in Aspergillus niger, was found to be essential for viability and establishment of polarisation [55]. In S. cerevisiae, Rho1 has also been shown to regulate the activity of TORC1 [56]. It has previously been suggested that the highly branched morphologies of C-variants resulted from pleiotropic effects of mutations in genes involved in carbon metabolism, membrane biosynthesis and/or wall biosynthesis [3]. It is therefore possible that mutations affecting the function of FvLRG1 may impact these pathways through the Rho1 signalling pathway, highlighting potential functional mechanisms responsible for C-variant morphology. Furthermore, cAMP and cGMP, molecules involved in other small GTPase signalling pathways, have previously been observed to alter branching morphology of F. venenatum A3/5 and C-variant isolates [17], further highlighting the importance of these signalling pathways in C-variant morphology.

LRG1 has previously been reported to be responsible for maintaining hyphal apical tip extension and restricting excessive branch formation in subapical hyphal regions in N. crassa [11, 48], and has been implicated in control of hyphal branching in other fungal species, including C. albicans [57] and Magnaporthe oryzae [58]. Most of the mutations identified in FvLRG1 were missense mutations, responsible for amino acid substitutions primarily affecting the RhoGAP domain. This domain is essential for the correct functioning of LRG1 [59], and these mutations could potentially result in alterations to the shape and structure of the active site and/or changes to key residues involved in the catalytic process. Previous work in N. crassa assessed the morphological effects of an LRG1 mutant strain, with respect to the wild-type morphology [48]. The morphological abnormalities in N. crassa LRG1 mutants were analogous to those observed previously in F. venentatum C-variant strains [4, 9, 10, 13], as well as to the colony morphology observed in this study. Moreover, the mutation responsible for the hyper branching phenotype in N. crassa, Y926H, was also localised to the RhoGAP domain [48].

Two of the Cv isolates identified in this study did not have missense mutations in FvLRG1. However, Cv2 and Cv8 shared the same SNP 22 bp upstream of the FvLRG1 start codon. In Cv2 this mutation correlated with a significant downregulation in FvLRG1 expression (Fig. 3). This mutation was also present in the Cv2 RNA-seq data (data not shown), occurring in the 5’ UTR, suggesting that this mutation may inhibit mRNA processing and/or translation in Cv2 and Cv8 [60,61,62].

The only Cv isolate that did not have a mutation in FvLRG1 was Cv3. Instead, Cv3 had a missense mutation in the ortholog of CLA4, a Cdc42-activated signal transducing kinase involved in cell growth, septin ring assembly, and cytokinesis [63,64,65]. Like LRG1, CLA4 has also been implicated in the control of hyphal branching in a number of different fungal species, including Trichoderma reesei [66], Ashbya gossypii [67], C. albicans [68], and Aspergillus flavus [69]. Moreover, crosstalk between the Cdc42 and Rho1 pathways has also previously been reported [70, 71]. It is possible then that the CLA4 mutation in Cv3 either directly altered hyphal morphology, or indirectly because of the downregulation of FvLRG1 observed in Cv3 due to the crosstalk between these signal transduction pathways.

Although mutations in FvLRG1 were observed in the majority of the Cv isolates, other mutations unique to these isolates were also identified that could potentially affect growth rates. For example, the Cv isolates from the second batch of samples had much more variable growth rates than the first (Figs. 1 and 4), which appeared to separate into two groups: those around 2–4 mm day^− 1 and those around 6 mm day^− 1. The differences in growth rate could be explained by epistatic effects of additional mutations in these isolates. Cv8 had one other missense mutation in the ortholog of ROM1, a Rho-GTP exchange factor, which acts as the opposing Rho1 activator to LRG1 [72]. Cv10 has three additional missense mutations, one in the ortholog of GUF1, a mitochondrial matrix GTPase important for translation under temperature and nutrient stress in S. cerevisiae [73], another in an ortholog of a c6 transcription factor, and another in an ortholog of BCD1 which is required for the accumulation of box C/D snoRNA [74]. Cv11 has two additional missense mutations, one in an ortholog of MET5, a sulphite reductase beta subunit involved in amino acid biosynthesis [75], and another in an ortholog of Oxidoreductase NAD-binding domain-containing protein. However, we did not investigate these mutations further, and future studies could explore the epistatic effects of these genes on C-variant growth rates.

Whilst we identified mutations unique to Cv isolates, a smaller number of mutations unique to pfWT isolates were also identified. It is possible that some of these mutations may confer beneficial effects. For example, pfWT1 had a significant increase in radial growth rate compared to WT and possessed a single missense mutation in jg6717 gene, the ortholog of Activator of stress genes 1, which is potentially responsible for the observed increase in radial growth rate. Future studies could provide greater understanding of the effects of these mutations on the growth rate of pfWT isolates.

Conclusions

This study represents a significant step towards understanding the genetic basis of C-variant morphology in the Quorn mycoprotein production strain F. venenatum A3/5. The mutation of FvLRG1 has been identified as the causal event for the appearance of C-variant morphology during commercial fermentation campaigns. By revealing the key gene involved in hyper branching, we anticipate that our findings will enable the development of targeted strategies to control or prevent the appearance of C-variants, ultimately enhancing the efficiency, productivity, and quality of Quorn fermentation campaigns. This research provides a foundation for future studies aimed at achieving more sustainable mycoprotein production methods.