De Novo Transcriptome Assembly and Functional Annotation in Five Species of Bats

High-throughput RNA sequencing is a powerful tool that allows us to perform gene prediction and analyze tissue-specific overexpression of genes, but also at species level comparisons can be performed, although in a more restricted manner. In the present study complete liver transcriptomes of five tropical bat species were De novo assembled and annotated. Highly expressed genes in the five species were involved in glycolysis and lipid metabolism pathways. Cross-species differential expression analysis was conducted using single copy orthologues shared across the five species. Between 22 and 29 orthologs were upregulated for each species. We detected upregulated expression in Artibeus jamaicensis genes related to fructose metabolism pathway. Such findings can be correlated with A. jamaicensis dietary habits, as it was the unique frugivorous species included. This is the first report of transcriptome assembly by RNA-seq in these species, except for A. jamaicensis and as far as our knowledge is the first cross-species comparisons of transcriptomes and gene expression in tropical bats.

The order Chiroptera is the second largest order of mammals and is divided into: two suborders: Yinpterochiroptera and Yangochiroptera 1 . Its diversity includes and estimated ~1,331 species distributed throughout the world, except for the polar regions and isolated islands. Bats present a wide diversity of feeding habits and may be carnivorous, frugivorous, hematophagous, insectivorous or nectarivorous 2 ; as a consequence, chiropters play a crucial roles in the maintenance of the ecosystem balance by providing important ecological services; two-thirds of bats species are insectivorous and as such are considered biological pests controls of agricultural importance 2 . Similar to birds and insects, nectarivores act as pollinators and contribute to genetic exchange between flowering plants by pollination 2 . Frugivorous species play a crucial role in the renovation of disturbed landscapes and in the maintenance of vegetation within the ecosystem due to seed dispersion 2 . The success of bats in a wide diversity of niches is a result of evolutionary processes that have resulted in the appearance of diverse adaptations such as echolocation and ability to fly, which are conspicuous characteristics of bats.
Molecular studies using next generation sequencing have made it possible to analyse the whole genome and transcriptome of bats species and have contributed to a better understanding of these evolutionary processes as well as to a new taxonomic classification. For instance, molecular research has led to the rearrangement of Chiroptera phylogeny, bringing a new classification dividing the order into the two suborders Yinpterochiroptera and Yangochiroptera 1,3 based on the observation that the Microchiroptera constitute a paraphyletic clade instead of a monophyletic group as previously thought 4 .
High-throughput sequencing is not only useful in phylogenetics, but also helps to reveal evolutionary and adaptative evolution in bats, such as the correlation between the increment in energy metabolism demand with the evolution of true self-powered flight, which is conspicuous in bats 5 . Whole-genome and transcriptome sequencing studies have also generate new findings regarding the bat immune system response 6 . Transcriptome analysis in the Australian flying fox (P. alecto) 7,8 has provided insights into the co-evolution of bats and viruses by contributing to understanding the evolution of the chiropteran immune system and the apparent immunity or greater resistance of bats to viral infection in comparison with other mammals 9 . The study of the whole genome and transcriptome can unveil knowledge about evolution insights as well as ecological interactions in

