Transcriptome analysis of the fungal pathogen Fusarium oxysporum f. sp. medicaginis during colonisation of resistant and susceptible Medicago truncatula hosts identifies differential pathogenicity profiles and novel candidate effectors

Quantification of Fusarium growth in resistant and susceptible M. truncatula accessions

We previously identified the reference M. truncatula accession A17 to display moderate to strong resistance to Fom [29] whilst the DZA315 accession was susceptible [5]. To examine differences in disease progression between these two hosts we infected both alongside each other and quantified Fom growth in root and shoot tissues over an infection time-course (Fig. 1). Within 14 days post inoculation (dpi) 97 % of DZA315 plants had visible disease symptoms of wilting, and chlorotic or necrotic leaves, while only 30 % of A17 plants were diseased and of these, only 6 % of their leaves on average had visible disease symptoms (Fig. 1a, b). By 21 dpi all A17 plants survived, although the number of diseased plants had risen to 40 %, these again displayed limited symptoms (on average 5 % of leaves were diseased per plant, Fig. 1a). In contrast, all DZA315 plants were dead by 21 dpi (Fig. 1c). At 21 dpi the Fom inoculated A17 seedlings were visibly smaller than mock or control inoculated seedlings, showing a 30 % reduction in shoot fresh weight (Fig. 1d). The limited disease progression observed in A17 suggests Fom is able to colonise A17 seedlings but that the molecular resistance mechanisms employed by A17 plants to control pathogen spread results in reduced growth.

To determine the extent of Fom colonisation in A17 and compare this to levels in susceptible accession DZA315, we quantified the relative amount of fungal biomass in root and shoot tissues of A17 and DZA315 seedlings by qRT-PCR over the early stages of infection between 1 and 7 dpi. Both A17 and DZA315 had similar levels of Fom biomass at 1, 2 and 4 dpi in roots, with some growth also detected in shoots (Fig. 1e). However, by 7 dpi the relative abundance of Fom had risen sharply and significantly in DZA315 root tissues to levels ~5-fold greater than levels in A17 which remained similar to those detected at its 4 dpi time-point.

In vitro and in planta Fusarium RNA-sequencing

We hypothesized that increased Fom colonisation of susceptible DZA315 seedlings may be associated with changes in its pathogenicity gene expression profile compared to colonisation of a resistant A17 host. To capture Fom genes expressed in resistant and susceptible plants and to compare and contrast their expression profiles we generated high coverage RNA-Seq data from infected root and shoot tissues of A17 and DZA315 seedlings at 2 dpi where we detected no difference in host colonisation levels and at the later time-point of 7 dpi when disease symptoms start to manifest in susceptible plants (Fig. 1). We confirmed that disease progressed as expected in seedlings used for RNA-sequencing by scoring the percentage of diseased plants and survival at 7-21 dpi on plants from the same experiment that were not harvested for RNA extractions (Additional file 1). We also collected data from Fom mycelia grown vegetatively (in vitro) in order to compare levels of induction upon host detection and use this as a feature to facilitate identification of effectors or other pathogenicity-associated genes.

Identification of differentially expressed Fusarium genes

