Accounting for cis-regulatory constraint prioritizes genes likely to affect species-specific traits

Measuring allele-specific expression in interspecies hybrids is a powerful way to detect cis-regulatory changes underlying adaptation. However, it remains difficult to identify genes most likely to explain species-specific traits. Here, we outline a simple strategy that leverages population-scale allele-specific RNA-seq data to identify genes that have constrained cis- regulation within species yet show divergence between species. Applying this strategy to data from human-chimpanzee hybrid cortical spheroids, we identify signatures of lineage-specific selection on genes related to cellular proliferation, speech, and glucose metabolism. We also highlight cis-regulatory divergence in CUX1 and EDNRB that may shape the unique trajectory of human brain development.


2
Background 27 28 Changes in gene expression are thought to play a major role in the evolution of complex traits 29 [1][2][3][4]. As a result, comparing gene expression between species can enable the identification of 30 molecular changes underlying phenotypic divergence. However, obtaining accurate 31 comparisons of gene expression between species is challenging due to confounding factors like 32 age, environmental effects, differential cell type abundances, differences in developmental 33 timing, and batch effects [5][6][7]. The use of interspecies hybrids overcomes these issues through 34 the measurement of allele-specific expression (ASE) [8,9]. In hybrids, the genomes of both 35 species share the same nucleus and are exposed to identical environments, so there are no 36 confounding environmental, batch, compositional, or developmental effects. This approach has 37 been successfully applied in many species and advanced our understanding of the evolution of 38 gene regulation and its role in establishing phenotypic variation [10][11][12][13]. Furthermore, the recent 39 development of human-chimpanzee allotetraploid "hybrid" cells and organoids enables detailed, 40 accurate quantification of differences in gene expression between humans and our closest living 41 relatives [8,9,14]. 42 43 Hybrids also enable the separation of cis and trans components of interspecies differences in 44 gene expression [8,9]. The cis-component is caused by differences in regulatory elements such 45 as promoters or enhancers that only affect the expression of a nearby gene or genes on the 46 same DNA molecule. The trans-component stems from changes in diffusible molecules such as 47 transcription factors that can regulate gene expression throughout the genome. In hybrids the 48 genomes of the two species are exposed to the same trans factors. As a result, allele-specific 49 differences in gene expression can only be explained by cis-regulatory differences. In addition, 50 using ASE to identify differentially expressed genes (referred to as AS-DEGs) enables 51 8 143 First, we tested whether the variance in ASE in the human population is a reasonable proxy for 144 constraint on gene expression. We compared the variance of the GTEx ASE distribution for 145 each gene to its probability of haploinsufficiency score (pHI), a measure of sensitivity to a 50% 146 reduction in gene dosage [24]. We observe a significant negative correlation (Spearman's rho = 147 -0.28, p < 10 -170 ) between the two measures (Supp. Fig. 1) indicating that the variance of the 148 population-scale ASE distribution provides a reasonable proxy for constraint on expression 149 levels. This correlation remains significant when only five (as opposed to ten) reads from each 150 allele for a gene in a sample are required for inclusion (Spearman's rho = -0.27, p < 10 -151 163 ). Furthermore, as pHI is computed from copy number variation, there is no circularity in 152 reasoning when comparing pHI to the variance of the GTEx ASE distribution [24]. While this 153 indicates that variation in ASE is sufficient to approximate evolutionary constraint on gene 154 expression, the strength of this relationship should not be interpreted as a quantitative estimate 155 of how well ASE represents constraint. Overall, the correlation with pHI indicates that within-156 species ASE variance contains useful information about evolutionary constraint on gene 157 expression. 158