Results and Discussion
RNA-seq library construction and sequencing. Total RNA was extracted from liver tissue obtained from three biological replicates of each of five bats species (Table 1). All specimens were collected from three locations at Yucatan State in the south-eastern Mexico (Fig. 1). The RNA integrity numbers (RIN values) considered for library preparation ranged from 6.5 and 8. The mean size of the fifteen cDNA libraries was 360 bp. HiSeq. 4000 multiplex sequencing generated a total of 403 million paired-end reads of 101 bp length. For each individual sequencing, the depth ranged between 30 to 69 million paired end reads (Table 1). After quality filtering with cutadapt 15 software, 99% of paired-end reads were conserved; these reads had final lengths between 30 and 101 bp and quality scores ≥ 30 (Table 1). A high proportion of the reads was retained after quality trimming, suggesting that although library construction was not performed using the traditional sample preparation kit manufactured by Illumina, but instead with the KAPA Stranded RNA-Seq with RiboErase kit from KapaBiosystems ® , we were able to obtain libraries of sufficient quality for an accurate sequencing, and therefore a high coverage of good-quality reads, suggesting that the use of an alternative kit for RNA-seq library construction in non-model organisms was successful.
De novo transcriptome assembly. Although several reference genomes and transcriptomes of bats are available, in the present study, transcripts were assembled by de novo instead of by genome reference mapping to decrease the possibility of losing novel or rare transcripts and species-specific sequences. De novo transcriptome assembly was performed with the Trinity 16 bioinformatics tool by applying two strategies: The first strategy involved transcript reconstruction for each biological replicate; this reconstruction yielded 15 assemblies in total (three assemblies per species). For the second strategy we merged the three replicates of each species, to obtain a more complete assembly, i.e., one assembly per species instead of one per individual. The concatenated assemblies showed higher N50 values and higher completeness. The subsequent analysis and results presented in this publication correspond to the concatenated assembly, statistics from each assembly are presented in Table 2.
M. keaysi showed the lowest number of assembled bases (413.9 Mb) and a total of 603,100 trinity transcripts, these numbers are higher than those reported for the assembled transcriptomes of Myotis ricketti 3,17 , which included 82,000-104,000 contigs. www.nature.com/scientificreports www.nature.com/scientificreports/ A. jamaicensis presented the highest number of assembled bases (682.4 Mb) and therefore a greater number of assembled transcripts (895,777), the number of contigs was significantly greater than in a previous study of A. jamaicensis that reported 214,707 contigs assembled from spleen with SOAPdenovo 12 .
According to N50 values, 50% of the assembled bases were incorporated in transcripts between 1,100 and 1,500 nucleotides in length. A preliminary assessment of de novo assembly quality was performed by obtaining the percentage of paired reads represented in the assembly. For A. jamaicensis, M. megalophylla, M. keaysi, N. laticaudatus and P. macrotis, the overall alignment rates were 91.82, 92.96, 93.76, 92.93 and 91.67% respectively. We interpreted these results as a first indication of a good quality assembly, based on the Trinity protocol 16 , for which it is established that 80% of read mapping should be considered indicative of good quality.

