Comprehensive annotation of the Chinese tree shrew genome by large-scale RNA sequencing and long-read isoform sequencing.

The Chinese tree shrew (Tupaia belangeri chinensis) is emerging as an important experimental animal in multiple fields of biomedical research. Comprehensive reference genome annotation for both mRNA and long non-coding RNA (lncRNA) is crucial for developing animal models using this species. In the current study, we collected a total of 234 high-quality RNA sequencing (RNA-seq) datasets and two long-read isoform sequencing (ISO-seq) datasets and improved the annotation of our previously assembled high-quality chromosome-level tree shrew genome. We obtained a total of 3 514 newly annotated coding genes and 50 576 lncRNA genes. We also characterized the tissue-specific expression patterns and alternative splicing patterns of mRNAs and lncRNAs and mapped the orthologous relationships among 11 mammalian species using the current annotated genome. We identified 144 tree shrew-specific gene families, including interleukin 6 (IL6) and STT3 oligosaccharyltransferase complex catalytic subunit B (STT3B), which underwent significant changes in size. Comparison of the overall expression patterns in tissues and pathways across four species (human, rhesus monkey, tree shrew, and mouse) indicated that tree shrews are more similar to primates than to mice at the tissue-transcriptome level. Notably, the newly annotated purine rich element binding protein A (PURA) gene and the STT3B gene family showed dysregulation upon viral infection. The updated version of the tree shrew genome annotation (KIZ version 3: TS_3.0) is available at http://www.treeshrewdb.org and provides an essential reference for basic and biomedical studies using tree shrew animal models.

INTRODUCTION

Suitable animal models are essential for expanding our knowledge regarding fundamental biological questions and for developing new drugs, vaccines, and therapeutics (McGonigle & Ruggeri, 2014; Robinson et al., 2019; Yao et al., 2015). An ideal animal model should possess many features, including high genetic similarity to humans, similar pathobiology and symptoms, efficacy with drug prediction and response, low cost, and low restriction (Bennett & Panicker, 2016; McGonigle & Ruggeri, 2014; Robinson et al., 2019; Yao et al., 2015). The Chinese tree shrew (Tupaia belangeri chinensis) is a small rat-sized (100–150 g) mammal with a short reproductive cycle (~6 weeks) (Yao, 2017; Zheng et al., 2014), and is widely distributed in Southeast Asia and South and Southwest China. In the past few decades, tree shrews have been widely used in a variety of biomedical studies, including research on viral infections (Amako et al., 2010; Li et al., 2018; Xu et al., 2007, 2020c), cancer (Ge et al., 2016; Lu et al., 2021), myopia (He et al., 2014; Levy et al., 2018; Phillips et al., 2000), visual cortex function (Fitzpatrick, 1996; Lee et al., 2016; Petry & Bickford, 2019), and neuroscience (Dimanico et al., 2021; Fan et al., 2018; Ni et al., 2018; Savier et al., 2021; Wei et al., 2017). Our research on tree shrews began with the genetic dissection of the Chinese tree shrew genome (Fan et al., 2013). We also aimed to promote the use of this animal in basic and biomedical research by continuing to update relevant genome information (Fan et al., 2014, 2019). Moreover, we developed two immortalized tree shrew cell lines for resource sharing (Gu et al., 2019b; Zhang et al., 2020b) and established the first genetic manipulation of tree shrews using spermatogonial stem cells to successfully generate transgenic offspring (Li et al., 2017). Compared with commonly used animal models such as rodents, tree shrews are phylogenetically closer to primates (Fan et al., 2013, 2019), and can more accurately mimic the physiological and pathological conditions of humans.

Accurate genome assembly and annotation are crucial for understanding tree shrew biology and for developing disease models using this animal. Indeed, creating an animal model of human disease using a tree shrew genome-based approach (Yao, 2017) is dependent on high-quality annotations of the tree shrew genome. Many attempts have been made to decipher the tree shrew genome in great detail and accuracy (Fan et al., 2013, 2019; Sanada et al., 2019). We successfully assembled the first high-quality genome of the Chinese tree shrew (KIZ version 1: TS_1.0) using high depth (~79X) short-read sequencing technology (Fan et al., 2013) and the first chromosome-level tree shrew genome (KIZ version 2: TS_2.0) using single-molecule real-time (SMRT) sequencing technology (Fan et al., 2019). The release of two versions of the tree shrew genome at www.treeshrewdb.org (Fan et al., 2014, 2019) has undoubtedly enhanced our knowledge on the usage of this species. Recently, Sanada et al. (2019) assembled a tree shrew genome using short reads for coding sequence (CDS) annotation. However, despite efforts to improve the annotation of the tree shrew genome, our understanding of the coding and non-coding genes of the tree shrew remains incomplete and unlikely to meet the growing needs of the research field.

RNA sequencing (RNA-seq) technology provides accurate and massive amounts of information regarding the direct transcription status of a genome (Cardoso-Moreira et al., 2019; Stark et al., 2019). The emergence of third-generation sequencing, which features long sequence reads that can cover the full-length of most transcripts (Gordon et al., 2016; Sharon et al., 2013), has greatly improved the accuracy of transcript structure annotation. Previous studies using next-generation sequencing (NGS) based on RNA-seq and long-read isoform sequencing (ISO-seq) have revealed the complexity and characteristics of eukaryotic transcriptomes (Chen et al., 2017). Using ortholog and de novo annotations (Garber et al., 2011; Yandell & Ence, 2012), transcriptome sequencing has been widely used to annotate the genomes of plants (Purugganan & Jackson, 2021; Wang et al., 2019), model animals (Ji et al., 2020; Nudelman et al., 2018; Zhang et al., 2020a), and livestock (Beiki et al., 2019; Foissac et al., 2019). In this study, we aimed to provide a more comprehensive tree shrew genome annotation using a wide range of transcriptome sequencing data. We collected high-quality RNA-seq datasets of tree shrews from publicly available sources (Supplementary Table S1), as well as two ISO-seq datasets and 139 RNA-seq datasets newly generated in this study. These transcriptome datasets included expression data of tree shrew cells and tissues under different conditions, including viral infection (Sanada et al., 2019; Yan et al., 2012), normal tissue (Fan et al., 2013; Han et al., 2020), and pathological tissue (Li et al., 2017; Lin et al., 2014; Tu et al., 2019; Wu et al., 2016b; Zhang et al., 2020b). Using a stringent pipeline, we obtained a total of 53 298 newly annotated coding transcripts and 115 562 newly annotated non-coding transcripts and produced a relatively complete and reliable tree shrew genome annotation (KIZ version 3: TS_3.0). Based on this comprehensive annotation, we further explored the spatial expression and alternative splicing patterns of the tree shrew transcripts and characterized the orthologous relationships among tree shrews and other species. We also compared expression similarity across species and provided a landscape of the innate immune response in tree shrews upon viral infection.

