Methodology:

25 individuals were sequenced from the 5 sets of sampling/experimental work detailed above. All analyses were performed by Novogene and comprised the following stages and application of software packages. Raw reads were quality controlled for error rate using Phred and GC content using the Illumina CASAVA v1.8 software. Reads were cleaned with the removal of Illumina kit adapter sequences and adapter contamination, where the level of uncertain nucleotides N > 10% and where low quality bases as defined by the Phred Q20 score constituted more than 50% of a read. De novo transcriptome assembly was performed using Trinity version r20140413p1 with parameters min_kmer_cov=2, min_glue=2, others were set to default, using the modules Inchworm, Chrysalis and Butterfly (Grabherr et al., 2011). Hierarchical clustering was performed using the Corset program in Trinity to remove read redundancy. The longest transcripts from each cluster were selected as unigenes. Annotation of the unigenes was performed using seven databases (NR, NT, KO, SwissProt, Pfam, Go and KOG). Blast searching against NT was performed using NBCI blast 2.2.28+ with an e-value threshold of 1e-5 (Altschul et al., 1997). Diamond v0.8.22 (Buchfink et al., 2015) was used to blast search the unigenes against NR, SwissProt and KOG. The e-value threshold for NR and SwissProt was 1e-5 and 1e-3 for KOG. Pfam (Finn et al., 2008) was screened using the hmmscan package in HMMER v3.1b1 with an e-value threshold of 0.01. GO annotation was based on the protein annotation results from NR and Pfam using Blast2GO vb2g4pipe_2.5 (Goetz et al., 2008; Young et al., 2010) with an e-value threshold of 1e-6. KEGG mapping was performed using KAAS (KEGG Automated Annotation Server) v.r140224 with an e-value threshold of 1e-8 (Mao et al., 2005; Moriya et al., 2007; Kanehisa et al., 2008). GO enrichment was performed using GOSeqtopGO vGOSeq 1.32.0, topGO-2.32.0) with a corrected p value of <0.05. KEGG enrichment was performed using KOBAS v3.0 with a corrected p value of <0.05. The de novo transcriptome was used as a reference assembly and the reads from each library mapped back to the transcriptome and quantified using Bowtie2 vbowtie2-2.2.2.2 and RSEM vRSEM-v1.3.0 (Li et al., 2011) with output referenced as FPKM (Fragments Per Kilobase of transcript sequence per Million base pairs sequenced). The threshold for expression was set at FPKM >3.0. Differential expression between the different sets of samples was calculated using DEGseq v1.12.0 (Wang et al., 2010) with normalization via TMM and FDR calculated using BH (Benjamini and Hochberg, 1995) with output threshold of log2fold change >1 and adjusted p value <0.005. Protein gene identifiers were extracted from the SwissProt annotations for the MI v. MO differential expression analysis for analysis via the STRING v11 program (https://string-db.org/) to visualize protein-protein interactions (Szklarczyk et al., 2019).