Transcripts filtering.
To reduce the number of potential spurious contigs, we filtered them based on the abundance and quality of the assembly. For abundance, we took into consideration isoforms expression level and  www.nature.com/scientificreports www.nature.com/scientificreports/ redundancy. Expression levels were detected using the RSEM method, and only highly expressed isoforms were retained for further analysis. After this first filtering step between 76 and 81% of the assembled transcripts were conserved (Table 3).
Redundancy was eliminated by clustering our de novo assembled transcripts with highly similar contigs using CD-HIT-EST at a nucleotide identity of 95%. Fewer than 3% of the Trinity transcripts were redundant and were therefore removed (Table 3).
Regarding quality, we retrieved bad assembled contigs such as chimaeras and incomplete contigs using TransRate software. The assemblies with the highest proportion of good quality-contigs were those obtained for A. jamaicensis and M. megalophylla, which contained less than 2% of poor-quality contigs. On the other hand, M. keaysi contigs were depleted drastically after quality contig assessment, with 40% of poor-quality contigs that were removed for this study ( Assembly completeness. The results for quantitative measures for determining assembly completeness with BUSCO 18 (Benchmarking Universal Single Copy Orthologs), showed the percentage of conserved orthologues among vertebrates, mammals and the Laurasiatheria superorder, represented in our assembled bat transcriptomes.
Completeness of the non-redundant contigs containing all the isoforms resulted in a high percentage of complete vertebrate orthologues (74 to 82%), but also a significantly higher percentage of putative paralogues, i.e., complete genes with more than one copy (Fig. 2a). Although these results might be considered as high for genome data, for transcriptomic data they seem logical due to the presence of multiple isoforms. To probe this, we performed BUSCO analysis of the strictly filtered transcripts, excluding weakly expressed transcripts and redundant contigs; in comparing these results, we observed that the percentage of duplicates decreased markedly (Fig. 2b). Of the 2,586 orthologues searched in the BUSCO set of vertebrates, between 60 and 65% were recovered completely; of the latter, less than 1.2% were putative paralogues, i.e., duplicates) 18 . Of the 4,104 single-copy orthologues in mammals, 56 to 59% were complete, and 12 to 19% were partially recovered or fragmented. Finally, approximately 49% of the 6,253 laurasiatherian orthologues were complete, and 12 to 19% were fragmented (Fig. 2b). The observed recovery of more than 60 percentage of the complete single-copy orthologues from vertebrates and 50% of those from the mammalian and laurasiatherian databases is indicative of good coverage and of high recovery of conserved orthologues for the five generated non-redundant transcripts.
These results lead us to the conclusion that additional filtering steps such as debug redundancy and weakly expressed isoforms are required to increase the percent recovery of single-copy orthologues percentage and reduce the presence of putative paralogues in the assembly if required.
By comparing the completeness of our assemblies against the published transcriptomes 10,11,17,19 of other bat species, we found that we were able to recover a higher percentage of complete vertebrate orthologues except in the case of the Rousettus aegyptiacus 10 transcriptome, which presented 93% of complete orthologues. Myotis ricketti 17 was the transcriptome with the lowest percentage of recovery, with only 51%, while in the case of the Myotis keaysi transcriptome we were able to recover 63% of the vertebrate orthologues, suggesting that the transcriptome generated in this study can be considered as a good quality reference for Myotis genus. For mammals and laurasiatherian orthologues, the characteristics of the database were similar (Supplementary Table S1).