To identify genes with a high potential for involvement in pathogenicity, we set out to identify Fom genes differentially expressed between vegetative (in vitro) and in planta growth conditions with the premise that genes involved in fungal pathogenicity would be switched on or more highly and rapidly expressed upon detection of a suitable host [6, 18, 19, 31–33]. Aligned read counts (generated from the maploci and genDEseq subprocesses within the BioKanga toolkit [http://sourceforge.net/projects/biokanga/files/]) were used with a normalisation step to identify differentially expressed genes (DEGs) between each dataset (growth condition, dpi, plant accession) with EdgeR [34]. DEGs were selected based on a ≥ 2-fold change and a False Discovery Rate (FDR) ≤ 0.05. Due to low proportion of Fom reads in the RNA-Seq datasets, to decrease the number of false positives based on mapping of only a few reads or mapping of multiple reads to only one location, we added an additional criteria of at least 25 % read coverage of the predicted gene model in each of the 3 in planta replicates. In the early 2 dpi dataset replicates this equated to on average 84−92 % of transcripts with a minimum of 5 mapped reads, with individual gene models covered on average by 85 reads. Similar additional DEG criteria have been used in other lowly in planta expressed pathogen or endophyte studies [35, 36]. Full details of DEGs and characteristics of their encoded proteins are listed in Additional files 2 and 3.

At both 2 and 7 dpi the number of DEGs induced/repressed in planta in DZA315 was greater than those detected in A17 suggesting a larger degree of transcriptional reprogramming in the susceptible plant interaction (Fig. 2, Table 2). This was particularly evident at 2 dpi even when Fom fungal biomass was comparable between susceptible and resistant plant roots (Fig. 1e) with 11 % more DEGs in the DZA315 up-regulated dataset compared to A17, and 30 % more genes with ≥ 95 % read coverage (Table 2). By 7 dpi the difference in number of DEGs between DZA315 and A17 increased by over 3.5-fold, likely associated with the greater fungal colonisation of DZA315 roots (Fig. 1e, Table 2). In DZA315 80 % of DEGs up-regulated at 2 dpi were also up-regulated at 7 dpi, while this was only observed for 57 % of DEGs in A17 (Figs. 2b, c). Considering the larger number of mapped reads and DEGs identified in DZA315 at 7 dpi, less overlap between the 2 and 7 dpi DZA315 datasets might be expected. This result suggests a fast, prolonged, and concerted expression of components of the Fom pathogenicity arsenal during infection of susceptible plants. With the exception of Fom_00898 (expression down at 2 dpi DZA315 but up at 7 dpi DZA315) there were no genes detected in both up- and down-regulated datasets of either genotype. Interestingly Fom_00898 encodes a protein with characteristics of a small secreted protein (SSP; protein length ≤ 300 amino acids, predicted to be secreted (SignalP) and containing ≤ one transmembrane domain in the N-terminal region [5]), contains a CFEM domain possibly associated with fungal pathogenicity [37], and is predicted as a putative effector by the fungal effector prediction software EffectorP [38].

Proportion of gene model coverage	A17 2 dpi vs IV		A17 7 dpi vs IV		DZA315 2 dpi vs IV		DZA315 7 dpi vs IV
	up	down	up	down	up	down	up	down
≥25 %	107	90	128	88	117	108	446	333
≥50 %	64	21	63	34	66	35	281	141
≥95 %	13	1	25	5	17	4	105	28

IV in vitro. Note, DEGs in ≥ 95 % are also captured within the ≥ 50 % and ≥ 25 % datasets. Likewise, DEGs in ≥ 50 % are captured within the ≥ 25 % dataset

We also examined RNA-Seq data from the same experiment, at the same level of sample read coverage obtained from the shoots of the Fom inoculated DZA315 or A17 plants but due to the relative low abundance of fungal transcripts in shoot tissues at our early sampling time-points only a handful of DEGs could be identified (data not presented).

Protein characteristics of differentially expressed genes with a focus on in planta up-regulated datasets

To assess differences in pathogenicity profiles of Fom during interactions with susceptible or resistant Medicago hosts, we focussed the remainder of our studies on DEGs up-regulated in planta versus in vitro. To interpret putative functions, we interrogated their predicted protein characteristics including gene ontology (GO) terms, Pfam domains and pathogenicity associated characteristics such as encoding SSPs or similarity to known fungal effectors. Although fungal effectors generally lack similarity to other proteins, they may share common protein motifs or characteristics such as a secretion signal, small size (generally less than 300 aas), increased number of paired cysteines, proximity to repetitive DNA and these were therefore included in our analyses. Fom protein characteristics were previously determined by [5] and supported by annotations listed in publicly available Fusarium spp. genome projects, PHI-base [39, 40] and GenBank. We also incorporated an assessment of putative chromosomal location based on our previous study [5]. Fom scaffolds representing putative CD chromosomes were predicted based on the criteria of having no match to designated core chromosomes of F. oxysporum f. sp. lycopersici or F. solani, whose genome assemblies contain well characterised core and CD/accessory chromosome sequences.

One fifth of the predicted proteins in up-regulated DEG datasets could not be assigned functional annotations and were annotated as uncharacterised. This included 24 (22 % of total in the dataset) and 27 (21 %) in A17 at 2 and 7 dpi respectively. In DZA315 infected plants there were a similar number of uncharacterised proteins at 2 dpi (22 DEGs; 19 % of dataset) but within the 7 dpi dataset the total number of uncharacterised proteins more than tripled to 78 (representing 17 % of the dataset). Nearly half of the latter were predicted as secreted and as putative effectors by EffectorP.

Of proteins in the up-regulated datasets that could be assigned GO terms, 14−18 % were assigned to biological processes, 16−18 % molecular functions, and 3−4 % cellular components. Thirty-six and 33 GO terms could be assigned respectively to A17 and DZA315 2 dpi datasets, and 45 and 69 respectively in the 7 dpi datasets. Of these, metabolic and catalytic activity represented the majority of classified proteins in all datasets, followed by catabolism, binding, hydrolase activity and biosynthesis in different percentages (Additional file 4). The number of encoded proteins represented by these GO terms were 3-6 times more in the DZA315 7 dpi dataset compared to A17 at the same time-point. For example, DZA315:A17 carbohydrate metabolism 22:7, lipid metabolism 13:2, hydrolase activity 39:18, catabolism 37:10, signal transduction 5:0 (Additional file 5).

For Fom genes up-regulated in planta, a significant over-representation (Fisher’s exact test p ≤ 0.05) of domains associated with degradation of proteins and sugars/carbohydrates (e.g. glycoside hydrolase, pectate lyase, protease), membrane transport (e.g. sugar, amino acid and Major Facilitator Superfamily (MFS) transporters, nucleobase cation symporter-1, permease) and oxidative processes (e.g. 3-beta hydroxysteroid dehydrogenase, oxidoreductase, cytochrome p450s) was observed, particularly in DZA315 (Fig. 3, Additional file 6). Encapsulated, these results suggest host cell wall and membrane degradation along with nutrient transport are initiated earlier upon infection of susceptible plants, evident as early as 2 dpi compared to the resistant host interaction. By 7 dpi these processes were even more apparent during Fom infection of susceptible plants, correlating with an increase in fungal biomass at this time point (Fig. 1e).

Searching for proteins involved in regulation of Fom pathogenicity-associated gene expression, we identified four proteins annotated as fungal transcription factors with one (Fom_15733) common to three of the datasets (2 dpi A17, 7 dpi A17 and DZA315). Fom_15733 is predicted CDC encoded and shares greatest similarity with proteins from F. oxysporum f. sp. raphani (Brassica pathogen, 90 %) and f. sp. pisi (pea pathogen, 87 %) but low similarity with other F. oxysporum ff. spp. (e.g. 26−36 % with ciceris, melonis, lycopersici, Fo5176, see [5]). This protein contains a fungal specific transcription factor domain (Pfam:PF04082) and a fungal Zn(2)-Cys(6) binuclear cluster domain (Pfam:PF00172). The other transcription factors identified in the up-regulated datasets were only found in the DZA315 7 dpi dataset (Fom_14279, 14162, 07202), are predicted to reside on core chromosomes and share most similarity to proteins from other F. oxysporum ff. spp.. The Fom homologue (Fom_08318) of Fol transcription factor SIX Gene Expression 1 (SGE1) whose expression is up-regulated during infection of tomato roots and is required for expression of most secreted Fol effectors [41, 42], was detected as expressed in vitro but not detected as significantly up-regulated in our in planta DEG datasets.

Next we applied several criteria to identify candidate effector and host-specific pathogenicity genes. These included characteristics such as predicted secretion, in planta up-regulated gene expression and chromosomal location, as F. oxysporum effectors have previously been identified on non-core or conditionally dispensable chromosomes (CD chromosomes). The majority of DEGs in all datasets were located on scaffolds predicted to form part of core Fom chromosomes [5] (Fig. 4a). Most DEGs predicted to lie on putative CD chromosomes were only identified within the up-regulated datasets (Figs. 4a, b). DEGs encoding SSPs were also predominantly identified within the in planta up-regulated datasets (Fig. 4c, Additional files 2 and 3). Interestingly the number of SSPs in A17 datasets didn’t differ much between the two sampled time-points but in DZA315 the number almost tripled between the 2 and 7 dpi up-regulated datasets (Fig. 4c). Further, DEGs within the susceptible DZA315 datasets also contained a larger proportion of genes expressed more highly in planta and these were enriched for SSPs (Additional file 7).

Overall, the enrichment during the early stages of infection of genes encoding SSPs, uncharacterised proteins and proteins with roles in protein/sugar degradation and their transport, oxidative stress and other pathogenicity associated processes indicates substantial changes in pathogen gene expression upon colonisation of a susceptible host.

Highly up-regulated genes unique to in planta up-regulated datasets from a susceptible or resistant host and links to pathogenicity

As we were most interested in identifying pathogenicity factors we focused firstly on the differences between genes expressed during infection of resistant or susceptible host accessions (indicating changes that Fom may undergo when it detects a susceptible or resistant host) and secondly on those that were expressed during both infection of a susceptible and a resistant host (indicating key roles in pathogen attack).

Firstly we compared DEGs unique to infection of each host accession. Under this analysis, 11 and 6 DEGs were uniquely significantly up-regulated in A17 at 2 and 7 dpi respectively, while 16 and 280 were unique to DZA315 at the same time-points (Fig. 5). Of those expressed during infection of A17 at 2 dpi, 4 of the most highly induced (32-13,000-fold) were co-localised at the genomic scale (Table 3) and another, Fom_08477, is a predicted SSP showing similarity to an uncharacterized protein from the rice pathogen Fusarium fujikori. At 7 dpi in A17 all six unique DEGs were core scaffold located with two induced >100-fold and encoding SSPs (Fom_05133, Fom_09362, encoding a hypothetical protein and a glycoside hydrolase respectively).

Dataset	Gene IDs	Chrom./ scaffolda	Pfam domain / best BlastP match / other characteristics
A17 2 dpi	Fom_11200, 11207, 11209, 11211	core (scaffold_15; 1.38 Mb)	GST, MFS transporter, FAD-binding, uncharacterised protein
DZA315 7 dpi	Fom_01388, 01394, 01397	core (scaffold_2; 3.35 Mb)	glycoside hydrolase, fasciclin, MFS transporter
	Fom_01935, 01936, 01938	core (scaffold_2; 3.35 Mb)	arca-like protein, MFS transporter, glycoside hydrolase
	Fom_03730, 03733, 03734	core (scaffold_4; 2.62 Mb)	Hypothetical, MFS transporter, synthase
	Fom_15853, 15855, 15856, 15860	core (scaffold_131; 0.03 Mb)	FAD-binding protein, ADH, FAD-binding protein, p450
	Fom_14230, 14231, 14232, 14241	CDC (scaffold_31; 0.21 Mb)	Oxidoreductase, MFS transporter, FAD-binding protein, oxidase

^apredicted core or CD chromosome (CDC) location and scaffold size based on [5]

The 16 unique DEGs from the DZA315 2 dpi dataset comprised of both core and non-core located genes with inductions ranging from 3.5 to >65,000-fold over in vitro conditions and encoding several FAD-binding proteins, hydrolases, peptidase, MFS transporter and p450 amongst others. Of the 280 significant DEGs in the DZA315 7 dpi unique dataset the majority were core scaffold located (243) and 43 were up-regulated in planta over 100-fold. Several were also co-localised (Table 3). Two Fom SSPs (Fom_08816, Fom_11981) with similarity to the phytotoxin cerato-platanin were also only identified in the DZA315 in planta 7 dpi up-regulated dataset. Some cerato-platanin members in other phytopathogens are implicated in disruption of host cell walls through expansin-like activity, chitin oligomer sequestration and the ability to induce plant cell necrosis [43, 44]. Other proteins only upregulated in DZA315 at 7 dpi included a xylosidase arabinosidase (Fom_12400) and an l-arabinitol 4-dehydrogenase (Fom_00399) with implied roles respectively in the hydrolysis of major hemicellulose component xylan to xylose and other sugars, and the catabolism of L-arabinose, an important constituent of plant cell wall polysaccharides.

Overlapping highly in planta up-regulated genes in susceptible and resistant host accessions and links to pathogenicity

We first assessed the overlapping DEGs for the presence of those encoding the well documented F. oxysporum Secreted In Xylem (SIX) proteins, with known roles in virulence and/or avirulence described for some in F. oxysporum f. sp. lycopersici, F. oxysporum f. sp. melonis or Fo5176 (Brassica infecting) [6, 13, 14, 17, 18, 45–49]. SIX proteins can broadly be defined as small, generally cysteine rich, and possessing a secretion signal [13, 49]. Of the 14 F. oxysporum f. sp. lycopersici SIX proteins described so far, all four that we previously identified in the Fom genome (Fom SIX1, SIX8, SIX9 and SIX13) [5] were significantly upregulated in all in planta datasets where they were amongst the top in planta induced DEGs. Fom_SIX13 was the most highly induced at levels 1 million times greater than in vitro (Table 4). qRT-PCR verification of Fom SIX gene expression over a 1 to 7 day infection time-course revealed all but Fom_SIX13 could be reliably detected under in vitro conditions suggesting expression of this SIX gene in Fom responds specifically to detection of a possible host (Fig. 6, Additional file 8). Correlating with our RNA-Seq data (Table 4), qRT-PCR examination of all four Fom SIX genes confirmed these genes were overall more highly expressed in DZA315 than A17. Interestingly Fom_SIX13 and Fom_SIX1 didn’t meet the statistical cut-off of the EffectorP [38] prediction algorithm (as also observed for putative Fol effectors SIX7 and SIX13).

After Fom_SIX13, the next three most highly in-planta induced DEGs are strong effector candidates. These encode Fom_16257 a Fom-unique protein with unknown function, Fom_16301 a SSP encoding a LysM-domain (~40 residue lysin motif implicated in chitin binding [50–52]), and Fom_16235 an uncharacterised protein with orthologs in other F. oxysporum ff. spp.. qRT-PCR verification of Fom_16257 expression confirmed increased expression in DZA315 over A17 and absence of expression under in vitro conditions (Fig. 6). Similar results were observed for Fom_16301. Interestingly, Fom_16235 was lowly expressed under in vitro conditions. Along with Fom_16301, another LysM-domain containing SSP gene (Fom_04092) was also significantly upregulated, but only at 7 dpi. The Fom_16301 encoded protein contains one LysM domain, sharing best similarity to a protein from F. oxysporum f. sp. pisi (95 %), other F. oxysporum ff. spp. and Colletotrichum species. Fom_04092 contains two LysM domains and is most similar to proteins from other F. oxysporum ff.spp. and the rice pathogen F. fujikuroi.

Other highly in planta up-regulated genes included three proteases/peptidases and genes involved in oxidative stress responses (Table 4). Of the latter type, Fom_15948 and Fom_15949 both encode peroxidases predicted to reside alongside each other. Fom_15948 is a predicted effector and SSP sharing 95 % identity to other F. oxysporum peroxidases/catalases or hypothetical proteins, as well as to proteins from F. solani, F. fujikori, and Colletotrichum and Verticillium species.

Several uncharacterised proteins with properties of effectors were also identified as highly up-regulated in-planta during infection of both resistant and susceptible M. truncatula accessions (Table 4). Of those predicted as SSPs, only two were not predicted as effectors by EffectorP and these encode a putative endoglucanase (Fom_02353) and a cupredoxin domain-containing protein (Fom_05481) sharing 94-96 % similarity with several other serine rich proteins from F. oxysporum ff. spp. The serine rich nature of this protein may imply possible involvement in post-translational regulatory processes or protein-protein interactions. The remaining novel uncharacterised effector-like proteins identified contain no other known functional protein domains.

Isolate sources and growth conditions

F. oxysporum f. sp. medicaginis (Weimer) W.C. Snyder & H.N. Hansen, (Fom-5190a, BRIP 5190a/IMI 172838, collection number 19911) was obtained from the Queensland Plant Pathology Herbarium (BRIP) and has been described previously [5]. Fom was maintained on sterile filter paper and grown on ½ strength potato dextrose agar at 22 °C in the dark. M. truncatula accessions A17 and DZA315 were at least second generation inbred lines derived from seed obtained from the South Australian Research and Development Institute (SARDI). M. truncatula seeds were germinated on damp filter paper, transplanted into 30 mm Jiffy-7 peat pellets and grown under a short day-light regime (8-h light/16-h dark) at 21 °C.

F. oxysporum disease assays

Fom disease assays of M. truncatula were conducted as described previously [5, 29]. Briefly, a 1 × 10⁶ spores mL^–1 spore suspension was used to inoculate the roots of two week old seedlings which had roots protruding from the peat pellets trimmed prior to inoculation. Seedlings in peat pellets were placed in a petri dish of spore suspension for 5 min, followed by addition of a further 1 mL of spore suspension to the base of the hypocotyl. Inoculated pellets were transferred to growth trays lined with a plastic sheet and a thin layer of damp vermiculite, covered with a clear plastic dome to maintain humidity, and incubated under a long-daylight regime (16-h light/8-h dark) at 28 °C. Plants were scored every 7 days over 4 weeks to assess disease progression.

RNA isolation

RNA for qRT-PCR and RNA-Seq experiments on Fom infected M. truncatula root tissue was conducted as described previously [5]. Briefly, root tissue was collected from 10 plants and pooled per replicate at 1, 2, 4 and 7 days post inoculation. For Fom in vitro samples, mycelia were grown in a petri dish containing one-half-strength potato dextrose broth for 7 days at 22 °C and mycelia harvested by filtering through Miracloth. RNA extraction was performed independently for each of 3 replicates using a Trizol extraction (Sigma-Aldrich, St. Louis, MO) followed by DNase treatment using TURBO DNase (Ambion). RNA samples were cleaned via RNeasy mini spin columns (Qiagen).

qRT-PCR

Following RNA isolation and DNase treatment, complementary DNA synthesis was performed with 1ug of input RNA followed by qRT-PCR performed using SsoFast EvaGreen Supermix (Bio-Rad) on a CFX384 (Bio-Rad) system as described previously [5]. Absolute gene expression levels relative to F. oxysporum housekeeping gene GPDA (Fom_05751) were used for each complementary DNA sample using the equation: relative ratio gene of interest/GPDA = (Egene^{-Ct gene})/(EGPDA^{-Ct GPDA}) where Ct is the cycle threshold value. Medicago root samples were verified for even abundance of plant input material using the M. truncatula B-tubulin reference gene [5, 83] which was found to be within ±1 Ct across all samples. Quantification of Fom biomass in M. truncatula tissues was conducted as described previously by determining the abundance of Fom_18S relative to M. truncatula_18S [5]. Primer sequences have been previously published [5] and/or are also listed in Additional file 9.

RNA-Seq library construction, Illumina sequencing and read-mapping

RNA integrity was confirmed using the Agilent 2100 Bioanalyser Plant Nano system (Agilent Biotechnologies). Stranded Illumina TruSeq libraries were generated from 1 μg of total RNA and sequenced (100 bp paired end reads) on an Illumina HiSeq platform by the Australian Genome Research Facility (AGRF). As per [5], RNA-Seq paired-end reads were trimmed for low-quality base-calls (≥ q30) and Illumina adapter sequences using Cutadapt v1.1 [84] (parameters: --quality-cutoff 30 --overlap 10 --times 3 –minimum-length 25) with reads trimmed to less than 25 bp discarded and remaining reads sorted into pairs and singleton reads. RNA-Seq reads were then mapped to the Fom-5190a genome assembly via Tophat2 (TopHat2 v2.0.9) (parameters: --b2-very-sensitive -r 80 --mate-std-dev 40 -i 20 -I 4000 -g 20 --report-secondary-alignments --report-discordant-pair-alignments -m 0 --min-coverage-intron 20 --microexon-search --library-type fr-firststrand) [30]. BAM files generated from Tophat2 output were sorted using SAMtools version 0.1.19 [85] and visualised using IGV (version 2.3) [86]. Trimmed sequencing data is available from the NCBI/GenBank database under BioProject number PRJNA294248 (http://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA294248).

Differential gene expression analysis

BAM files generated by TopHat2 [30] were used by BioKanga maploci [http://sourceforge.net/projects/biokanga/files/] to assign aligned reads to known loci. This step was followed by the genDEseq subprocess of BioKanga to generate counts files for differential expression. Read counts from different biological replicates and samples were combined for each gene and the resulting count matrix was normalized and analysed for differential expression using the EdgeR [34] package in R version 3.2.2 http://www.r-project.org/ using a false discovery rate cutoff of 0.05 and a minimum fold change of log2 ≥ 1. Identification of unique or overlapping genes within the DEG datasets and the generation of Venn diagrams was determined using Draw Venn Diagram http://bioinformatics.psb.ugent.be/webtools/Venn/ (accessed 12-01-16). Gene coverage was calculated from TopHat2 mapped RNA-Seq reads using BEDTools (v2.21.0) coverageBed [87].

Protein characteristics

Protein characteristics were previously reported by [5] with GO term counts determined via CateGOrizer [88] (accessed 15-01-16). Statistical examination of over-represented protein functional attributes based on the number of proteins with specific Pfam domains were compared between DEG datasets and those from the whole annotated Fom genome using Fisher’s exact test with a significance threshold of p ≤ 0.05.