RNA sequencing of kidney and liver transcriptome obtained from wild cynomolgus macaque (Macaca fascicularis) originating from Peninsular Malaysia

Abstract

Objective

Using high-throughput RNA sequencing technology, this study aimed to sequence the transcriptome of kidney and liver tissues harvested from Peninsular Malaysia cynomolgus macaque (Macaca fascicularis). M. fascicularis are significant nonhuman primate models in the biomedical field, owing to the macaque’s biological similarities with humans. The additional transcriptomic dataset will supplement the previously described Peninsular Malaysia M. fascicularis transcriptomes obtained in a past endeavour.

Results

A total of 75,350,240 sequence reads were obtained via Hi-seq 2500 sequencing technology. A total of 5473 significant differentially expressed genes were called. Gene ontology functional categorisation showed that cellular process, catalytic activity, and cell part categories had the highest number of expressed genes, while the metabolic pathways category possessed the highest number of expressed genes in the KEGG pathway analysis. The additional sequence dataset will further enrich existing M. fascicularis transcriptome assemblies, and provide a dataset for further downstream studies.

Introduction

Cynomolgus macaques (Macaca fascicularis) are nonhuman primate (NHP) models significant to biomedicine due to their close evolutionary relationship with humans. The cynomolgus macaque’s recapitulation of human physiology, genetics, and behaviour is advantageous as translational models for various studies in the biomedical field, including drug development and safety testing [1, 2]. Cynomolgus macaque individuals from different geographical locations have been shown to exhibit genetic characteristics that varies between geographical origins [3]. Therefore, it is vital that the genomes and transcriptomes of cynomolgus macaque NHP models originating from different geographical locations are sequenced as a reference for future biomedical research design and implementation.

This study describes the RNA sequencing (RNA-seq) of two tissues—kidney and liver—harvested from wild Peninsular Malaysian cynomolgus macaques, and is an extension of a previous study [4] whereby the lymph node, spleen, and thymus transcriptomes of wild Peninsular Malaysian M. fascicularis were also sequenced with RNA-seq technology. An additional 75,350,240 sequence reads were obtained from this study, supplementing the previous wild Peninsular Malaysian M. fascicularis dataset for downstream applications. Furthermore, additional Malaysian cynomolgus macaque RNA-seq datasets will further furnish the cynomolgus macaque genome and transcriptome annotations and also provide a valuable asset for biomedical studies involving cynomolgus macaques.

Main text

Materials and methods

Samples

Two male conflict macaques that appeared to belong in the same family group were captured by the Department of Wildlife and National Parks (DWNP) from the state of Selangor. Visual examination showed that the macaques were free of disease. Euthanisation and harvesting of organs were also performed by authorised and qualified veterinarians of DWNP. In brief, the macaques were sedated intramuscularly using general anaesthesia (combination of ketamine, 5–10 mg/kg and xylazine, 0.2–0.4 mg/kg) before a lethal dosage of Dolethal^® were given intravenously. Harvested kidney and liver tissues were stored in RNAlater RNA Stabilization Agent (Qiagen). Total RNA extraction was carried out with RNeasy Mini Kit (Qiagen). A modified protocol was employed to remove genomic DNA contamination by incorporating Epicentre’s DNase I solution. RNA extract integrity was determined using 2100 Bioanalyzer (Agilent Technologies, Inc., United States of America) via RNA Pico chip. Samples with RIN > 7.0 were selected for library construction.

High-throughput sequencing, data analysis, and validation