159
In addition, it is important to determine the extent to which sample heterogeneity may affect 160 ASE variance and thereby impact our results. For example, GTEx contains data from a wide 161 range of adult tissues and donors have variable sex and ancestry. To determine to what extent 162 these factors affect our results, we divided the GTEx dataset along lines of sex (male/female), 163 ancestry (African descent/European descent), and brain/non-brain tissues. In all cases, the 164 results obtained were very highly correlated (spearman's rho > 0.97, p < 10 -170 , Figure 2 A-C). 165 This indicates that sample heterogeneity is unlikely to have major effects on our results and that 166 aggregating across all the GTEx samples provides a reasonable measure of constraint on gene 167 expression within the human population.

169
In addition, the wide range of adult tissues in GTEx may not accurately reflect gene expression 170 constraint in cortical spheroids, which mimic fetal development [25,26]. To test this, we used an 171 ASE dataset generated from neural progenitors, neurons, and fetal cortical wall which 172 resembles the developmental stage of cortical spheroids more closely than GTEx samples [27]. 173 There is a strong correlation between the rankings generated by comparing to the fetal cortex 174 dataset and GTEx (Spearman's rho = 0.92, p < 10 -170 for each of day 50, day 100, and day 150 175 after initial differentiation of cortical spheroids; Fig 2D). While it is likely that the ASE variance 176 for some genes changes over development, these genes appear to be somewhat rare. Overall, 177 these results (Fig 2) indicate that even though gene expression levels vary considerably across 178 samples, ASE variance is robust to sample heterogeneity. This suggests that within-and 179 between-species ASE values can be meaningfully compared even when they are not based on 180 the same tissues or developmental timepoints. We therefore focused our analysis on the full 181 GTEx data in order to maximize statistical power but verified results for specific genes in the 182 fetal cortex data when appropriate. To highlight differences between our method and the traditional method for ranking genes (Fig  194   1A vs. Fig 1C), we computed the difference in ranks between the two methods and used 195 GSEAPY to identify gene sets enriched near the top or bottom of the resulting sorted list 196 (referred to as the difference in ranks list, Fig 3A) [28]. A top-ranking gene would be one that is 197 lowly ranked by the traditional method but highly ranked by the population comparison method.

Applying the Sign Test to Constrained Differentially Expressed Genes 255 256
While our method is designed to prioritize genes that are relevant to phenotypic differences 257 between species, it does not, on its own, imply selection on changes in gene expression. On the 258 other hand, the sign test is designed to detect lineage-specific selection based on systematic 259 up-or down-regulation of genes that deviates from neutral expectation [39,40]. One important 260 requirement of the test is that the expression divergence of every gene being tested should be 261 driven by independent genetic differences. For example, a single mutation in a trans-acting 262 factor could cause a whole pathway to be down-regulated, but this would only count as one 263 genetic difference. For this reason, hybrids are ideal for applying the sign test to gene 264 expression, since cis-regulation of genes is typically independent (except in the case of some 265 neighboring genes that share cis-elements, which can be accounted for by using a minimum 266 distance threshold between genes). For example, in the human-chimpanzee hybrid, if a 267 particular pathway contains significantly more genes with higher expression from the human 268 allele than the chimpanzee allele (significance determined by a binomial test), that would 269 provide evidence of selection on that pathway in the human and/or chimpanzee lineage. Similar 270 to all tests of selection, the sign test cannot discern what the cause of the selective pressure is. 271 For example, many cis-regulatory changes could be compensating for a change in a single 272 trans-acting factor. Alternatively, changes in gene expression might cause a phenotypic change 273 that is being selected for. In either case, the sign test provides evidence of lineage-specific 274 selection on cis-regulation. 275 276 While our method ( Fig 1C) does not require comparison to the traditional method (Fig 1A), it is 277 still of interest to determine if any of the categories that contain genes with strong differences in 278 ranking between our method and the traditional method show signatures of lineage-specific 279 selection. We therefore applied the sign test to identify gene sets in which highly ranked genes 280 (i.e. those with a difference in ranks between the population comparison method and the 281 traditional method above the cutoff found by GSEAPY) show a systematic bias in directionality. to within-human GTEx ASE. Fetal cortical ASE is not shown due to insufficient data. P-value is 304 the same as in C. E) ASE of EDNRB across timepoints in cortical spheroids. Expression from 305 the human allele is consistently higher than expression from the chimpanzee allele. F) ASE of 306 EDN1 across timepoints in cortical spheroids. Expression from the chimpanzee allele is 307 consistently higher than expression from the human allele except at D200. 308 309 Box 2: The Cytokine Activity category is enriched near the top of the difference in ranks list with 310 15 out of 16 genes driving the enrichment displaying chimpanzee-biased expression (Fig. 4A,  311 B). It is dominated by genes that generally promote neural stem cell proliferation, which is 312 unexpected considering the directionality bias and the higher proliferative capacity of human 313 neural stem cells. EDN1 had the strongest chimpanzee bias (mean log 2 fold-change = -0.92, 314 FDR < 0.025 at Day 150, Fig 4B). As EDN1 primarily signals through EDNRB in the brain [41], 315 we also investigated the expression of EDNRB. Surprisingly, EDNRB is one of the most strongly 316 human-biased genes across all timepoints in both hybrid and parental cortical spheroids (mean 317 log2 fold-change = 4.46, FDR < 0.005 at Day 150 in hybrids, log2 fold-change = 2.85, FDR < 318 0.0005 in parental samples). Human-chimpanzee ASE generally exceeds ASE found in human 319 populations for EDNRB (Fig 4C; Mann Whitney U Test comparing EDNRB distribution to human 320 population distribution p = 0.00027 for EDNRB), although this is not the case for EDN1 (Fig 4D). 321 The human-biased EDNRB expression and chimpanzee-biased EDN1 expression is generally 322 consistent across timepoints in hybrid cortical spheroids (Fig 4E, F). 323