Coding regions identification and functional annotation. Coding domain sequences (cds) and
annotation were identified with TransDecoder for each species after two filtering steps. We first selected the single best open reading frame (ORF) per transcript, then, only transcripts more than 200 bp in length were retained. The total number of protein coding transcripts among the final non-redundant transcripts was 57,919 for A.jamaicensis, 35,289 for M. megalophylla, 29,461 for M. keaysi, 46,116 for N. laticaudatus and 41,152 for P. macrotis. Between 38 and 45% of the coding sequences were complete for all five species, approximately 30% were 5′ partial, 9% were 3′ partial, and 15 to 23% were internal (Fig. 3). Protein-coding transcripts were compared against the Uniprot reference database using a BLAST search with an e-value cut-off of 10 −3 . For the five species ~60% of the coding transcripts matched with this search criteria, and between 80% and 85% of these showed full-length or nearly full-length recovery, respectively. Following the Trinotate script and GOseq pipeline, Gene Ontology terms were assigned to the assembled transcriptomes. GOseq results are classified in three categories:  www.nature.com/scientificreports www.nature.com/scientificreports/   open reading frames in each species identified using the TransDecoder pipeline. "Complete ORF" refers to sequences in which the first codon and the stop codon are present. "5 prime" refers to sequences that contain the start codon but lack the stop codon; "3 primer" refers to partial ORF sequences that lack the start codon. Finally, "internal ORF" refers to sequences that lack both the start and stop codons. (2019) 9:6222 | https://doi.org/10.1038/s41598-019-42560-9 www.nature.com/scientificreports www.nature.com/scientificreports/ Myotis lucifugus presented the highest percentage of homology against all query species; as expected, the species with greatest homology was Myotis keaysi with 68%, whereas A. jamaicensis, M. megalophylla, N. laticaudatus and P. macrotis displayed only 25-31% homology (Fig. 5). In contrast, the de novo assembled transcriptomes had lower percentages of homologous transcripts (<18%) against the Yinpterochiroptera database, including Hipposideros armiger (F. Rhinolophidae), which was formerly classified within the suborder Microchiroptera; according to the recently proposed molecular phylogeny 1,3 , H. armiger was reclassified within the suborder Yinpterochiroptera along with the family Pteropodidae.
Detection of putative orthologues and orthology grouping of proteins was performed with Orthofinder 20 for the five transcriptomes using the BLAST all-v-all algorithm. A total of 209,937 transcripts were assigned to 19,205 orthogroups. An orthogroup is defined as a set of genes descended from a single gene of the last common ancestor within species groups 20 . Among the five species, 5,848 (36.4%) of the orthogroups were shared (Fig. 6a); this will be discussed further in the section on transcript abundance. A total of 72 genes were classified as species-specific, these were grouped in 15 inferred orthogroups, of which only 10 were annotated. Three orthogroups classified as species-specific from M. keaysi caught our special attention because these groups, according to our functional annotation, corresponded to the retroviral genes GAG, POL and ENV, respectively. The Gag p24 gene protein forms the inner protein layer of the nucleocapsid during virus replication, specifically during the assembly, maturation and infection stages of retroviruses. These findings might be a first indication that an active viral cycle was occurring at the moment of capture of the M. keaysi individuals used in our study. However, it is retroviral genetic material can be inserted into the host genome (endogenous retroviruses), and despite the assumption that these endo-retroviruses are not functional, in some cases they have been shown to be expressed. Further research is needed to determine whether the observed expression is due to an exogenous retrovirus or to a functional endogenous retrovirus.
In addition to detecting orthologues, Orthofinder infers a species tree based on single-copy orthogroups 20,21 . The tree was rooted with Yinpterochiroptera species as an outgroup (Fig. 6b). In the orthologues phylogram, a clear separation between the suborders Yinpterochiroptera and Yangochiroptera can be observed. As expected,  www.nature.com/scientificreports www.nature.com/scientificreports/ the query species were classified within Yangochiroptera; Myotis keaysi constitutes a monophyletic group with Eptesicus fuscus and Myotis lucifugus, and all three belong to the family Vespertillionidae (Fig. 6b). The four remaining species constituted one clade with two clear divisions. Artibeus jamaicensis (family Phyllostomidae) and Mormoops megalophylla (family Mormoopidae) formed the first clade, consistent with recent molecular classification that groups them within the superfamily Noctilionoidea 1,3 . However, this was not the case for Peropteryx macrotis (family Emballonuridae) and Nyctinomops laticaudatus (family Molossidae); this clade is inconsistent with recent molecular phylogeny 1 but agrees with two different phylogenetic reconstructions, one that was inferred from the sequences of the cytochrome B (cytb) molecular marker 3,22 and another that was inferred from transcriptomic data 18 . It must be emphasized that this phylogram was not constructed under an evolutionary model but is only based on the substitution rates of single-copy orthologues.

