Widespread naturally variable human exons aid genetic interpretation

Domain: neuroscience
Published: 2025-12-16
Editorial outcome: major_revision

Abstract
Most mammalian genes undergo alternative splicing. The splicing of some exons has been acquired or lost in specific mammalian lineages, but differences in splicing within the human population are poorly understood. Using GTEx tissue transcriptomes from 838 individuals, we identified 57,271 “naturally variable exons” (NVEs) – exons which are included in mRNAs in some individuals but entirely excluded from others (or vice versa). NVEs impact three quarters of protein-coding genes, occur at all population frequencies, and are often absent from reference annotations. NVEs are more abundant in genes depleted of genetic loss-of-function mutations and aid in the interpretation of causal genetic variants. Genetic variants modulate the splicing of many NVEs, and 5’ untranslated region and coding-region NVEs are often associated with increased and decreased gene expression, respectively. Together, our findings characterize abundant splicing variation in the human population, with implications for a range of human genetic analyses.

Introduction
Nearly all human genes undergo alternative splicing, in which distinct mRNAs derived from different combinations of exons and splice sites are produced from the same gene1. Small changes in mRNA sequence can exert large effects, including the addition/deletion of protein domains2,3, or production of functional versus nonfunctional isoforms, and dysregulation of splicing contributes to human disease4,5. Understanding variation in alternative splicing is therefore important in understanding physiological and disease mechanisms.

Evolutionarily conserved alternative exons tend to preserve protein reading frame6, and often exhibit tissue-specific regulation7. In addition, comparison across mammalian species has identified thousands of cases of complete loss or gain of the splicing of exons8, with more recently evolved exons more likely to occur in 5’ UTRs, where they are often associated with increased gene expression8. Gain or loss of the splicing of an exon in the coding sequence (CDS) may also influence expression by changing the reading frame, often triggering nonsense-mediated mRNA decay (NMD)9,10. However, the evolution of alternative splicing within the human species remains poorly understood.

Most alternative exons are generally assumed to be spliced in all individuals, perhaps with different inclusion levels, a pattern which we call ‘canonical’ alternative splicing. Some such splicing differences are driven by genetic variation. Large-scale genetic association studies have associated changes in the usage of many alternatively spliced exons with single-nucleotide polymorphisms (SNPs), designated splicing quantitative trait loci (sQTLs)11,12. However, these variants typically have modest effect sizes, and their molecular mechanisms are often unclear. Ultra-rare genetic variants can sometimes lead to inclusion of previously unseen, “cryptic” exons or splice sites13. Some recent studies suggest that many exons are missing from transcriptome reference annotations, occurring only in a subset of transcriptomic datasets11,14,15. The availability of tissue transcriptomes from nearly 1000 individuals in the Genotype-Tissue Expression (GTEx) project presents an opportunity to characterize human variation at the exon level, asking how often the presence/absence of individual exons or splice sites varies between the transcriptomes of different individuals.

One key challenge in the study of unannotated transcripts is assessing which reflect technical artifacts or other types of noise16,17, which are reproducible, and which have phenotypic effects18. A cross-species analysis has suggested that many unannotated splicing events in longer-lived organisms may represent “transcript drift”, and are often nonfunctional until they begin to splice in at higher levels19. However, another recent study found that lowly used exons that trigger NMD broadly impact gene expression in human lymphoblastoid cell lines, and can help interpret complex trait loci10.

Here, we define “naturally variable exons” (NVEs) as exons whose splicing is variable within the human population (i.e., not present in every individual), focusing on alternative “skipped” or “cassette” exons, and alternative splice sites, which extend or truncate exons. We introduce a statistical approach that estimates the population frequency with which NVEs are meaningfully included (at levels above 5%) in a human cohort, generating a quantitative atlas of splicing variation across tissues. Using data from GTEx, we find that NVEs impact most human genes, are sometimes associated with changes in gene expression, and can aid in the interpretation of non-coding genetic variants.

Results
Generating a catalog of human NVEs
In order to generate a catalog of high-quality NVEs, we obtained mapped RNA-seq reads from 14,000 GTEx samples representing 49 tissues from 838 individuals11. To perform quality control and ensure accurate estimates of NVEs in the dataset, we devised a straightforward Bayesian approach. The method de-noises estimates by sharing signals across samples and enables estimation of the uncertainty in exon usage of each NVE.

To illustrate the method, consider an NVE of the alternative splice site variety, used in some individuals as an alternative to a canonical or “cognate” splice site, used in all people. The proportional use of the NVE within an individual – the “percent spliced in” (PSI or Ψ) – can be estimated from the proportion of RNA-seq reads that span the NVE and cognate exon junctions (EJs) (Supplementary Fig. 1A)20. NVEs tend to have fewer total junction reads per sample, making some estimates of Ψ less certain (Supplementary Fig. 1B). We therefore de-noised these estimates using a beta-binomial model, employing a mixture of 3 beta distributions for reasons discussed below. The workflow is described in Supplementary Fig. 1C, Methods and Supplementary Note Sections A-D. The method fits a smooth function representing the NVE’s underlying Ψ distribution across the population from the observed splice junction reads for an NVE across individuals in a tissue. We chose a mixture of 3 betas model because of its flexibility to model many different distributions and because many NVEs are associated with sQTLs (see below), which have 3 possible genotypes that may each generate distinct levels of splicing. The mixture of betas model can be used to assess uncertainty in Ψ values, and to estimate the proportion of individuals that splice the NVE at any given level. Here, we defined NVEs as alternative exons or splice sites that have estimated Ψ values of at least 5% in between 1% and 99% of the population in a particular tissue (see “Methods”); exons with Ψ ≥ 5% in over 99% of individuals are considered canonical alternative exons, and exons that don’t meet either of these criteria may be considered rare or cryptic exons. Using this definition, we identified 414,141 NVE-tissue pairs across all 49 tissues in the GTEx dataset (Supplementary Data 1), representing 57,271 unique NVEs. Like canonical alternative exons, NVEs are roughly balanced between the alternative splice site and skipped exon types (Supplementary Data 2).

Using these data, we define an interpretable summary statistic, exon frequency (EF), that describes the frequency of splicing of an NVE in a population, somewhat analogous to allele frequency (AF). The EF of an NVE in a tissue is the estimated proportion of individuals that splice the NVE at a threshold Ψ level or above. We typically employ a 5% Ψ threshold here because previous studies have suggested that this level of inclusion is near the lowest level where sequence conservation is commonly observed19. Distributions of EF values using different Ψ cutoffs are shown in Supplementary Fig. 1D

Methods
Ethics statement
No primary data were generated for this study. Person-related data were obtained through authorized access from primary data controllers.