324
Changes in EDNRB expression appear to be human-derived with respect to gorillas and 325 macaques (although orangutans may have independently acquired similar expression to 326 humans in early-stage brain organoids) (Supp. Fig. 3 A, B). We were unable to confidently 327 determine whether changes in EDN1 expression were human-or chimpanzee-derived using 328 currently available data (data not shown). Next, we analyzed single-cell RNA-seq data 329 generated from 1 month, 2 month, and 4 month-old human and chimpanzee brain organoids to 330 identify the cell types driving increased EDNRB expression. We found that a previously 331 identified radial glial cell (RGC) cluster was characterized by high EDNRB expression with non-332 zero expression in over 50% of cells [42]. Furthermore, this cluster had higher expression than 333 any chimpanzee cluster (Mann-Whitney U test, p < 10 -16 , Supp. Fig. 4-5). Overall, our results 334 suggest that a subpopulation of human radial glia have much higher EDNRB expression than 335 chimpanzee radial glia. Thus far, we have primarily focused on the differences between our new method and the 348 traditional method for ranking genes. Having established their differences, we now turn to 349 analysis of results from our new approach (Fig. 1 and Methods). To examine human-and 350 chimpanzee-biased genes separately, we sorted the list so that highly ranked genes with Spastic Dysarthria by binomial test), suggesting lineage-specific selection on genes with 361 constrained expression in these gene sets [40]. Notably, the bias in gene expression found in 362 the hybrids for these genes generally matched the bias in expression found in the parental 363 cortical spheroids (Supp. Fig. 6A-C). with spastic dysarthria and show systematic human bias, which may be connected to the human 374 capacity for speech ( Fig 5D). Finally, six genes that were ranked very highly by our method and 375 all have DESeq2 FDR < 0.1 drive the "Other glycan degradation" enrichment. Interestingly, loss Here we presented a method incorporating population-scale ASE data as a proxy for constraint 392 on expression. This ranking method helps reveal candidate genes and signatures of selection 393 that may explain phenotypic differences between humans and chimpanzees. The test is based 394 on the logic of comparing within-to between-species variation, similar in spirit to the  Kreitman-Aguade test although it is based on variation in ASE rather than protein-coding 396 sequences [51]. Although the phenotypic consequences of these differences remain to be 397 determined, our finding of polygenic lineage-specific selection on several gene sets suggests 398 that these changes must have some phenotypic effects in order to be under natural selection. 399 Collectively, these findings more than double the number of known cases of lineage-specific 400 polygenic selection on gene expression between humans and chimpanzees (the two previous 401 examples being Hedgehog signaling and astrocyte-related genes) [8,9]. 402 403 Importantly, the strong correlation between the GTEx brain vs. non-brain rankings ( Fig. 2A) and 404 GTEx vs. fetal cortex rankings (Fig. 2D) suggests that comparison to the GTEx ASE distribution 405 can be meaningfully compared with between-species ASE measured in diverse cell types and 406 organoids and that comparison to the GTEx population distribution will be useful for other cell 407 types and organoids. The method can also be applied to any species with sufficient gene 408 expression data, e.g. comparing ASE in Arabidopsis interspecies hybrids to ASE within A. 409 thaliana [52]. 410

