Downloaded from genome.cshlp.org on September 26, 2022 - Published by Cold Spring Harbor Laboratory Press Research Transcriptome innovations in primates revealed by single-molecule long-read sequencing Luis Ferrández-Peral,1,12 Xiaoyu Zhan,2,12 Marina Alvarez-Estape,1 Cristina Chiva,3,4 Paula Esteller-Cucala,1 Raquel García-Pérez,1 Eva Julià,3 Esther Lizano,1,5 Òscar Fornas,3,4 Eduard Sabidó,3,4 Qiye Li,2,6 Tomàs Marquès-Bonet,1,4,5,7,8,12 David Juan,1,12 and Guojie Zhang9,10,11,12 1 Institute of Evolutionary Biology (UPF-CSIC), PRBB, 08003 Barcelona, Spain; 2BGI-Shenzhen, Shenzhen 518083, China; 3Center for Genomic Regulation (CRG), Barcelona Institute of Science and Technology (BIST), 08003 Barcelona, Spain; 4Universitat Pompeu Fabra (UPF), 08003 Barcelona, Spain; 5Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Cerdanyola del Vallès, 08193 Barcelona, Spain; 6College of Life Sciences, University of Chinese Academy of Sciences, Beijing 100049, China; 7Institució Catalana de Recerca i Estudis Avançats (ICREA), 08010 Barcelona, Spain; 8CNAG-CRG, Center for Genomic Regulation (CRG), Barcelona Institute of Science and Technology (BIST), 08028 Barcelona, Spain; 9State Key Laboratory of Genetic Resources and Evolution, Kunming Institute of Zoology, Chinese Academy of Sciences, Kunming 650223, China; 10Section for Ecology and Evolution, Department of Biology, University of Copenhagen, DK-2100 Copenhagen 2200, Denmark; 11Evolutionary and Organismal Biology Research Center, School of Medicine, Zhejiang University, Hangzhou 310058, China Transcriptomic diversity greatly contributes to the fundamentals of disease, lineage-specific biology, and environmental adaptation. However, much of the actual isoform repertoire contributing to shaping primate evolution remains unknown. Here, we combined deep long- and short-read sequencing complemented with mass spectrometry proteomics in a panel of lymphoblastoid cell lines (LCLs) from human, three other great apes, and rhesus macaque, producing the largest full-length isoform catalog in primates to date. Around half of the captured isoforms are not annotated in their reference genomes, significantly expanding the gene models in primates. Furthermore, our comparative analyses unveil hundreds of transcrip- tomic innovations and isoform usage changes related to immune function and immunological disorders. The confluence of these evolutionary innovations with signals of positive selection and their limited impact in the proteome points to changes in alternative splicing in genes involved in immune response as an important target of recent regulatory divergence in primates. [Supplemental material is available for this article.] The vast complexity of eukaryotic transcriptomes arises from the Most previous comparative studies are based on short-read co-occurrence of multiple RNA processing events, leading to the RNA-seq, which cannot capture the actual isoform landscapes in generation of different isoforms encoded by the same gene. different species. The recent emergence of full-length isoform Particularly, the fine-tuning of alternative splicing (AS) is a critical sequencing overcomes this limitation and has contributed to the mechanism underlying disease and phenotypic evolution (Scotti refinement of isoform repertoires in model and nonmodel and Swanson 2016; Bush et al. 2017). Hence, previous research organisms (Byrne et al. 2019). Thus, single-molecule sequencing on the conservation of RNA processing events in humans and in combination with short-read RNA-seq and proteomic data ap- nonhuman primates (NHPs) has revealed the functional impor- pears as a powerful tool to disentangle the recent primate evolu- tance of transcriptomic diversity (Blekhman et al. 2010; Mazin tion of transcriptomes with isoform resolution, especially with et al. 2013, 2018; Zhang et al. 2017; Dougherty et al. 2018; the most recent high-quality genome assemblies for NHPs Kronenberg et al. 2018; Xiong et al. 2018; Mittleman et al. (Kronenberg et al. 2018; He et al. 2019; Warren et al. 2020; Mao 2021). On the other hand, the potential of increasingly complex et al. 2021), allowing the accurate definition of orthologous re- transcriptomes to produce alternative protein isoforms is under in- gions in primates. tense debate as most of the predicted proteins remain undetected In this study, we inspect the recent evolutionary dynamics of (Sheynkman et al. 2013; Abascal et al. 2015; Liu et al. 2017; Tress splicing programs using high-quality full-length transcriptomes et al. 2017; Carlyle et al. 2018; Agosto et al. 2019; Lau et al. obtained by integrating multispecies deep Iso-Seq and RNA-seq 2019), and alternative transcripts can carry out other regulatory data for a panel of lymphoblastoid cell lines (LCLs) from humans functions (Mignone et al. 2002). and NHPs (Fig. 1). These full-length transcriptomes substantially expand the current isoform repertoires in primates, and additional proteomic experiments (MS/MS) provide evidence of protein 12 These authors contributed equally to this work. translation from novel transcripts (Fig. 1). We leveraged our Corresponding authors: tomas.marques@upf.edu, david.juan@upf.edu, guojiezhang@zju.edu.cn Article published online before print. Article, supplemental material, and publi- © 2022 Ferrández-Peral et al. This article, published in Genome Research, is cation date are at https://www.genome.org/cgi/doi/10.1101/gr.276395.121. available under a Creative Commons License (Attribution 4.0 International), Freely available online through the Genome Research Open Access option. as described at http://creativecommons.org/licenses/by/4.0/. 1448 Genome Research 32:1448–1462 Published by Cold Spring Harbor Laboratory Press; ISSN 1088-9051/22; www.genome.org www.genome.org Downloaded from genome.cshlp.org on September 26, 2022 - Published by Cold Spring Harbor Laboratory Press Transcriptome innovations in primates Figure 1. Characterization of transcriptomes and proteomes of human (Homo sapiens), chimpanzee (Pan troglodytes), gorilla (Gorilla gorilla), orangutan (Pongo pygmaeus), and rhesus macaque (Macaca mulatta) lymphoblastoid cell lines (LCLs). We performed multispecies full-length isoform sequencing (PacBio Iso-Seq), RNA-seq (Illumina), and tandem mass spectrometry experiments (MS/MS). (Mya) Million years ago. transcriptomes to reconstruct transcript expression gains and loss- set, which comprises full-length isoforms from 10 human cell lines es in the primate lineage, linking these expression differences to derived from multiple tissues (Supplemental Fig. S4). The AS genetic changes in splice sites and novel exonization events. events newly disclosed by deep long-read sequencing (73% of all Innate immune system–related genes and cell proliferation genes detected events, 31,544 out of 43,179) unveil a high level of iso- are preferential targets of evolutionary innovations. Moreover, form diversity in primates that is missed by current annotations, the convergent accumulation of these innovations and their oc- especially for NHPs (Supplemental Fig. S5). These events are most- currence in genes under positive selection in primates highlight ly associated with canonical splice sites (GT-AG, GC-AG, and the strength of the evolutionary pressures suffered by the immune AT-AC), and thus, they agree with the previous knowledge on response and other cell type–specific functions. RNA splicing mechanisms. Moreover, the contribution of different AS modes is similar across all primate species, and in line with pre- vious observations, exon skipping (SE) and retained introns (RI) Results are the most prevalent mechanisms (Fig. 2B; Supplemental Table Isoform sequencing largely improves human and other S6; Zhang et al. 2017). primate annotations To define high-quality transcriptomes, we combined high-depth Mass spectrometry detects missing proteins in NHP long reads with short reads in our panel of LCLs derived from reference proteomes humans and NHPs (Supplemental Table S1), resulting in the larg- Because most Iso-Seq transcripts are predicted to be protein cod- est full-length isoform catalog in primates to date (Supplemental ing, we performed high-throughput mass spectrometry experi- Table S2). Briefly, Iso-Seq subreads for two human, two chimpan- ments to search for evidence of peptide translation. To do so, we zee, two gorilla, one (two isogenic cultures) orangutan, and two predicted the ORF of all Iso-Seq transcripts in all species and built rhesus macaque LCLs were processed to produce an average of a comprehensive protein database that constitutes our search over 2.3 million circular consensus sequences (CCS) per species space for peptides across human and NHP LCLs, allowing cross- (range from 1.9 to 2.9 million CCS). A strict quality control, includ- species peptide identifications for genes detected by Iso-Seq (Sup- ing base correction and splice junction adjustment with RNA-seq, plemental Methods). Within this search space, we detect a total was conducted to refine the transcript models (Methods) (Supple- of 49,428 distinct tryptic peptides (FDR < 5% and Mascot ion mental Methods). As a result, we defined a total of 148,610 iso- score > 20, excluding contaminants), from which over 24,900 are forms (on average, 29,722 per species) in 42,814 genes (on confidently detected, on average, in each species (reporter ion in- average, 8563 genes per species), which were further assessed by tensity signal [abundance] per sample greater than 50 in any sam- SQANTI (for the data production in each processing step, isoform ple of a given species, and species-wise median ratio of sample length distributions, and saturation analyses, see Supplemental abundance to pool abundance greater than 0.6) (see Supplemental Tables S3, S4; Supplemental Figs. S1, S2; Tardaguila et al. 2018). Table S7; Supplemental Data S2; Supplemental Methods). These Our refined transcriptomes reveal a considerable fraction of peptides confidently detected in our mass spectrometry experi- novel isoforms (unannotated in each reference transcriptome; ments are more frequently annotated in reference proteomes 49%, 73,161 out of 148,610) derived from unreported combina- from any species (>99% are annotated, 34,214 out of 34,266 de- tions of annotated splice sites, novel splice sites, fusion events, tected peptides) compared with the total set of tryptic peptides pre- and intergenic transcription (Methods) (Fig. 2A; Supplemental dicted from Iso-Seq isoforms (∼88% are annotated, 885,709 out of Fig. S3; Supplemental Table S5; Supplemental Data S1; Supplemen- 1,009,703 predicted peptides). tal Methods). We also captured novel human isoforms (∼19% of A number of Iso-Seq-derived peptides detected in each species our human Iso-Seq transcripts, 9762 out of 50,769) compared are novel (i.e., not annotated) in the corresponding UniProt refer- with the Universal Human Reference RNA (UHRR) Iso-Seq data ence proteomes (Supplemental Methods), including 22 for human Genome Research 1449 www.genome.org Downloaded from genome.cshlp.org on September 26, 2022 - Published by Cold Spring Harbor Laboratory Press Ferra ́ n dez-Peral et al. A B C D Figure 2. Isoform sequencing and mass spectrometry proteomics largely improve reference annotations. (A) Percentage of known versus novel tran- scripts (top) and junctions (bottom) captured by Iso-Seq (Ensembl V91). (B) Percentage of alternative splicing (AS) events detected by Iso-Seq (left) and percentage of splice site combinations associated with each AS class (right). (SE) Skipping exons, (RI) retained introns, (A5SS) alternative 5′ splice sites, (A3SS) alternative 3′ splice sites, (MEX) mutually exclusive exons. (C ) Number of annotated genes with novel detected peptides in each species for UniProt reference proteomes, NCBI RefSeq CDS annotations, or both. UniProt annotations are based on hg38, panTro5, gorGor4, ponAbe2, and rheMac8. NCBI RefSeq CDS annotations are based on hg38, panTro6, gorGor6, ponAbe3, and rheMac10. Peptides mapping to multiple isoforms are in- cluded. (D) Example of an unannotated gene (RIOX1 locus shown in gorGor6 assembly) supported by Iso-Seq, RNA-seq, and mass spectrometry in gorilla samples. [Minimum–maximum] intervals for the number of mapped RNA-seq reads are represented for each gorilla LCL. (hg38), 320 for chimpanzee (panTro5), 882 for gorilla (gorGor4), incompleteness in primate reference proteomes, reflecting the dif- 1590 for orangutan (ponAbe2), and 207 for rhesus macaque ferences in the previous efforts made by the community to com- (rheMac8) (for the number of annotated genes with novel detected plete the gene repertoire and study the proteome of each species. peptides, see Fig. 2C). To explore these results in the context of the These improvements in proteome annotations provide a van- most recent primate assemblies, we then searched the detected tage point to study the expression of proteins in primates. An peptides in the new CDS annotations (NCBI RefSeq CDS), which example of an unannotated gene with multiple supportive pep- provide not only curated proteins (NP) but also in silico predic- tides detected in NHP is RIOX1 (Fig. 2D), which is classified as a tions (XP). On average, nearly 78% of our detected tryptic peptides pseudogene in Ensembl (gorGor4), and it is absent in RefSeq in NHPs are only present in predicted entries (XP) from NHP (gorGor6) despite its importance in cancer research (Bundred RefSeq proteomes, highlighting the value of our multispecies et al. 2018). In other cases, we expand the exon and CDS annota- mass spectrometry experiments in primates to confirm such pro- tions, as in the gene SON in orangutan (ponAbe3), a key splicing tein predictions for the first time. When considering both curated factor in cell cycle progression involved in the response to the hep- (NP) and predicted (XP) proteins, we detected 25 novel peptides in atitis B virus (Supplemental Fig. S8; Sun et al. 2001; Ahn et al. 2011; hg38, 98 in panTro6, 94 in gorGor6, 225 in ponAbe3, and 75 in Kim et al. 2016; Tokita et al. 2016). Our results show that integrat- rheMac10 RefSeq CDS sets corresponding to unannotated and an- ing cross-species Iso-Seq with mass spectrometry significantly im- notated genes (Fig. 2C; Supplemental Data S3, S4). Based on these proves proteomic annotations even in the most recent NHP novel peptides, we only observed unannotated protein-coding assemblies. genes in NHPs, with all their human orthologs being actually an- notated in the hg38 gene models. The structural classification of Iso-Seq isoforms from which these novel peptides were predicted Emergence of species-specific genomic regions under and the number of species in which the peptides were simultane- active transcription ously detected are found in Supplemental Figures S6 and S7, re- To study the recent transcriptome evolution in primates, we first spectively. Our results clearly illustrate the distinct levels of focused on the genomic presence/absence of the expressed 1450 Genome Research www.genome.org Downloaded from genome.cshlp.org on September 26, 2022 - Published by Cold Spring Harbor Laboratory Press Transcriptome innovations in primates transcripts in the different species. To do so, we inspected orthol- that are consistently expressed or unexpressed based on the ogous transcripts for the set of isoforms detected in at least one RNA-seq from all biological replicates across species (Methods) species’ Iso-Seq, as long as they are transcribed from genes with (Fig. 3A,B; Supplemental Fig. S14; Supplemental Data S9; Supple- one-to-one orthology in the five primates (Supplemental Fig. S9; mental Methods). The number of high-confidence expressed tran- Supplemental Methods). The genomic coordinates of 49,797 tran- scripts (H: 14,783; C: 14,602; G: 14,438; O: 14,014; M: 13,880) and script models (out of 55,888, 89.10%) from 7858 genes were the number of species in which they are simultaneously expressed effectively mapped to the five primates’ assemblies to define (for an evaluation of transcriptome complexity across the five spe- orthologous transcripts (Supplemental Fig. S10; Supplemental cies, see Supplemental Fig. S15) are balanced across primates, with Methods), highlighting the high genomic conservation of tran- a high percentage of novel transcript models in each reference script structures (Supplemental Data S5, S6). In addition to this transcriptome being also unannotated in the remaining species high conservation, we found 61 species-specific exons without (63% on average) (Supplemental Fig. S16; Supplemental Methods). an orthologous counterpart in the most recent assemblies of the Based on these criteria, we defined 1215 high-confidence spe- other primates (Methods) (Supplemental Data S7; Supplemental cies-specific transcripts (only expressed in a single species; H: 332; Methods). A significant fraction of them (72%, 44 out of 61) C: 201; G: 167; O: 128; M: 387 transcripts) produced by genes with have not been annotated before, which is possibly owing to their one-to-one orthology in the five primates. The vast majority of lower expression levels (Supplemental Fig. S11). them (83%, 1009 out of 1215) arise from species-specific splice Compared with conserved exons, these species-specific exons junctions, with SE being the main mechanism involved in spe- are more frequently located in repetitive regions (79% of them cies-specific isoform expression gains (Supplemental Fig. S17). To overlap repeat elements, 48 out of 61), particularly SINEs (mainly inspect the usage of our set of species-specific splice junctions Alu elements) and LINEs (Supplemental Fig. S12), transposable el- (H: 266; C: 156; G: 140; O: 116; M: 337 junctions) across human ements that have been previously associated with the birth of new tissues, we searched RNA-seq junction reads in GTEx samples exons (Makałowski et al. 1994; Sorek et al. 2002; Cordaux and and found that human-specific junctions are widely used, being Batzer 2009; Schwartz et al. 2009; Avgan et al. 2019). We then supported in a median of 29 human tissues. In contrast, species- wondered if these exonizations might impact the protein function specific junctions from NHPs are largely unsupported in human by inspecting the combinations of protein domains encoded in tissues (median = 0 supportive GTEx tissues), suggesting that our the transcript (Supplemental Methods). In line with previous ob- results are extensible to other cell types (Supplemental Figs. S18– servations (Corvelo and Eyras 2008), we found that the incorpora- S20; Supplemental Methods). Furthermore, we confirmed that tion of these exons in the transcript usually results in annotated our detections of species-specific gains are not driven by increased protein domain combinations with likely conserved functions, gene expression levels in the corresponding species, because 783 with only ∼7% of exonization events (four out of 61) deriving in out of 868 genes expressing species-specific transcripts (90.21%) novel protein domain combinations. In most cases (85% of exoni- are not significantly up-regulated in the corresponding species zations, 52 out of 61), the new exon is partially (N = 35) or (DESeq2 Wald test, Benjamini–Hochberg adjusted P-value < 0.1 completely (N = 17) included in untranslated regions (UTRs), being and log2FoldChange > 0; see Supplemental Methods) (for the dis- significantly depleted in coding sequences compared with con- tribution of total gene expression levels and number of projected served exons (Fisher’s exact test P = 0.02391, OR = 0.51). Thus, transcript models in genes only producing transcripts expressed our findings are consistent with the previously observed relation- in all five primates and genes that express species-specific tran- ship between exon age and coding potential (Zhang and Chasin scripts, see Supplemental Figs. S21, S22). This set of transcripts is 2006; Merkin et al. 2015). As an example, in our LCLs, the exoni- enriched in unannotated isoforms (OR = 3.86, Fisher’s exact test zation of three repeat elements in human MRNIP results in a longer P = 1.4 × 10−117), demonstrating the suitability of our strategy to 3′ UTR with predicted small interfering RNA (siRNA) binding sites uncover isoform gains otherwise missed by reference transcrip- (Supplemental Fig. S13). tomes (Supplemental Figs. S23, S24). Importantly, the genes expressing species-specific transcripts show an overrepresentation in innate immune system functions Isoform innovations revealed by single-molecule sequencing compared with the background of LCLs-expressed genes under The high conservation of the gene structures led us to analyze the evaluation (Methods) (Fig. 3C; Supplemental Table S8) and are evolution of transcript expression in primates. We focused on the more often involved in immune system–related diseases set of orthologous transcript models in the five genomes, and (e.g., OAS2, OAS3, RNASEL, IRF7, and JAK1; Fisher’s exact test quantified transcript expression levels using high-depth RNA-seq P = 0.016, OR = 1.50). Given the involvement of these genes in data from 15 LCLs (three LCLs per species, including the cell lines the innate immune response, we then asked if they were enriched used to produce long-read transcriptomes) (see Supplemental Ta- in any functional subgroup according to a curated classification of ble S1) to obtain accurate expression estimates (Methods) (Sup- innate immunity genes (Deschamps et al. 2016). In particular, plemental Methods; Supplemental Data S8). Consistent with immune sensors (Fisher’s exact test P = 0.019, OR = 1.96), effectors widespread transcriptome conservation, a high number of tran- (P = 0.007, OR = 2.5), and secondary receptors (P = 0.027, OR = scripts are expressed in all five species (N = 18,935 transcripts), 3.31) are more likely to express species-specific isoforms in any pri- whereas 54% (22,198 out of 41,133) of all projected transcripts mate species. Furthermore, genes expressing species-specific tran- are only expressed in some species. Thus, to further understand scripts in any primate species tend to be highly expressed in the the patterns of transcript expression gains and losses in the pri- human spleen, which actively participates in the immune response mate lineage, we reconstructed the evolutionary dynamics under (GTEx expression data across all human tissues, hypergeometric test Wagner parsimony (Csű ös 2010). We reasoned that at least part −log10 corrected P = 3.14, FC = 1.95; Methods). We also found of the transcripts not detected in all species can result from inter- that these genes producing species-specific isoforms are enriched individual transcript expression variability. To account for this var- in genes involved in cell proliferation, consistent with iability, we established stringent criteria to keep the transcripts an enrichment in LCL-specific functions (Fisher’s exact test P = Genome Research 1451 www.genome.org Downloaded from genome.cshlp.org on September 26, 2022 - Published by Cold Spring Harbor Laboratory Press Ferra ́ n dez-Peral et al. A B C D Figure 3. Transcript expression gains and losses in the primate lineage. (A) Dynamics of transcript expression gains and losses in five primate species (inferred by Wagner parsimony). The phylogenetic tree is built based on the known phylogeny. Transcript expression gains and losses in macaque versus great apes are shown in macaque’s branch, although the direction (gain or loss in macaque or great apes) cannot be determined. (B) Overlap of genes expressing highly conserved (blue) and species-specific transcripts (orange). (C) Functional enrichment for genes displaying species-specific isoform inno- vations (N = 868 genes). Overrepresentation analysis (ORA) was conducted using the Gene Ontology Biological Process database (affinity propagation clus- tering). False-discovery rates (FDRs) were adjusted using the Benjamini–Hochberg method (Benjamini et al. 2001). The color legend indicates the number of overlapping genes between the gene set under evaluation and the pathway gene set. (D) Example of a gene, CYLD, expressing transcripts in all five spe- cies (conserved; blue shadow), exclusive expression in human (human-specific; orange shadow), and exclusive expression in rhesus macaque (macaque- specific; green shadow). Gray shadows indicate species-specific splice junctions. Untranslated regions (UTRs) and coding sequences (CDSs) are colored in blue and orange, respectively. 0.002, OR = 1.61). In contrast, we observed that this set of genes is transcripts encoding a different ORF than the isoforms expressed depleted in housekeeping genes (Fisher’s exact test P = 5.46 × 10−6, in all species (Supplemental Fig. S26). Our results indicate that im- OR = 0.64), showing that the essential role of housekeeping genes mune system– and cell proliferation–related genes are preferential- is also reflected in the conservation of their global splicing ly targeted by recent transcriptomic species-specific innovations structures. within the 7858 orthologous genes expressed in LCLs. An illustrative example of transcriptomic innovation is CYLD, a tumor-suppressor deubiquitinase involved in NF-kB acti- vation (Sun 2010; Mathis et al. 2015), whose splicing patterns im- The impact of cis-regulatory genetic changes in isoform pact the regulation of immune cells activation (Hövelmeyer et al. innovation 2007; Srokowski et al. 2009). In LCLs, we observe species-specific The efficient inclusion of exons into mature transcripts is mainly expression in CYLD isoforms (Fig. 3D; Supplemental Fig. S25). In determined by splicing signals located in the exon boundaries. this case, the inclusion of a noncoding exon generates two hu- Thus, interspecies genetic changes affecting the orthologous splice man-specific transcripts, whereas an alternative acceptor site sites might contribute to events of transcript diversification, in- produces a macaque-specific transcript. Another example is inter- cluding the emergence of species-specific isoforms. In the context feron induced protein 44 like (IFI44L), which is involved in the im- of our transcriptomes, we examined the sequence conservation of mune response to infection and autoimmune disorders (Herberg orthologous splice sites (terminal intron dinucleotides). Overall, et al. 2016; Zhao et al. 2016; de Oyarzabal et al. 2019; Bian et al. we observed widespread sequence conservation of splice sites 2020; Busse et al. 2020; Hadjadj et al. 2020; Mick et al. 2020; across primates, reflecting strong selective constraints in these re- Shaath et al. 2020; Lempainen et al. 2021; Li et al. 2021; gions (Supplemental Fig. S27). To explore high-impact interspecies Unterman et al. 2022) and whose isoform expression has also changes in splice sites, we then classified all orthologous junctions been found to be clinically relevant (Feenstra et al. 2014; into canonical and noncanonical, associated with high and low Haralambieva et al. 2017). IFI44L expresses three human-specific splicing efficiencies, respectively. After filtering out lowly 1452 Genome Research www.genome.org Downloaded from genome.cshlp.org on September 26, 2022 - Published by Cold Spring Harbor Laboratory Press Transcriptome innovations in primates expressed isoforms and genes (Supplemental Methods), we ob- can produce identical proteins, we also evaluated the conservation served 82,527 splice junctions present in our transcript models, of the ORF and of the combination of protein domains encoded by from which 1318 show a change in canonicity in any species com- the dominant transcripts (rank 1) in different species. We found pared with their orthologous positions. This scenario affects 16% of that all species share the same dominant transcript for 55.44% of genes (N = 1008 genes), which are enriched in T cell (FDR = 0.019) the genes (1814 out of 3272), a trend that is more evident as bigger and B cell (FDR = 0.019) activation as well as JAK/STAT signaling fold changes of expression between rank 1 and rank 2 transcripts (FDR = 0.022) pathways compared with the background of genes are considered (Fig. 4B), suggesting that changes in the dominant under evaluation (Methods) (Supplemental Table S9), showing transcript are less common. This conservation is higher for the pro- that cis-regulatory changes in orthologous splice sites are primarily tein-coding features, with 78.33% and 90.01% of genes sharing associated with genes involved in immune system functions. identical ORF or combinations of protein domains in their domi- We then investigated the contribution of these genetic nant transcripts in all species, respectively, consistent with chang- changes to the expression of species-specific isoforms. We ob- es in the dominant transcript being almost always functionally served that 69% of the high-confidence species-specific isoforms conservative. Our results highlight the greater constraints associat- (833 out of 1215 isoforms) harbor conserved canonicity of their ed with protein structure and function conservation, whereas dis- splice sites in all species. In contrast, the expression of 184 spe- tinct preferences in isoform UTR or even subtle ORF changes are cies-specific transcripts is driven by genetic differences producing more frequent in the dominant transcripts across the primate canonical splice sites only in the species in which they are ex- phylogeny. pressed. The enrichment in immune system functions holds for We investigated other changes of relative transcript usage species-specific transcripts regardless of if they occur with highly across primates. To do so, we conducted pairwise comparisons to conserved orthologous splice sites or not, showing the global rele- find interspecies differential isoform usage (DIU) (Supplemental vance of recent transcriptome changes in the evolution of the im- Methods) after excluding lowly expressed genes and isoforms mune system (Supplemental Fig. S28). (Soneson et al. 2016), as well as genes expressing highly similar An example of interspecies genetic differences in the splice transcripts (variability in internal exon boundaries) (Supplemental sites is RRP15, a gene involved in cell proliferation and apoptosis Methods). In a scenario of broadly conserved IU across species, we (Wu et al. 2016). In RRP15, a high-frequency variant only present detected over 800 transcripts showing species-specific usage in modern humans (Kuhlwilm and Boeckx 2019) disrupts an ac- changes (Supplemental Data S10; for a comparison with a less con- ceptor site of a coding exon (from AG to TG), and a downstream servative strategy to retrieve species-specific IU changes, see Sup- splice site is used instead. This mutation results in a human-specif- plemental Fig. S33), a significant fraction of which are missed by ic exon shortening representing a deletion of three amino acids reference annotations (Fig. 4C; Supplemental Fig. S34). LCL-specif- (disordered residues involved in protein binding), whereas the ic genes display species-specific DIU changes more frequently reading frame is preserved compared with other primates (Fisher’s exact test P = 0.015, OR = 1.69), although we did not (Supplemental Fig. S29; Zhang and Kurgan 2019; Zhao et al. observe any pathway enrichment for the genes displaying spe- 2021). In other scenarios, the genetic change in splice sites pro- cies-specific DIU when considering all evaluated genes as the back- motes the transcription of a new exon, as we observe in gorilla ground (Methods). An example of human-specific DIU is the LCLs for NAP1L4 gene (Supplemental Fig. S30). interleukin 4 receptor (IL4R), a gene under positive selection in primates (van der Lee et al. 2017), which plays a major role in the immune response (Nelms et al. 1999). We detected four up-reg- Evolutionary dynamics of isoform usage patterns ulated IL4R transcripts in human LCLs, two of them encoding in the primate lineage shorter IL4R proteins (Fig. 4D). After studying the evolution of switching on/off isoform trans- cription, we aimed to investigate the evolution of isoform tran- scription levels. To disentangle isoform-specific transcription Local differential splicing dynamics in transcribed regions regulation from changes in total gene expression, we assessed Although full-length isoform expression levels provide an accurate the relative transcript usage across the five species. To do so, we view of major usage changes in clearly different transcripts, other computed sample-wise fractions of isoform usage (IU) in our quan- local splicing changes are better detected by total exon expression tified orthologous transcripts (Supplemental Methods) (for the dis- levels that accumulate the signal of all transcripts including the tribution of IU values across LCLs and species, see Supplemental exon. We assessed the differential exon usage (DEU) across the pri- Fig. S15). Principal component analysis and hierarchical clustering mate phylogeny to obtain a comprehensive view of isoform tran- indicate that IU values cluster by species and that differences across scription evolution (Supplemental Methods). We first leveraged species increase with evolutionary distance, reflecting the known our projected Iso-Seq models to define 201,191 orthologous exon- phylogeny (Fig. 4A; Supplemental Methods). As expected, IU val- ic parts (i.e., transcribed regions) from 7766 one-to-one ortholo- ues are higher in annotated transcripts than in novel isoforms gous genes. Then, RNA-seq read counts in these exonic parts (Wilcoxon rank-sum test P < 2.2 × 10−308), indicating that isoforms were obtained across all samples to detect interspecies differences displaying important contributions to the expression of their in relative exon usage (Supplemental Methods). We observed near- genes are better represented in reference annotations (also see Sup- ly 4000 exonic parts displaying species-specific local usage chang- plemental Figs. S31, S32). es (Supplemental Data S11), with skipping exons (SE) being a To investigate the usage conservation of transcripts with a prevalent AS event among them (for a description of features in major contribution to the total gene abundances, we then focused SE showing species-specific usage changes, see Fig. 5A). As in the on the isoforms with the highest IU across primates (rank 1 and species-specific exonic structures (see section “Emergence of spe- rank 2 transcripts) after filtering out genes with very similar tran- cies-specific genomic regions under active transcription”), we ob- scripts (internal exon extensions leading to quantification ambi- served significant enrichment of Alu elements in human-specific guities) (see Supplemental Methods). Because different isoforms up-regulated exonic parts compared with down-regulated or Genome Research 1453 www.genome.org Downloaded from genome.cshlp.org on September 26, 2022 - Published by Cold Spring Harbor Laboratory Press Ferra ́ n dez-Peral et al. A B C D Figure 4. Isoform usage (IU) dynamics across primates. (A) Principal component analysis (PCA) for isoform fraction values of 43,484 transcripts from one- to-one orthologous genes. A hierarchical clustering for the Euclidean distances of IU Spearman’s correlations across samples is also shown. IU values were calculated from transcripts per million (TPM) after batch effect correction (comBat) and trimmed mean of M-values (TMM) normalization (Supplemental Methods). (B) Percentage of rank 1 isoform (green), ORF (orange), and protein domain combination (purple) shared by all primates by different thresholds of expression fold change between rank 1 and rank 2 isoforms. The number of genes considered for each threshold is indicated in the inner bar plot. (C) Number of isoforms showing species-specific usage changes. Up-regulated (dark shadow) and down-regulated (light shadow) isoforms are classified into novel (orange) and known (blue) isoforms based on the reference transcriptome in the corresponding species (NHP novel assemblies). (D) IL4R human- specific IU. Human-specific up-regulated isoforms are colored in red; nondifferentially used isoforms, in gray (left). IU values across species are shown in the heatmap (right), where colors represent row-scaled IU values (legend). regions of conserved usage (Supplemental Fig. S35). This suggests analyses (FDR ≤ 0.05; Methods) (Supplemental Table S10). An ex- that Alu insertions contribute to increasing orthologous exon in- ample of this is CENPQ, a centromere protein involved in autoim- clusion frequencies, especially for the case of SE. This is consistent mune disease and mitosis progression (Song et al. 2013; Bancroft with previous reports on the Alu-mediated modulation of exon us- et al. 2015), whose species-specific DEU is driven by a genetic age rates (Payer et al. 2019). change in the primate splice sites. In macaque, two junctions har- We found that genes with species-specific DEU are enriched bor GT–AG splice sites, leading to the recognition of two 5′ UTR ex- in functions related to the regulation of immunity, such as inter- ons and the preferred usage of the isoforms that include them. On feron signaling and innate immune response (IRF2, IRF3, IRF5, the contrary, great apes, in which these splice sites are not canon- IRF7, IRF9, IFI35, IFIT1, ISG20, IFNGR2, HLA-DQB1, CIITA, etc.), ical, only express the transcript that skips these exons (Fig. 5B). We as well as proliferative (DNA repair, centrosome localization, also noticed that 13.89% of genes displaying species-specific DEU etc.), housekeeping (ncRNA processing, mitochondrial gene ex- (221 out of 1591) co-occur with gene expression up- or down-reg- pression, etc.), and viral life-cycle control functions compared ulation in the same species. Six genes down-regulated in human with a background set consisting of all genes evaluated in DEU LCLs (SYNRG, DUSP14, TADA2A, ACACA, MYO19, and DDX52) 1454 Genome Research www.genome.org Downloaded from genome.cshlp.org on September 26, 2022 - Published by Cold Spring Harbor Laboratory Press Transcriptome innovations in primates A B C D Figure 5. Exon usage dynamics across primates. (A) Usage coefficients for skipping exons displaying species-specific changes (row-wise scaled). Overlap with coding sequence (CDS) and Pfam protein domains (UCSC), divisibility by three, average percentage-splice-in (PSI), and average phastCons scores are shown for each exon (see Supplemental Methods). (B) Example of gene (CENPQ) with different splice sites in macaque and great apes. Exon usage coef- ficients are illustrated at the top. Hominidae- and macaque-specific Iso-Seq transcripts are color-coded (bottom). (C ) Chr17q12 region with high density of human-specific DEU. Hi-C interaction intensity is indicated in the color legend. Gene orientation, DNase hypersensitive sites (DHSs), and topologically as- sociating domains (TADs) are also displayed. (D) Average conservation score in primates for exonic parts with conserved usage (N = 95,453), species-spe- cific up-regulation (N = 941) or down-regulation (N = 1340), and exonic parts showing other usage changes (N = 277). “Other changes” correspond to usage differences between groups of species (e.g., any two species show differential usage vs. the remaining three species). Exonic parts longer than five bases (all included as passed in 1000 Genomes hg38 strict mask) are shown. Statistical significance of the difference across groups was assessed by the Wilcoxon test, and shown P-values are adjusted by the Holm method. are comprised in the same topologically associating domain in forms (defined by APPRIS) (Rodriguez et al. 2018) compared with Chr17q12 (Rao et al. 2015) and display human-specific exon usage human-specific down-regulated or regions of conserved usage changes (Fig. 5C). Genomic rearrangements in this locus have (Supplemental Fig. S38). Taken together, these findings suggest been extensively studied given their evolutionary relevance and that species-specific gains in exon usage are associated with fast- implications in many diseases, including neurodevelopmental dis- evolving regions whose changes in usage lead to subtle or regula- orders (Mefford et al. 2007; Cardone et al. 2008; Moreno-De-Luca tory changes in the main function of these genes. et al. 2010; Nagamani et al. 2010). To date, this is the first report of We reasoned that genes with splicing changes detected across human-specific splicing changes occurring in multiple genes en- primates might show characteristic selective pressures in their cod- coded in this region. ing sequences. To investigate this possibility, we assessed the ratios We found that exonic parts experiencing species-specific us- of nonsynonymous to synonymous substitutions (dN/dS) in hu- age up-regulation display lower sequence conservation in primates man populations for the genes under evaluation (van der Lee than exonic parts down-regulated in a given species (Fig. 5D). This et al. 2017). To do so, we combined our results from isoform higher variation in species-specific up-regulated regions is also re- gains and losses and DEU and DIU analyses to define gene classes flected in the nucleotide diversity within human populations according to their splicing patterns: those showing any species- (Supplemental Methods; Supplemental Fig. S36). Furthermore, specific up-regulation signal in humans (N = 278 genes) or NHPs when focusing on transcribed regions with human-specific usage (N = 1239 genes), those displaying multiple events of species-spe- changes, we noticed up-regulated exonic parts are depleted in cific up-regulation in several species (N = 287 genes), those experi- CDSs (Supplemental Fig. S37) and are less present in principal iso- encing other splicing changes consistent with differences between Genome Research 1455 www.genome.org Downloaded from genome.cshlp.org on September 26, 2022 - Published by Cold Spring Harbor Laboratory Press Ferra ́ n dez-Peral et al. groups of species (N = 1174 genes), and genes with fully conserved mune system–related splicing programs in the context of recent splicing patterns (N = 4532 genes) (Supplemental Methods). We human evolution (Sams et al. 2016; Rotival et al. 2019). LCLs dis- observed significantly higher dN/dS in human populations for play transcriptional profiles that are very similar to those of plas- the gene sets showing splicing changes involving NHPs than for mablasts and early plasma cells (Mrozek-Gorska et al. 2019), genes with fully conserved splicing, whereas the genes with only allowing us to study the transcription of immune system–related human-specific up-regulation show intermediate dN/dS (Supple- genes. The full-length isoform characterization of our multispecies mental Fig. S39). These results suggest a link between dynamic LCLs panel showed a significant number of isoform gains and loss- splicing evolution in primates and accelerated or less constrained es in a context of widespread transcript expression conservation. coding changes in human populations, with intermediate selec- We found that genes involved in the regulation of innate immuni- tion constraints in genes with only human-specific up-regulation. ty are enriched in species-specific isoforms even when controlling To further investigate the occurrence of interspecies splicing for the composition of LCLs transcriptomes. For example, interfer- changes in genes under adaptive selection, we then asked if the on-response genes, which show early expression divergence in ac- above-defined gene classes showing splicing changes were en- tivated leukocytes from different primates (Hawash et al. 2021), riched in signals of positive selection in primates (van der Lee often produce species-specific isoforms in primate LCLs. This is et al. 2017) in comparison to genes with conserved splicing pat- also true for cytokine and chemokine genes, which display major terns. We found that only genes showing independent events of transcriptional differences in stimulated fibroblasts and phago- species-specific up-regulation in multiple primates are signifi- cytes from different mammals (Hagai et al. 2018). The concurrence cantly enriched in signals of positive selection (Fisher’s exact of global gene expression and splicing differences across closely re- test, Bonferroni-adjusted P = 0.014, OR = 2.90), indicating that lated species underlines the fast, multilayered transcriptome evo- they might be important targets of adaptive evolution. This set lution of genes involved in the immune response. of genes accumulating recurrent species-specific changes in the We also detected high-impact genetic changes located in splicing programs of different primates is enriched in genes associ- orthologous splice sites, which are clearly involved in shaping ated with the regulation of the immune response (Gene Ontology the splicing preferences of immune system–related genes in differ- Biological Process overrepresentation analysis [ORA], FDR = 0.03, ent primate species. Particularly, SE events are enriched in species- enrichment ratio = 2.17), reinforcing the link between evolution- specific transcript gains, posing this type of AS as a prevalent ary changes in splicing and the immune function. mechanism in the recent evolution of primate transcriptomes. Taking advantage of the most recent NHP genome assemblies, we have also characterized the emergence of new exons, many of Discussion them absent in reference annotations, and confirmed the associa- To date, our knowledge of the evolution of primate transcriptomes tion between the exonization process and the insertion of trans- is very limited. Here, we have generated the most extensive catalog posable elements such as SINEs and LINEs (Makałowski et al. of full-length isoforms across the primate lineage complemented 1994; Sorek et al. 2002; Cordaux and Batzer 2009; Schwartz et al. with matching Illumina RNA-seq data. The combination of both 2009; Avgan et al. 2019). allowed us to create highly curated transcriptomes consisting of The patterns of IU recapitulate the phylogenetic structure in a high proportion of novel transcripts from annotated genes, as primates, with the most abundant transcripts of each gene fre- well as a number of antisense and intergenic transcripts. Further- quently encoding the same ORF or showing conservative changes more, new AS events have been identified, and we have observed across species. This highlights that recent isoform changes usually a conserved composition of AS classes in great apes and rhesus ma- preserve the functional integrity of the protein encoded by the caque. This resource greatly expands the available gene annota- gene, representing a source of subtle regulatory fine-tuning for im- tions for model and nonmodel primates by providing a more portant genes. Here, we report a number of interesting species-spe- accurate and complete view of genes and isoforms, which are cru- cific changes in isoform preferences, involving a significant cial for further genomic and functional studies. fraction of isoforms absent in reference transcriptomes. Further In this study, we have taken advantage of our curated isoform analyses will be required to understand the functional implica- catalog to predict the corresponding ORF for coding transcripts, tions of the splicing changes affecting coding and noncoding representing the majority of the sequenced transcriptomes. This regions. In this direction, it is remarkable that our set of differen- ORF collection has been fundamental to build proteomic databas- tially used exons is enriched in UTRs, which harbor binding motifs es assisting coupled mass spectrometry analysis. In this way, we de- for RNA-binding proteins and small RNA molecules (e.g., tected evidence of translation for Iso-Seq transcripts identifying miRNAs), playing fundamental roles in the regulation of isoform peptides that are not present in reference proteomes. We detected stability and localization (Mayr 2019). fewer novelties at the proteome than at the transcriptome level, re- The enrichment in signals of positive selection in the coding flecting a combination of the differences in experimental sensitiv- sequences of genes accumulating independent species-specific ities, in the completeness of proteome and transcriptome gains of isoform expression, IU, or exon usage in multiple primates annotations, and in the proteome and transcriptome complexities suggests a link between recurrent splicing changes and adaptive (Ezkurdia et al. 2015). All considered, our analyses show that there evolution in primates. The immune response is highly dynamic, is still room for improvement of proteome annotations by integrat- and global gene expression differences have been previously re- ing isoform resolution sequencing data. ported between leukocytes of apes and African and Asian monkeys From an evolutionary point of view, many immunity-related (Hawash et al. 2021). These changes have also been seen when genes have undergone strong selective pressures in response to comparing human, chimpanzee, and rhesus macaque primary constantly evolving host–pathogen interactions, leading to the ac- monocytes (Barreiro et al. 2010) across multiple time points after celerated evolution of these genes (Barreiro and Quintana-Murci stimulation. Considering the functional specialization and dyna- 2010; Barreiro et al. 2010). However, very few studies have ad- mism of immune cells, characterizing interspecies isoform expres- dressed how these rapid evolutionary rates are reflected in the im- sion changes in different immune cell types and multiple time 1456 Genome Research www.genome.org Downloaded from genome.cshlp.org on September 26, 2022 - Published by Cold Spring Harbor Laboratory Press Transcriptome innovations in primates points after infection will be required to further understand this re- formation for the panel of LCLs is provided in Supplemental Table lationship. In the same vein, comparative analyses based on iso- S1. form sequencing in other tissues and cell types will help to Cells were cultured in RMPI-1640 with L-Glu (Gibco) supple- dissect changes in splicing preferences in primate-accelerated mented with 15% FBS (Gibco) and 1% P/S (10000 U/mL; Gibco). genes that are not expressed in LCLs, such as those involved in ner- Multiple cell pellets totaling 12 × 106 cells per cell line were stored vous system function and neurodevelopment (Dorus et al. 2004). at −80°C until further extraction after two washes with DPBS. The evaluation of interspecies splicing differences across develop- mental stages in primates is particularly interesting, because fast- RNA extraction evolving regions in apes have been shown to locate nearby regula- A pellet consisting of 2 × 106 cells was washed once and resuspend- tors of embryonic development (Kostka et al. 2018). ed in DPBS to use as starting material to extract DNA using the The Iso-Seq transcriptomes produced here are limited to a sin- DNeasy blood and tissue kit (Qiagen) following the manufacturer’s gle cell type and two biological replicates per species. In this regard, instructions. RNA was extracted from a 3 × 106 cell pellet, previous- we foresee that comprehensive studies based on long reads will fur- ly stored in QIAzol lysis reagent, using the miRNeasy mini kit ther extend our knowledge on the transcriptome complexity by as- (Qiagen) and following the manufacturer’s instructions. sessing multiple tissues from numerous species and individuals. The increasing availability of long-read transcriptomes will enable RNA-seq data production and processing the identification of tissue-specific and cell type–specific isoforms, their inter-individual variability, and evolutionary conservation. Illumina TruSeq stranded RNA-seq experiments were conducted at In our study, we addressed the unbalanced number of high-quality Macrogen for the following LCLs: GG05, CRL-1850, GM19238, CH507, CH170, CH322, EB185, and R94011, whereas the RNA- isoforms captured per species by leveraging deep short-read RNA- seq for the remaining LCLs (except for GM12878) was generated seq data (three biological replicates per species) to perform differ- by García-Pérez et al. (2021). RNA-seq data from GM12878 were re- ential expression and usage analyses based on our full-length tran- trieved from the NCBI Sequence Read Archive (SRA; https://www script models, revealing a similar number of expressed transcripts .ncbi.nlm.nih.gov/sra) under accession numbers SRR998197 and in all five species. However, the lack of saturation in Iso-Seq iso- SRR998198. RNA-seq reads were filtered with SOAPnuke (version form discovery rates implies that a number of species-specific tran- 1.5.6, parameters: -Q 2 -l 20 -q 0.5 -3 -E 60 -5 1) (Chen et al. scripts might be missed from our transcriptomes. Future advances 2018) after adapter removal. Specifically, 3′ adapter sequence of in sequencing platforms will allow direct comparisons of saturated fastq1 file (-f in SOAPnuke) was AGATCGGAAGAGCACACGTCT long-read transcriptomes, which will be essential for the detection GAACTCCAGTCAC, and 5′ adapter sequence of fastq2 file (-r) of interspecies differences affecting low-abundance transcripts was AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCT that are hard to infer from short-read RNA-seq alone. CGGTGGTCGCCGTATCATT. Filtered reads were then aligned In summary, this study expands our current knowledge on to each primate genome with STAR (version 2.7.0a, parameters: transcript and protein diversity in humans and NHPs, providing ‐‐readStrand Reverse ‐‐quantMode TranscriptomeSAM ‐‐outSAM a comprehensive evaluation of the dynamics of transcriptome evo- type BAM SortedByCoordinate) (Dobin et al. 2013). lution in the primate lineage as reflected in B cell–derived LCLs. The combination of high-depth long- and short-read transcrip- Isoform sequencing data production tomics and high-throughput mass spectrometry provides an excel- Two LCLs per species were used for isoform sequencing experi- lent framework to characterize many evolutionary novelties ments (PacBio Iso-Seq; human: GM12878, GM19150; chimpan- missed by reference annotations. The discovery of these recent zee: CH114, CH391; gorilla: DIAN, OMOYE; orangutan: two changes highlights the prevalence of splicing evolution in pri- isogenic cultures of PPY6; macaque: R05040, R02027). In the first mates. Our results show that genes involved in the innate immune run, three libraries of 1–2 kb, 2–3 kb, and 3–6 kb with poly(A) en- response, in cell proliferation, and with cell type–specific expres- richment (Oligo(dT)) were generated using SMRTbell template sion profiles have been targets of numerous isoform changes. prep kits and sequenced on a Pacific Biosciences (PacBio) sequel The strong adaptive demands of the immune response and cell system at Novogene. A total of 26 Gb (subread data) was generated proliferation at the transcriptional level are satisfied mostly by in- per species, with an average length and N50 of 1980 bp and corporating fast regulatory and subtle innovations without drasti- 2399 bp, respectively. In a second run, poly(A)-enriched libraries cally impacting the viability of the response. Moreover, this set of (Oligo(dT)) were produced without size selection using SMRTbell evolutionary changes constitutes the basis for future functional as- template prep kits and sequenced on a PacBio sequel system. In says for studying their role in interspecies phenotypic divergences this run, a total of 56 Gb (subread data) was generated per species, in the immune response and cell proliferation, areas of major inter- with an average length and N50 of 2162 bp and 2647 bp, est in evolutionary and biomedical research. respectively. Iso-Seq data processing Methods SMRT Link (version 6.0) was used to generate full-length polished transcript sequences from PacBio data by combining subreads LCLs origin and culture from the two LCLs per species. These transcript sequences were LCLs were kindly provided by Dr. Antoine Blancher (chimpanzee corrected with Illumina RNA-seq data using LoRDEC (Salmela LCLs CH507, CH170, and CH322), Dr. Aurora Ruiz (gorilla LCL and Rivals 2014). GMAP (Wu and Watanabe 2005) was used to GG05), Dr. Chris Tyler-Smith (orangutan LCL CRL-1850), and map the sequences to the corresponding reference genomes after Dr. Gaby Dioxiadis (macaque LCL R94011). GM19238 (human) poly(A) trimming. TranscriptClean (Wyman and Mortazavi and EB185(JC) (orangutan) were purchased from Coriell Institute 2019) was used to refine transcript splice junctions and reduce er- and Sigma-Aldrich, respectively. The remaining cell lines used in rors derived from Iso-Seq. Cupcake ToFU scripts were used to this article are described by García-Pérez et al. (2021). Detailed in- remove redundancy from high-quality, uniquely mapped Genome Research 1457 www.genome.org Downloaded from genome.cshlp.org on September 26, 2022 - Published by Cold Spring Harbor Laboratory Press Ferra ́ n dez-Peral et al. transcript sequences. Nonredundant transcripts were provided to Peptides were separated in a 90-min linear gradient from 5% to SQANTI for quality control and artifact filtering. Detailed informa- 90% acetonitrile in 10 mM ammonium bicarbonate (pH 8) at a tion on Iso-Seq data processing and artifact filtering is described in 0.8 mL/min flow rate over an Agilent Zorbax Extend-C18 (4.6 the Supplemental Methods. A comparison of our pipeline with an mm ID, 220 mm in length, 5-m particles, 80 Å pore size). alternative strategy involving independent processing of Iso-Seq Fractions were collected every minute into a total of 96 fractions, subreads from each LCL can be found in Supplemental Figures which were consolidated into 24, of which only 12 nonadjacent S40 through S42. samples were analyzed. Fraction volume was reduced to 300 µL ap- proximately and acidified by adding formic acid up to 10% of the final volume. Fractions were desalted with a MicroSpin C18 col- Protein extraction and precipitation umn (The Nest Group) and dried by vacuum centrifugation. Two million cell pellets were thawed on ice, washed within 1 mL of cold DPBS, and further centrifuged at 2500g 5 min at 4°C. After re- Species-specific exon gains (genomic structure) moving the supernatants, 450 µL of cold RIPA buffer (Thermo Fisher Scientific) supplemented with 10 µL/mL Halt protease and LiftOver and Liftoff (Shumate and Salzberg 2020) were used to phosphatase inhibitor cocktail (100; Thermo Fisher Scientific) identify conserved and species-specific exons at sequence level was added to disrupt cell pellets. After incubating 5 min with the (Iso-Seq exons derived from one-to-one orthologous genes). An tubes swirling on ice, the suspension was sonicated for 30 sec exon was defined as conserved when successfully mapped from with 50% pulse a total of three times with 10 sec of rest using NHP to human and then to all four NHPs by liftOver. An exon the Branson 250 sonicator. Later, tubes were further incubated was defined as species-specific when NHP exons failed to be for 30 min on ice swirling. Finally, tubes were centrifuged at mapped to the human reference genome or human exons failed 14,000g for 15 min at 4°C to pellet the cell debris, and supernatants to be mapped to all four NHPs by liftOver and achieved <50% were transferred into new protein LoBind tubes (Eppendorf). Liftoff alignment coverage in the remaining four assemblies. Extracted proteins were quantified using a Pierce BCA protein as- Additional information can be found in the Supplemental say kit (Thermo Fisher Scientific) and stored at −80°C until further Methods. analysis. To precipitate the proteins, six volumes of overnight chilled Transcript expression calculation and reconstruction acetone (−20°C) were added, and then tubes were inverted three of transcript gains and losses times and incubated overnight at −20°C. The next day, tubes RNA-seq of the corresponding species was used to quantify the were centrifuged at 16,000g for 10 min at 4°C, and acetone was re- projected transcript models in each genome by kallisto (‐‐rf-strand- moved carefully without disturbing the protein pellet. ed argument, default settings) (Bray et al. 2016). After filtering orthologous transcripts whose splicing structure is not supported Mass spectrometry sample preparation by RNA-seq and mono-exonic transcripts, Count (Csű ös 2010) was used to reconstruct the tree of transcript expression gains Three biological replicates for each species were processed. Samples and losses in primates considering the known phylogenetic struc- (12–120 µg) were solubilized in 100 mM triethylammonium bi- ture and the resulting binary matrix of transcript expression across carbonate (TEAB) at a concentration of 1 µg/µL, reduced with species. Two data sets were generated with different criteria: In the tris(2-carboxyethyl)phosphine (0.12–1.20 µmol, 60 min, 37°C), less conservative criteria, a transcript was defined as expressed in and alkylated in the dark with iodoacetamide (0.23–2.25 µmol, each species if TPM > 0 in at least two LCLs, whereas the strict 30 min, 25°C). The resulting protein extract was first diluted to 2 criteria required TPM > 0 in all LCLs from the same species M urea with 200 mM ammonium bicarbonate for digestion with (Supplemental Methods). Results of isoform expression across endoproteinase LysC (1:10 w:w, o/n, 37°C; Wako 129-02541) LCLs and functional enrichment of genes expressing species-spe- and then diluted twofold with 200 mM ammonium bicarbonate cific transcripts using TPM = 0.1 as an expression cutoff are shown for trypsin digestion (1:10 w:w, 8 h, 37°C; Promega V5113). in Supplemental Figures S43 and S44. Splice junction usage and After digestion, the peptide mixes were acidified with tri- isoform expression conservation estimated from Iso-Seq or RNA- fluoroacetic acid (TFA) to have a final concentration of 1% TFA, de- seq are illustrated in Supplemental Figure S45. Isoform expression salted with a MicroSpin C18 column (The Nest Group), and dried conservation was also analyzed in the context of independent Iso- by vacuum centrifugation. Peptide mixes were also solubilized Seq processing for each LCL (without combining subreads from with 100 mM TEAB at a concentration of 1 µg/µL. each species) (Supplemental Fig. S46). Principal component anal- Samples were labeled to generate three different TMT mixes; ysis of transcript TPM values before and after batch effect correc- in each TMT mix, one sample of each species was included with tion (Supplemental Methods) is depicted in Supplemental Figure a common pool consisting of a mix of all the samples. Labels S47. were used in a way that, in each TMT mix, each species uses a dif- ferent TMT reagent (see Supplemental Table S11). TMT-6 label re- agents were resuspended adding 41 µL of anhydrous acetonitrile dN/dS in human populations and positive selection in primates (ACN) to each tube. Two microliters of the resuspended TMT label Precomputed ratios of nonsynonymous to synonymous substitu- reagent was added to 10 µg of each sample. Additionally, 10 µL of tions (dN/dS) in human populations and genes undergoing positive anhydrous ACN and 8 µL of TEAB 100 mM were added to achieve selection in the primate lineage were retrieved from van der Lee approximately a final concentration of 30% of acetonitrile. Tubes et al. (2017). Briefly, dN/dS ratios were calculated using human var- were incubated for 1 h at room temperature, and then 2 µL of 5% iation data from 60,706 exomes collected by the Exome hydroxylamine was added and incubated to quench the reaction Aggregation Consortium (ExAc, release 0.3.1) (Lek et al. 2016), (15 min, 25°C). Samples were combined in three mixes as shown and positively selected genes were obtained through multiple in Supplemental Table S11, desalted with a MicroSpin C18 column alignments of orthologous protein sequences in nine primate spe- (The Nest Group), and dried by vacuum centrifugation. cies (human, chimpanzee, gorilla, orangutan, gibbon, macaque, TMT mixes were fractionated using basic pH reversed-phase baboon, vervet, and marmoset). Comparisons across gene classes fractionation (Isasa et al. 2015) in an Agilent 1200 system. (see Supplemental Methods, “Classification of genes according 1458 Genome Research www.genome.org Downloaded from genome.cshlp.org on September 26, 2022 - Published by Cold Spring Harbor Laboratory Press Transcriptome innovations in primates to their splicing and usage patterns”) were restricted to genes with CNGB Sequence Archive (CNSA; https://db.cngb.org/cnsa/) of a precomputed dN/dS and to background genes for positive selec- CNGBdb under accession number CNP0002604. Mass spectrome- tion analysis from van der Lee et al. (2017). try proteomics data have been submitted to the ProteomeXchange Consortium via the PRIDE (Perez-Riverol et al. 2022; https://www .ebi.ac.uk/pride/archive/) partner repository with the data set Functional enrichment analysis identifier PXD031213. In-house scripts for transcriptome data Functional enrichment was performed using ORA (minimum processing are available at GitHub (https://github.com/lfp1/ number of IDs in the category = 10, maximum number of IDs in isoseq_primate_LCLs/tree/XiaoyuZhan) and as a Supplemental the category = 500, BH FDR method, significance level FDR ≤ Code file. 0.05) in WebGestalt (Liao et al. 2019), based on Gene Ontology Biological Processes (The Gene Ontology Consortium et al. 2000; The Gene Ontology Consortium 2021) and PANTHER pathways Competing interest statement (human data) (Mi et al. 2021). For each functional enrichment, The authors declare no competing interests. we adjusted the background gene set to the orthologous genes ex- pressed in LCLs (detected by Iso-Seq) that were under evaluation (background gene sets are different for species-specific expression Acknowledgments gains, DIU, and DEU; see the different filtering criteria used in each analysis in Supplemental Methods). In addition, tissue-spe- We thank China National GeneBank for providing the computa- cific expression of genes producing species-specific transcripts tional resources. We thank Dr. Antoine Blancher, Dr. Aurora was tested using tissueEnrich (Jain and Tuteja 2019) against Ruiz, Dr. Chris Tyler-Smith, and Dr. Gaby Dioxiadis for providing GTEx RNA-seq expression data in 29 human tissues (The GTEx NHP LCLs, and Arianna Bruguera Reyes for the technical support Consortium 2015) with default settings (significance-level BH-ad- in cell line cultures. This work was supported by the Strategic justed P-value < 0.01). Background gene set was restricted to the list Priority Research Program of the Chinese Academy of Sciences of genes under evaluation in isoform expression gains/losses (XDB31020000); National Key R&D Program of China (China’s analysis. Ministry of Science and Technology [MoST]) grant 2018YFC We tested the enrichment and depletion in particular gene 1406901; the International Partnership Program of the Chinese lists related to the innate immune response, immune disease, Academy of Sciences (no. 152453KYSB20170002); the Carlsberg LCL-specific expression, and housekeeping genes. A manually cu- Foundation (CF16-0663); the Villum Foundation (no. 25900) to rated classification of innate immunity genes (eight subgroups: ac- G.Z.; and the La Caixa Foundation (ID 100010434) Fellowship cessory molecule, adaptor, effector, regulator, secondary receptor, Code LCF/BQ/DE16/11570011 (L.F.-P.). The Center for Genomic sensor, signal transducer, and transcription genes) was retrieved Regulation (CRG) / Universitat Pompeu Fabra (UPF) Proteomics from Deschamps et al. (2016). Immune disease genes correspond Unit is part of the Spanish Infrastructure for Omics Technologies to those provided by the International Union of Immunological (National Map of Unique Scientific and Technical Infrastructures Societies (IUIS; December 2020) (Tangye et al. 2020). Genes [ICTS] OmicsTech) and a member of the ProteoRed PRB3 Consor- with LCL-specific expression patterns were downloaded from tium, which is supported by grant PT17/0019 of the PE I + D + i Ryaboshapkina and Hammar (2019), and housekeeping genes 2013–2016 from the Instituto de Salud Carlos III (ISCIII), European were retrieved from the Housekeeping and Reference Transcript Regional Development Fund (ERDF), and “Secretaria d’Universi- Atlas (Hounkpe et al. 2021). tats i Recerca del Departament d’Economia i Coneixement de la Generalitat de Catalunya” (2017SGR595). T.M.-B. is supported by funding from the European Research Council (ERC) under Statistical analyses the European Union’s Horizon 2020 research and innovation pro- Gene set enrichment/depletion against customized gene lists (in- gramme (grant agreement no. 864203), BFU2017-86471-P nate immunity genes, immune disease-associated genes, genes (MINECO/FEDER, UE); “Unidad de Excelencia María de Maeztu,” with LCL-specific expression, and housekeeping genes) was tested funded by the Agencia Estatal de Investigación (AEI) (CEX2018- using the Fisher’s exact test (two-sided). For enrichment analyses 000792-M); Howard Hughes International Early Career; National involving multiple comparisons (Alu, CDS, and APPRIS principal Institutes of Health 1R01HG010898-01A1; and Secretaria d’Uni- transcript overlap with exonic parts, and positively selected genes versitats i Recerca and Centres de Recerca de Catalunya (CERCA) across splicing gene classes), the Fisher’s exact test (two-sided) with Programme del Departament d’Economia i Coneixement de la Bonferroni P-value correction was applied. IU differences between Generalitat de Catalunya (GRC 2017 SGR 880). annotated and novel isoforms were assessed using the Wilcoxon Author contributions: G.Z., T.M.-B., and D.J. designed the rank-sum test (two-sided). Pi diversity and phastCons score differ- study; L.F.-P. and X.Z. performed most of the data analyses; C.C., ences between different groups of exonic parts were assessed using Q.L., and P.E.-C. contributed to the data analysis; C.C., M.A.-E., the Wilcoxon test (two-sided) with P-value adjustment (Holm E.L., R.G.-P., E.J., E.S., and O.F. performed the experimental method). dN/dS differences between gene classes were evaluated work; L.F.-P., X.Z., and D.J. wrote the manuscript with input using the Dwass–Steel–Critchlow–Fligner all-pairs test (single- from all authors; and all authors read and approved the step P-value correction). Isoform/exon usage and global gene ex- manuscript. pression differences were tested using the statistical methods de- scribed in the DEXSeq and DESeq2 packages, respectively (Anders et al. 2012; Love et al. 2014). References Abascal F, Ezkurdia I, Rodriguez-Rivas J, Rodriguez JM, del Pozo A, Vázquez J, Valencia A, Tress ML. 2015. Alternatively spliced homologous exons have ancient origins and are highly expressed at the protein level. Data access PLoS Comput Biol 11: e1004325. doi:10.1371/journal.pcbi.1004325 Agosto LM, Gazzara MR, Radens CM, Sidoli S, Baeza J, Garcia BA, Lynch KW. The transcriptomic data generated in this study have been submit- 2019. Deep profiling and custom databases improve detection of proteo- ted to the NCBI BioProject database (https://www.ncbi.nlm.nih forms generated by alternative splicing. Genome Res 29: 2046–2055. .gov/bioproject/) under accession number PRJNA797478 and doi:10.1101/gr.248435.119 Genome Research 1459 www.genome.org Downloaded from genome.cshlp.org on September 26, 2022 - Published by Cold Spring Harbor Laboratory Press Ferra ́ n dez-Peral et al. Ahn E-Y, DeKelver RC, Lo M-C, Nguyen TA, Matsuura S, Boyapati A, Pandit vous system genes in the origin of Homo sapiens. Cell 119: 1027–1040. S, Fu X-D, Zhang D-E. 2011. SON controls cell-cycle progression by co- doi:10.1016/j.cell.2004.11.040 ordinated regulation of RNA splicing. Mol Cell 42: 185–198. doi:10 Dougherty ML, Underwood JG, Nelson BJ, Tseng E, Munson KM, Penn O, .1016/j.molcel.2011.03.014 Nowakowski TJ, Pollen AA, Eichler EE. 2018. Transcriptional fates of hu- Anders S, Reyes A, Huber W. 2012. Detecting differential usage of exons man-specific segmental duplications in brain. Genome Res 28: 1566– from RNA-seq data. Genome Res 22: 2008–2017. doi:10.1101/ 1576. doi:10.1101/gr.237610.118 gr.133744.111 Ezkurdia I, Rodriguez JM, Carrillo-de Santa Pau E, Vázquez J, Valencia A, Avgan N, Wang JI, Fernandez-Chamorro J, Weatheritt RJ. 2019. Tress ML. 2015. Most highly expressed protein-coding genes have a sin- Multilayered control of exon acquisition permits the emergence of nov- gle dominant isoform. J Proteome Res 14: 1880–1887. doi:10.1021/ el forms of regulatory control. Genome Biol 20: 141. doi:10.1186/ pr501286b s13059-019-1757-5 Feenstra B, Pasternak B, Geller F, Carstensen L, Wang T, Huang F, Eitson JL, Bancroft J, Auckland P, Samora CP, McAinsh AD. 2015. Chromosome con- Hollegaard MV, Svanström H, Vestergaard M, et al. 2014. Common var- gression is promoted by CENP-Q- and CENP-E-dependent pathways. J iants associated with general and MMR vaccine–related febrile seizures. Cell Sci 128: 171–184. doi:10.1242/jcs.163659 Nat Genet 46: 1274–1282. doi:10.1038/ng.3129 Barreiro LB, Quintana-Murci L. 2010. From evolutionary genetics to human García-Pérez R, Esteller-Cucala P, Mas G, Lobón I, Di Carlo V, Riera M, immunology: how selection shapes host defence genes. Nat Rev Genet Kuhlwilm M, Navarro A, Blancher A, Di Croce L, et al. 2021. 11: 17–30. doi:10.1038/nrg2698 Epigenomic profiling of primate lymphoblastoid cell lines reveals the Barreiro LB, Marioni JC, Blekhman R, Stephens M, Gilad Y. 2010. Functional evolutionary patterns of epigenetic activities in gene regulatory archi- comparison of innate immune signaling pathways in primates. PLoS tectures. Nat Commun 12: 3116. doi:10.1038/s41467-021-23397-1 Genet 6: e1001249. doi:10.1371/journal.pgen.1001249 The Gene Ontology Consortium. 2021. The Gene Ontology resource: en- Benjamini Y, Drai D, Elmer G, Kafkafi N, Golani I. 2001. Controlling the riching a GOld mine. Nucleic Acids Res 49: D325–D334. doi:10.1093/ false discovery rate in behavior genetics research. Behav Brain Res 125: nar/gkaa1113 279–284. doi:10.1016/S0166-4328(01)00297-2 The Gene Ontology Consortium, Ashburner M, Ball CA, Blake JA, Botstein Bian P, Ye C, Zheng X, Luo C, Yang J, Li M, Wang Y, Yang J, Zhou Y, Zhang F, D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, et al. 2000. et al. 2020. RIPK3 promotes JEV replication in neurons via downregula- Gene Ontology: tool for the unification of biology. Nat Genet 25: 25– tion of IFI44L. Front Microbiol 11: 368. doi:10.3389/fmicb.2020.00368 29. doi:10.1038/75556 Blekhman R, Marioni JC, Zumbo P, Stephens M, Gilad Y. 2010. Sex-specific The GTEx Consortium. 2015. The genotype-tissue expression (GTEx) pilot and lineage-specific alternative splicing in primates. Genome Res 20: analysis: multitissue gene regulation in humans. Science 348: 648–660. 180–189. doi:10.1101/gr.099226.109 doi:10.1126/science.1262110 Bray NL, Pimentel H, Melsted P, Pachter L. 2016. Near-optimal probabilistic Hadjadj J, Yatim N, Barnabei L, Corneau A, Boussier J, Smith N, Péré H, RNA-seq quantification. Nat Biotechnol 34: 525–527. doi:10.1038/nbt Charbit B, Bondet V, Chenevier-Gobeaux C, et al. 2020. Impaired type .3519 I interferon activity and inflammatory responses in severe COVID-19 Bundred JR, Hendrix E, Coleman ML. 2018. The emerging roles of ribosomal patients. Science 369: 718–724. doi:10.1126/science.abc6027 histidyl hydroxylases in cell biology, physiology and disease. Cell Mol Hagai T, Chen X, Miragaia RJ, Rostom R, Gomes T, Kunowska N, Henriksson Life Sci 75: 4093–4105. doi:10.1007/s00018-018-2903-z J, Park J-E, Proserpio V, Donati G, et al. 2018. Gene expression variability Bush SJ, Chen L, Tovar-Corona JM, Urrutia AO. 2017. Alternative splicing across cells and species shapes innate immunity. Nature 563: 197–202. and the evolution of phenotypic novelty. Philos Trans R Soc Lond B doi:10.1038/s41586-018-0657-2 Biol Sci 372: 20150474. doi:10.1098/rstb.2015.0474 Haralambieva IH, Ovsyannikova IG, Kennedy RB, Larrabee BR, Busse DC, Habgood-Coote D, Clare S, Brandt C, Bassano I, Kaforou M, Zimmermann MT, Grill DE, Schaid DJ, Poland GA. 2017. Genome- Herberg J, Levin M, Eléouët J-F, Kellam P, et al. 2020. Interferon-induced wide associations of CD46 and IFI44L genetic variants with neutralizing protein 44 and interferon-induced protein 44-like restrict replication of antibody response to measles vaccine. Hum Genet 136: 421–435. doi:10 respiratory syncytial virus. J Virol 94: e00297-20. doi:10.1128/JVI .1007/s00439-017-1768-9 .00297-20 Hawash MBF, Sanz-Remón J, Grenier J-C, Kohn J, Yotova V, Johnson Z, Byrne A, Cole C, Volden R, Vollmers C. 2019. Realizing the potential of full- Lanford RE, Brinkworth JF, Barreiro LB. 2021. Primate innate immune length transcriptome sequencing. Philos Trans R Soc Lond B Biol Sci 374: responses to bacterial and viral pathogens reveals an evolutionary 20190097. doi:10.1098/rstb.2019.0097 trade-off between strength and specificity. Proc Natl Acad Sci 118: Cardone MF, Jiang Z, D’Addabbo P, Archidiacono N, Rocchi M, Eichler EE, e2015855118. doi:10.1073/pnas.2015855118 Ventura M. 2008. Hominoid chromosomal rearrangements on 17q map He Y, Luo X, Zhou B, Hu T, Meng X, Audano PA, Kronenberg ZN, Eichler EE, to complex regions of segmental duplication. Genome Biol 9: R28. doi:10 Jin J, Guo Y, et al. 2019. Long-read assembly of the Chinese rhesus ma- .1186/gb-2008-9-2-r28 caque genome and identification of ape-specific structural variants. Nat Carlyle BC, Kitchen RR, Zhang J, Wilson RS, Lam TT, Rozowsky JS, Williams Commun 10: 4233. doi:10.1038/s41467-019-12174-w KR, Sestan N, Gerstein MB, Nairn AC. 2018. Isoform-level interpretation Herberg JA, Kaforou M, Wright VJ, Shailes H, Eleftherohorinou H, Hoggart of high-throughput proteomics data enabled by deep integration with CJ, Cebey-López M, Carter MJ, Janes VA, Gormley S, et al. 2016. RNA-seq. J Proteome Res 17: 3431–3444. doi:10.1021/acs.jproteome Diagnostic test accuracy of a two-transcript host RNA signature for dis- .8b00310 criminating bacterial vs. viral infection in febrile children. JAMA 316: Chen Y, Chen Y, Shi C, Huang Z, Zhang Y, Li S, Li Y, Ye J, Yu C, Li Z, et al. 835–845. doi:10.1001/jama.2016.11236 2018. SOAPnuke: a MapReduce acceleration-supported software for in- Hounkpe BW, Chenou F, de Lima F, De Paula EV. 2021. HRT atlas v1.0 da- tegrated quality control and preprocessing of high-throughput sequenc- tabase: redefining human and mouse housekeeping genes and candi- ing data. Gigascience 7: 1–6. doi:10.1093/gigascience/gix120 date reference transcripts by mining massive RNA-seq datasets. Nucleic Cordaux R, Batzer MA. 2009. The impact of retrotransposons on human ge- Acids Res 49: D947–D955. doi:10.1093/nar/gkaa609 nome evolution. Nat Rev Genet 10: 691–703. doi:10.1038/nrg2640 Hövelmeyer N, Wunderlich FT, Massoumi R, Jakobsen CG, Song J, Wörns Corvelo A, Eyras E. 2008. Exon creation and establishment in human genes. MA, Merkwirth C, Kovalenko A, Aumailley M, Strand D, et al. 2007. Genome Biol 9: R141. doi:10.1186/gb-2008-9-9-r141 Regulation of B cell homeostasis and activation by the tumor suppressor Csű ös M. 2010. Count: evolutionary analysis of phylogenetic profiles with gene CYLD. J Exp Med 204: 2615–2627. doi:10.1084/jem.20070318 parsimony and likelihood. Bioinformatics 26: 1910–1912. doi:10.1093/ Isasa M, Rose CM, Elsasser S, Navarrete-Perea J, Paulo JA, Finley DJ, Gygi SP. bioinformatics/btq315 2015. Multiplexed, proteome-wide protein expression profiling: yeast de Oyarzabal E, García-García L, Rangel-Escareño C, Ferreyra-Reyes L, deubiquitylating enzyme knockout strains. J Proteome Res 14: 5306– Orozco L, Herrera MT, Carranza C, Sada E, Juárez E, Ponce-de-León A, 5317. doi:10.1021/acs.jproteome.5b00802 et al. 2019. Expression of USP18 and IL2RA is increased in individuals Jain A, Tuteja G. 2019. TissueEnrich: tissue-specific gene enrichment anal- receiving latent tuberculosis treatment with isoniazid. J Immunol Res ysis. Bioinformatics 35: 1966–1967. doi:10.1093/bioinformatics/bty890 2019: 1297131. doi:10.1155/2019/1297131 Kim J-H, Shinde DN, Reijnders MRF, Hauser NS, Belmonte RL, Wilson GR, Deschamps M, Laval G, Fagny M, Itan Y, Abel L, Casanova J-L, Patin E, Bosch DGM, Bubulya PA, Shashi V, Petrovski S, et al. 2016. De novo mu- Quintana-Murci L. 2016. Genomic signatures of selective pressures tations in SON disrupt RNA splicing of genes essential for brain develop- and introgression from archaic hominins at human innate immunity ment and metabolism, causing an intellectual-disability syndrome. Am J genes. Am J Hum Genet 98: 5–21. doi:10.1016/j.ajhg.2015.11.014 Hum Genet 99: 711–719. doi:10.1016/j.ajhg.2016.06.029 Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Kostka D, Holloway AK, Pollard KS. 2018. Developmental loci harbor clus- Chaisson M, Gingeras TR. 2013. STAR: ultrafast universal RNA-seq align- ters of accelerated regions that evolved independently in ape lineages. er. Bioinformatics 29: 15–21. doi:10.1093/bioinformatics/bts635 Mol Biol Evol 35: 2034–2045. doi:10.1093/molbev/msy109 Dorus S, Vallender EJ, Evans PD, Anderson JR, Gilbert SL, Mahowald M, Kronenberg ZN, Fiddes IT, Gordon D, Murali S, Cantsilieris S, Meyerson OS, Wyckoff GJ, Malcom CM, Lahn BT. 2004. Accelerated evolution of ner- Underwood JG, Nelson BJ, Chaisson MJP, Dougherty ML, et al. 2018. 1460 Genome Research www.genome.org Downloaded from genome.cshlp.org on September 26, 2022 - Published by Cold Spring Harbor Laboratory Press Transcriptome innovations in primates High-resolution comparative analysis of great ape genomes. Science 360: Nagamani SCS, Erez A, Shen J, Li C, Roeder E, Cox S, Karaviti L, Pearson M, eaar6343. doi:10.1126/science.aar6343 Kang S-HL, Sahoo T, et al. 2010. Clinical spectrum associated with recur- Kuhlwilm M, Boeckx C. 2019. A catalog of single nucleotide changes distin- rent genomic rearrangements in Chromosome 17q12. Eur J Hum Genet guishing modern humans from archaic hominins. Sci Rep 9: 8463. 18: 278–284. doi:10.1038/ejhg.2009.174 doi:10.1038/s41598-019-44877-x Nelms K, Keegan AD, Zamorano J, Ryan JJ, Paul WE. 1999. The IL-4 receptor: Lau E, Han Y, Williams DR, Thomas CT, Shrestha R, Wu JC, Lam MPY. 2019. signaling mechanisms and biologic functions. Annu Rev Immunol 17: Splice-junction-based mapping of alternative isoforms in the human 701–738. doi:10.1146/annurev.immunol.17.1.701 proteome. Cell Rep 29: 3751–3765.e5. doi:10.1016/j.celrep.2019.11.026 Payer LM, Steranka JP, Ardeljan D, Walker J, Fitzgerald KC, Calabresi PA, Lek M, Karczewski KJ, Minikel EV, Samocha KE, Banks E, Fennell T, Cooper TA, Burns KH. 2019. Alu insertion variants alter mRNA splicing. O’Donnell-Luria AH, Ware JS, Hill AJ, Cummings BB, et al. 2016. Nucleic Acids Res 47: 421–431. doi:10.1093/nar/gky1086 Analysis of protein-coding genetic variation in 60,706 humans. Nature Perez-Riverol Y, Bai J, Bandla C, García-Seisdedos D, Hewapathirana S, 536: 285–291. doi:10.1038/nature19057 Kamatchinathan S, Kundu DJ, Prakash A, Frericks-Zipper A, Lempainen J, Korhonen LS, Kantojärvi K, Heinonen S, Toivonen L, Räty P, Eisenacher M, et al. 2022. The PRIDE database resources in 2022: a Ramilo O, Mejias A, Laine A-P, Vuorinen T, et al. 2021. Associations be- hub for mass spectrometry-based proteomics evidences. Nucleic Acids tween IFI44L gene variants and rates of respiratory tract infections dur- Res 50: D543–D552. doi:10.1093/nar/gkab1038 ing early childhood. J Infect Dis 223: 157–165. doi:10.1093/infdis/ Rao SSP, Huntley MH, Durand NC, Stamenova EK, Bochkov ID, Robinson jiaa341 JT, Sanborn AL, Machol I, Omer AD, Lander ES, et al. 2015. A 3D map Li Y, Zhang J, Wang C, Qiao W, Li Y, Tan J. 2021. IFI44L expression is regu- of the human genome at kilobase resolution reveals principles of chro- lated by IRF-1 and HIV-1. FEBS Open Bio 11: 105–113. doi:10.1002/ matin looping. Cell 162: 687–688. doi:10.1016/j.cell.2015.07.024 2211-5463.13030 Rodriguez JM, Rodriguez-Rivas J, Di Domenico T, Vázquez J, Valencia A, Liao Y, Wang J, Jaehnig EJ, Shi Z, Zhang B. 2019. WebGestalt 2019: gene set Tress ML. 2018. APPRIS 2017: principal isoforms for multiple gene analysis toolkit with revamped UIs and APIs. Nucleic Acids Res 47: sets. Nucleic Acids Res 46: D213–D217. doi:10.1093/nar/gkx997 W199–W205. doi:10.1093/nar/gkz401 Rotival M, Quach H, Quintana-Murci L. 2019. Defining the genetic and evo- Liu Y, Gonzàlez-Porta M, Santos S, Brazma A, Marioni JC, Aebersold R, lutionary architecture of alternative splicing in response to infection. Venkitaraman AR, Wickramasinghe VO. 2017. Impact of alternative Nat Commun 10: 1671. doi:10.1038/s41467-019-09689-7 splicing on the human proteome. Cell Rep 20: 1229–1241. doi:10 Ryaboshapkina M, Hammar M. 2019. Tissue-specific genes as an underuti- .1016/j.celrep.2017.07.025 lized resource in drug discovery. Sci Rep 9: 7233. doi:10.1038/s41598- Love MI, Huber W, Anders S. 2014. Moderated estimation of fold change 019-43829-9 and dispersion for RNA-seq data with DESeq2. Genome Biol 15: 550. Salmela L, Rivals E. 2014. LoRDEC: accurate and efficient long read error cor- doi:10.1186/s13059-014-0550-8 rection. Bioinformatics 30: 3506–3514. doi:10.1093/bioinformatics/ Makałowski W, Mitchell GA, Labuda D. 1994. Alu sequences in the coding btu538 regions of mRNA: a source of protein variability. Trends Genet 10: 188– Sams AJ, Dumaine A, Nédélec Y, Yotova V, Alfieri C, Tanner JE, Messer PW, 193. doi:10.1016/0168-9525(94)90254-2 Barreiro LB. 2016. Adaptively introgressed Neandertal haplotype at the Mao Y, Catacchio CR, Hillier LW, Porubsky D, Li R, Sulovari A, Fernandes JD, OAS locus functionally impacts innate immune responses in humans. Montinaro F, Gordon DS, Storer JM, et al. 2021. A high-quality bonobo Genome Biol 17: 246. doi:10.1186/s13059-016-1098-6 genome refines the analysis of hominid evolution. Nature 594: 77–81. Schwartz S, Gal-Mark N, Kfir N, Oren R, Kim E, Ast G. 2009. Alu exonization doi:10.1038/s41586-021-03519-x events reveal features required for precise recognition of exons by the Mathis BJ, Lai Y, Qu C, Janicki JS, Cui T. 2015. CYLD-mediated signaling splicing machinery. PLoS Comput Biol 5: e1000300. doi:10.1371/jour and diseases. Curr Drug Targets 16: 284–294. doi:10.2174/ nal.pcbi.1000300 1389450115666141024152421 Scotti MM, Swanson MS. 2016. RNA mis-splicing in disease. Nat Rev Genet Mayr C. 2019. What are 3′ UTRs doing? Cold Spring Harb Perspect Biol 11: 17: 19–32. doi:10.1038/nrg.2015.3 a034728. doi:10.1101/cshperspect.a034728 Shaath H, Vishnubalaji R, Elkord E, Alajez NM. 2020. Single-cell transcrip- Mazin P, Xiong J, Liu X, Yan Z, Zhang X, Li M, He L, Somel M, Yuan Y, tome analysis highlights a role for neutrophils and inflammatory mac- Phoebe Chen Y-P, et al. 2013. Widespread splicing changes in human rophages in the pathogenesis of severe COVID-19. Cells 9: 2374. doi:10 brain development and aging. Mol Syst Biol 9: 633. doi:10.1038/msb .3390/cells9112374 .2012.67 Sheynkman GM, Shortreed MR, Frey BL, Smith LM. 2013. Discovery and Mazin PV, Jiang X, Fu N, Han D, Guo M, Gelfand MS, Khaitovich P. 2018. mass spectrometric analysis of novel splice-junction peptides using Conservation, evolution, and regulation of splicing during prefrontal RNA-Seq. Mol Cell Proteomics 12: 2341–2353. doi:10.1074/mcp.O113 cortex development in humans, chimpanzees, and macaques. RNA .028142 24: 585–596. doi:10.1261/rna.064931.117 Shumate A, Salzberg SL. 2020. Liftoff: accurate mapping of gene annota- Mefford HC, Clauin S, Sharp AJ, Moller RS, Ullmann R, Kapur R, Pinkel D, tions. Bioinformatics 37: 1639–1643. doi:10.1093/bioinformatics/ Cooper GM, Ventura M, Ropers HH, et al. 2007. Recurrent reciprocal ge- btaa1016 nomic rearrangements of 17q12 are associated with renal disease, diabe- Soneson C, Matthes KL, Nowicka M, Law CW, Robinson MD. 2016. Isoform tes, and epilepsy. Am J Hum Genet 81: 1057–1069. doi:10.1086/522591 prefiltering improves performance of count-based methods for analysis Merkin JJ, Chen P, Alexis MS, Hautaniemi SK, Burge CB. 2015. Origins and of differential transcript usage. Genome Biol 17: 12. doi:10.1186/s13059- impacts of new mammalian exons. Cell Rep 10: 1992–2005. doi:10 015-0862-3 .1016/j.celrep.2015.02.058 Song G, Hu C, Zhu H, Wang L, Zhang F, Li Y, Wu L. 2013. New centromere Mi H, Ebert D, Muruganujan A, Mills C, Albou L-P, Mushayamaha T, autoantigens identified in systemic sclerosis using centromere protein Thomas PD. 2021. PANTHER version 16: a revised family classification, microarrays. J Rheumatol 40: 461–468. doi:10.3899/jrheum.120264 tree-based classification tool, enhancer regions and extensive API. Sorek R, Ast G, Graur D. 2002. Alu-containing exons are alternatively Nucleic Acids Res 49: D394–D403. doi:10.1093/nar/gkaa1106 spliced. Genome Res 12: 1060–1067. doi:10.1101/gr.229302 Mick E, Kamm J, Pisco AO, Ratnasiri K, Babik JM, Castañeda G, DeRisi JL, Srokowski CC, Masri J, Hövelmeyer N, Krembel AK, Tertilt C, Strand D, Detweiler AM, Hao SL, Kangelaris KN, et al. 2020. Upper airway gene ex- Mahnke K, Massoumi R, Waisman A, Schild H. 2009. Naturally occur- pression reveals suppressed immune responses to SARS-CoV-2 com- ring short splice variant of CYLD positively regulates dendritic cell func- pared with other respiratory viruses. Nat Commun 11: 5854. tion. Blood 113: 5891–5895. doi:10.1182/blood-2008-08-175489 doi:10.1038/s41467-020-19587-y Sun S-C. 2010. CYLD: a tumor suppressor deubiquitinase regulating NF-κB Mignone F, Gissi C, Liuni S, Pesole G. 2002. Untranslated regions of mRNAs. activation and diverse biological processes. Cell Death Differ 17: 25–34. Genome Biol 3: reviews0004. doi:10.1186/gb-2002-3-3-reviews0004 doi:10.1038/cdd.2009.43 Mittleman BE, Pott S, Warland S, Barr K, Cuevas C, Gilad Y. 2021. Sun CT, Lo WY, Wang IH, Lo YH, Shiou SR, Lai CK, Ting LP. 2001. Divergence in alternative polyadenylation contributes to gene regulato- Transcription repression of human hepatitis B virus genes by negative ry differences between humans and chimpanzees. eLife 10: e62548. regulatory element-binding protein/SON. J Biol Chem 276: 24059– doi:10.7554/eLife.62548 24067. doi:10.1074/jbc.M101330200 Moreno-De-Luca D, SGENE Consortium, Mulle JG, Simons Simplex Tangye SG, Al-Herz W, Bousfiha A, Chatila T, Cunningham-Rundles C, Collection Genetics Consortium, Kaminsky EB, Sanders SJ, GeneSTAR, Etzioni A, Franco JL, Holland SM, Klein C, Morio T, et al. 2020. Myers SM, Adam MP, Pakula AT, et al. 2010. Deletion 17q12 is a recur- Human inborn errors of immunity: 2019 update on the classification rent copy number variant that confers high risk of autism and schizo- from the International Union of Immunological Societies Expert phrenia. Am J Hum Genet 87: 618–630. doi:10.1016/j.ajhg.2010.10.004 Committee. J Clin Immunol 40: 24–64. doi:10.1007/s10875-019- Mrozek-Gorska P, Buschle A, Pich D, Schwarzmayr T, Fechtner R, Scialdone 00737-x A, Hammerschmidt W. 2019. Epstein–Barr virus reprograms human B Tardaguila M, de la Fuente L, Marti C, Pereira C, Pardo-Palacios FJ, Del Risco lymphocytes immediately in the prelatent phase of infection. Proc H, Ferrell M, Mellado M, Macchietto M, Verheggen K, et al. 2018. Natl Acad Sci 116: 16046–16055. doi:10.1073/pnas.1901314116 SQANTI: extensive characterization of long-read transcript sequences Genome Research 1461 www.genome.org Downloaded from genome.cshlp.org on September 26, 2022 - Published by Cold Spring Harbor Laboratory Press Ferra ́ n dez-Peral et al. for quality control in full-length transcriptome identification and quan- Wyman D, Mortazavi A. 2019. TranscriptClean: variant-aware tification. Genome Res 28: 396–411. doi:10.1101/gr.222976.117 correction of indels, mismatches and splice junctions in long-read Tokita MJ, Braxton AA, Shao Y, Lewis AM, Vincent M, Küry S, Besnard T, transcripts. Bioinformatics 35: 340–342. doi:10.1093/bioinformatics/ Isidor B, Latypova X, Bézieau S, et al. 2016. De novo truncating variants bty483 in SON cause intellectual disability, congenital malformations, and Xiong J, Jiang X, Ditsiou A, Gao Y, Sun J, Lowenstein ED, Huang S, failure to thrive. Am J Hum Genet 99: 720–727. doi:10.1016/j.ajhg Khaitovich P. 2018. Predominant patterns of splicing evolution on hu- .2016.06.035 man, chimpanzee and macaque evolutionary lineages. Hum Mol Genet Tress ML, Abascal F, Valencia A. 2017. Alternative splicing may not be the 27: 1474–1485. doi:10.1093/hmg/ddy058 key to proteome complexity. Trends Biochem Sci 42: 98–110. doi:10 Zhang XH-F, Chasin LA. 2006. Comparison of multiple vertebrate genomes .1016/j.tibs.2016.08.008 reveals the birth and evolution of human exons. Proc Natl Acad Sci 103: Unterman A, Sumida TS, Nouri N, Yan X, Zhao AY, Gasque V, Schupp JC, 13427–13432. doi:10.1073/pnas.0603042103 Asashima H, Liu Y, Cosme C Jr, et al. 2022. Single-cell multi-omics re- Zhang J, Kurgan L. 2019. SCRIBER: accurate and partner type-specific pre- veals dyssynchrony of the innate and adaptive immune system in pro- diction of protein-binding residues from proteins sequences. gressive COVID-19. Nat Commun 13: 440. doi:10.1038/s41467-021- Bioinformatics 35: i343–i353. doi:10.1093/bioinformatics/btz324 27716-4 Zhang S-J, Wang C, Yan S, Fu A, Luan X, Li Y, Sunny Shen Q, Zhong X, Chen van der Lee R, Wiel L, van Dam TJP, Huynen MA. 2017. Genome-scale J-Y, Wang X, et al. 2017. Isoform evolution in primates through inde- detection of positive selection in nine primates predicts human-virus pendent combination of alternative RNA processing events. Mol Biol evolutionary conflicts. Nucleic Acids Res 45: 10634–10648. doi:10 Evol 34: 2453–2468. doi:10.1093/molbev/msx212 .1093/nar/gkx704 Zhao M, Zhou Y, Zhu B, Wan M, Jiang T, Tan Q, Liu Y, Jiang J, Luo S, Tan Y, Warren WC, Harris RA, Haukness M, Fiddes IT, Murali SC, Fernandes J, et al. 2016. IFI44L promoter methylation as a blood biomarker for sys- Dishuck PC, Storer JM, Raveendran M, Hillier LW, et al. 2020. temic lupus erythematosus. Ann Rheum Dis 75: 1998–2006. doi:10 Sequence diversity analyses of an improved rhesus macaque genome en- .1136/annrheumdis-2015-208410 hance its biomedical utility. Science 370: eabc6617. doi:10.1126/science Zhao B, Katuwawala A, Oldfield CJ, Dunker AK, Faraggi E, Gsponer J, .abc6617 Kloczkowski A, Malhis N, Mirdita M, Obradovic Z, et al. 2021. Wu TD, Watanabe CK. 2005. GMAP: a genomic mapping and alignment DescribePROT: database of amino acid-level protein structure and function program for mRNA and EST sequences. Bioinformatics 21: 1859–1875. predictions. Nucleic Acids Res 49: D298–D308. doi:10.1093/nar/gkaa931 doi:10.1093/bioinformatics/bti310 Wu T, Ren M-X, Chen G-P, Jin Z-M, Wang G. 2016. Rrp15 affects cell cycle, proliferation, and apoptosis in NIH3T3 cells. FEBS Open Bio 6: 1085– 1092. doi:10.1002/2211-5463.12128 Received November 16, 2021; accepted in revised form July 12, 2022. 1462 Genome Research www.genome.org Downloaded from genome.cshlp.org on September 26, 2022 - Published by Cold Spring Harbor Laboratory Press Transcriptome innovations in primates revealed by single-molecule long-read sequencing Luis Ferrández-Peral, Xiaoyu Zhan, Marina Alvarez-Estape, et al. Genome Res. 2022 32: 1448-1462 originally published online July 15, 2022 Access the most recent version at doi:10.1101/gr.276395.121 Supplemental http://genome.cshlp.org/content/suppl/2022/08/10/gr.276395.121.DC1 Material References This article cites 105 articles, 24 of which can be accessed free at: http://genome.cshlp.org/content/32/8/1448.full.html#ref-list-1 Open Access Freely available online through the Genome Research Open Access option. Creative This article, published in Genome Research, is available under a Creative Commons Commons License (Attribution 4.0 International), as described at License http://creativecommons.org/licenses/by/4.0/. Email Alerting Receive free email alerts when new articles cite this article - sign up in the box at the Service top right corner of the article or click here. To subscribe to Genome Research go to: https://genome.cshlp.org/subscriptions © 2022 Ferrández-Peral et al.; Published by Cold Spring Harbor Laboratory Press