Datasets
GTEx release v8. From the GTEx download data portal, we downloaded splicing “phenotype” files, which were originally used as inputs of the original sQTL study in GTEx (v8 phg38). We accessed it using the GTEx data browser, and you can also use the requester pay google cloud bucket (gs://gtex-resources/GTEx_Analysis_v8_QTLs). Protected data, such as genotypes and metadata, were accessed via dbGaP (study accession: phs000424.v8.p2).

GTEx fine-mapped sQTLs and eQTLs. Baberia et al. fine-mapped sQTLs and eQTLS in GTEx using the DAP-G fine-mapping method. We accessed this dataset using their Zenodo link. Note that fine-mapping was performed only on European ancestry samples from GTEx.

gnomAD. We accessed LOEUF scores from gnomAD v4.0 using the gnomAD data browser (linked here), which computed the LOEUF across nearly all genes in the genome.

GENCODE exon annotations. We accessed comprehensive transcript annotations from GENCODE v44 using their data browser, for use as our reference set of exons.

SpliceAI scores. SpliceAI scores have been computed across the genome (all possible SNP mutations). To compare the differences in SMS and spliceAI scores, we used a dataset that released all variants with SpliceAI scores above a nominal threshold (>0.1), which are available in the Illumina data browser (link: https://basespace.illumina.com/s/otSPW8hnhaZR). All other variants were set to 0.

Fine-mapped BioBank (GWAS) data
UKBB (96 traits): https://github.com/mkanai/finemapping-insights, including ASGR1 exon 4
Biobank Japan: https://pheweb.jp/downloads
FinnGen: https://www.finngen.fi/en/researchers/data_available

Data preprocessing
GTEx LeafCutter phenotypes. LeafCutter intron clusters were required to include only 2 or 3 introns, so we could classify them as alternative 3’ or 5’ splice sites or skipped exons, respectively. We mapped to genes and strands by scoring every splice site in the cluster, and picking the strand that had MaxEnt scores > 0 for all splice sites, excluding any clusters that did not score above zero across all splice sites on either strand as potential artefacts. Only protein-coding genes were included. Lastly, exons were required to be ≤ 500 bp in length and alternative splice sites were required to be within 500 bp of each other.

Generating estimates of Ψ values of NVEs
Many NVEs have fairly low inclusion levels, so to estimate their usage accurately, we built a model to assess their abundance in an RNA-sequencing sample. We considered a potential NVE in relation to a more frequently observed “cognate” junction in transcripts from the same gene.

PSI (Ψ) estimates the fraction of a gene’s transcripts that contain the exon or splice site of interest and is a widely used statistic to quantify splicing. It can be estimated from any RNA-seq dataset with sufficient read coverage of the alternative region:

ψ NVE = EJ reads NVE / (EJ reads NVE + EJ reads cognate)

We consider the exon or splice junction with fewer reads across all individuals a candidate NVE, and the exon or splice junction with more reads as the cognate.

Conditional on Ψ, the number of NVE reads in an RNA-seq dataset can be modeled as a binomial distribution, as defined by LeafCutter splicing33 clusters:

NVE reads | ψ ~ Binomial(EJ reads NVE + EJ reads cognate, Ψ NVE)

Note that we have to frame the exon junction reads slightly differently in alternative splice sites as compared to skipped exons to model the NVE (Supp. Fig. 1A).

Mixture of Betas binomial model
Partial pooling involves using information across samples to estimate effect sizes of individual samples by fitting a population distribution to the observed data and letting that inform effect size estimates for individuals54,55.

For proof of concept, the NVE population distribution can be modeled with a single beta distribution with parameters α and β:

ψ NVE ~ Beta(α, β)

Where α can be described as the pseudocount on the splicing of the NVEs splicing event, and β can be described as the pseudocount on the splicing of the cognate splicing event, respectively. In other words, in the absence of read count data, α and β provide an initial estimate of the splicing of the NVE.

A single beta distribution can capture several common situations (but not all – see Supplementary Note). One example of a distribution that a single beta cannot handle is a trimodal distribution, which would be expected for an NVE whose splicing differs between three subsets of the population. Because trimodal distributions are interesting biologically, potentially representing cases where a single variant (that may be absent or present in heterozygous or homozygous form) drives the inclusion of an NVE, we decided to fit a three-component mixture of betas model to all exons in the dataset. We compare this result to a single beta model in the Supplementary Note Sections E-F, and find that the summary statistic, EF, is fairly robust in both three- and one-component models.

For a three-component beta distribution, with parameters α and β, the NVE population distribution can be modeled with nine parameters:

ψ NVE ~ Σ i=1 to 3 γ i Beta(α i, β i)

Where γ i are weights that represent the fraction of the distribution coming from each individual beta distribution, requiring that all γ i ≥ 0, and that they sum to 1. Note that the triple-beta distribution can model a two-beta or single beta distribution by setting one or two of the γ parameters to zero. We can then estimate the population distribution of a given tissue sample:

PDF NVE = Σ i=1 to 3 γ i * PDF Beta(α i, β i)

The task is to predict the population distribution of a splicing event in a given tissue sample. We employ maximum likelihood estimation with the expectation maximization algorithm. The algorithm is initialized with uniform weights, α=(1,2,3), and β=(3,2,1) as an unbiased starting point able to capture symmetric and asymmetric distributions. It then iteratively performs the E-step to analytically find optimal weights γ, and the M-step to determine the α, β using the L-BFGS-B optimization algorithm. The method is represented graphically in Fig. 1A.

Alternative splice site NVEs
We searched for the optimal α, β, and γ parameters that best fit the data, using the expectation maximization framework (see Supplementary Note Sections B-C).

Skipped exon NVEs
Skipped exons are represented by three intron clusters: two introns supporting the inclusion of the exon, and one longer intron excluding the exon. We include events where both inclusion introns have numbers of reads that are not drastically different (see Supplementary Fig. 11), keeping one of the two inclusion introns to simplify the analysis. We then followed the same approach as for alternative splice sites to find the optimal parameters.

Calculation of exon frequency (EF) at a given Ψ threshold
We estimated the Ψ of all NVEs of length ≤ 500 (for skipped exons) or inter-splice site distance ≤ 500 bp (see Supp. Fig. 16), as this size capture the vast majority (>99%) of known internal exons and alternative splice sites in the human transcriptome. After obtaining the population parameters for all splicing events, we calculated the EF for a tissue at particular Ψ thresholds by using the population distribution of junction reads across individuals in the tissue. For each of four cutoff values c=1%, 5%, 10%, or 20%, we computed EF as the percentage of individuals with Ψ ≥ c from the CDF of the population distribution that was generated from the fitted model, and provide this value in a summary table (Supplementary Table 1). For a tissue with n individuals (almost always >100), we included events where the EF 5% was between 1% and 99%, excluding poorly spliced exons with EF below 1% and “canonical alternative” exons with EF > 99%.

Following ref. 56, to obtain an estimated Ψ value for an NVE in a single individual, we computed the mean posterior on Ψ and standard deviation of the estimate in each individual in the tissue by using the α and β parameters estimated from the 3-beta mixture binomial fit described above. To obtain a posterior mean on Ψ, a posterior on Ψ is computed using component priors is first estimated, and updated weights using the posterior estimates of α and β parameters are estimated. The weighted sum of these parameters for all components is the posterior mean on Ψ. That is

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this manuscript, Jacobs et al. discuss naturally variable exons (NVEs) – exons that vary in splicing across individuals – by leveraging the large-scale tissue data made available by the GTEx project. The authors introduce a statistical approach that captures uncertainty in exon inclusion and estimate 56,415 NVEs. Further, the study defines exon frequency (EF) as a metric to describe the prevalence of NVEs in a population and demonstrate that most NVEs generally have low EFs and occur in coding regions. The authors characterize the functional importance of NVEs by assessing gene tolerance to LoF, as well as effects on gene expression and GWAS interpretation. They consider the contribution of both genetic and non-genetic variables to NVE presence. Moreover, the study tackles the challenging problem of variant interpretation and yields valuable insights by means of their logistic regression model SMS. Overall, the manuscript presents novel approaches to understanding splicing regulation at the individual level and I believe that this study would be a significant contribution to the field. Please see the questions below for further discussion.\n1. How does the distribution of average total EJ reads per individual in Extended Data Fig. 1B compare to that of alternatively spliced non-NVEs. Also, please define top EF NVEs before mentioning for the first time in Extended Data Fig. 1B.\n2. For highly covered splice junctions, an absolute log 10 ratio just under 1 could still represent a significant bias in coverage for the two splice junctions flanking a skipped exon. How was this threshold chosen, and could this introduce some noise downstream?\n3. Given that exon skipping is generally more prevalent than alternative splice site usage, is the observed balance between the two types of splicing events expected among NVEs?\n4. NVEs tend to occur at low frequencies. Is it be more appropriate to represent the default EF value of an NVE with a measure of central tendency (e.g. median) rather than the top EF across all tissues, which may be an outlier?\n5. What is the exon length distribution of NVEs? Are high-frequency NVEs longer/shorter than low-frequency NVEs, and do they generally maintain/disrupt the reading frame?\n6. For NVEalt_ss, 5’ and 3’ splice site motif scores are weaker than their cognates. Is this pattern similar or different for skipped exon NVEs?\n7. Are there any similarities in biological pathways or functions among NVE-containing genes that share similar LEOUF scores?\n8. The correlation shown in Extended Data Fig. 6 seems to be driven by splicing outliers. What is the correlation coefficient? Is it significant? Also, do the authors observe this correlation in the heart and blood?\n9. In relation to the question above, can the authors provide support for their speculation of NMD being a potential mode of gene regulation by perhaps confirming the presence of a premature stop codon when the splicing of the ASGR1 NVE is altered?\n10. Regarding the contribution of genetics to the presence of NVEs, do individuals that share an NVE also share similar genetic variants? What is the distribution of sQTLs (Fig. 5A) around and within the NVEs?\n11. Most brain regions do not show virtually any correlation between age and NVEs. Interestingly, the substantia nigra has a relatively higher negative correlation. Can the authors comment or speculate on this pattern seen in the brain?\n12. Page 22: “The variant alters predicted affinity for two of these RBPs (HNRNPK, PCBP1) by 5-fold and 4-fold, respectively, highlighting a potential molecular mechanism for this variant”. How were the RBP affinities quantified? Are these inferred from the SMS model or are they obtained from other external tools?\n13. In relation to the previous question, do RBPs largely play a ubiquitous role in splicing of NVEs and non-NVEs (e.g. other AS exons), or do NVEs have stronger affinities for certain RBPs?\nMinor comments:\n14. Page 10: “The RNA-seq read depths available in GTEx are more than adequate to distinguish these moderate ψ values from absence of splicing, and NVEs possess splice site motifs very similar to those of other exons (Extended Data Fig. 2D), even for low-EF NVEs (Extended Data Fig. 2E).” Did the authors rather mean to reference Extended Data Figures 2E and 2F respectively?\n15. Please add a legend for Supplementary Figure 5", "reviewer_2": "The authors make the interesting observation that a substantial number of exons and splice sites are used (i.e., spliced) in only a subset of individuals. They denote such events as NVEs. They show this formally using a beta-binomial model that learns a per tissue, cross-individual Beta prior on the PSI, which is used to both denoise the per individual estimates and better estimate the frequency of the NVE across individuals. They investigate various properties of these NVEs (sequence motifs) and the genes (conservation) they are found in, and relate them to know splicing QTLs. They show fine-mapped GWAS variants are somewhat enriched in NVEs (slightly more so than synonymous variants). By classifying NVEs as likely NMD-causing vs not, they find genetically regulated nmdNVEs tend to result in the expected effect on expression (as shown by eQTL data). Finally they develop a sQTL predictor, \"SMS\", which, when including NVE splice sites in its features, performs substantially better than spliceAI.\nOverall, this is a interesting study with a number of interesting findings about interindividual variation in splicing in paricular and gene regulation in general. I have two major concerns.\n1. The beta binomial model assumes a beta distribution for the underlying (unknown) distribution of true PSIs across individuals. As they mention, this is somewhat flexible but will not model all distributions. Even one of the examples shown in the manuscript (CAST, fig 1c) appears to be tri-modal with modes close to 0, 0.6, and 1, a distribution that a single beta will not be able to capture. Mixtures of Gaussians are known to be able to approximate any smooth distribution, so using a mixture of betas would be a tractable and flexible alternative to the single beta. By summing over the possible mixture assignments it is still possible to use LBFGS to fit this model. Given that the beta-binomial model is a key component of the analysis, I feel it should be ascertained whether the results are consistent using such an approach. It would additionally allow analysis of how often the distribution is unimodal/bi/trimodal etc.\n2. The SMS part of the paper feels somewhat underdeveloped. There is no comparison to methods beyond SpliceAI (e.g. MTSplice or Pangolin), no mention of the sample size available for the logistic regression, and no assessment of which features are most important for the prediction. Unless I missed it it also wasn't specified whether an SMS is tissue specific or not? I'd almost be tempted to suggest cutting it, but the bump in AUPRC adding the NVE SSs is remarkable and would be a shame to miss.\nIt's not a major concern but the finding that NVEs are _more_ common in constrained genes warrants more discussion. Are the constrained genes more highly expressed making NVEs easily to detect? It might be worthwhile to ask if the NVE exons with highly constrained themselves are constrained as a class, compared to constitutive or canonical alternative exons (since testing individual exons would be underpowered).\nMinor comments:\n- \"This exon appears completely absent or very lowly included in most people, but occurs at a Ψ of 5-10% in about 10% of individuals, with credible intervals on Ψ well separated from zero\". This phrasing gives the impression there is bimodality, whereas really the data looks like a continuum to me.\n- \"This splicing variation does occur against\" cut \"does\"\n- Some amount of inter-individual variation in splicing for events differentially spliced between cell-types could be due to variable cell-type proportions. Cell type proportion estimates could be used to control for this possibility.\n- the legend in 2b is confusing.\n- several of the comparisons in \"Most NVEs occur in coding regions and have low EF\" would benefit from comparisons to a control set of exons.\n- label 3a&b with which direction is more constrained\n- what's going on with substantia nigra in 5b??\n- psi_NVE should be labeled differently in eqns 1 and 2 since one is an estimate and one is a parameter\n- \"Where (cid:1) can be described as the prior\"... alpha is a prior pseudocount, not a prior.\n- \"...using the L-BFGS-B algorithm (see Supplementary Note), using a beta prior with parameters α = β = 1\" not sure what this means.\n- \"keeping one of the two inclusion introns to simplify the analysis\" this is throwing out a lot of data. I think you could just take the average - the binomial likelihood in practice works fine for non integers.\n- \" we computed EF as the percentage of individuals with Ψ ≥ c from the modeled CDF of the population\" assuming the CDF means the CDF of the fit beta(alpha,beta) (which should be clarified) then this implies getting the distribution correct is v important for defining NVEs, and further motivates the mixture of betas direction\n- \"90% of above\" -> or\n- \"sQTL that we also eQTLs\"", "reviewer_3": "Jacobs et al. characterize inter-individual variation on splicing using the GTEx dataset and report many “naturally variable exons” unannotated in GENCODE. Jacobs et al. analyzes transcriptomes from the GTEx project to identify exons that are present in some individuals but completely absent in others that the authors call “naturally variable exons” (NVEs). The authors show that NVEs are frequently missing from reference annotations (gencode) and that they are more common in genes that are less prone to loss-of-function mutations. Many NVEs are influenced by genetic variants and could help interpret causal genetic variants.\nOverall, we find this work relevant as it reveals studies splicing variation in humans from a different perspective and offers significant implications for genetic research. We have some concerns/suggestions that we would like the authors to consider.\nMajor comments:\nSome of the conclusions are not driven by statistically significant results, only by trends in the data, most likely due to power issues. However, these conclusions should be toned down and statistics should be added to some plots. In particular:\nIn Figure 2B (page 10) the authors mention: “NVEs with low EF are predominantly unannotated, while around half of high-EF NVEs are annotated”. However, in the plot the proportions do not seem different, all bars seem to be around 60%. The text could be improved by adding the exact percentages in the text or in the plotting area instead and if the percentages are not different enough, tone down the message of finding differences between low EF and high EF NVEs.\nThe authors report a correlation for Extended Data Fig. 6. However, a correlation coefficient, and a p-value (independently of it being significant or not) should be added to the plot to help readers interpret the correlation.\nIn Figure 3c, authors should explain what they refer to as “Enrichment”, similar to the explanation in the supplementary figure. Importantly, they report “The observed enrichments of causal variants in NVEs are above the level of synonymous variants, particularly at high PIP values, but below those of missense annotations”. According to the figure, the trend is not “particularly” but “only” observed at high PIP values, and the confidence intervals completely overlap, so it is not a significant result. The authors should tone down the results and state that “tend to be above” instead of “are above”. Furthermore, the sentence “Unannotated splicing events are enriched above synonymous at high PIPs” should be removed from the figure captions.\nFigure 6A would benefit from having the (AUPRC) of the three curves, especially to support the claim: “When given access only to reference gene annotations, where the distance feature uses GENCODE splice sites (SMS-GENCODE), the method’s performance improves moderately over SpliceAI”. If not on the plot, the AUPRC should be stated in the text.\nThe Analysis shown in figure 5B should be redone, correcting for batch effects to confirm for age being the main explanatory factor for more NVEs. Batch effects such as ischemic time and hardy scale are a very important source of variation in the GTEx dataset.\nIn the first section of the results, the authors need to describe the methodology / cutoffs of actually calling NVEs in greater detail so that readers do not have to read the methods text to get a general idea of what cutoffs were used.\nIn the introduction, the authors mention evolution of new exons and draw some connections to why NVEs might represent evolutionarily new exons. This connection could be explored by examining sequence conservation characteristics of NVEs, likely as part of the “Most NVEs occur in coding regions and have low EF” section.\nIn the sentence “NVEs rarely occur in 3' UTRs, likely because they contain very few introns relative to other transcript regions”, the authors should elaborate to mention NMD, as this is the reason that 3’ UTRs do not contain introns.\nIn the sentence “Such a pattern could occur if NVEs arise more commonly in more constrained genes but often experience selection favoring lower ψ values, driving EFs lower, although other scenarios are possible”, the authors should elaborate on other scenarios. To me, what they propose seems possible but as it is the only scenario presented it is hard to assess. In connection with this, I think the hypothesis that the authors present is really interesting, but it would be strengthened through complementary analyses. For instance, the authors could examine the distribution of the metrics used to call NVEs of exons stratified by high or low LOEUF genes. Maybe this isn’t the best way to get an answer but it is what I can think of.\nRegarding the analysis of the relationship of age versus number of NVEs, it would be interesting to see if the effect of age on NVEs could be correlated to the somatic mutation rate in the different GTEx tissues (https://doi.org/10.1186/s13073-020-00746-1) with the hypothesis being that the greater the mutation rate, the more variable splicing will be due to the genetic effects of somatic mutations.\nMinor comments:\nSupplementary Figure 5E figure caption is missing.\nAuthors should be consistent in the naming of the supplementary note figures. Sometimes they are referred to as “Supplementary Fig. 3” and sometimes as “Supplementary Note Fig. 5”.\nSometimes the authors use “an NVE” and some other “a NVE”. For example: “that describes the frequency of splicing of an NVE in a population, somewhat analogous to allele frequency (AF). The EF of a NVE…” in page 5.\nMention what NMD stands for the first time you mention it\nThe title of method section “sQTL that we also eQTLs” is unclear\nReferences are missing in the paragraph: “The identification of variants that impact splicing remains challenging…” in page 18\nIn page 21: “When training on the top splicing signal (top PIP), i.e. the best splicing association for each unique variant, we obtained qualitatively similar results and the same ranking of algorithms (data not shown)”. Even if the data is not shown, the (AUPRC) could be reported in the text.\nRegarding the statement “recent studies suggest that many exons are missing from transcriptome reference annotations, occurring only in a subset of transcriptomic datasets11,14,15 “. I don’t think GTEx is the paper for this statement.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "1\nResponses to Reviewer Comments\nNote: Figures made specifically for the reviewers are labeled and numbered\nsequentially as Reviewer Figures. Other figures present in the main text or supplement\nare numbered as in the revised text.\nReviewer 1\nR1 Major comments\n1. How does the distribution of average total EJ reads per individual in Extended Data\nFig. 1B compared to that of alternatively spliced non-NVEs. Also, please define top EF\nNVEs before mentioning for the first time in Extended Data Fig. 1B.\nThe mean and median EJ read count for NVEs (13 and 43, respectively) are somewhat\nlower than those for non-NVE alternative exons (16 and 66, respectively), as expected\n(Review Figure 1). We have explained top EF NVEs when first mentioned, as\nrequested. We also add the new non-NVE curve to Extended Data Fig. 1B.\nReview Figure 1. Distribution of non-alternative NVE reads compared to NVE reads. The\naverage number of reads across all individuals was calculated by taking the total reads\nspanning all junctions (including cognate junctions) divided by the number of people in the\ntissue. The average number of reads are shown for both NVE events and non-NVE alternative\nsplicing events (100% EF). The p-value of KS test between two distributions is shown at top.\n2. For highly covered splice junctions, an absolute log 10 ratio just under 1 could still\nrepresent a significant bias in coverage for the two splice junctions flanking a skipped\nexon. How was this threshold chosen, and could this introduce some noise\ndownstream?\nThe counts of the two splice junctions flanking a skipped exon can vary for a number of\nreasons. For example: 1) sampling variation (more prominent when total read counts\nare low); 2) the RNA secondary structure potential of the region may differentially\nimpact priming in these two regions during library prep; 3) the base composition of\n2\nthese regions may differ, resulting in differences in the amplification and sequencing\nsteps; or 4) differences in mappability, e.g., if the skipped exon has very similar\nsequence to the downstream exon, then upstream junctions and exon skipping\njunctions may be so similar as to prevent some upstream (but not downstream) reads\nfrom being uniquely mappable (and conversely if the skipped exon is similar in\nsequence to the upstream exon). Thus, applying a stringent threshold on this ratio will\nexclude many true skipped exons. We assessed several mean log ratio thresholds. We\nchose a threshold of 1 based on the distributions and scatterplots in Supplementary\nFigure 1 and because the correlation between log a and b counts was still high, while\nretaining a large portion of the data. Other downstream filters including filtering on\nmaxent score, filtering on protein coding genes, and setting a maximum exon size also\nhelp to reduce misclassification of skipped exons. We have added a command line\noption to modify the log ratio threshold, if desired.\nunskipped junction of skipped exons (a and b referred to above) in an arbitrarily\nselected tissue: Ovary. Blue line indicates the range around 0 that passed our filter\n(absolute mean log ratio < 1). B) The read counts of introns a and b for NVEs that\npassed filtering in arbitrarily selected tissue, Ovary (Pearson and spearman correlation\n0.97, 0.95 respectively). C) Selection of log 10 threshold of 1. Different log 10 ratio\ncutoffs [0.1, 0.2, 0.5, 1, 2, 3] were assessed, varying on stringency. C) Shown is\ndistribution of reads same as A and different assessed log 10 ratio cutoffs, plotted\ncutoffs as lines to indicate where in the data that filter would throw away anything\nbeyond that line. D) Proportion of exons left after the filter E) correlation between log of\na counts and log of b counts for all introns that pass filter, colored by stringency of filter.\nF) sum across all reads for a and b counts across all individuals in the tissue.\n4\n3. Given that exon skipping is generally more prevalent than alternative splice site\nusage, is the observed balance between the two types of splicing events expected\namong NVEs?\nThe observed balance between skipped exons and alternative splice sites is similar for\nboth NVEs and non-NVE alternative exons, with roughly twice as many skipped exons\nas alternative splice sites of either type (Review Figure 2).\nReview Figure 2 Proportions of splicing event types for NVEs (purple) and non NVEs\n(grey). Showing top EF NVEs.\n4. NVEs tend to occur at low frequencies. Is it more appropriate to represent the default\nEF value of an NVE with a measure of central tendency (e.g. median) rather than the\ntop EF across all tissues, which may be an outlier?\nSince many low EF exons are tissue specific (Extended Data Fig. 5D), the median EF\nwould be 0 for many such exons, collapsing the variability we observe when looking at\ntissue-specific values. This tendency toward very low EFs when using the median is\napparent in Review Figure 3, which compares median to top EF and other possible EF\ndefinitions. Top EF has a similar distribution to that defined at the 90th or 95th\npercentiles, indicating that it does not behave like an outlier, which we mention in the\ntext.\ntissues. No. of NVEs plotted in each bin, with the EF taken for each NVE using the\nA) median across all tissues that contain the NVE. B) the 90th percentile of EF across tissues.\nC) 95th percentile of EF across tissues. D) the maximum of EF across tissues. Note that D is\nthe EF spectrum shown in the main text.\n5. What is the exon length distribution of NVEs? Are high-frequency NVEs\nlonger/shorter than low-frequency NVEs, and do they generally maintain/disrupt the\nreading frame?\nAn interesting question. We observed that the length distribution of NVEs varies\nsomewhat between the types of splice sites (Extended Data Fig. 6A-D). The smallest\nexons on average are A3SS, driven by NAGNAGs, which are 3 bp in length and are\ncommon (PMID: 22235189, 15516930). We find that the differences in length\ndistribution of A5SS NVE lengths are impacted by position within the gene, as\nalternative 5SS NVEs located in 5' UTRs tend to be longer than those located in other\ngene regions sites (Extended Data Fig. 6E) This trend could potentially be related to\nthe frequent presence of nearby alternative promoters in 5' UTRs, which could drive use\nof different alternative 5SS at distances a few hundred bases apart, compared to typical\nalternative 5SS located closer together, though we did not investigate this further.\nThese observations are now mentioned in the Results.\nA-D) shown split by NVE splicing event: the distribution of exon sizes across the EF spectrum.\nKolmogorov–Smirnov (KS) tests were between low and high EF NVEs and p-values are shown\nabove. In C), denoting where a likely “NAGNAG” signature is seen, which a higher percentage\n(~20%) of the data having 3bp 3’ splice sites.\nE) shown split by whether 5p splice site NVE is in the 5’UTR or not: the distribution of exon\nsizes in those two bins. Shown on right is a schematic detailing the two types of NVEs shown in\nthe plot.\nWe find that higher frequency NVEs are less likely to have NMD potential than low\nfrequency NVEs for all NVE types for alternative splice sites, but less so for skipped\nexons (Review Figure 4). This trend likely reflects a combination of factors. For\nexample, NMD-sensitive exons may more often have individual-specific PSI values\nbelow 5% than NMD-insensitive exons because of the action of NMD, thus reducing\ntheir measured EFs. Furthermore, frame-preserving exons may more often retain\nprotein function, and therefore natural selection may more often allow higher PSI and\nEF values for such exons than for frame-disrupting exons.\nThe proportion of NMD potential events across the EF spectrum, grouped by NVE type. Exons\nhad NMD potential if they caused a frameshift, or they had a stop codon in their frame. See\nmethod for detailed description of how we defined exons with NMD potential. Z-test of\nproportions were performed between low and high EF NVEs and p-values are shown above.\n6. For NVE alt_ss, 5’ and 3’ splice site motif scores are weaker than their cognates. Is\nthis pattern similar or different for skipped exon NVEs?\nExtended Data Fig. 5 shows that the NVE 5' and 3' SS tend to be weaker than cognate\n5' and 3'SS. Former Supplemental Figure 5 (now 10) shows that “included NVEs” (PSI >\n5% in some people, absent in others) tend to have weaker SS than their flanking\ncognate exons, as do “excluded NVEs” (PSI < 95% in some people, fully included in\nothers). We apologize that this figure legend was cropped in PDF conversion.\nSupplementary Figure 9. Splice site scores of NVEs and cognate exons. Panels C, D only.\nC-D) Skipped exon scores, stratified on whether an NVE is excluded or included.\n8\n7. Are there any similarities in biological pathways or functions among NVE-containing\ngenes that share similar LEOUF scores?\nWe analyzed sets of NVE-containing constrained genes (LOEUF NVE decile 10%) and\nNVE-containing unconstrained genes (LOEUF NVE decile 90%). We performed GO\nanalyses between NVE-containing genes and all constrained genes (LOEUF decile\n10%) and found the following pathways enriched in NVE-containing constrained genes:\nhttps://cbl-gorilla.cs.technion.ac.il/GOrilla/nhr9a6a2/GOResults.html\nWe performed GO analyses between NVE-containing unconstrained genes and all\nunconstrained genes (LOEUF decile 90%) and found the following pathways enriched in\nNVE-containing unconstrained genes: https://cbl-\ngorilla.cs.technion.ac.il/GOrilla/gk6afl36/GOResults.html\nFor both sets, the level of enrichment for particular subsets of genes was relatively\nmodest, in the ~1.1-fold to at most about 1.3-fold range. We suspect that NVEs can\narise in all genes at some frequency and did not feel that any category was sufficiently\nenriched to draw a strong biological conclusion.\n8. The correlation shown in Extended Data Fig. 6 seems to be driven by splicing\noutliers. What is the correlation coefficient? Is it significant? Also, do the authors\nobserve this correlation in the heart and blood?\nIt seems that some information was cut off this figure in the initial submission in PDF\nconversion, and we apologize for that. We used Spearman correlations, which are\nrelatively insensitive to outliers. P-values are now shown in the plot (Extended Data Fig\n6F). This analysis was not tissue-specific. We wanted to compute a global NVE rate, as\nmany NVEs are tissue specific, and so a tissue specific NVE right might underestimate\nthe number of NVEs observed across a gene. We reflected this by editing the x-axis to\nsay “global NVE rate” as supposed to the previous “NVE rate” and describe this term in\nthe legend.\nExtended Data Figure 6F. Size distributions of NVEs and LOEUF score analysis\nSpearman correlation between global NVE rate (number of NVEs per gene considering all\ntissues, divided by the number of exons in reference transcriptome, -log pval > 90) and the\nbecause some NVEs are specific to particular tissues.\n9\n9. In relation to the question above, can the authors provide support for their speculation\nof NMD being a potential mode of gene regulation by perhaps confirming the presence\nof a premature stop codon when the splicing of the ASGR1 NVE is altered?\nThis exon was characterized in another study. In the example studied, a frameshift\ncreated by a deletion variant near the exon boundary resulted in reduced expression at\nmRNA and protein levels, consistent with NMD (PMID: 27192541). This paper is now\ncited in the main text.\n10. Regarding the contribution of genetics to the presence of NVEs, do individuals that\nshare an NVE also share similar genetic variants? What is the distribution of sQTLs\n(Fig. 5A) around and within the NVEs?\nThe presence of sQTLs associated with a majority of NVEs (Figure 5A) indicates that\npeople who share an NVE very often share a genetic variant associated with the\nsplicing of the NVE. sQTL variants often occur very close to the involved splice sites\n(Extended Data Fig. 8A). Nearly half of the high PIP variants in GTEx are within 60 bp\nof the splice site. These observations are consistent with extensive literature indicating\nthat splicing regulatory elements occur more frequently and have higher activity when\nlocated near splice sites (PMID 25784052, 18369186).\nExtended Data Figure 8A. Details on SMS.\nPlotting the distance (in bp) between the variant in given QTL to the nearest splice site in the\nsplicing phenotype associated with the variant. Separated by which s-QTL PIP range the variant\nis in. Inset plot shows the range from 0-50 bp and demonstrates that half of all high PIP (90% or\nabove) s-QTL variants tested are within 50 bp of the nearest splice site. K.S. test between\nhighest and lowest PIP range is shown.\n10\n11. Most brain regions do not show virtually any correlation between age and NVEs.\nInterestingly, the substantia nigra has a relatively higher negative correlation. Can the\nauthors comment or speculate on this pattern seen in the brain?\nWe completed a new analysis for effects of age using a generalized linear model to\nbetter control for certain technical factors (see response to reviewer 3 comment below).\nFor some brain regions, NVE counts increase with age, while for others NVE counts\ndecrease. Unfortunately, most brain regions are quite poorly sampled at younger ages\n(see Review Figure 8 below), reducing our confidence in these trends. Because of the\npotential for technical confounders uncovered in these new analyses and the relatively\nmodest effect sizes and variability in direction, we are now less confident in the aging\nresults and decided to remove them from the manuscript. It will be interesting to revisit\nthis question in future, particularly using datasets with better sampling of younger\nindividuals.\n12. Page 22: “The variant alters predicted affinity for two of these RBPs (HNRNPK,\nPCBP1) by 5-fold and 4-fold, respectively, highlighting a potential molecular mechanism\nfor this variant”. How were the RBP affinities quantified? Are these inferred from the\nSMS model or are they obtained from other external tools?\nWe used RBPamp, a method developed in our lab to estimate affinities from RNA Bind-\nn-Seq data generated at different protein concentrations:\nhttps://www.biorxiv.org/content/10.1101/2022.11.08.515616v1. We have clarified this in\nthe Methods section.\n13. In relation to the previous question, do RBPs largely play a ubiquitous role in\nsplicing of NVEs and non-NVEs (e.g. other AS exons), or do NVEs have stronger\naffinities for certain RBPs?\nTo address this question, we leveraged KATMAP\n(https://www.biorxiv.org/content/10.1101/2024.06.25.600605v2), a computational\nmethod that learns distance-dependent splicing factor activity from the combination of a\nbinding model derived from in vitro data and splicing changes observed in\nknockdown/RNA-seq data. It also calculates the “splicing impact” of each factor on\nevery exon, which is an estimate of the extent to which a change in the given factor’s\nactivity would directly perturb the splicing of the exon. We analyzed all NVE and non-\nNVE skipped exons (EF above 0.99) in GTEx, also separating into groups by top EF.\nWe additionally compared to “simulated novel exons” generated by sampling a random\nintronic region of similar size within the same gene as each NVE. For the 18 splicing\nfactors analyzed by KATMAP, following differed significantly between simulated exons\nand low EF NVEs: PCBP1, HNRNPK, HNRNPC, SRSF1, PUF60, and U2AF2. As a\ngeneral rule, positive-acting factors such as PCBP1, PUF60 and U2AF2 tended to have\nhigher splicing impacts in all subsets of NVEs than in simulated exons, almost as high\nas non-NVE skipped exons, while splicing repressors like HNRNPK had lower splicing\nimpacts in NVEs, roughly comparable to non-NVE skipped exons (Reviewer Fig. 5).\nThese observations may simply reflect the necessity to possess or acquire positive\n11\nelements such as a strong polypyrimidine tract and to avoid negative elements such as\nC-rich HNRNPK binding sites in order to have sufficient levels of splicing to be classified\nas an NVE.\nReview Figure 5. Several RBPs have slightly different affinities for NVEs than non-NVEs,\nand are above the baseline of exons simulated from intronic sequences.\nA-B) Exons were simulated by randomly sampling intronic sequences above 500 bp in length,\nsize matched to the skipped exon NVE set (11k exons). The predicted impacts of RBP model\n(using KATMAP, add citation) were calculated for regions proximal surrounding splice sites (50\non each end). Plotted is, for each exon, the sum of all impacts for a given RBP (shown in x-\naxis), and binned by whether the exon is a non-NVE (100% top EF in the dataset, ~3k exons),\nNVE (colored by top EF), or simulated exon (grey).\nC-D) Exons were simulated by randomly sampling intronic sequences above 500 bp in length,\nsize matched to the skipped exon NVE set (11k exons). The predicted impacts of RBP model\n(using KATMAP, add citation) were calculated for regions proximal surrounding splice sites (50\non each end). Plotted is, for each exon, the sum of all impacts for a given RBP (shown in x-\naxis), and binned by whether the exon is a non-NVE (100% top EF in the dataset, ~3k exons),\nNVE (colored by top EF), or simulated exon (grey).\nMinor comments:\n14. Page 10: “The RNA-seq read depths available in GTEx are more than adequate to\ndistinguish these moderate ψ values from absence of splicing, and NVEs possess\nsplice site motifs very similar to those of other exons (Extended Data Fig. 2D), even for\nlow-EF NVEs (Extended Data Fig. 2E).” Did the authors rather mean to reference\nExtended Data Figures 2E and 2F respectively?\n15. Please add a legend for Supplementary Figure 5\nThis may have been cutoff in PDF conversion. We apologize. We have ensured that the\nlegend is visible in the revised manuscript.\nReviewer 2\nR2 Major comments\n1. The beta binomial model assumes a beta distribution for the underlying (unknown)\ndistribution of true PSIs across individuals. As they mention, this is somewhat flexible\nbut will not model all distributions. Even one of the examples shown in the manuscript\n(CAST, Main Text Fig. 1C) appears to be tri-modal with modes close to 0, 0.6, and 1, a\ndistribution that a single beta will not be able to capture. Mixtures of Gaussians are\nknown to be able to approximate any smooth distribution, so using a mixture of betas\nwould be a tractable and flexible alternative to the single beta. By summing over the\npossible mixture assignments it is still possible to use LBFGS to fit this model. Given\nthat the beta-binomial model is a key component of the analysis, I feel it should be\nascertained whether the results are consistent using such an approach. It would\nadditionally allow analysis of how often the distribution is unimodal/bi/trimodal etc.\nWe implemented a 3-component mixture of betas model and applied it across all exons\nin GTEx. Briefly, we use LBFGS and expectation maximization to approximate the most\nlikely parameters that can be described using a mixture of betas given the data. We\nchose a 3-component model because of its flexibility to model many different\ndistributions and because of our finding that many NVEs are associated with sQTLs,\nwhich have 3 possible genotypes that might each generate distinct splicing patterns.\nAfter fitting the mixture to all exons, we compared the EFs fit with this new 3-component\nmodel and the previous (1-component) model. In many tissues, the single component\nand new mixture models give similar EFs for the NVEs reported previously (Reviewer\nFigure 6). We have added a new Supplementary Note Section E: that provides more\ndetails, including simulations showing that a 3-component mixture model can handle\nmany different distributions of PSI values. The released code has also been updated so\nthat users can choose whether to use a 1-component, 2-component or 3-component\nmodel to fit observed RNAseq data. The major difference between our results with a\nsingle beta and those with the mixture model is that the resulting distribution of EFs is\nsomewhat shifted toward higher values (Reviewer Figure 7). This reflects a\ncombination of slight changes in EFs of previously defined NVEs and changes in the\nNVEs identified, especially an increase in the proportion of NVEs with higher PSI values\nthat were previously classified as non-NVE constitutive exons.\nestimated using the single beta model and the new 3-component mixture models are shown for\ndifferent subsets of NVEs (described above the plot), together with associated R2 values. B)\nsame as A, but for top EF NVEs. NVEs from single beta were merged with 3 betas and EFs\nwere compared. The color of the point represents the type of splicing event.\nBecause the PDFs are more accurate with a mixture of betas in cases where the PSI\ndistribution has multiple modes, including the CAST exon shown in Fig. 1C, and appear\ncomparably accurate in other cases, we chose to redo all of our analyses using the mixture of 3\nbetas model. In most cases, EF estimates don’t change substantially, even when the fit is\nbetter. For example, the CAST exon goes from EF 5% of 60% with the single component model\nto EF 5% of 50% with the mixture model.\n1 beta) versus new Figure 2B (right, using mixture of 3 betas). Star indicates significance\n(p=10^-255 by performing a simple normal test for proportions of proportions of low and high\nEFs in reference annotations and are connected by grey lines.\n2. The SMS part of the paper feels somewhat underdeveloped. There is no comparison\nto methods beyond Splice AI (e.g. MTSplice or Pangolin), no mention of the sample size\navailable for the logistic regression, and no assessment of which features are most\nimportant for the prediction. Unless I missed it also wasn't specified whether an SMS is\ntissue specific or not? I'd almost be tempted to suggest cutting it, but the bump in\nAUPRC adding the NVE SSs is remarkable and would be a shame to miss.\nWe agree that the SMS results are strong and add significantly to paper, so we have\nchosen to keep and expand them in the revised manuscript. SMS is not tissue-specific.\nThe sample size has been added to Methods (Methods: SMS algorithm). We describe\nSMS in more detail in Extended Data Fig 8. The coefficients for all features are shown\nin new Extended Data Fig. 9.\nA) Plotting the distance (in bp) between the variant in given QTL to the nearest splice site in the\nsplicing phenotype associated with the variant. Separated by which s-QTL PIP range the variant\nis in. Inset plot shows the range from 0-50bp and demonstrates that half of all high PIP (90% or\nabove) s-QTL variants tested are within 50 bp of the nearest splice site. K.S. test between\nhighest and lowest PIP range is shown.B) Model set-up and C) number of variants scored,\n16\nincluding training and test sizes. D) Performance of model in held-out test sets, using either all\nfeatures without distance to splice site (far left), or only distance to splice sites as features:\nusing GENCODE splice sites (middle) or using distance to GTEx splice site as a single feature\nmodel (right). Model was tested using a held-out test set of 20% of negative and positive set\nvariants AUC shown on top of all plots for respective models.\nExtended Data Fig. 9. Model performance with varying features\nA) Coefficients in SMS-full model (y-axis), split by gene region annotation (x-ais).\nB) Other features that are split by a different set of gene region annotations: “in any exon” and\n“not ss”\nWe compared SMS to Panoglin as well as Splice AI and both are included in the\nupdated model comparison shown in Figure 5B:\nhelps interpret GWAS.\nB) Comparison of model performance by precision recall curves (PRC) on a held out test set.\nPRC of full logistic regression model trained on GTEx distances, PRC of full logistic regression\ntrained on GENCODE distances only, and PRC of logistic regression trained only on SpliceAI\nscore is shown, well as Pangolin, and the maximum FM, and AM scores.\nR2 Minor comments\n3. It's not a major concern but the finding that NVEs are _more_ common in constrained\ngenes warrants more discussion. Are the constrained genes more highly expressed\nmaking NVEs easy to detect? It might be worthwhile to ask if the NVE exons with highly\nconstrained themselves are constrained as a class, compared to constitutive or\ncanonical alternative exons (since testing individual exons would be underpowered).\nGood question. Yes, highly constrained genes are more highly expressed (PMID:\n32461654), and this may contribute to NVE detection. However, the fraction of genes\ncontaining NVEs was still lower in a set of less constrained genes that was sampled to\nmatch the expression distribution of highly constrained genes, as shown below.\nFigure 3B and Extended Data Figure 6G. (Shown side-by-side to save space.)\n3B) B) Percentage of genes with NVEs in constrained and unconstrained genes, matched on\ngene expression, in Whole Blood. For each gene, the median TPM across all individuals was\n18\ntaken to assess average expression in blood. Genes were separated as either being\nconstrained (bottom quintile of LOEUF scores in blood) or unconstrained (top quintile of\nLOEUF scores in blood). Constrained genes were then filtered to have a median expression of over 1\nTPM, and then we found matching unconstrained genes by taking genes that had the same\nmedian expression to 3 decile points, n = 178. Significant by proportion normal z-test in scipy,\np=10^-11. Extended Data Fig. 6G) gene expression distributions of two sets of genes shown in\nFig. 3B in gene expression matched genes with NVEs.\nWe also observed that NVEs in highly constrained genes have lower EF and PSI\nvalues. One interpretation of these observations is that highly constrained genes may\nhave features – such as their exon-intron organization, sub-nuclear localization or other\nproperties – that favor the emergence of new exons, relative to the features of lowly\nconstrained genes. But that more stringent selection on the expression levels of\nconstrained genes are less tolerant of higher PSI values (which would impact the levels\nof canonical mRNA and proteins more), so the PSI and EF values of these NVEs are\nless free to drift to higher values. Another possibility is discussed in the text.\n4. \"This exon appears completely absent or very lowly included in most people, but\noccurs at a Ψ of 5-10% in about 10% of individuals, with credible intervals on Ψ well\nseparated from zero\". This phrasing gives the impression there is bimodality, whereas\nreally the data looks like a continuum to me.\nWe have changed the text to avoid implying bimodality. The new sentence reads:\n“This exon has variable inclusion across individuals, ranging from completely absent or\nvery lowly included up to a PSI of 5% or higher in about 10% of individuals.”\n5. \"This splicing variation does occur against\" cut \"does\"\nFixed\n6. Some amount of inter-individual variation in splicing for events differentially spliced\nbetween cell-types could be due to variable cell-type proportions. Cell type proportion\nestimates could be used to control for this possibility.\nYes, this is a possibility. We considered a few papers on this topic (citations listed\nbelow) but found a lack of consistency between different estimates for the same\ntissues/samples. Therefore, although this might be productively explored in the future,\nwe feel that the estimation methods are not sufficiently mature to perform a definitive\nanalysis at this stage. We have acknowledged the potential contribution of variability in\ncell-type proportions in the text. Citations for cell type proportion estimates:\nhttps://www.biorxiv.org/content/10.1101/2022.05.17.492324v1.full.pdf\nhttps://www.nature.com/articles/s41467-022-28655-4\nhttps://pmc.ncbi.nlm.nih.gov/articles/PMC10888823/\n7. the legend in 2b is confusing.\nWe rewrote the legend and added percentages to improve clarity.\n19\n8. Several of the comparisons in \"Most NVEs occur in coding regions and have low EF\"\nwould benefit from comparisons to a control set of exons.\nAgreed. We added a control set of annotated GENCODE exons for the analyses shown\nin Figure 2 (new figure below). The PSI values of NVEs tend to be lower than control\nexons (Fig. 2C). The proportions of NVEs in coding regions is above that of control\nexons at low EF values, and below that of control exons at high EF values (Fig. 2D).\nThese observations are discussed in the text.\nNew Figure 2. The distribution of NVEs across genes, gene regions, and their relative\nusage in individuals.\nA) Histogram of the number of NVEs in a given gene, including those protein coding genes with\nno NVEs (gray) at least one NVE (purple). B) The top EF spectrum in GTEx. For each NVE, the\ntissue with the most usage of the NVE was taken, indicating the highest possible EF in GTEx.\nWe use this color scheme throughout the paper with lower EFs shown in blue and higher EFs\nshown in red. Lighter boxes indicate that the NVEs are not on GENCODE reference\nannotations. Star indicates significance (p=10^-255 by performing a simple normal test for\nproportions of proportions of low and high EFs in reference annotations and are connected by\ngrey lines. C) Average ψ of NVE in individuals who splice in ≥ 5% ψ. 1000 randomly sampled\nNVEs were selected, repeated 10 bootstraps. We approximated the conditional expectation of\nΨ ≥ 5% using the empirical mean of observed Ψ values above that value, within each EF bin.\nUncertainty was represented by the 2.5th and 97.5th percentiles. D) Proportion of genes in\nMANE select/RefSeq transcripts in different gene regions (right) and Proportion of NVEs in\ngiven gene regions, split by exon frequency (left). Stars indicate to significance (p < 10^-24 and\np = 10^-24) by performing a simple normal test for proportions of CDS proportions between\ngroups connected by grey lines.\n20\n9. label 3a&b with which direction is more constrained\nDone\n10. what's going on with substantia nigra in 5b??\nWe looked into this tissue and observed that there are very few samples available in\nGTEx for younger ages. We also investigated some technical factors raised by another\nreviewer (see response to Review 3 comment 4 below) and found that controlling for\nthese factors generally reduced the magnitudes of associations between age and NVE\ncount. These observations reduced our confidence in the aging-related results generally\nand for substantia nigra in particular, so we removed these results from the revised\nmanuscript.\n11. psi_NVE should be labeled differently in eqns 1 and 2 since one is an estimate and\none is a parameter\nAgreed. Fixed.\n12. \"Where 𝛼 can be described as the prior\"... alpha is a prior pseudocount, not a prior.\nThe revised text reads: “Where 𝛼 can be described as a prior pseudocount.\n13. - \"...using the L-BFGS-B algorithm (see Supplementary Note), using a beta prior\nwith parameters α = β = 1\" not sure what this means.\nThis has been changed because we are using a mixture model now, as discussed\nabove.\n14. - \"keeping one of the two inclusion introns to simplify the analysis\" this is throwing\nout a lot of data. I think you could just take the average - the binomial likelihood in\npractice works fine for non integers.\nWe consider the counts of the upstream and downstream inclusion junctions as mostly\nredundant in that they both reflect the same underlying inclusion event. To analyze\nskipped exons and alternative splice sites on the same footing, we use the count of just\none of these junctions. See Supplementary Figure 1 for comparison of both intron\nreads.\n15. we computed EF as the percentage of individuals with Ψ ≥ c from the modeled CDF\nof the population\" assuming the CDF means the CDF of the fit beta(alpha,beta) (which\nshould be clarified) then this implies getting the distribution correct is v important for\ndefining NVEs, and further motivates the mixture of betas direction\ngenerated from the fitted model” (In Supplementary Note) and Supplementary Figure\n6. We agree that this motivates use of a mixture of betas, which we have implemented.\n16. \"90% of above\" -> or\nFixed.\n17. \"sQTL that we also eQTLs\"\nFixed.\nReviewer 3\nR3 Major comments\n1. In Figure 2B (page 10) the authors mention: “NVEs with low EF are predominantly\nunannotated, while around half of high-EF NVEs are annotated”. However, in the plot\nthe proportions do not seem different, all bars seem to be around 60%. The text could\nbe improved by adding the exact percentages in the text or in the plotting area instead\nand if the percentages are not different enough, tone down the message of finding\ndifferences between low EF and high EF NVEs.\nWe have modified the text to clarify that: “NVEs with low EF are predominantly (67%)\nunannotated, and around half of high-EF NVEs are annotated.”\n2. The authors report a correlation for Extended Data Fig. 6. However, a correlation\ncoefficient, and a p-value (independently of it being significant or not) should be added\nto the plot to help readers interpret the correlation.\nWe moved the correlation coefficient and p-values below the plot, as these values seem\nto have been cut off in the PDF conversion process. The updated figure is shown below:\nExtended Data Figure 6F.\nSpearman correlation between global NVE rate (number of NVEs per gene considering all\ntissues, divided by the number of exons in reference transcriptome, highly significant:\n22\n-log10pval > 90) and the LOEUF score of the gene. We compute this metric globally and not in\na tissue specific manner because many NVEs are tissue specific.\n3. In Figure 3c, authors should explain what they refer to as “Enrichment”, like the\nexplanation in the supplementary figure. Importantly, they report “The observed\nenrichments of causal variants in NVEs are above the level of synonymous variants,\nparticularly at high PIP values, but below those of missense annotations”. According to\nthe figure, the trend is not “particularly” but “only” observed at high PIP values, and the\nconfidence intervals completely overlap, so it is not a significant result. The authors\nshould tone down the results and state that “tend to be above” instead of “are above”.\nFurthermore, the sentence “Unannotated splicing events are enriched above\nsynonymous at high PIPs” should be removed from the figure captions.”\nAgreed. We have modified the text as requested.\n4. Figure 6A would benefit from having the (AUPRC) of the three curves, especially to\nsupport the claim: “When given access only to reference gene annotations, where the\ndistance feature uses GENCODE splice sites (SMS-GENCODE), the method’s\nperformance improves moderately over SpliceAI”. If not on the plot, the AUPRC should\nbe stated in the text.\nWe added the AUPRC to Fig. 5 in the text and to the legend of Fig. 5B.\n4. The Analysis shown in figure 5B should be redone, correcting for batch effects to\nconfirm for age being the main explanatory factor for more NVEs. Batch effects such as\nischemic time and hardy scale are a very important source of variation in the GTEx\ndataset.\nWe perform a generalized linear model (with standardized effect sizes for continuous\nvariables) associating the number of low-to-medium EF (≤ 50%) NVEs with age,\nischemic time, Hardy scale, sex, and BMI. We find that effect sizes of the age coefficient\nacross tissues are modest (ranging from about +0.035 to –0.035), with 12 tissues\nshowing a positive association and 8 tissues showing a negative association. We show\nall model coefficients in Reviewer Fig. 7. While some associations with age remain, the\nreduced magnitudes of these effects and the variability in direction – as well as the\npaucity of samples from younger individuals for some tissues – reduced our confidence\nin the aging results. We therefore decided to remove the aging results from the main\nfigure. This change also freed up space for us to expand discussion of the SMS results.\nWe show more relevant results in Reviewer Fig. 8-9.\nWe counted the number of mid to low EF NVEs in the dataset (50% EF or below) and\nperformed a generalized linear regression on age, sex, BMI, Hardy Scale, and ischemic time of\neach individual in GTEx. We plot here the tissues where after the GLM, age is still significantly\nassociated with age (p-value of coefficient is less than BF corrected 0.05). Plotted are the age\ncoefficients and error bars, sorted and colored by decreasing effect of age. Positive values\nmeans that age is positively associated with NVEs, and negative values mean age (more NVE\nover lifespan) is negatively associated with NVEs (less NVE over lifespan).\noccurrence of NVEs.\nA) Showing all features in GLM for Uterus, tissue with highest positive relationship between\nNVEs and age.\nB) Pearson correlation between age of individual and number of reads in the sample (-0.13).\nC) Tissues with selected marginally significant (BF corrected alpha 0.05) spearman correlations\nbetween age and low-to-medium EF (≤50%) NVEs. y-axis is the number of NVEs vs. age of an\nindividual on x. Each dot is an individual in the tissue sample. Note that the y-axis differs among\nplots.\nModel coefficients of age (left) and all other features (right) for each tissue in GTEx. Boxed in\nblack are marginally significant (bonferoni corrected p-value).\n5. In the first section of the results, the authors need to describe the methodology /\ncutoffs of actually calling NVEs in greater detail so that readers do not have to read the\nmethods text to get a general idea of what cutoffs were used.\nAgreed. We have added more details.\n6. In the introduction, the authors mention evolution of new exons and draw some\nconnections to why NVEs might represent evolutionarily new exons. This connection\ncould be explored by examining sequence conservation characteristics of NVEs, likely\nas part of the “Most NVEs occur in coding regions and have low EF” section.\nlineage, and should not be highly conserved across mammals. We assessed\nconservation by calculating PhastCons scores in 25 bp regions around A3SS NVEs\nlocated in annotated introns, and compared to the 25 bp around the canonical splice\njunction across the NVE intron (Reviewer Figure 10).\nReviewer Figure 10. Conservation upstream and downstream of alternative 3SS NVEs.\nConservation near a set of alternative 3’ss NVEs. PhastCons score was calculated within 25 bp\nof both NVE and constant 5'SS (which splices to cognate and non-cognate 3SS). NVEs were\nfiltered to those occurring in MANE select introns of at least 100 bp in length, flanked by coding\nexons. Mean fraction conserved at each position is shown, considering PhastCons score >0.95\nas conserved.\nNVE alternative 3p splice sites located in introns were lowly conserved, though had\nslightly higher average conservation than background intronic regions. These\nobservations suggest that most NVEs are newly evolved (and therefore not conserved),\nas we hypothesized, but that a small fraction may represent more ancient alternative\nsplice sites, which have conserved regulatory regions, and the presence of these\nancient sites raises the average conservation above background levels.\n7. In the sentence “NVEs rarely occur in 3' UTRs, likely because they contain very few\nintrons relative to other transcript regions”, the authors should elaborate to mention\nNMD, as this is the reason that 3’ UTRs do not contain introns.\nWe have added this to the text.\n8. In the sentence “Such a pattern could occur if NVEs arise more commonly in more\nconstrained genes but often experience selection favoring lower ψ values, driving EFs\nlower, although other scenarios are possible”, the authors should elaborate on other\nscenarios. To me, what they propose seems possible but as it is the only scenario\npresented it is hard to assess. In connection with this, I think the hypothesis that the\nauthors present is interesting, but it would be strengthened through complementary\nanalyses. For instance, the authors could examine the distribution of the metrics used to\ncall NVEs of exons stratified by high or low LOEUF genes. Maybe this isn’t the best way\nto get an answer, but it is what I can think of.\nWe agree that this is an interesting direction. We have added another possible scenario\nto the text and commented on an approach to explore these possibilities:\n“Such a pattern could occur if NVEs arise more commonly in more constrained genes\n(e.g., due to more favorable nuclear environment for splicing) but often experience\nselection favoring lower ψ values, driving EFs lower. NVEs in less-constrained genes\ntend to have higher EFs, some of which might become canonical alternative exons over\nevolutionary time. This suggests the alternative possibility that NVEs may arise at\nsimilar frequencies in all genes, but their EFs may evolve more rapidly toward 1\n(canonical alternative exon) or 0 (ψ below threshold in all people) in less constrained\ngenes due to relaxed selection, making them more evolutionarily transient in this subset\nof genes. Distinguishing among these possibilities would be difficult with current data\nbut might be enabled by comparing NVEs across human populations to help assess\ntheir evolutionary ages.”\nTo further explore this topic, we assessed the correlation between EF values of\nNVEs with the LOEUF scores of their host genes separately for NVEs in different gene\nregions. Overall, we observed small but significant positive Spearman correlations for\nNVEs in all regions studied, consistent with the trend observed above, but a somewhat\nhigher correlation for coding region NVEs that have potential to trigger NMD (Figure\n4D). The increased correlation for NMD-associated exons likely reflects a combination\nof factors. For example, UTR exons generally preserve the full reading frame and\nencoded protein, and frame-preserving exons may often retain protein function, while\nexons with NMD potential typically disrupt protein function. Thus, since constrained\ngenes are typically under stronger pressure to maintain protein levels in a narrow range,\npressure for NMD-inducing exons to have lower PSI and EF values may be stronger in\nthis subset of genes, contributing to the stronger correlation for this subset of exons.\nNew Figure 4D.\n28\n4D) NVEs were filtered into whether they were in coding (CDS) and have potential to\ncause NMD (w/ potential for NMD, orange) or no potential (CDS no NMD, grey) or non-\ntranslated regions of the gene (UTR, green). Spearman correlations were performed on\n1000 sampled NVE per gene region (without replacement), to control for the fact that\ntotal coding region (CDS) NVEs are larger in number than total untranslated region\n(UTR) NVEs. This process was repeated for 100 bootstraps, and the range of estimates\nare plotted as box plots. Shown is the average p-value for the correlations performed in\nbootstrapping.\nR3 Minor comments\n9. Supplementary Figure 5E figure caption is missing.\nThis somehow got cut off in PDF conversion but is now included.\n10. Authors should be consistent in the naming of the supplementary note figures.\nSometimes they are referred to as “Supplementary Fig. 3” and sometimes as\n“Supplementary Note Fig. 5”.\nFixed\n11. Sometimes the authors use “an NVE” and some other “a NVE”. For example: “that\ndescribes the frequency of splicing of an NVE in a population, somewhat analogous to\nallele frequency (AF). The EF of a NVE…” in page 5.\nFixed\n12. Mention what NMD stands for the first time you mention it\nFixed\n13. The title of method section “sQTL that we also eQTLs” is unclear\nApologies. This was a typo. It has been corrected to “sQTLs that were also eQTLs”\n14. References are missing in the paragraph: “The identification of variants that impact\nsplicing remains challenging…” in page 18\nWe added a reference and clarified that the statement is about difficulty in finding causal\nvariants that impact splicing.\n15. In page 21: “When training on the top splicing signal (top PIP), i.e. the best splicing\nassociation for each unique variant, we obtained qualitatively similar results and the\nsame ranking of algorithms (data not shown)”. Even if the data is not shown, the\n(AUPRC) could be reported in the text.\nDone\n29\n16. Regarding the statement “recent studies suggest that many exons are missing from\ntranscriptome reference annotations, occurring only in a subset of transcriptomic\ndatasets11,14,15 “. I don’t think GTEx is the paper for this statement.\nAgreed. This has been fixed.\nREVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author)\nIn this manuscript, Jacobs et al. discuss naturally variable exons (NVEs)\nobserved in the GTEx dataset. The authors introduce a statistical approach\nthat captures uncertainty in exon inclusion. They define exon frequency (EF)\nas a metric to describe the prevalence of NVEs in a population and\ndemonstrate that most NVEs generally have low EFs and occur in coding\nregions. The authors characterize the functional importance of NVEs by\nassessing gene tolerance to LoF. Moreover, the study tackles the\nchallenging problem of variant interpretation and yields valuable insights by\nmeans of their logistic regression model SMS. The manuscript presents\nnovel approaches to understanding splicing regulation at the individual\nlevel and I believe that this study will be a significant contribution to the\nfield.(Remarks on code availability)\nThank you for your comments\nReviewer #2 (Remarks to the Author)\nThe authors have done a thorough job addressing my concerns and those\nof the other reviewers. I recommend the manuscript be accepted.(Remarks\non code availability)\nThank you for your comments\nReviewer 3:\nWe want to thank the authors for taking the time to address our concerns and\ngenerally agree that most of ours have been resolved. In particular, we are satisfied\nwith the author’s decision to remove the section based on aging. Otherwise, we have\na few minor issues:\nIn Figs. 1b-c, the error bars are difficult to see / interpret. Could the authors use a\nsmoothed CI estimate instead?\nWe thought about using a smoothed CI, but there is quite a bit of variability in CI\nbetween individual samples that would be lost if a smoothed estimate were used, so\nwe prefer to show the individual CIs and let readers zoom in if needed.\nThe authors should mention the significant result from the test done for Fig. 3a in\nthe text, as it seems statistically unsupported just by reading the text at the moment.\nGood catch. We added a p value to the text that was previously shown in the legend.\nGiven that a large focus of the analyses of this work pertains to genetics, I am\ncurious if the splice motif and / or motif score analyses are done using the human\nreference genome or genomes that are somehow personalized using the genotypes\navailable in GTEx? I was unable to find details in the methods, which should be\nadded.\nFor the SMS analysis, we used variant information to calculate the splice site scores\nand the change in splice site score between the reference and variant. We did not use\ngenotype to assign NVEs, as EF is a population-level rather than an individual-specific\nmeasure.\nRelated to the splice motif analyses, for the NVEs such as in CAST where some\nindividuals constitutively use the NVE and in others it’s completely absent, we are\nwondering if there are genetic differences in the splice sites for the NVEs? This is\nnot a request but a curiosity.\nThere are variants that impact the usage of the CAST exon, and this is mentioned.\nIn response to our comment about the alternative hypotheses of Psi distributions\nand mutational constraints, I definitely think the supporting text that was added\nhelps, but I don’t understand what the authors mean by “favorable nuclear\nenvironment for splicing”; particularly so when it’s not accompanied by a citation.\nWe added a citation to this sentence and clarified our reasoning: we added that\nconstrained genes tend to be more highly expressed, and highly expressed gene loci\ntend to be located near nuclear speckles, where the splicing machinery is more\nconcentrated, so this may create a nuclear environment that is more conducive to the\nemergence of new splicing patterns than that of lower-expressed genes, including\nmany less-constrained genes.\nThe authors clearly did a lot of work to perform analyses for the reviewers and we are\na bit confused about their decision to exclude the “reviewer figures” from the work. If\nit were us, we would have found places for the most meaningful reviewer figures in\nmain or supplementary text to support the work as a whole. Perhaps the authors\nshould think about this as well.\nWe thank the reviewers for their insightful comments and suggested analyses. We did\nadd several new analyses based on reviewer comments to the final manuscript, but\nchose not to include a few others that we felt did not impact the main conclusions and\nmight be distracting.\nReviewer 4:\n* https://hub.docker.com/r/hnjacobs/python-11-image-hj--> The Docker code\nlinked here is provided with no documentation, and was last updated 11 months\nago. It's impossible to comment on the utility of the code when everything is\ninside the Docker container.\nThis is the model from the initial submission, and will not be used in the final release.\nWe apologize for the confusion.\nAs for the code linked under the \"Code availability\" section,\n* https://github.com/atgu/mixture_betas --> I get a 404 error\nWe have since made it public. Apologies for this error!\n* https://github.com/jacobs-hannah-mit/post_process_POVS --> Very little and\nconfusing documentation (ie references to other webpages without links to\nthem) and inclusion of a library dependency that says that you must contact the\nlab to get access to the tool\nWe do not currently have the bandwidth to fully optimize this code and documentation\nfor easy use. This is something we are actively working on, and hope to make\navailable in coming months. Interested parties are encouraged to contact the authors.\n* https://github.com/jacobs-hannah-mit/SMS --> Only has documentation for the\ntemplate used to generate the repo; no documentation related to the actual\nanalyses performed\nWe have updated this repo."}, {"round": 2, "reviewer_comments": {"reviewer_1": "In this manuscript, Jacobs et al. discuss naturally variable exons (NVEs) observed in the GTEx dataset. The authors introduce a statistical approach that captures uncertainty in exon inclusion. They define exon frequency (EF) as a metric to describe the prevalence of NVEs in a population and demonstrate that most NVEs generally have low EFs and occur in coding regions. The authors characterize the functional importance of NVEs by assessing gene tolerance to LoF. Moreover, the study tackles the challenging problem of variant interpretation and yields valuable insights by means of their logistic regression model SMS. The manuscript presents novel approaches to understanding splicing regulation at the individual level and I believe that this study will be a significant contribution to the field.", "reviewer_2": "The authors have done a thorough job addressing my concerns and those of the other reviewers. I recommend the manuscript be accepted.", "reviewer_3": "We want to thank the authors for taking the time to address our concerns and generally agree that most of ours have been resolved. In particular, we are satisfied with the author’s decision to remove the section based on aging. Otherwise, we have a few minor issues:\nIn Figs. 1b-c, the error bars are difficult to see / interpret. Could the authors use a smoothed CI estimate instead?\nThe authors should mention the significant result from the test done for Fig. 3a in the text, as it seems statistically unsupported just by reading the text at the moment.\nGiven that a large focus of the analyses of this work pertains to genetics, I am curious if the splice motif and / or motif score analyses are done using the human reference genome or genomes that are somehow personalized using the genotypes available in GTEx? I was unable to find details in the methods, which should be added.\nRelated to the splice motif analyses, for the NVEs such as in CAST where some individuals constitutively use the NVE and in others it’s completely absent, we are wondering if there are genetic differences in the splice sites for the NVEs? This is not a request but a curiosity.\nIn response to our comment about the alternative hypotheses of Psi distributions and mutational constraints, I definitely think the supporting text that was added helps, but I don’t understand what the authors mean by “favorable nuclear environment for splicing”; particularly so when it’s not accompanied by a citation.\nThe authors clearly did a lot of work to perform analyses for the reviewers and we are a bit confused about their decision to exclude the “reviewer figures” from the work. If it were us, we would have found places for the most meaningful reviewer figures in main or supplementary text to support the work as a whole. Perhaps the authors should think about this as well.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": ""}]

Decision letter