Conclusion: 411
We outlined a strategy that uses allele-specific expression data from interspecies hybrids and 412 population-scale studies to prioritize genes that are more likely to impact species-specific traits 413 and applied this method to data from human-chimpanzee cortical spheroids. Our findings 414 provide opportunities for targeted follow-up experiments and increase our understanding of how 415 polygenic selection has shaped human and chimpanzee evolution. Overall, we anticipate that 416 our method will become a useful tool for identifying functionally significant gene expression 417 changes between species, and will contribute to our understanding of how gene expression 418 drives phenotypic diversification. SNPs and indels to correct for mapping bias. Reads for every sample were aligned to both the 426 human and chimpanzee genomes and the log2 fold-change from both alignments was 427 compared. Any genes with log2 fold-change that differed by greater than 1 were removed. We 428 used the ASE log2 fold-change (log2 FC) values available in the supplemental tables of Agoglia 429 et al. Although this dataset is restricted to hybrids from two humans and two chimpanzees, 430 previous work has shown that interspecies differences dominate over differences between 431 populations within a species so we expect that our results generalize well. It also contains 432 multiple independently derived hybrid lines and independent differentiations, reducing 433 confounding by technical differences. Throughout, for hybrids the mean log2 FC between human 434 genome mapped and chimpanzee genome mapped reads is stated as well as the highest p- of the chimpanzee allele). However, we did not require 10 counts from each allele and did not 482 normalize the medians. We did not require 10 counts from each allele because we expect 483 extreme differences in expression to be relatively common in between species comparisons. 484 We compared the log2(ASE Ratio) interspecies distribution to the population distribution using 485 the Mann-Whitney U Test (a nonparametric test robust to the distribution of data) and used the 486 resulting p-values to rank genes as described below. with DESeq2 FDR less than or equal to 0.1 and greater than 0.1. We then concatenated the list 493 so that genes with FDR less than or equal to 0.1 were always ranked higher than genes with 494 FDR greater than 0.1 (referred to as the MWU ranking). This consensus gene ranking was then 495 used in GSEAPY preranked with the rankings used as the score that GSEAPY uses to sort the 496 list. All results highlighted in the text replicated when using 0.05 as a cutoff instead of 0.1. We 497 used REVIGO in conjunction with a custom python script to generate the plots shown in figure 2  498 [62]. 499 500 To compare to the traditional method, we also ranked genes by the DESeq2 derived FDR and 501 used that in GSEAPY preranked (referred to as the DESeq2 ranking). To highlight differences 502 between the two methods, we computed the difference in ranks between the two methods by 503 subtracting the DESeq2 ranking from the MWU ranking and sorting the list on those rankings for 504 use in GSEAPY. In this context, highly ranked genes are likely those that show relatively mild 505 gene expression changes but have more constrained expression. 506 507 We next performed the expression sign test. First, we generated a list of all gene sets across all 508 tested ontologies that were nominally enriched at an FDR of 0.25 (using the FDR from GSEAPY 509 preranked) in at least one time point and that had greater than 10 genes driving the enrichment. 510 To avoid testing the same gene set multiple times, we only tested each gene set at the 511 timepoint that had the lowest GSEAPY FDR. We used the average of the DESeq2 log2 FC from 512 mapping to the human allele and from mapping to the chimpanzee allele as input for the 513 binomial test to identify gene sets with significantly more human-biased or chimpanzee-biased 514 changes than expected by chance. This log2 FC was generated by comparing the reads from 515 each species' allele in the cortical spheroid data. For example, if a gene set had 8 human-516 biased and 3 chimpanzee-biased genes, then the binomial test was used with k = 8, n = 11, and 517 p = 0.5. We considered any gene set with Benjamini-Hochberg corrected FDR < 0.1 to be 518