Preparation of M. fascicularis kidney and liver RNA-seq libraries followed methods previously described in Ee Uli et al. [4]. Raw sequencing reads were submitted to NCBI Short Read Archive under accessions SRX2499144-SRX2499147. For data analyses, we used the same bioinformatic workflow and tools as depicted in our earlier publication [4]. Briefly, the assembly of raw sequence reads, normalisation, empirical analysis of differentially expressed genes (EDGE test), and gene ontology (GO) and KEGG pathway annotations were performed in CLC Genomics Workbench (CLCGW) version 7.5.1 (Qiagen Denmark), while functional classification of expressed genes in the kidney and liver transcriptomes was performed using Panther Database version 11.1 released on 2016-10-24 (http://pantherdb.org/) and DAVID Bioinformatics Resources version 6.8 (https://david-d.ncifcrf.gov/home.jsp). To better present the tissue-specific genes, we combined gene expression data from this study and our previous study [4] to create a Venn diagram representation of tissue-specific genes using the web tool Venn Diagrams (http://bioinformatics.psb.ugent.be/webtools/Venn). Validation of the RNA-seq dataset was performed with NanoString nCounter Elements XT (NanoString Technologies Inc., Seattle, WA, USA) gene expression assay. The methodology was the same as previously described in Ee Uli et al. [4]. Additional file 1 lists the genes utilised for the RNA-seq validation and their respective probe pair sequences.

Results

Sequence reads trimming and RNA-seq mapping

The number of raw kidney and liver sequence reads obtained were 44,452,083 and 30,898,157 respectively, each read possessing a sequence length of 70 bp. After trimming, the number of kidney and liver sequence reads were 41,854,334 and 28,943,585 singly. Summing up to a total of 70,788,919 good quality reads retained for assembly, with the average lengths of the reads ranging from 58.8 to 60.1 bp. Post-trimming, 93.94% of the reads were retained and were considered suitable for RNA-seq mapping. The overall mapping percentage of the sequence reads for each RNA-seq library ranged from 65.00 to 68.59%. Additional file 2 summarises the statistics of the sequence reads trimming and mapping.

Normalisation and differential gene expression analysis

A total of 5473 significant differentially expressed genes were called. Additional file 3 lists the differentially expressed genes identified from the EDGE test with their respective annotated gene names and descriptions. Within the Kidney vs. Liver tissue comparison, the top three most upregulated genes in liver were PLA2G2A, SAA1, and ORM1, while the three most downregulated genes include UMOD, TINAG, and DHDH.

Tissue-specific genes

Using gene expression data obtained from our previous study [4] and this study, a Venn diagram representation of tissue-specific genes in kidney, liver, lymph node, spleen, and thymus tissues was generated (Fig. 1). Lymph node showed the highest number of tissue-specific genes (2156 genes) among the five tissues compared; while kidney and liver had 310 and 154 tissue-specific genes respectively. The three most highly expressed tissue-specific genes in the kidney include MIOX, AQP2, and CDH16, while the top three tissue-specific genes expressed in the liver include APCS, APOA2, and LOC102121341. A complete list of tissue-specific genes can be found in Additional file 4.

Fig. 1

Venn diagram of the number of expressed genes (normalised expression value > 1) in their respective tissues

Functional annotation and classification

In the kidney, 14,001, 7046, and 5635 expressed genes were categorised to biological process, molecular function, and cellular component domains respectively, while 12,004, 6059, and 4935 expressed genes in the liver were categorised to BP, MF, and CC domains respectively. The distribution of the expressed genes into different level 2 GO categories are shown in Table 1. The BP category with the highest number of expressed genes for both tissues is cellular process (kidney, 4039; liver, 3496). For the MF domain, the catalytic activity category contains the highest number of expressed genes for both tissues (kidney, 2709; liver, 2419). As for CC, both tissues show the highest number of expressed genes in the cell part category (kidney, 2185; liver, 1954).

Table 1 Categorisation of expressed genes to gene ontology categories in kidney and liver tissues

In the kidney and liver, a total of 4406 and 3960 genes were assigned to KEGG pathway categories respectively. The metabolic pathways category contains the highest number of expressed genes for both kidney and liver tissues—903 and 839 genes respectively. The distribution of the expressed genes into different KEGG pathway categories are shown in Table 2.

Table 2 Categorisation of expressed genes to KEGG pathway categories in kidney and liver tissues

Validation of RNA-seq differential gene expression

Thirteen genes were selected to validate the RNA-seq differential gene expression dataset via NanoString gene expression assay. Normalised fold change values between the RNA-seq and NanoString platforms show concordance in magnitude directionality (Additional file 5). Our previous study also employed the NanoString gene expression assay to validate the RNA-seq dataset, and has shown similar concordance between the two platforms [4].

Discussion

In our previous endeavour [4], the lymph node, spleen, and thymus transcriptomes of the Malaysian cynomolgus macaque were sequenced via Illumina Hi-seq 2500 technology. The dataset presented in this study provides additional transcriptomic sequences from the kidney and liver of the Malaysian cynomolgus macaque, with a sum of 75,350,240 sequence reads obtained from four tissue libraries. In Ee Uli et al. [4], the transcriptome sequences were derived from immune-related tissues, whereas the kidney and liver transcriptomes generated in this study were obtained from metabolism-related tissues, and to a certain extent, these tissues also pertain to immune functions as well. Cynomolgus macaques are utilised as NHP models for liver-expressed drug metabolising enzyme studies, for example, the pharmacokinetics of the cytochrome P450 family of enzymes have been extensively studied in cynomolgus macaques [2, 5, 6].

It was observed that the percentages of the processed sequence reads that mapped back to the reference genome for all four libraries were below 70%. To rule out sequence contamination, we performed a BLAST alignment of the unmapped reads to NCBI’s nucleotide (nt) database and classified the hits into their specific taxonomic classification. The majority of the unmapped reads were assigned to Cercopithecidae and other primate family taxa, while the remaining reads were unassigned or did not have any BLASTN hits (Additional file 6). This suggests the absence of contaminating sequences, and also suggests the presence of novel transcripts in the kidney and liver dataset of the Malaysian cynomolgus macaque, which are yet to be described and curated in the M. fascicularis reference genome.

Comparisons were made between the kidney and liver tissues due to the similarities in metabolic function between the two tissues. However, the specific metabolic functions of both kidney and liver tissues are ultimately different. For instance, the KEGG pathway enrichment analysis highlights 58 liver genes that were specifically enriched in the biosynthesis of amino acids pathway, while 55 kidney genes were involved in the Inositol phosphate metabolism pathway.

Based on the differential gene expression analysis of kidney and liver tissues, the top three most upregulated genes in liver were PLA2G2A, SAA1, and ORM1, while the three most downregulated genes in kidney include UMOD, TINAG, and DHDH. The PLA2G2A gene codes for phospholipase A2 group IIA, and is suggested to be involved in the metabolism of phospholipids, inflammatory responses, and antimicrobial defence [7]. SAA1 codes for the acute-phase serum amyloid A1 protein, and is involved in inflammation response [8]. ORM1 codes for orosomucoid 1, which is suggested to induce the conversion of monocytes to M2b macrophages in response to type 2 cytokines in cancer patients [9]. UMOD codes for uromodulin, a urinary protein which inhibits the crystallisation of renal fluids that causes renal stone formation, and also expedites transepithelial migration of neutrophils in inflammation responses [10]. Tubulointerstitial nephritis antigen is a glycoprotein coded by TINAG, which is a nephritogenic antigen implicated in tubular homeostasis and cell survival in the kidney [11]. DHDH codes for the enzyme dihydrodiol dehydrogenase, and is suggested to be involved in detoxification processes in the kidney and liver [12].

In the GO analysis, both kidney and liver transcriptomes were observed to have the majority of genes represented by the same categories in their respective GO domains. In the BP domain, the highly represented category was cellular process, while the category with the highest gene representation in the MF domain was catalytic activity. As for CC domain, the highest number of expressed genes was represented in the cell part category. The representation of functional GO categories in the M. fascicularis transcriptome were observed to be similar to that of other chordates, including that of the silver carp [13], pig [14], black grouse [15], and pufferfish [16], and suggest a fair representation of the M. fascicularis transcriptome with that of other chordates.

The Malaysian cynomolgus macaque, M. f. fascicularis, is a unique subspecies that is suggested to have high nucleotide diversity and is genetically unique compared to other South East Asian M. fascicularis populations [17,18,19]. Genetic heterogeneity is suggested to contribute to varied responses in cynomolgus macaques towards drugs and pathogens [20], making research into the Malaysian population of cynomolgus macaque invaluable to the biomedicine field. The transcriptome dataset of the Peninsular Malaysian cynomolgus macaque obtained from this study can be compared with cynomolgus macaques originating from other geographical locations including Mauritius, the Philippines, Indochina, Indonesia, and China. Such comparisons can be utilised for the detection of single nucleotide polymorphisms in the macaque transcriptome, to infer phylogeny utilising a larger set of genic sequences, and also to determine the nucleotide diversity of the Peninsular Malaysian cynomolgus macaque. Furthermore, the transcriptomic data obtained from this study may be utilised to design genic molecular markers for rapid assessment of wild populations of cynomolgus macaques in Peninsular Malaysia in a relatively shorter period of time and at a lower cost.

Conclusion

Via the Illumina Hi-seq 2500 platform, 75,350,240 sequence reads in total were obtained from kidney and liver transcriptomes of the wild Peninsular Malaysian cynomolgus macaque. Additional sequence reads from the Malaysian cynomolgus macaque individuals will potentially further supplement the present genomic and transcriptomic data of this significant NHP primate model.

Limitations

• Only two biological replicates were available for library construction and data analyses.

• RNA sequencing of additional tissues will provide more robust transcriptomic dataset of wild Peninsular Malaysia M. fascicularis.