Transcript Abundance and Differential Expression Analysis. For estimation of the level of expression
in each species, we first performed transcript quantification for each biological replicate and then used the ExN50 statistic to retrieve expression levels of the isoform level in each species. To reduce the number of comparisons, we focused only on transcripts that were shared among the five species.
These values were obtained from the ExN50 statistics calculated from theTMM matrix. For A. jamaicensis, the orthologous transcript with the highest expression values corresponded to the gene ALDOB which is expressed in the liver and encodes the aldolase B enzyme, a member of the glycolytic protein family (Pfam PFF00274). The fructose-biphosphate aldolase type B glycolytic enzyme catalyses the condensation of fructose 1-6-bisphosphate or fructose-1-phosphate 22,23 ; the high level of expression of this enzyme in A. jamaicensis is clearly correlated with the species' feeding habits as a frugivorous bat, as fructose when absorbed is metabolized in the liver by ALDOB, producing substrate for fatty acid synthesis, all this process takes place in the liver 24 . On the other hand, for the other four insectivorous bats, the most highly expressed transcripts were annotated as serum albumin protein, which is the main component of blood plasma and is related to fatty acid binding, metabolites, hormones and bilirubin. This protein was also highly expressed in A. jamaicensis. Other transcripts that were highly expressed in the five species corresponded to the apolipoprotein E gene (APOE). APOE is a plasma protein that is secreted primarily in hepatic tissues with a role in lipid transportation, and is also related to innate and adaptative immune response 23,25 . Now then, making a comparison between transcripts of different species with high evolutionary distances, is challenging for three main factors: (1) because we are studying wild, non-model organisms it is difficult to establish a control group that can be used to perform traditional differential expression analysis; (2) all individuals were captured under aleatory circumstances, making it impossible to know their feeding and metabolic status at the time of capture and making a comparison related to liver metabolic degradation of nutriments is not possible; (3) the lack of reference genomes for the studied species is an important obstacle to conducting a proper and reliable comparison of cross-species transcripts expression, due to the potential loss of transcript information. The ultimate (but not available) solution would be the comparison of a reference genome of at least a closely related species of each target study object; herein lies the importance of generating new available genomic and transcriptomic resources for non-model species such as bats. Therefore, to ensure that our data are as uniform as possible, we proposed, as an alternative solution, to conduct the DE analysis using only in the orthologous transcripts that were shared among the five species (Fig. 7). We retrieved 3,052 single-copy orthologues from the output of Orthofinder analysis; of these, only approximately 1,000 genes were considered differentially expressed according to the four statistical methods employed (EdgeR, DeSeq, DEseq. 2 and NOIseq). Volcano plots of each of the ten comparisons can be found in Supplementary Fig. S2. A total of 27 genes were up regulated in Artibeus jamaicensis against the other four species (Supplementary Table S3). Transcript that encodes for Myo-inositol oxygenase enzyme (MIOX) was the one with higher level of expression among the A. jamaicensis single-copy orthologues and was significantly downregulated (p < 0.01) in the other four species (Fig. 7). MIOX enzyme is www.nature.com/scientificreports www.nature.com/scientificreports/ involved the first step of myo-inositol (MI) degradation into D-glucuronate, which is essential for the pentose phosphate cycle 26 . Upregulation of MIOX might be related to A. jamaicensis diet, as MI has been found in high concentrations in fruits and seeds 26,27 . In humans this enzyme is restricted to kidney expression 22 , contrary as expected, transcript abundance in the three biological replicates of A. jamaicensis indicates an overexpression of a MIOX like transcript in liver, further research has to be done to confirm that liver MIOX expression actually occurs in other frugivorous bat species.
Glucose-6-phoshphate dehydrogenase (G6PD) enzyme was also upregulated in A. jamaicensis and downregulated in the studied insectivorous species. It has been found that G6PD upregulation is related with insulin resistance metabolic syndrome 28 , which has been previously reported in organisms with high fructose levels of consumption 29 . However, more studies need to be done for probe insulin resistance in frugivorous bats as a result of an evolutionary adaptation related to feeding habits.
The biological process Gene Ontology terms involved within the upregulated genes in A. jamaicensis according to GO enrichment analysis were metabolic process (GO:008152), regulation of phosphate metabolic process (GO:0019220) and immune system (GO: 0071840), among others. For molecular function catalytic activities such as transferase activity and oxidoreductase activity were involved in the upregulated transcripts. The overrepresented pathways corresponded to the pentose phosphate pathway (P02762) and pyruvate metabolism (P02772). Most of the upregulated transcripts in A. jamaicensis are correlated with species feeding habits, in Fig. 7 it is appreciated how Artibeus replicates are grouped in an independent cluster separated from the insectivorous species.