significant. 519
Finally, we ranked the genes using a signed version of the MWU ranking. More specifically, 520 genes were effectively ranked by the log10(MWU p-value) multiplied by the sign of mean 521 DESeq2 log2 fold-change so that top ranked genes with negative L2FC are at the bottom of the 522 list and top ranked genes with positive L2FC are at the top of the list. This ranking was then 523 used in GSEAPY preranked. All statistical tests (Mann Whitney U Test, Binomial Test, 524 correlations) were performed in python using the implementation in scipy. 525 For enrichment testing, we tested gene sets from the Gene Ontology Cellular Component and 526 Molecular Function categories, the Human Phenotype Ontology, KEGG, and REACTOME using 527 the same version as in Gokhman et. al [8,[33][34][35][36][37]. Regardless of which ranking was used, we 528 used GSEAPY preranked with the following arguments: processes=4, permutation_num=1000, 529 seed=6, min_size = 10, max_size = 300 to test for enrichment [28]. Following the authors 530 suggestion, we considered any category with an FDR below 0.25 to be nominally enriched [28]. 531 We required that ASE data be available from at least 50 individuals in GTEx for a gene to be 532 included in the ranking. We conducted gene set enrichment analysis with the four different 533 rankings described above. Single cell data from human and chimpanzee organoids and associated metadata were 538 downloaded from E-MTAB-7552 [42]. We used SCANPY to read in the counts matrix and filter 539 the data so that only data from 1 month, 2 month, and 4 month old organoids remained [63]. We 540 used a two-sided Mann Whitney U Test to compare EDNRB log2(counts per million) between 541 the "RGC early 2" cluster and all chimpanzee clusters with and without cells with 0 EDNRB 542 counts included. Mean counts by cell type and tissue for fetal human expression were 543 downloaded from GSE156793 [64]. We compared cerebrum counts to all other organs except 544 eye and cerebellum due to their similarity to the cerebrum and used a binomial test to determine 545 if more tissues exhibited higher expression than brain than expected by chance. (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE127898), GSE106245 558 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE153076), and phs000755.v2.p1 560 (https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000755.v2.p1), 561 phs002493.v1.p1 (link unavailable, data was received directly from the Stein lab). The GTEx 562 data was downloaded from https://www.gtexportal.org/home/datasets (although the specific file 563 that was downloaded is no longer available). The software used to perform the analysis in this 564 paper is available at: https://github.com/astarr97/Z_Score. There is no restriction on the use of 565 this code. Acknowledgements: We thank Jason Stein and Nil Aygun for providing processed fetal cortex, 574 neuron, and neural progenitor ASE calls. We thank members of the Fraser lab for helpful 575 discussion. Illustrations were created with BioRender. 576 Description of Additional file 1: File name is Additional file 1.csv. Enrichment testing for 577 difference between traditional and population-based method. This file contains the output from 578 GSEAPY for enrichment testing as well as the dataset that was used for the enrichment 579

analysis. 580
Description of Additional file 2: File name is Additional file 2.csv. Sign test on enriched terms 581 for difference between traditional and population-based method. This file contains the results for 582 the sign test applied to a subset (see Methods) of the significant terms from Additional file 1. 583