Mao-Sen Ye, Jin-Yan Zhang, Dan-Dan Yu, Min Xu, Ling Xu, Long-Bao Lv, Qi-Yun Zhu, Yu Fan,Yong-Gang Yao1,2,,＊

1 Key Laboratory of Animal Models and Human Disease Mechanisms of the Chinese Academy of Sciences & Yunnan Province, KIZ/CUHK Joint Laboratory of Bioresources and Molecular Research in Common Diseases, Kunming Institute of Zoology, Chinese Academy of Sciences, Kunming, Yunnan 650204, China

2 Kunming College of Life Science, University of Chinese Academy of Sciences, Kunming, Yunnan 650204, China

3 National Resource Center for Non-Human Primates, National Research Facility for Phenotypic & Genetic Analysis of Model Animals(Primate Facility), Kunming Institute of Zoology, Chinese Academy of Sciences, Kunming, Yunnan 650107, China

4 State Key Laboratory of Veterinary Etiological Biology, Lanzhou Veterinary Research Institute, Chinese Academy of Agricultural Sciences, Lanzhou, Gansu 730046, China

ABSTRACT

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 chromosomelevel tree shrew genome.We obtained a total of 3 514 newly annotated coding genes and 50 576 lncRNA genes.We also characterized the tissuespecific 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.

Keywords: Tree shrew; Genome annotation;Transcriptome; Gene family; Virus infection

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 (Tupaiabelangerichinensis) 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) shortread 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 nextgeneration sequencing (NGS) based on RNA-seq and longread isoform sequencing (ISO-seq) have revealed the complexity and characteristics of eukaryotic transcriptomes(Chen et al., 2017).Using ortholog anddenovoannotations(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 highquality 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?UltraTMRNA 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 BluePippinTMSize 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?UltraTMRNA Library Prep Kit for Illumina?.The quality of each library was assessed using the Agilent Bioanalyzer 2 100 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 highconfidence 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 highconfidence 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_T YPE=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 theMammaliaandEukaryotaBUSCO 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 aPadjust＜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.ThePadjustvalues were calculated using the BH method, as described above.Genes were identified as dysregulated upon viral infection ifPadjust＜0.05 and |log2foldchange|＞1 were met.KEGG and GO enrichment analyses were performed using the R package ClusterProfiler (Yu et al.,2012), withP-values adjusted by the BH method.A pathway with aPadjust＜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),includingHomosapiens(GRCh38.p13),Pantroglodytes(Pan_tro_3.0),Gorillagorillagorilla(gorGor4),Macacamulatta(Mmul_10),Rattusnorvegicus(Rnor_6.0),Musmusculus(GRCm39),Susscrofa(Sscrofa11.1),Bostaurus(ARSUCD1.2),Canislupusfamiliaris(CamFam3.1), andOryctolaguscuniculus(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 11 059 one-to-one orthologous genes and 3 291 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

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 36 381 non-redundant non-chimeric full-length transcripts located at 12 366 loci (hereafter referred as ISO-seq transcripts).The mean length of the ISO-seq transcripts was 2 371 nt and the longest ISO-seq transcript (death effector domain containing [DEDD] gene) was 9 417 nt.A total of 10 96 8 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% (1 317/146 347) of exons and 2.6% (319/12 366) 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 403 792 transcripts located in 98 142 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 403 792 newly obtained transcripts,115 562 (＞200 nt) were located in 56 401 loci and were predicted to be lncRNA transcripts.Among these lncRNA transcripts, 50 576 were antisense, 60 118 were intergenic,and 4 868 were bidirectional.We predicted 53 298 coding transcripts located in 19 242 loci.In addition, 234 932 transcripts located in 93 426 loci had inconsistent prediction regarding the characteristics of coding and lncRNA transcripts.

論文針對“機(jī)械工程測試技術(shù)”課程多媒體教學(xué)資源的建設(shè)進(jìn)行研究，總結(jié)和歸納了本課程在實(shí)際教學(xué)過程中的多媒體教學(xué)資源制作軟件和方法。通過改進(jìn)傳統(tǒng)教學(xué)資源，旨在增強(qiáng)教學(xué)的靈活性和生動(dòng)性，提升學(xué)生學(xué)習(xí)興趣和積極性，讓學(xué)生更為直觀地理解和掌握課程知識(shí)點(diǎn)。

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 cisregulatory 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.

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 shrewTLRgene family members (Yu et al., 2016) and sixIL7transcripts (Yu et al.,2014) reported in our previous studies in TS_3.0.The fourIL7transcripts 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.

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% forEukaryotaBUSCOs (255 genes) and from 81.52% to 92.80% forMammaliaBUSCOs (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).

Figure 2 Characteristics of tree shrew TS_3.0 transcripts

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 testisspecific genes were involved in spermatogenesis(Supplementary Figure S2), as also reported in rats (Ji et al.,2020).

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.

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

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).

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 180LIRE1genes 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 thatIL6(Supplementary Figure S3) andSTT3B(Supplementary Figure S4) were significantly expanded, with 34 of 39STT3Bgene family members being newly annotated and four of 13IL6gene family members being refined in the TS_3.0 genome annotation, respectively.TheIL6family containsIL6, 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 theseIL6family members,and 12 of the 13IL6copies in tree shrews were not detected in other species.We found that each member of theIL6family was grouped into a single clade in the ML tree of theIL6family genes, confirming the close relationship of the expanded copies of eachIL6gene family member (Figure 4B).IL6LI7appeared to be the ancestral gene of the tree shrewIL6gene copies.The tree shrewIL6family members were all located on pseudochromosome 6 (Figure 4C), with consistent exon numbers for each family member (Supplementary Figure S5).This suggests that theIL6gene copies were most likely generated from tandem duplication and segmental duplication.We constructed an ML tree for the tree shrewSTT3Bgene family members, together with those of the other mammals,and theSTT3Aparalog.TheSTT3AandSTT3Bcopies showed a gene-specific clustering pattern (Figure 4D).In the clade for the expandedSTT3Bcopies from the tree shrew,STT3BLI27diverged first and appeared to be the ancestral gene of the tree shrewSTT3Bfamily.Intriguingly, all 39 copies ofSTT3Bin the tree shrew were distributed on 15 pseudochromosomes and one unplaced contig (Figure 4E).Of note, the tree shrewSTT3BLI27had 16 exons, while the other copies ofSTT3Bcontained no more than four exons(Supplementary Figure S6), suggesting that expansion of the tree shrewSTT3Bwas 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).

Figure 4 Orthologous relationships and gene family size changes among different species

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.

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 neurorelated 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.

Figure 5 Expression similarities among different species

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).AsPURAis 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

Among the expanded gene family members, 19 of the 39 copies of theSTT3Bgene family in the tree shrew were upregulated 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 19STT3Bgene copies were not located on the same pseudochromosome, we speculated that they were upregulated by the same transcription regulation system.The ancestral gene of theSTT3Bgene copies,STT3BLI27, was not dysregulated upon viral infection, suggesting that the 19STT3Bgene 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).

Table 1 Comparisons of five tree shrew genome annotations

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, theIL6andSTT3Bfamilies 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 expandedIL6gene family played a role in this process is an interesting and important question.Cloning all 13IL6copies 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 theSTT3Bgene family, 19 were up-regulated upon SeV infection, whereas the other copies, including ancestralSTT3BLI27, 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 theSTT3Bgene 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 shrewSTT3Bgene 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 asPURAandIFI35.ThePURAgene 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,IFI35is 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 thatinvivooverexpression of bothPURAandIFI35may 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

Supplementary data to this article can be found online.

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.