De Novo Transcriptome Assembly and Annotation of Liver and Brain Tissues of Common Brushtail Possums (Trichosurus vulpecula) in New Zealand: Transcriptome Diversity after Decades of Population Control.

The common brushtail possum (Trichosurus vulpecula), introduced from Australia in the mid-nineteenth century, is an invasive species in New Zealand where it is widespread and forms the largest self-sustained reservoir of bovine tuberculosis (Mycobacterium bovis) among wild populations. Conservation and agricultural authorities regularly apply a series of population control measures to suppress brushtail possum populations. The evolutionary consequence of more than half a century of intensive population control operations on the species' genomic diversity and population structure is hindered by a paucity of available genomic resources. This study is the first to characterise the functional content and diversity of brushtail possum liver and brain cerebral cortex transcriptomes. Raw sequences from hepatic cells and cerebral cortex were assembled into 58,001 and 64,735 transcripts respectively. Functional annotation and polymorphism assignment of the assembled transcripts demonstrated a considerable level of variation in the core metabolic pathways that represent potential targets for selection pressure exerted by chemical toxicants. This study suggests that the brushtail possum population in New Zealand harbours considerable variation in metabolic pathways that could potentially promote the development of tolerance against chemical toxicants.


Introduction
New Zealand separated from other landmasses approximately 85 million years ago and has remained largely isolated ever since [1]. Ecological isolation created distinctive evolutionary dynamics in the archipelago rarely seen elsewhere on the planet. By the time the first humans colonised in approximately 1300 AD [2,3], no native land mammals, other than three species of bats, were present [4]. The negative impacts of two consecutive waves of human colonisation triggered a series of fundamental alterations in the biological diversity of New Zealand. The introduction of three species of rats (Rattus spp.), three species of mustelids (Mustela spp.), and the common brushtail possum (Trichosurus vulpecula), in combination with extensive human hunting and habitat destruction, triggered a spiral of catastrophic events that drove fifty-nine bird species extinct and severely threatens the survival of other native species, including endemic marine mammals, native bat species, frogs, and much of the invertebrate fauna [5,6].
The common brushtail possum (hereafter 'brushtail possum') was introduced into New Zealand from localities in Tasmania and the eastern Australian mainland in the 1850s to establish a profitable fur industry [7,8]. Soon after its introduction, the species established itself in the wild, proliferated, and spread across all habitats in the archipelago. In these new habitats, brushtail possums reached population densities far exceeding those in their native Australian range [4]. The phenomenal success of brushtail possums in New Zealand has been linked to the high abundance of palatable and nutritious food sources, and the presence of fewer parasites, predators, and competitors [9]. The scale of devastation upon New Zealand's endemic ecosystems prompted the brushtail possum's designation as both an agricultural pest and a major threat to New Zealand's unique flora and fauna [9,10]. Wild populations of brushtail possums also harbour a self-sustained reservoir of bovine tuberculosis (Mycobacterium bovis) that threatens New Zealand's dairy and livestock industries [4]. As a result, the brushtail possum has been subject to decades of targeted population control operations.
Conservation and agricultural authorities in New Zealand regularly apply a series of population eradication and control measures to mitigate the negative impacts of brushtail possums on native ecosystems. Since 1957, aerial dispersal of cereal baits containing the chemical toxicant sodium fluoroacetate (or 'compound 1080 ) has constituted the cornerstone of management plans to control brushtail possums and other invasive mammals in New Zealand [10,11]. Synthetic compound 1080 is a water-soluble salt modelled on the naturally occurring plant toxin fluoroacetate [12]. Several plant species across the globe, such as Gastrolobium spp. and Oxyylobium spp., use this toxicant as a deterrent against herbivory, including species within the possum's native Australian home range [13][14][15].
Brushtail possum populations in New Zealand face a combination of stressors and human-induced selection pressures that differ from those in their native Australia. The speed of local adaptation typically found in many invasive species suggests that the genetic recipe for adaptation predominately arises from standing genetic variations present in the founder populations [16,17]. Triggs and Green [18] and Sarre et al. [19] reported that current populations of the brushtail possum in New Zealand originated from two distinct subspecies, Trichosurus vulpecula fuliginosus from Tasmania and T. vulpecula vulpecula from mainland Australia, and that extensive hybridisation took place following their introduction to New Zealand. Genetic variations among the founder populations propagated in the non-native range and most likely became reshuffled into new genotypes that do not occur in Australia.
Recent advances in DNA sequencing technology and high-performance computing have facilitated the comprehensive study of functional regions in eukaryote genomes. Variability in the genome has made it possible to address some fundamental questions in evolutionary genetics, such as the prediction of pathways associated with fruit colour [20], intraspecific variation in morphology and behaviour [21], susceptibility to pathogens [22][23][24][25], responses to chemicals [23], adaptation to environmental changes [26], species diagnostics markers in endangered shrimps [27], sex determination [28], and changes in the expression profiles of an organ during embryonic development [29].
Understanding the brushtail possum's evolutionary dynamics in New Zealand is hindered by a paucity of available genomic resources for this species [30]. In particular, the evolutionary consequences of more than half a century of intensive population control measures on its genetic diversity and population structure merit in-depth investigation. To address these fundamental questions, genomic resources are required.
In the current study, we characterise the functional content and genetic diversity of brushtail possum liver and brain cerebral cortex transcriptomes, tissues of major importance in responding to chemical toxicants and in modulating complex behaviour [31][32][33][34]. The core metabolic pathways in these tissue types are among the primary targets of selective pressures imposed by anthropogenic chemicals [35][36][37]. This study serves as a stepping-stone for a better understanding of the brushtail possum's evolutionary dynamics in New Zealand and will contribute towards the optimisation of a comprehensive population management plan for the species, in line with the New Zealand Predator-Free 2050 horizon objectives.

Animal Ethics Statement
Animal ethics approval for this study was obtained from the Lincoln University Animal Ethics Committee (AEC 586). Approved standard operation procedures were used for both animal husbandry and euthanasia. All necessary steps were taken to minimise the pain and suffering of the subject animals.

Specimen Collection
In total, four brushtail possums (two males and two females) were live trapped on the east and west coasts of New Zealand's two main islands (Banks Peninsula and Buller, South Island, and Taranaki and Hawkes Bay, North Island). The animals were transferred to the wild animal husbandry facilities (Johnstone Memorial Laboratory) at Lincoln University. On the day of specimen collection, sex, weight, and health status of each animal were recorded by Lincoln University's wildlife animal facility staff. The brushtail possums were humanely euthanised by intracardiac administration of 300 mg of Pentobarbital. The dissection and collection of the liver and cerebral cortex biopsies were performed after the death of the animals was confirmed on site. Approximately 1 mm 3 tissue from liver and brain cerebral cortex were excised, transferred to 10:1 volume of RNAlater solution (Qiagen-Hilden) and stored at −80 • C until RNA was extracted within one week after specimen collection. Liver cells are assumed to be more homogeneous in their gene expression profiles, but see Reference [38], whereas brain cell expression profiles often differ considerably between different brain regions [39,40]. To account for variations in the gene expression profiles in the possums' cerebral cortex, brain biopsies from three individuals were excised from different sections of the proximal cerebral cortex in the left and right cortical lobes. The remaining biopsy was taken from deeper white matter in the cerebral cortex. A wildlife histologist examined all brain biopsies to validate their origin.

Nucleic Acid Extraction, Genomic Library Preparation, and Sequencing
Total RNA was purified from liver and brain tissue using Qiagen's RNeasy Mini Kit (Hilden, Germany). The extracted RNA was eluted in 50-100 µl of water. The quality of the purified RNA was measured using Agilent's Bioanalyzer (Santa Clara, United States) to obtain an RNA Integration Number score before proceeding with complementary DNA preparation. Nugen's Ovation RNA-Sequencing System (Redwood City, United States), automated on an Apollo324 liquid handler (Takara Bio, Kyoto, Japan), was used to create complementary DNA using a single primer isothermal amplification method and 100 ng of template RNA. The resulting complementary DNA was quantified on a Nanodrop spectrophotometer (ThermoFisher Scientific, Massachusetts, United States) and sheared to 200 bp fragments using the Covaris M220 ultrasonicator (Woburn, Unites States). Illumina-compatible sequencing libraries were prepared using KAPA Biosystem's LTP library preparation kit (Pleasanton, United States). Individual sequencing adapters with sample-specific barcodes (Bioo Scientific, Austin, United States) were ligated to the end-repaired fragments, the resulting libraries were cleaned using AMPure beads (Beckman-Coulter, Brea, United States) and amplified for ten PCR cycles with KAPA's HiFi polymerase. The concentration of each library was quantified separately using qPCR on the Quantstudio 5 (Applied Biosystems, Foster City, United States) and the fragment sizes were verified using the Agilent Bioanalyzer. All sequencing libraries were normalised to 2 nM concentration before pooling for multiplex sequencing on two 1 × 50 flow cells on an Illumina HiSeq2000 platform (San Diego, United States) at the Arizona State University genomics core facility.
A combination of SAMtools v1.9 mpileup command [53] and Varscan2 v2.3.7 [54] were used to identify single nucleotide polymorphic sites (SNPs) in each tissue type. Default settings were used for variant calling, an exception being the p-values in Varscan2, which were lowered to 0.001 to increase prediction accuracy. Functional impacts of the polymorphic sites on protein sequences were predicted using the KisSplice2reftranscriptome v1.3.3 pipeline [55].
Two datasets were created and investigated in parallel to unravel variability in the functional content of each tissue type. The first consisted of ten transcripts within each tissue's reference transcriptome that showed the highest number of polymorphic sites, as reported in the Varscan2 output files. In the second dataset, the abundance of each transcript was estimated in RSEM v1.3.1 [56], and 5% of the transcripts with the highest normalised level of expression (estimated as per million mapped reads, TPM) were sub-selected for further functional analysis.
The enrichments of the ten most polymorphic transcripts in each tissue type, for particular biological functions in terms of Gene Ontology categories (GO terms) [57], were investigated using the g: Profiler online server [58] and the annotated genome of the grey short-tailed opossum, Monodelphis domestica, was used as the closest high-quality annotated gene set.
The tissue-specific expression profiles and intra-individual variations in the functional content of polymorphic sites were statistically tested using the WEGO v2.0 online tool [59] and visualised using a combination of the REVIGO online server [60] and CirGO [61]. The signal of pervasive selection (the ratio of non-synonymous to synonymous amino acid substitution, dN/dS ratio) in one particular liver transcript with an exceptionally high number of polymorphic sites was tested using a fixed-effects likelihood model [62] implemented in the online server Datamonkey [63].
The functional content of the assembled transcripts and two generated datasets, in terms of core metabolic pathways, was reconstructed by mapping transcripts against the KEGG database using the KEGG mapper online server [64] and results were summarized in GAEV [65].

Results
The histological reports confirmed that three of the brain cerebral cortex biopsies originated from the proximal cortex, while the remainder originated from deeper within the cerebral cortex white matter (Supplementary Materials Figures S1-S4).
The RNA sequencing run yielded 183,517,606 and 178,725,338 raw sequences from liver and brain cerebral cortex, respectively. Liver tissue sequences were assembled into 58,001 transcripts (candidate genes and isoforms), with a mean contig length of 813.37 and a GC content of 42.57%. From the cerebral cortex brain tissue, Trinity assembled sequences with a total of 64,735 transcripts, a mean contig length of 738.13, and a GC content of 42.68% (Table 1). CDHIT subsumed all transcripts into 57,221 and 63,679 unique clusters with more than 98% sequence identity for liver and brain cerebral cortex, respectively. Table 1. Summary statistics of the reference transcriptome assembly of the New Zealand brushtail possum's liver and brain cerebral cortex tissue as estimated in QUAST v.4.0 [66].

Assembly Name
Liver Cells Brain Cells The semi-global BWA-MEM aligner performed better, in terms of mapping quality and percentage of aligned sequences, compared to two other aligners, and was selected for downstream analysis. The overall alignment rate was estimated at 90.93% for liver sequences and 79.80% for brain cerebral cortex sequences. The alignment rate in both tissue types lies within the recommended range in the Trinity pipeline for a complete reference transcriptome assembly (https://github.com/trinityrnaseq/ trinityrnaseq/wiki/RNA-Seq-Read-Representation-by-Trinity-Assembly). From a total of 9226 searched BUSCOs, 3290 and 2446 of the core mammalian's complete BUSCOs were identified in the hepatic and brain cerebral cortex transcriptomes, respectively (Supplementary Materials Table S1). Tissue-specific transcriptomes, especially in non-model organisms, are very unlikely to produce complete mammalian BUSCOs and that should be interpreted within the biological context of the study, rather than as a direct metric for estimating the completeness of the assembled transcriptomes.
Varscan2 identified 35,300 (one SNP for every 1339 bp) and 23,805 (one SNP for every 2007 bp) variable sites in liver and brain cerebral cortex transcripts, respectively. The number of polymorphic sites with transitions was estimated to be three times higher than the number of transversions in both tissue types. The polymorphism landscape in brushtail possum cerebral cortex was more heterogeneous compared to hepatic cells. Only 25% (5939 out of 23,805) of the polymorphic sites in cerebral cortex tissues were successfully genotyped across all four sampled tissues compared to 70% shared polymorphisms (24,924 out of 35,225) among hepatic tissues.
KisSplice identified 34,329 (22,887 synonymous and 11,442 synonymous) and 4695 (3244 synonymous and 1451 non-synonymous) codon substitutions in hepatic and cerebral cortex tissues, respectively. The reported high levels of heterogeneity in gene expression profiles among different sections of the cerebral cortex provide the most likely explanation for the observed pattern of codon substitutions.
In the hepatic cells, non-synonymous mutations of the codons GTT to ATT, ACG to ATG (substitution of amino acid valine to isoleucine) and GTC to ATC (substitution of amino acid threonine to methionine) had the highest frequency. In the brain cerebral cortex, non-synonymous mutation of the codons GTG to CTG (substitution of amino acid valine to leucine) AAA to AGA (substitution of amino acid lysin to arginine) and AGT to AAT (substitution of amino acid serine to asparagine) were prevalent. The highest rate of synonymous amino acid substitutions was observed for aspartate (mutation in codon GAC to GAT or GAT to GAC) and asparagine (mutations in the codon AAC to AAT) in the hepatic cells and for amino acids aspartate (mutation in the codon GAC to GAT) and proline (mutation in the codon CCA to CCG) in the brain cortical cells (Supplementary Materials Tables  S4 and S5) A transcript in the liver was of particular interest because of an exceptionally large number of polymorphic sites (127 SNPs). The BLASTX search of this transcript showed high sequence identity (88% similarity) to various isoforms of multidrug resistance proteins reported in koala, Phascolarctos cinereus (NCBI reference XP_020843128.1), and common wombat, Vombatus ursinus (NCBI Reference Sequence: XP_027703707.1). In all four brushtail possums, alignment of the protein-coding sequences of this transcript demonstrated 19 candidate sites that were under pervasive purifying selection (dN/dS < 1, p-value < 0.05). The majority of the non-synonymous amino acid substitutions in these transcripts (6 out of 30) were in the form of valine to isoleucine substitutions. The synonymous substitution for amino acids alanine and isoleucine (11 out of 45) constituted a considerable number of substitution patterns in the same transcript.
Searching all assembled transcripts against multiple databases showed a considerable number of uncharacterised transcripts within both tissue types (Supplementary Materials Figures S5 and S6).
The analysis of the functional content of both transcriptomes in terms of Gene Ontology (GO) illustrated that GO terms associated with biochemical pathways involved in cellular components (cells, organelle and organelle parts, membrane, and extracellular regions), molecular functions (catalytic activity, binding, molecular function regulators, and molecular transducer activity), and biological processes (response to stimulus, cellular process, biological regulation, metabolic process, and biogenesis) dominated the functional diversity of both tissue types ( Figure 1; Figure 2).
The hepatic cells were more homogeneous in their gene ontology profiles across individuals. Only two metabolic pathways involved in the biological process, anatomical structure developments and multicellular organism developments, showed statistically significant differences in the number of associated polymorphic sites ( 2 test p-value < 0.05) across the four individuals ( Figure 3). The brain biopsies taken from different sections of the cerebral cortex demonstrated more heterogenous profiles in their polymorphism landscape. Twelve metabolic pathways involved in cellular components (neurone part, extracellular organelle-complex, and extracellular space), molecular functions (catalytic, hydrolase, transferase, and protein-binding activities), and a biological process linked to transposition showed a significant level of variation across different parts of the cerebral cortex ( 2 test p-value < 0.05) ( Figure 4).
Comparing the liver and brain datasets consisting of the ten most polymorphic sites against the closely related grey short-tailed opossum illustrated significant enrichment for various metabolic pathways ( Figure 5; Figure 6).
The dataset of 5% highly expressed transcripts illustrated that the expression profile between two tissue types differs significantly ( 2 test p-value < 0.05) in the metabolic pathways linked to 58 GO terms (Supplementary Materials Figure S7).
The metabolic pathways analysis of the transcripts against the KEGG database indicated that the liver transcriptome was linked to 394 core metabolic pathways and 23 complete metabolic modules. Brain cerebral cortex transcriptome was involved in 379 core metabolic pathways and eight complete metabolic modules (Supplementary Materials Tables S2 and S3).
In the liver tissue, the most polymorphic transcripts were involved in eight KEGG metabolic pathways including lipid metabolism, nucleotide metabolism, amino acid metabolism, cofactor and vitamin metabolism, xenobiotics biodegradation and metabolism, and the endocrine system. Mapping the ten most polymorphic transcripts in the brain cerebral cortex against the KEGG database indicated their involvement in 14 core biological pathways, including core metabolic pathways, biosynthesis of secondary metabolites, metabolism of cofactors and vitamins, signal transduction, signalling molecules and interaction, transport and catabolism, the cellular community, and the excretory system.     The dataset consisting of 5% of the most abundant transcripts of hepatic cells was linked to 354 biological pathways and eight complete metabolic modules ( Table 2). The highly expressed transcripts in the brain cerebral cortex were associated with 342 biological pathways and four complete metabolic modules (Table 3). A total of 22 genes in the liver transcriptomes and ninegenes in the cerebral cortex were directly linked to the metabolic pathways that involved drug and xenobiotic substances metabolism (Table 4; Table 5).    Table 3. Kyoto Encyclopaedia of Genes and Genomes (KEGG) metabolic pathway names and number of instances identified in the top 5% of the highly expressed genes within the brushtail possum cerebral cortex transcriptome.  Table 4. Name and description of the putative genes involved in drug and xenobiotic substance metabolic pathways in brushtail possum liver cells as predicted using the KEGG database.

Discussion
The potential of pests, pathogens, and invasive species to evade conventional population control measures represents a major threat across the globe [67]. Resistance against antibiotics and anthropogenic toxicants occurs in all domains of life, from bacteria to arthropods and vertebrates [68][69][70][71][72]. The genomic basis of developing resistance in non-rodent vertebrates continues to be underrepresented in research and suffers from a lack of functional genomic studies. The current study is the first to characterise the functional content and diversity of the liver and brain cerebral cortex transcriptomes of a non-rodent mammalian pest, the marsupial common brushtail possum, in New Zealand.
The functional content and polymorphism landscape of liver and brain cerebral cortex transcriptomes in brushtail possums are involved in a diverse array of biological processes. The functional content of these tissue types plays a critical part in the detoxification of xenobiotic substances [73,74] and modulation of complex behaviours [75,76] that potentially promote the development of physiological and behavioural tolerance against chemical compounds.
Enzymes produced from cytochrome P450 gene groups (such as CYP1A2, CYP2E1, CYP3A4, CYP3A5, CYP2B6, CYP2C8, and CYP2C19) dominate the biochemical pathways involved in drug and xenobiotic substance metabolism in brushtail possum hepatic cells. Polymorphisms in this group of heme-containing proteins influence the speed and efficiency of chemical compound breakdown in the liver and other organs [77,78]. Similarly, the CYP2C group of genes has been reported to show high levels of expression in liver tissues of the herbivorous koala (Phascolarctos cinereus) [79].
The detoxification of the harmful component of compound 1080, fluoroacetate, takes place in the hepatocytes [14,[80][81][82]. Polymorphism in the liver transcripts for the metabolic pathways involved in response to drugs and xenobiotic substances is not unexpected. Nonetheless, an exceptionally high number of polymorphic sites associated with these pathways, particularly a contig with high similarity to the multidrug resistance protein family reported from other marsupials, suggest the potential of standing genomic variation within New Zealand populations that may allow alternative responses to population control chemicals.
Multidrug resistance proteins confer increased drug resistance through a decrease in drug accumulations resulting from increased drug efflux [83,84]. Members of this ATP-binding cassette protein super-family function as active transporters of various organic substances, anionic conjugates, and xenobiotic substances, including glutathione, glucuronide, and sulfate conjugates [85,86]. Metabolic pathways that involve glutathione and a unique glutathione transferase are closely linked to the detoxification of compound 1080 and other similar-acting toxicants in hepatic cells [14,[80][81][82]. Moreover, the estimated ratio of synonymous to non-synonymous amino acid substitutions in this transcript indicates that parts of the multi-drug resistance metabolic pathways in brushtail possum hepatic cells may be subject to a significant level of purifying selection against less beneficial genotypes. The pre-adaptations of founder populations with natural exposure to this toxicant in Australia, and the fitness implication of the observed evolutionary dynamics, are yet to be fully understood and require comparative functional genomic studies between a larger sample of brushtail possums in their native range in Australia and invaded habitat in New Zealand.
In the current study, several enriched metabolic pathways were identified in the brushtail possum cerebral cortex transcriptome that regulate a wide range of core neurobiological processes, such as organic substance transport, response and regulation of biological processes, cell communications, neurochemical signalling, and locomotion. The biochemical pathways involved in drug and xenobiotic substance metabolism in the cerebral cortex are largely different from those identified in liver cells. The products of several genes in the brushtail possum cerebral cortex (Table 5) are associated with the development of resistance against platinum-based anti-malignancy drugs in humans. The complex dynamics of resistance involve a combination of decreased drug influx, increased drug efflux, intracellular detoxification by glutathione, decreased binding, increased DNA repair, decreased mismatch repair, and defective apoptosis [87,88]. However, the physiological significance of these metabolic pathways in the brushtail possum cerebral cortex in response to the chemicals remains to be investigated.
The delicate interplay between brain cerebral cortex polymorphisms and modulation of complex behaviours is an active area of research in neurobiology [89,90]. While addressing these associations lies beyond the scope of the current study, our results suggest a considerable variation in the transcriptomic repertoire of thebrushtail possum cerebral cortex that is necessary to modulate complex behaviours. Ross et al. [91] suggested that a single exposure to a sub-lethal dose of chemical toxicants in brushtail possums triggers the development of long-lasting bait-shyness behaviours in wild and captive populations. To what extent bait-shyness represents learned behaviour within an individual subsequent to sub-lethal exposure, and its heritability across multiple generations, requires more neurogenomic research.
Instances of developing tolerance against 1080 have already been reported in some wild and laboratory-populations with natural exposure to fluoroacetate [82], including wild populations of European rabbit (Oryctolagus cuniculus) in Australia [71], brushtail possums in south-western Australia [92], and laboratory brown rats (Rattus norvegicus) [93]. The current study did not aim to identify the development of practical resistance against toxins in New Zealand's brushtail possum population. Rather, we demonstrated that these populations harbour considerable levels of functional variation in the metabolic pathways that could potentially promote the development of tolerance against chemicals, such as compound 1080 and other similar-acting toxicants. To determine the sensitivity of the current populations to long-term toxin exposure requires formal toxicology trials to measure the lethal concentrations in the laboratory and over time. The power of genetic data to conclusively predict the development of tolerance against toxins remains limited. The small scale of this study only allows for some preliminary conclusions that are consistent with what we would expect to see from a species that has been under intense anthropogenic selection pressures, including both trapping and toxins, for generations.
As conservation efforts to eradicate invasive mammals intensify towards the objectives of the New Zealand Predator-Free 2050 horizon, we strongly advocate combined and multifaceted population control strategies that simultaneously target different aspects of an invasive species' life history, behaviour, and physiology in time and space, and thus prevent or delay the potential development of resistance against any specific population control measure [67,94]. The ecological and economic impacts of brushtail possums potentially developing tolerance against population control measures necessitate constant monitoring of current populations and their response to management plans.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/11/4/436/s1, Figures S1-S4: The histological images and reports of the four brain cerebral cortex biopsies, Table S1: The gene content of the liver and brain transcriptomes estimated in the BUSCO v3, Figures S5 and S6: The number of transcripts that was annotated by producing significant matches against different databases, Figure S7: The Gene Ontology terms of the top 5% highly expressed transcripts that differ significantly between the liver and brain cerebral cortex cells. Table S1. The gene content of the brushtail possum's liver and brain cerebral cortex transcriptomes compared to the core-mammalian BUSCOs. Tables S2 and S3: Definition and number of KEGG metabolic pathways in the liver and brain cerebral cortex transcriptomes. Tables S4 and S5: Liver and brain synonymous/non-synonymous amino acid substitution. The raw sequences generated for this study were submitted to the Sequence Read Archive (SRA) at the NCBI database under Bio Project submission numbers PRJNA623153 and PRJNA623060.