MATERIALS AND METHODS

Animals and tissue collection

Nine adult Chinese tree shrews were purchased from the Experimental Animal Center of the Kunming Institute of Zoology, Chinese Academy of Sciences. Animals were anesthetized with pentobarbital and intracardially perfused with phosphate-buffered saline (PBS). Eight tissues (small intestine, liver, heart, kidney, spleen, ovary, brain, and testis) from four animals were collected and snap-frozen in liquid nitrogen. The remaining animals were used for the isolation of tree shrew primary renal cells (TSPRCs) for viral infection assays. All animal experiments were approved by the Institutional Review Board of the Kunming Institute of Zoology, Chinese Academy of Sciences.

ISO-seq for tree shrew tissues

Tissues from two adult Chinese tree shrews were used for ISO-seq (Supplementary Table S2). RNA extraction, library construction, and sequencing were performed by Annoroad Gene Technology (China). In brief, total RNA from each sample was isolated using a NEBNext® UltraTM RNA Library Prep Kit for Illumina® (Catalog # 7530; New England Biolabs Inc., USA) and processed following the manufacturer’s protocols. RNA degradation and contamination were monitored by 1% agarose gels and RNA purity was checked using a NanoPhotometer® spectrophotometer (IMPLEN, USA). RNA integrity was assessed using a Qubit® RNA Assay Kit with a Qubit® 2.0 Fluorometer (Life Technologies, USA) and an RNA Nano 6000 Assay Kit with the Bioanalyzer 2100 system (Agilent Technologies, USA). Equal amounts of RNA from each tissue of the two tree shrews were pooled as one mixed RNA sample for ISO-seq. Two ISO-seq libraries (<4 kb and >4 kb) were prepared according to the Isoform Sequencing protocol (ISO-seqTM) using a Clontech SMARTer PCR cDNA Synthesis Kit and the BluePippinTM Size Selection System (Sage Science, USA) protocol as described by PacBio (Menlo Park, USA). SMRT sequencing was performed on the Pacific Bioscience Sequel System using two SMRT cells.

The ISO-seq data were processed using IsoSeq v3.4.0 (https://github.com/PacificBiosciences/IsoSeq). Only sequence reads containing both 5' and 3' adaptors were retained to cover the entire transcript. We used LoRDEC (Salmela & Rivals, 2014) to correct errors in the SMRT reads by referring to the RNA-seq data. Subsequently, the corrected SMRT reads were aligned to the tree shrew reference genome TS_2.0 (Fan et al., 2019) using GMAP (Wu et al., 2016a) to locate the position of the predicted genes on the pseudochromosomes.

Compilation of publicly available tree shrew RNA-seq datasets

To ensure a robust and complete annotation of the tree shrew genome, we obtained all publicly available RNA-seq datasets of tree shrews from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/), DNA Data Bank of Japan (DDBJ, https://www.ddbj.nig.ac.jp/index-e.html), and China National Center for Bioinformation (BIGD, https://bigd.big.ac.cn/). These transcriptome sequencing datasets were originally obtained by sequencing normal and pathological tissues or cells and represent a wide spectrum of biological and pathological conditions (Supplementary Table S1).

We used the following strategy for quality control (QC) of the RNA-seq data and filtered those data that did not meet requirements. Briefly, raw sequencing reads were processed by Trimmomatic (v0.38) (Bolger et al., 2014) to trim adaptor and low-quality sequences, with the parameters “LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36”. After read filtering, the quality of the clean reads was assessed by FastQC (https://sourceforge.net/projects/fastqc.mirror/). Datasets that passed QC (Q30>20) were aligned to the chromosome-level tree shrew genome (KIZ version 2: TS_2.0; https://www.treeshrewdb.org/download) using STAR (v2.6.0c) (Dobin et al., 2013). We discarded those datasets that failed to pass QC, i.e., mapping ratio to genome below 75%. After QC, a total of 91 publicly available RNA-seq datasets were retained for analysis (Supplementary Table S1).

RNA-seq of tree shrew tissues and cells with or without viral infection

To enhance the credibility of the genome annotation, we generated new RNA-seq data from tissues and/or cells of nine tree shrews (Supplementary Table S2) for further analysis with the publicly available RNA-seq data.

To explore the changes in gene expression in tree shrew cells in response to viral infection, we performed RNA-seq of virus-infected TSPRCs. We used the same procedure to isolate and culture TSPRCs and challenge cells with or without virus as described in our previous studies (Xu et al., 2016, 2020a; Yu et al., 2014). Briefly, TSPRCs were infected with a DNA virus (herpes simplex virus type 1, HSV-1; multiplicity of infection (MOI)=10) and RNA viruses (Sendai virus (SeV, 20 hemagglutinating units/mL), encephalomyocarditis virus (EMCV, MOI=2), and Newcastle disease virus (NDV, MOI=1)), respectively, for the indicated times before harvesting for RNA-seq (Supplementary Table S3). RNA-seq of tree shrew tissues and infected cells was performed by Annoroad Gene Technology (China). Approximately 1 μg of total RNA from each sample was used to construct the RNA-seq libraries (500–1 000 bp) with a NEBNext® UltraTM RNA Library Prep Kit for Illumina®. The quality of each library was assessed using the Agilent Bioanalyzer 2100 system. Libraries were sequenced on the Illumina NovaSeq platform and 150 bp paired-end reads were generated. We followed the same QC procedures as described above for processing the publicly available RNA-seq data. We also included transcriptome datasets of lung tissues from influenza A virus (IAV)-infected tree shrews and of hepatitis C virus (HCV)-infected tree shrew primary hepatocytes in the current analyses (authors’ unpublished data).

Evaluation of coding ability of transcripts and annotation of coding genes

We assembled RNA transcripts based on the RNA-seq reads from publicly available datasets (Supplementary Table S1) and the new data generated in this study (Supplementary Tables S2, S3) using StringTie (v2.1.1) in a reference-guided manner (-G) (Pertea et al., 2015). The assembled RNA-seq transcript models were merged with the SMRT models by StringTie merge (--merge option) to obtain a transcript model covering all transcriptome datasets. To ensure the credibility of the transcripts, we only retained transcripts with high-confidence expression levels (fragments per kilobase per million (FPKM)>0.5 in at least one sample) during the merging of the RNA-seq and SMRT transcripts. We used Gffcompare (https://ccb.jhu.edu/software/stringtie/gffcompare.shtml) (Pertea & Pertea, 2020) to compare the assembled transcript models with the transcript models of the reference genome TS_2.0 (reference model). We treated those transcripts with no matches to the reference models (class code “=”) as newly identified transcripts.

We evaluated the coding potential of the newly identified transcripts by incorporating the results predicted using the following approaches to achieve a more reliable annotation: (1) The Coding Potential Assessment Tool (CPAT) (Wang et al., 2013) was applied to evaluate the transcript coding ability using a logistic regression model. The hexamer frequency table was built using “make_hexamer_tab.py” script, and the logit model was built using “make_logitModel.py” script. (2) We used the Coding Potential Calculator 2 (CPC2) (Kang et al., 2017; Kong et al., 2007) to evaluate the coding ability of the transcripts employing a novel discriminative model based on four sequence-intrinsic features. (3) We used TransDecoder (https://github.com/TransDecoder/) to predict the high-confidence open reading frames (ORF) of each transcript. (4) Pfam (El-Gebali et al., 2019), which contains a comprehensive archive of protein domains, and UniProtKB/Swiss-Prot (Boutet et al., 2007), which contains a comprehensive archive of protein sequences from multiple species, were used to identify potentially translated ORFs. Except for CPAT, we ran all other programs with their default parameters and integrated the prediction results with the following procedures. First, we integrated the prediction results from both CPAT and CPC and only transcripts that met the coding cut-off of both approaches (CPAT, coding potential>0.4; CPC, designated as “coding”) were subjected to further analyses. Second, the ORF of each transcript was predicted using TransDecoder (https://github.com/TransDecoder/) and only transcripts with at least one high-confidence ORF were retained. Third, we scanned the potentially translated ORFs against the Pfam (http://rfam.xfam.org/) (El-Gebali et al., 2019) and UniProtKB/Swiss-Prot databases (Boutet et al., 2007). Those transcripts with at least one predicted Pfam domain or high protein sequence identity (E-value>1e-5) with at least one known protein were defined as coding transcripts. We selected the longest transcript from each gene locus as the representative transcript of the gene. We BLASTed the representative transcripts against the UniProtKB/Swiss-Prot and UniProtKB/Trembl databases (Boutet et al., 2007) using blastall (https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=Download). The best hit gene name of each gene locus was designated as the name of the newly annotated gene. For multiple gene loci with the same best hit, we modified the gene name by adding “LI (like)+number” to avoid gene name redundancy. We used eggnog-mapper (Huerta-Cepas et al., 2017) to BLAST the translated ORFs for each transcript against the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (Kanehisa et al., 2017) and Gene Ontology (GO) database (http://geneontology.org/) for gene function annotation. The best KEGG pathway and GO term hits were designated for each gene. Fourth, lncRNAs were identified and annotated from the transcripts using the following criteria: (1) transcripts are >200 nucleotides (nt) long and meet the non-coding cut-offs of CPAT (coding potential<0.4) and CPC (designated as “noncoding”); (2) transcripts have a predicted ORF<100 nt; and (3) transcripts have a low similarity (E-value>1e-5) with the tRNA family in the Rfam database ( El-Gebali et al., 2019) and UniProtKB/Swiss-Pro database (Boutet et al., 2007). We defined those transcripts showing inconsistent prediction of coding RNAs and lncRNAs based on the above approaches as biased transcripts.

After coding potential evaluation and gene annotation, we merged the TS_2.0 genome annotation file (Fan et al., 2019) with the newly annotated transcripts to generate the TS_3.0 genome annotation. We verified the accuracy of the TS_3.0 transcripts by comparing the annotated transcripts with those characterized by molecular cloning. In total, 30 transcripts reported in our previous studies (Gu et al., 2019a; Luo et al., 2018; Yao et al., 2019; Yu et al., 2014, 2016) (Supplementary Table S4) were selected for comparison using blastall (https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=Download). An E-value<1e-5 was used as a cutoff to define whether a transcript obtained by molecular cloning was included in the TS_3.0 genome annotation.

Benchmarking universal single-copy orthologs (BUSCO) analysis

We used BUSCO (Seppey et al., 2019) to evaluate the completeness of the TS_3.0 genome annotation using the Mammalia and Eukaryota BUSCO datasets (Simão et al., 2015; Waterhouse et al., 2018). To ensure that the sequences of each species originated from the genome annotation, we selected protein-coding genes according to the information provided in the gene transfer format (GTF) file for each species, and the longest transcript of each gene was selected for consideration. We compared the TS_3.0 annotation with the previous tree shrew genome annotations (Supplementary Table S6). We used Gffread (Pertea & Pertea, 2020) to extract the sequences from the reference genome according to the annotation information. BUSCO evaluation was run in transcriptome mode (-m tran).

Alternative splicing prediction, differential gene expression, and dimension reduction analyses

We used Kallisto (Bray et al., 2016) to quantify the expression level of each transcript. Briefly, the clean reads obtained from each RNA-seq dataset were mapped to the transcriptome constructed using the TS_3.0 annotation. We used SUPPA2 (Trincado et al., 2018) to predict alternative splicing events, including skipped exons (SE), alternative 5' splice sites (A5), alternative 3' splice sites (A3), mutually exclusive exons (MXE), and retained introns (RI). The splicing level of each gene in each RNA-seq dataset was quantified using the Percent Spliced-In (PSI) index. The PSI values of each gene were calculated based on the transcript model and expression level (transcripts per million, TPM) of all transcripts for the gene (Trincado et al., 2018).

To characterize tissue-specific gene/transcript expression levels and alternative splicing events, we identified specifically expressed genes by applying the Wilcoxon rank-sum test and Dunn test in the R (http://www.R-project.org/) package Seurat (Butler et al., 2018). P-values were adjusted (Padjust) by the Benjamini-Hochberg (BH) method. Genes/transcripts with a Padjust<0.05 were defined as differentially/specifically expressed in a given condition. Using Seurat (Butler et al., 2018), we performed uniform manifold approximation and projection (UMAP) (McInnes et al., 2020) for each tree shrew tissue based on TPM of mRNA and lncRNA at the gene and transcript levels, respectively. Those genes/transcripts expressed in all samples of a given tissue and with a gene/transcript |log2 fold-change|>0.5 when compared with other tissues were regarded as tissue-specific genes/transcripts.

We used the R package DESeq2 (Love et al., 2014) to identify differentially expressed genes (DEGs) under virus infection conditions. The Padjust values were calculated using the BH method, as described above. Genes were identified as dysregulated upon viral infection if Padjust<0.05 and |log2 fold-change|>1 were met. KEGG and GO enrichment analyses were performed using the R package ClusterProfiler (Yu et al., 2012), with P-values adjusted by the BH method. A pathway with a Padjust<0.05 was defined as significantly enriched.

Gene family analyses

We obtained protein sequences of multiple mammals from the Ensembl database (https://asia.ensembl.org/index.html), including Homo sapiens (GRCh38.p13), Pan troglodytes (Pan_tro_3.0), Gorilla gorilla gorilla (gorGor4), Macaca mulatta (Mmul_10), Rattus norvegicus (Rnor_6.0), Mus musculus (GRCm39), Sus scrofa (Sscrofa11.1), Bos taurus (ARS-UCD1.2), Canis lupus familiaris (CamFam3.1), and Oryctolagus cuniculus (OryCun2.0) (Supplementary Table S5). The orthologous relationships among species were calculated using OrthoFinder (Emms & Kelly, 2019). We used CAFÉ (De Bie et al., 2006; Mendes et al., 2020) to detect gene family size changes, including expansion and contraction, based on the orthogroups and phylogenetic tree constructed by OrthoFinder (Emms & Kelly, 2019). The phylogenetic tree was constructed using all protein-coding genes of the genomes.

Two gene families showed expansion in this study, i.e., STT3 oligosaccharyltransferase complex catalytic subunit B (STT3B) and subunit A (STT3A) and the interleukin 6 (IL6) gene family, which were featured for their potential roles in viral infection. We constructed gene trees of these gene families using the maximum-likelihood (ML) method (K2+G model) with 1 000 bootstraps. Trees were based on protein sequence alignment and constructed using MEGA (Kumar et al., 2018).

Tissue expression pattern and pathway gene similarity across species

We retrieved expression data from five tissues (liver, brain, kidney, testis, and heart) of mice (https://www.ebi.ac.uk/arrayexpress/E-MTAB-6798), rhesus macaques (https://www.ebi.ac.uk/arrayexpress/E-MTAB-6813), and humans (https://www.ebi.ac.uk/arrayexpress/E-MTAB-6814) from ArrayExpress (https://www.ebi.ac.uk/arrayexpress/) (Cardoso-Moreira et al., 2019). Principal component analysis of gene expression levels (TPM of each gene) for each tissue was constructed using the R package FactoMineR based on 11059 one-to-one orthologous genes and 3291 highly variable genes (defined by a coefficient of variation>0.1 TPM). Pathway gene information was retrieved from the KEGG database (Kanehisa et al., 2017). Protein sequence identities were calculated by BLASTing the tree shrew genes against human homologs and mouse genes against human homologs, respectively. Comparisons of protein sequence identity between mice and humans and between tree shrews and humans were performed using the Wilcoxon rank-sum test adopted in the R package. Here, P<0.05 was considered statistically significant.

RESULTS

Identification of high-quality tree shrew transcripts

To refine and annotate the chromosome-level tree shrew reference genome TS_2.0 (Fan et al., 2019), we adopted a stringent pipeline to integrate the related RNA-seq and SMRT datasets (Figure 1A). We collected and curated tree shrew transcriptome data across different tissues and cells with different viral infections. A total of 234 tree shrew transcriptome datasets were used in this study (Supplementary Tables S1–S3). These datasets covered a wide range of biological and pathological conditions, including cells infected with different viruses (Figure 1B) and normal and pathological tissue expression (Figure 1C). We used these datasets of diverse conditions to ensure that we captured the diversity and quantity of the transcripts, especially those with low abundancy in normal conditions. After QC, we discarded four samples with a mapping ratio below 75% from further analysis. The remaining datasets had a mean mapping ratio of 93.26% (Figure 1C, D), reflecting relatively high completeness of the reference genome and high quality of the RNA-seq datasets used. Based on the transcriptome datasets, we assembled a transcript model with FPKM>0.5 for each sample using StringTie (Pertea et al., 2015) and obtained 423 965 transcripts for the tree shrew (hereafter referred as RNA-seq transcripts).

Figure 1

Reference-guided transcriptome assembly of tree shrew TS_3.0 genome annotation

A: Integrative pipeline for tree shrew genome annotation using publicly available and newly generated transcriptome datasets. B, C: Number of RNA-seq datasets of virus-infected tissues/cells (B) and normal tissues (C) analyzed in this study. D, E: Mapping ratio of RNA-seq data of virus-infected tissues/cells (D) and normal tissues (E) relative to reference tree shrew genome TS_2.0 (Fan et al., 2019). Sample information in (D) is listed in Supplementary Table S1.

To verify the accuracy of the tree shrew RNA-seq transcripts, we constructed two ISO-seq full-length transcriptome libraries based on pooled RNA samples from eight tissues of two tree shrews. After QC and error correction, we obtained 36381 non-redundant non-chimeric full-length transcripts located at 12366 loci (hereafter referred as ISO-seq transcripts). The mean length of the ISO-seq transcripts was 2371 nt and the longest ISO-seq transcript (death effector domain containing [DEDD] gene) was 9417 nt. A total of 10968 transcripts were matched between the RNA-seq and ISO-seq transcripts, accounting for 30.15% of the ISO-seq transcripts. Nearly all ISO-seq transcripts were captured by the RNA-seq transcripts, and only 0.9% (1317/146347) of exons and 2.6% (319/12366) of loci captured by ISO-seq were missed by the RNA-seq transcripts. We combined the RNA-seq and ISO-seq transcripts and obtained a total of 403792 transcripts located in 98142 loci in the tree shrew genome. The compiled tree shrew transcripts were deposited in the Tree shrew Database (http://www.treeshrewdb.org/download.html).

Expanded list of tree shrew coding and long non-coding transcripts

Both lncRNAs and mRNAs share similar biogenesis pathways and are involved in multiple biological processes (Jiang et al., 2019; Quinn & Chang, 2016), but exert their functions in different manners (Dahariya et al., 2019; Melé et al., 2017). We categorized transcripts into high-confidence coding transcripts (mRNAs), lncRNA transcripts, and biased transcripts. Among the 403792 newly obtained transcripts, 115562 (>200 nt) were located in 56401 loci and were predicted to be lncRNA transcripts. Among these lncRNA transcripts, 50576 were antisense, 60118 were intergenic, and 4868 were bidirectional. We predicted 53298 coding transcripts located in 19242 loci. In addition, 234932 transcripts located in 93426 loci had inconsistent prediction regarding the characteristics of coding and lncRNA transcripts.

Overall, the expression levels of tree shrew lncRNAs were significantly lower than the expression levels of mRNAs at the gene and transcript levels (Figure 2A), consistent with previous reports on the human transcriptome (Iyer et al., 2015; Jiang et al., 2019). The exon number per lncRNA transcript was also significantly lower than that of mRNA in the tree shrew (Figure 2B), as observed in humans (Jiang et al., 2019). Moreover, the length of lncRNA was significantly shorter than that of mRNA in the tree shrew (Figure 2C). It has been reported that lncRNAs can regulate mRNA expression in a cis-regulatory manner (Jiang et al., 2019; Ørom et al., 2010; Ponjavic et al., 2009). Here, we calculated the expression correlation (Pearson correlation coefficient) between 10 000 mRNAs and their closest lncRNAs in the same genomic region across all RNA-seq datasets and compared the expression correlation between 10 000 randomly selected pairs of mRNAs and lncRNAs. Results showed that the expression correlation between closely located mRNA-lncRNA pairs was significantly stronger (Wilcoxon tests, P<2.2e-16) than that between randomly chosen pairs (Figure 2D). This provides additional evidence for the good accuracy of the lncRNA and mRNA annotations in the tree shrew genome.

Figure 2

Characteristics of tree shrew TS_3.0 transcripts

A: Expression level of mRNAs is greater than that of lncRNAs at gene and transcript levels. B: mRNA transcripts have a higher number of exons than lncRNA transcripts. C: Average length of mRNA transcripts is longer than that of lncRNA transcripts. Density plot was drawn based on kernel density and statistical analysis was performed with Wilcoxon rank-sum test. D: Tree shrew lncRNAs exert cis-regulatory function on expression of proximal mRNAs. Close, mRNA and lncRNA pairs neighboring each other on pseudochromosome of tree shrew. Random, randomly selected mRNA and lncRNA pairs distant from each other in the genome. E: Percentages of alternative splicing events in reference genome annotations of tree shrew (TS_3.0), mouse (GRCm39), and human (GRCh38.p13). SE, skipped exon; A5, alternative 5’ splice site; A3, alternative 3’ splice site; MXE, mutually exclusive exons; RI, retained intron. F: Expression level of newly annotated coding genes is lower than that of previously annotated genes at gene and transcript levels. G: BUSCO evaluation of different tree shrew genome annotations showing that current version (TS_3.0) is superior. H: Pathway enrichment analysis of newly annotated genes showing enrichment in 11 pathways (P<0.05). Values in A, B and F are presented as a boxplot, and statistical analyses were performed by Wilcoxon rank-sum test.

We further compared 30 transcripts obtained by molecular cloning in our previous studies (Gu et al., 2019a; Luo et al., 2018; Yao et al., 2019; Yu et al., 2014, 2016) (Supplementary Table S4) with those predicted in the TS_3.0 genome annotation. All 30 transcripts showed very good alignment with the currently annotated transcripts (blastall, E-value<1e-5), and 47 additional transcripts were identified in these gene loci according to the TS_3.0 genome annotation, suggesting high accuracy and completeness of the transcript annotation. For instance, we observed all tree shrewTLR gene family members (Yu et al., 2016) and six IL7 transcripts (Yu et al., 2014) reported in our previous studies in TS_3.0. The four IL7 transcripts showed a complete sequence match with the corresponding transcripts in TS_3.0 (Supplementary Figure S1). Among the five alternative splicing events in the tree shrew transcriptome, SE was the most common type of alternative splicing event in the TS_3.0 transcripts (Figure 2E). This pattern is consistent with that of humans and mice (Figure 2E).

Compared to the TS_2.0 tree shrew genome (Table 1), we found 6 126 coding transcripts (including 207 single-exon transcripts) located in 3 514 loci, none of which had been previously annotated and thus represented newly annotated genes. We profiled the gene expression patterns of these newly annotated genes across all RNA-seq datasets and found that the expression levels of the genes were significantly lower than those of the previously annotated ones (Figure 2F). The low abundancies of these newly annotated genes may be the reason for missing annotation in our previous studies (Fan et al., 2013, 2019). The newly annotated genes were enriched (BH adjusted, Padjust<0.05) in immune-related KEGG pathways, such as “Pattern recognition receptors” and “Inflammatory bowel disease (ko05321)” (Figure 2H), partly due to the bias of the RNA-seq datasets of tree shrew cells with viral infection. Combined with the previously annotated genes (Fan et al., 2019), 27 082 coding genes were finally annotated in the tree shrew genome.

Table 1

Comparisons of five tree shrew genome annotations

Parameters	TupChi_1.0 (NCBI)	TupaiaBase	TS_1.0	TS_2.0	TS_3.0
Tree shrew genome annotations TS_1.0 (Fan et al., 2013), TS_2.0 (Fan et al., 2019), and TS_3.0 (this study) were established in our studies. Tupchi_1.0, NCBI tree shrew annotation (https://www.ncbi.nlm.nih.gov/assembly/GCF_000334495.1/). TupaiaBase was reported by (Sanada et al., 2019). BUSCO: Benchmarking with Universal Single-Copy Orthologs. –: Not available.
Coding genes
Total number of coding genes	23 527	19 230	22 121	23 568	27 082
Transcript per coding gene	1.59	1	1	1	2.17
Annotated coding genes	23 537	12 612	20 225	20 811	25 127
Average mRNA length	48 104		33 712	40 114	41 239
Average CDS length	1 682	1 419	1 404	1 527	1 684
Average exon number	8.34	7.68	7.54	8.86	9.32
Average exon length	229	185	186	172	181
Average intron length	6 003	3 411	4937	4907	4 863
Complete BUSCOs (Eukaryota 255 genes)	216 (84.7%)	195 (76.5%)	221 (86.7%)	235 (92.2%)	250 (98.0%)
Complete BUSCOs (Mammalia 9 224 genes)	7 884 (85.5%)	6 080 (65.9%)	7 568 (82.0%)	7 519 (81.5%)	8 559 (92.8%)
Non-coding genes
Total number of lncRNA genes	3 718	–	–	–	56 401
Transcripts per lncRNA gene	5 179	–	–	–	2.05
Average lncRNA length	914	–	–	–	823
Average exon number	3.54	–	–	–	3.06
Average intron length	17 614	–	–	–	4 658

We further compared the gene, transcript, and lncRNA numbers between TS_3.0 and the previously reported versions of tree shrew genome annotation and found remarkable improvement (Table 1). Evaluation of gene completeness of the TS_3.0 annotation relative to the TS_2.0 annotation by BUSCO (Figure 2G; Supplementary Table S6) showed that the ratio of complete BUSCOs increased from 92.16% to 98.04% for Eukaryota BUSCOs (255 genes) and from 81.52% to 92.80% for Mammalia BUSCOs (9 224 genes). Compared with the NCBI TupChi_1.0 (https://www.ncbi.nlm.nih.gov/assembly/GCF_000334495.1) and TupaiaBase (Sanada et al., 2019), the current TS_3.0 annotation showed better completeness and better quality (Table 1).

Tissue expression and alternative splicing profiles of TS_3.0

To characterize the tissue-specific expression and alternative splicing patterns of each gene, we analyzed the transcriptome datasets from each tissue using UMAP and calculated the expression correlations. Both mRNAs and lncRNAs showed clear tissue-specificity in the context of UMAP (Figure 3A). The correlation matrix showed that the testis had the most unique tissue expression pattern compared with other tissues (Figure 3B). Notably, the testis possessed the most unique mRNAs and lncRNAs at the gene and transcript level (Figure 3C), consistent with the patterns reported in humans (Djureinovic et al., 2014). Many of the tree shrew testis-specific genes were involved in spermatogenesis (Supplementary Figure S2), as also reported in rats (Ji et al., 2020).

Figure 3

Tissue expression and alternative splicing profiles of tree shrew TS_3.0 transcripts

A: Tissue expression profiles of mRNAs and lncRNAs annotated in TS_3.0 at gene level (left panel) and transcript level (right panel). UMAP_1, UMAP dimension 1; UMAP_2, UMAP dimension 2. Detailed information on RNA-seq datasets of 13 tree shrew tissues is listed in Supplementary Tables S1, S2. B: Expression correlation matrix among different tree shrew tissues based on expression levels of mRNAs and lncRNAs. C: Tissue-specific expression patterns of mRNAs and lncRNAs at gene level (left) and transcript level (right). D: Comparison of PSI across 13 tree shrew tissues. P-value was calculated based on Dunn test and adjusted by Benjamini-Hochberg method. *: P<0.05;***: P<0.0005. Different colors indicate differentP values in triangle map. E: UMAP was constructed based on PSI of each gene in 13 tree shrew tissues.

We also characterized the specificity and intensity of RNA alternative splicing. Alternative splicing intensity of genes differed significantly across the 13 tree shrew tissues under study (Kruskal-Wallis rank-sum test, P<2.2e-16) (Figure 3D). Brain-related tissues (including the brain, hippocampus, and cortex) showed the highest PSI, whereas heart tissue had the lowest PSI. Furthermore, UMAP using PSI showed that the alternative splicing pattern for each gene also presented tissue specificity (Figure 3E). Collectively, we found that alternative splicing intensity and gene specificity were meticulously regulated across different tree shrew tissues (Figure 3E), which may account for the different functions of the respective tissues and organs.

Orthologous relationships of tree shrew genes with other mammals

We identified a total of 272 814 genes in the 11 mammals under study (Supplementary Table S5). Among them, 95.1% (259 313 genes) could be assigned to an orthogroup (i.e., a set of genes from multiple species descended from a single gene from the last common ancestor of that set of species) (Emms & Kelly, 2019). In total, 25 249 tree shrew genes could be assigned to 12 485 orthogroups. Among these orthogroups, 191 were paralogs and appeared to be tree shrew specific. We identified 17 299 orthogroups (including 14 549 one-to-one orthologs) shared between humans and tree shrews (Supplementary Table S7), which is a substantial improvement compared with the 12 840 one-to-one orthologs in TS_2.0 (Fan et al., 2019). Based on the current comprehensive orthologous relationships among species, we constructed a phylogenetic tree using the STAG algorithm (Emms & Kelly, 2019). We confirmed that the tree shrew is phylogenetically closer to primates than to rodents (Figure 4A), as described in our previous study based on 2 117 single-copy one-to-one orthologs (Fan et al., 2013).

Figure 4

Orthologous relationships and gene family size changes among different species

A: Phylogenetic tree of 11 mammals using orthogroups. Numbers on tree branches refer to numbers of gene family expansion (+red) and contraction (-blue), respectively. B: Maximum-likelihood (ML) trees of IL6 gene family. Coding sequence of the longest transcript for each gene in each species was used to construct ML tree. Values on tree branches refer to support of 1 000 bootstraps. The tree shrewIL6 gene family had 13 copies, labeled in red in the tree. C: Locations of 13 tree shrew IL6 gene copies on pseudochromosome 6 (chr6). D: ML trees of STT3B gene family and STT3A. Tree shrew STT3B gene family had 39 copies. E: Locations of 39 tree shrew STT3B gene copies on 15 pseudochromosomes and one unplaced contig. Pseudochromosomes and unplaced contig were defined in reference tree shrew genome TS_2.0 (Fan et al., 2019). F: Pathway enrichment of genes from tree shrew-specific gene families with significant expansion.

Gene gain and loss are important evolutionary processes that allow organisms to adapt to their environment (Page, 1998). Here, we analyzed changes in gene family size across 11 mammal species, including tree shrews (Supplementary Table S8), to validate and refine the previously characterized gene expansion and contraction events. We identified 120 gene families showing rapid expansion and 22 gene families showing rapid contraction (Figure 4A). The gene family exhibiting the greatest expansion was long-interspersed element-1 (LINE-1) retrotransposable element ORF (LIRE1), with 180 LIRE1 genes identified in the TS_3.0 genome annotation. The gene family exhibiting the greatest contraction was immunoglobulin heavy variable 3–35 (Supplementary Table S9). The guanylate binding protein (GBP) gene family was found to have rapidly contracted, consistent with the findings of our previous study (Gu et al., 2019a). We found that IL6 (Supplementary Figure S3) and STT3B (Supplementary Figure S4) were significantly expanded, with 34 of 39 STT3B gene family members being newly annotated and four of 13 IL6 gene family members being refined in the TS_3.0 genome annotation, respectively. The IL6 family contains IL6, cardiotrophin like cytokine factor 1 (CLCF1), cardiotrophin 1 (CTF1), ciliary neurotrophic factor (CNTF), interleukin 11 (IL11), interleukin 27 (IL27), LIF interleukin 6 family cytokine (LIF), and oncostatin M (OSM) (Rose-John, 2018). The tree shrew contained all these IL6 family members, and 12 of the 13 IL6 copies in tree shrews were not detected in other species. We found that each member of the IL6 family was grouped into a single clade in the ML tree of the IL6 family genes, confirming the close relationship of the expanded copies of each IL6 gene family member (Figure 4B). IL6LI7 appeared to be the ancestral gene of the tree shrew IL6 gene copies. The tree shrew IL6 family members were all located on pseudochromosome 6 (Figure 4C), with consistent exon numbers for each family member (Supplementary Figure S5). This suggests that the IL6 gene copies were most likely generated from tandem duplication and segmental duplication. We constructed an ML tree for the tree shrew STT3B gene family members, together with those of the other mammals, and the STT3A paralog. The STT3A and STT3B copies showed a gene-specific clustering pattern (Figure 4D). In the clade for the expanded STT3B copies from the tree shrew, STT3BLI27 diverged first and appeared to be the ancestral gene of the tree shrew STT3B family. Intriguingly, all 39 copies of STT3B in the tree shrew were distributed on 15 pseudochromosomes and one unplaced contig (Figure 4E). Of note, the tree shrew STT3BLI27 had 16 exons, while the other copies of STT3B contained no more than four exons (Supplementary Figure S6), suggesting that expansion of the tree shrew STT3B was most likely caused by retrotransposon activity.

To further dissect the potential evolutionary roles of the tree shrew gene family size changes, we conducted enrichment analysis using the canonical genes of each rapidly changing gene family. Results showed that gene families that have undergone rapid size change were enriched in the “immune response to tumor cell”, “regulation of cytokine production”, and “regulation of DNA metabolic process” pathways (Figure 4F).

Expression similarity across different species

To study the mRNA expression patterns of tissues and related pathways across humans, rhesus monkeys, mice, and Chinese tree shrews, we retrieved tissue RNA-seq data of mice, monkeys, and humans (Cardoso-Moreira et al., 2019), and compared their clustering patterns via principal component (PC) analysis. The species clustering patterns based on expression data from brain, liver, testis, kidney, and heart tissues showed distant divergence of mice from primates and tree shrews in the second PC, whereas humans, monkeys, and tree shrews were mainly separated by the first PC (Figure 5A). However, these clustering patterns should be considered with caution as the first and second PCs only contributed to a proportion of expression variance.

Figure 5

Expression similarities among different species

A: Tissue expression similarities among humans, rhesus monkeys, tree shrews, and mice. Expression patterns in five tree shrew tissues more closely resembled that of primates than that of mice. B: Comparisons of protein sequence identity of genes in KEGG pathways between tree shrews and humans and between mice and humans. C: Expression patterns of genes in brain-related pathways in tree shrews, rhesus monkeys, humans, and mice. Brain-related pathways included “Alzheimer’s disease”, “Parkinson disease”, “Neuroactive ligand-receptor interaction”, “Pathways of neurodegeneration-multiple diseases”, and “Axon guidance”.

We further determined the gene identity of tree shrews to humans using pathway analysis. For each human KEGG pathway, we compared protein sequence identities between mice and humans and between tree shrews and humans. Genes in 13 pathways showed greater protein sequence identity between tree shrews and humans than between mice and humans (Figure 5B). These 13 pathways included neuro-related pathways such as “Axon guidance”, “Parkinson disease”, and “Alzheimer disease”. Furthermore, proteins belonging to the “pathway in cancer” also showed higher identity between tree shrews and humans than between mice and humans (Figure 5B), suggesting that the tree shrew could be used to create valid cancer and neurodegenerative animal models. We also profiled the expression patterns of five pathways related to the brain (Figure 5C) and found that mice had a more distant clustering pattern than tree shrews with primates. Collectively, these results suggest that tree shrews are more genetically similar to primates than to mice at the transcriptomic level.

Changes in newly identified genes upon viral infection

To profile the transcriptome patterns of host immune responses to viral infection using the TS_3.0 genome annotation, we focused on the differential expression of the newly identified genes upon viral infection. Results showed that the TS_3.0 genome annotation had better accuracy and resolution for gene identification in tree shrew cells with or without viral infection. Some of the newly annotated genes were significantly dysregulated upon infection with HBV (67 genes), SeV (99 genes), NDV (48 genes), EMCV (seven genes), and ZIKA (178 genes) (Figure 6A). Of the 3 779 DEGs identified from the RNA-seq datasets by comparing infected and uninfected cells, only purine rich element binding protein A (PURA) and interferon induced protein 35 (IFI35) were significantly dysregulated in cells with all virus infections. Both genes are reported to have a pro-viral effect in other species (Das et al., 2014; Gounder et al., 2018; Krachmarov et al., 1996). As PURA is a newly annotated gene in TS_3.0, more studies should be carried out to characterize its role in viral infections in tree shrews.

Figure 6

Changes in expression of genes in virus-infected tree shrew tissues and cells

A: Changes in newly annotated genes upon viral infection. HBV, hepatitis B virus; NDV, Newcastle disease virus; EMCV, encephalomyocarditis virus; HSV-1, herpes simplex virus type 1; SeV Sendai virus; ZIKA, Zika virus. Genes were identified as differentially expressed genes (DEGs) upon virus infection if Padjust<0.05 and |log2 fold-change|>1. B: Plot of DEGs upon virus infection.PURA and IFI35 were dysregulated in cells and tissues infected with different viruses. Horizontal bar on left represents number of DEGs in each RNA-seq dataset. Dots and lines represent subsets of DEGs. Vertical histogram represents number of DEGs in each subset. C: Changes in expression of 19 gene copies of STT3B gene family upon SeV infection. Results are mean±standard deviation (SD). *: P<0.05;**: P<0.005;***, P<0.0005.P values were calculated using DESeq2.

Among the expanded gene family members, 19 of the 39 copies of the STT3B gene family in the tree shrew were up-regulated upon SeV infection. The oligosaccharyltransferase complex is known to be an essential host factor for dengue virus (DENV) replication (Lin et al., 2017). Considering that the 19 STT3B gene copies were not located on the same pseudochromosome, we speculated that they were up-regulated by the same transcription regulation system. The ancestral gene of the STT3B gene copies, STT3BLI27, was not dysregulated upon viral infection, suggesting that the 19 STT3B gene copies may have acquired the association with pro-viral function at a later stage of gene family expansion. However, further studies are required to confirm this speculation.

DISCUSSION

Comprehensive tree shrew genome annotation is crucial for developing animal models and for studying basic scientific questions (Yao, 2017). In this study, we annotated the Chinese tree shrew genome by integrating diverse RNA-seq datasets and newly generated ISO-seq datasets. We obtained a total of 27 082 coding genes (including 3 514 previously unannotated coding genes in TS_2.0 (Fan et al., 2019)) and 56 401 lncRNAs. Evaluation of the completeness of multiple tree shrew genome annotations using BUSCO (Seppey et al., 2019; Simão et al., 2015; Waterhouse et al., 2018) indicated that the current TS_3.0 annotation showed remarkable improvement in terms of completeness, which was achieved by incorporating diverse RNA-seq datasets that covered a wide range of biological and pathological conditions. The newly updated tree shrew TS_3.0 genome annotation can be downloaded from the Tree shrew Database (http://www.treeshrewdb.org/download.html).

Compared with the previous tree shrew genome annotation (Table 1), TS_3.0 provides a complete list of lncRNAs, which could help in the interpretation of the roles of lncRNAs in tree shrew biology and disease. The lack of lncRNA conservation among species is a considerable obstacle for functional annotation (Iyer et al., 2015). Here, we compared multiple characteristics between tree shrew mRNAs and lncRNAs and found smaller average exon number and shorter transcript length in lncRNAs than in mRNAs. We confirmed that the identified tree shrew lncRNAs may exert a cis-regulatory role on mRNA expression (Figure 3D). The overall characteristics of the tree shrew lncRNAs versus mRNAs resembled that of human lncRNAs versus mRNAs reported in previous studies (Iyer et al., 2015; Jiang et al., 2019; Ørom et al., 2010; Ponjavic et al., 2009).

Alternative splicing plays a key role in transcript processing and biological functions (Baralle & Giudice, 2017; Ule & Blencowe, 2019). By generating multiple transcripts from a particular gene, alternative splicing events can dramatically increase the diversity and complexity of the transcriptome, and can impact mRNA stability, localization, and translation (Baralle & Giudice, 2017; Ule & Blencowe, 2019). We previously showed that alternative splicing events in STING have played an important role in the innate immunity response of tree shrews against DNA and RNA viral infections (Xu et al., 2020a). Our updated annotation of tree shrew transcripts, especially from SMRT reads, provided an accurate model to characterize alternative splicing events in tree shrews. Using the TS_3.0 transcripts, we quantified and characterized the alternative splicing events and found a tissue-specific pattern of splicing intensity. The high occurrence of alternative splicing events in the tree shrew brain-related tissues is consistent with that found in humans (Rodriguez et al., 2020), suggesting that alternative splicing constitutes a straightforward strategy for enacting diverse functions such as tissue formation (Baralle & Giudice, 2017). It would be worth performing functional characterization of important genes that exhibit alternative splicing in different tree shrew tissues and/or in response to viral infection in the TS_3.0 genome annotation, as exemplified by the elegant functional assay for the STING isoform described in our recent study (Xu et al., 2020a).

The updated tree shrew TS_3.0 genome annotation could also provide insightful information regarding cross-species comparisons to initiate genome-based methods for creating animal models of human disease (Yao, 2017). We systematically characterized the orthologous relationships among experimental animals, including mice, monkeys, and tree shrews, using the newly updated tree shrew genome annotation. Orthologous comparison confirmed the closer relationship between primates and tree shrews than between primates and mice (Figure 4A), suggesting that tree shrews would be better model animals for biomedical research. We also compared the tissue expression patterns and related genes in particular pathways across four species, which again showed that tree shrews are closer to primates than to mice at the transcriptomic level (Figure 5C).

Gene expansion and contraction play key roles in environment adaptation (Yim et al., 2014). We re-appraised gene expansion and contraction events using the TS_3.0 transcripts and confirmed the gene families highlighted in our previous study (Fan et al., 2013). Among the 144 gene families that experienced size changes in the tree shrew, the IL6 and STT3B families may have particular biological implications. Notably, IL6 is thought to be actively involved in the cytokine storms observed in COVID-19 patients (Mehta et al., 2020; Vabret et al., 2020; Zhou et al., 2020) and therapy with the IL-6-receptor antagonist tocilizumab is considered a promising treatment for COVID-19 patients (Fu et al., 2020; Jones & Hunter, 2021). In SARS-CoV-2-infected tree shrews (Xu et al., 2020c; Zhao et al., 2020), different individuals demonstrated different susceptibility to SARS-CoV-2 and showed different viral loads after infection, though none of the infected tree shrews showed severe symptoms. Whether the expanded IL6 gene family played a role in this process is an interesting and important question. Cloning all 13 IL6 copies and characterizing the respective roles of each gene copy could help clarify why this gene family underwent expansion in the tree shrew. Among the 39 gene copies of the STT3B gene family, 19 were up-regulated upon SeV infection, whereas the other copies, including ancestral STT3BLI27, showed no such effect. The STT3B protein is a part of the oligosaccharyltransferase complex in humans (Lu et al., 2019), and is reported to play a pro-viral role in Dengue virus and HSV-1 infections (Lin et al., 2017; Lu et al., 2019). Expansion of the STT3B gene family may indicate a new immune response mechanism for tree shrews to counteract or facilitate these viral infections. However, more studies are required to characterize the function of the tree shrew STT3B gene family and to confirm the above speculation.

An important update of the TS_3.0 genome annotation was the inclusion of newly generated RNA-seq data from tree shrew cells and tissues challenged with different viruses. The inclusion of these datasets offers the chance to identify genes that are up-regulated or down-regulated upon viral infection for further study. Indeed, previously reported tree shrew genes that show altered expression upon viral infection (Gu et al., 2019a, 2021; Xu et al., 2016, 2020a, 2020b) could be confirmed. We identified several important targets showing a universal regulator effect, such as PURA and IFI35. The PURA gene encodes Pur-alpha, which has a repeated nucleic acid binding domain (Daniel & Johnson, 2018), and is reported to be regulated by transcription start sites Ⅰ and Ⅱ (Wortman et al., 2010). PURA is known to activate the John Cunningham virus in the glial cells of many acquired immunodeficiency syndrome patients (Krachmarov et al., 1996). In addition, IFI35 is an interferon-stimulated gene that negatively regulates RIG-I antiviral signals to support vesicular stomatitis viral replication (Das et al., 2014) and enhances H5N1 influenza disease symptoms (Gounder et al., 2018). We speculate that in vivo overexpression of both PURA and IFI35 may create tree shrew models more permissive to different viruses, including HCV and HBV, which have no feasible animal models at present.

In summary, we generated an improved tree shrew genome annotation using comprehensive RNA-seq and ISO-seq datasets. The updated version of the tree shrew genome annotation (TS_3.0) fixed some of the issues with previous versions, such as TS_1.0 (Fan et al., 2013) and TS_2.0 (Fan et al., 2019). Detailed annotation of the genes, gene families, and alternative splicing events in the tree shrew genome, as well as cross-comparison of expression patterns among different tissues and species, further illuminated the unique and common genetic features of tree shrews and provided further evidence of the considerable potential of tree shrews in biomedical research.

DATA AVAILABILITY

The TS_3.0 genome annotation data and newly generated RNA-seq and ISO-seq data are available from the Tree shrew Database (http://www.treeshrewdb.org/download/). Related data were also deposited in GSA (accession No. PRJCA006366).

Supplementary data to this article can be found online.

Click here for additional data file.

COMPETING INTERESTS

The authors declare that they have no competing interests.

AUTHORS’ CONTRIBUTIONS

Y.G.Y. and M.S.Y. conceived and designed the experiments. L.B.L. provided living tree shrews and tissues. D.D.Y. and L.X. isolated tree shrew primary cells and performed viral infection and RNA extraction. Q.Y.Z. provided the transcriptome data of lung tissues from IAV-infected tree shrews. M.S.Y., J.Y.Z., M.X., and Y.F. collected transcriptome data and performed genome annotation and transcriptome analyses. M.S.Y. and Y.G.Y. wrote the manuscript. All authors read and approved the final version of the manuscript.