Conclusions
In this study, we generated whole transcriptomes from the liver tissue of five species of tropical bats classified into five different families: A. jamaicensis (F. Phyllostomidae), Mormoops megalophylla (F. Mormoopidae), Myotis keaysi (F. Vespertilionidae), Nyctinomops laticaudatus (F. Molossidae) and Peropteryx macrotis (F. Emballonuridae). Also, a comparative analysis was performed to obtain overexpressed genes of each species, this comparison was performed exclusively on single-copy orthologues shared across the five species, with this cross-species differential gene expression analysis we were able to identify species-specific upregulated genes related with dietary habits in bats, however this analysis was restricted to shared orthologous transcripts due to the lack of a closest reference genome or transcriptome database. Here relies the importance of generating new genomic database for non-model species, especially for wild species with high ecological and evolutionary importance such as bats. Handling and sacrifice methods applied to the specimens were approved at the document CB-CCBA-I-2017-006, which was issued by the Bioethics Committee of the Faculty of Veterinary Medicine and Zootechnics at Autonomous National University of Yucatan (Universidad Nacional Autónoma de Yucatán, UADY).

Methods
sampling collection. In this study we collected three biological replicates from each of five bat species belonging to five families (Table 1). Bats were collected from various locations in Yucatan State in south-eastern Mexico (Fig. 1). A. jamaicensis, M. megalophylla and N. laticaudatus were collected at Calcehtok cave; free-ranging bats were captured using two mist-nets located at the cave entrance. M. keaysi individuals were captured inside Hoctún cave using mist-nets, and P. macrotis were collected in their roosting at Hobonil cave using a sweep net.
The specimens were transported alive to the Arbovirology Laboratory at the University of Yucatan (Universidad Autónoma de Yucatán, UADY). There, they were sacrificed by lethal cardiac puncture and dis- De novo transcriptome assembly. The raw read quality of each paired-end library was examined using the bioinformatics tool FastQC v 0.11.5 31 . Adapters were removed using cutadapt v. 1.12 15 for paired-end reads (R1 and R2), in addition poly A/T tails, ambiguities (N), sites with PHRED scores lower than 20 and reads below 30 bp in length were removed. For this purpose, the following parameters were used: "-a <adapter sequence> -A<adapter sequence>-o<output_read_1.fastq> -p<output_read_2.fastq> -q 20 -b "A{101}" -B "T{101}" -trim-n -minimum-length 30 -o<output_read_1.fastq> -p<output_read_2.fastq>"; for adapter trimming we added an "A" at the beginning of each adapter sequence. The edited reads were re-examined on FastQC v 0.11.5 31 to verify their final quality.
Transcriptome de novo assembly was performed for each species using Trinity v.2.4.0 16 with previous normalization of the edited reads. The following parameters were used: "-seqType fq-SS_lib_type RF-left<input_ file>-right<input_file> -CPU 62-max_memory 400 G". Assembly statistics were computed using the script TrinityStats.pl contained in the Trinity 16 package. The proportion of reads mapped to the assembly was assessed with Bowtie2 32 .
Assembly filtering. To reduce the probability of obtaining of spurious transcripts and attenuate transcript redundancy, the contigs were filtered using three methods: First, weakly expressed isoforms were removed based on their expression values 16 . TPM values were obtained by the RSEM 33 method using Trinity script aling_and_ estimate_abundance.pl; then, weakly expressed isoforms were removed using the Trinity script filter_low_expr_ transcripts.pl with "-highest_iso_only" parameter. Second, a set of non-redundant representative transcripts was generated using the CD-Hit 34 package with an identity threshold of 95%. The parameters used were "cd-hit-est -c 0.95 -n 10 -M 60000 -T 8".
Finally, bad contigs, i.e., misassembled or incomplete contigs, were filtered out from non-redundant assemblies based on read mapping metrics. For this purpose, we used the bioinformatics tool TransRate 35 , which is designed for analysis of the quality of de novo transcriptome assemblies. Downstream analysis was performed on the final filtered contigs. ( Homology search and Gene orthology prediction. Analysis of homology between de novo assembled transcripts and the bat refseq database was performed by pairwise comparison using BLAST to provide full-length transcript analysis. Bat databases were constructed with the option -makeblastdb included in the BLAST software. For the BLASTn and BLASTp analysis, we built a database of the refseq genome sequences and coding DNA sequences from seven species: Rousettus aegyptiacus (BioProject: PRJNA309421), Pteropus vampyrus (BioProject: PRJNA275879), Pteropus alecto (BioProject: PRJNA232518), Eptesicus fuscus (BioProject: PRJNA232522), Hipposideros armiger (BioProject: PRJNA357596) and Myotis lucifugus (BioProject: PRJNA208947).
Orthologous groups (orthogroups) of protein sequences amongst the five species were identified with the Orthofinder v.2.1.2 20 bioinformatics tool, using default parameters. Using reciprocal best-hits by the BLAST all-v-all algorithm, Orthofinder determined the number of shared putative orthologues between the five species as well as species-specific transcripts. Graphic representation for Orthofinder output was performed using the VennDiagram 38

Functional Annotation. Candidate coding regions with a minimum cut-off of 200 amino acids and open
reading frames (OFRs) were predicted with TransDecoder 16 pipeline. Once the candidate peptides we obtained, we performed homology searches against known proteins using BLASTp using UniProtKB/Swiss-Prot database and against common protein domains using the Pfam 43 database. The BLAST and Pfam search outputs were integrated into coding regions prediction. Functional annotation was conducted using the Trinotate 16 utility (http:// trinotate.sourceforge.net/), in which homology search was performed with BLASTx, whilst for predicted coding region homology search we used BLASTp. Additionally, protein domains were identified with HMMER v.3.1 44 against the Pfam database. Signal peptide predictions was performed using signalP v.4 45 . Transmembrane regions were predicted using the tmHMM v.2 46 server and ribosomal RNA genes were detected with RNAMMER v.1.2 47 . Annotation outputs were loaded into a Trinotate SQLite Database.
Transcript Abundance and Differential Expression Analysis. To quantify transcript abundance we applied the alignment-based methods contained in the Trinity 16 package, by mapping the reads of each biological replicate against the respective assembled transcriptome. This was obtained with the align_and_estimate_abundance Perl script. In this analysis, we used RSEM 33 as the abundance estimation method and chose bowtie2 48 for the alignment. When the transcript abundance for each biological replicate had been obtained, we built a Gene Expression Matrix using the abundance_estimates_to_matrix.pl script to generate a normalized expression values matrix that was used to obtain the expression level of each transcript by ExN50 analysis.
Differential expression analysis (DE) was performed only on transcripts for which orthologues were present in the five species; such sequences were obtained from Orthofinder analysis. For DE, quantification of orthologous transcripts was performed using the Salmon 49 tool, a quasi-map index was built with each species assembled transcriptome, using a k value of 29, as the shortest length of the reads was 30 bp. As with RSEM, we aligned each biological replicate to its transcriptome. Each quantification file was edited by replacing the transcript ID generated by Trinity, for its respective Single Gene Orthogroup name. This was necessary because each assembly was generated independently and without any reference, therefore the Trinity ID headers were assigned randomly to each species. By assigning the orthologue ID we propose that a more accurate differential expression analysis in non-related species can be performed.
The edited quantification files generated by Salmon software were imported to R using the tximport 50 package contained in the Bioconductor 51 library. The counts matrix filed imported was submitted to the web tool IDEAmex 52 (Integrated Differential Expression Analysis MultiEXperiment, http://zazil.ibt.unam.mx/idea) which performed differential expression analysis, taking into consideration the biological replicates with four methods: EdgeR, NOISeq, DESeq and DESeq. 2. In this analysis, we used only genes that were differentially expressed according to the four methods. Finally, GeneOntology (GO) and KEGG enrichment were obtained from the PANTHER web tool, using Sus scrofa as a reference.