The availability of high-coverage genomes from our closest extinct relatives, the Neanderthals and Denisovans, constitutes a significant advance in the range of questions one can ask about the deep history of our species (Meyer et al., 2012; Prüfer et al., 2014; Prüfer et al., 2017; Mafessoni et al., 2020). One of the main themes emerging from this progress is interbreeding. In recent years, a fairly large number of admixture events between Neanderthals, Denisovans and Sapiens populations have been postulated. A recent review (Bergström et al., 2021) considers that at least four such events are supported by strong evidence.

While it is important to ask whether our species benefited from these admixture events (so-called adaptive introgression, where alleles inherited from other hominins rose to high frequency as a result of positive selection after gene flow), it is also worth examining regions of the genomes that are depleted of alleles resulting from gene flow from other hominins (Sankararaman et al., 2016; Vernot et al., 2016; Chen et al., 2020; Skov et al., 2020; Rinker et al., 2020). Such regions are called introgression deserts (sometimes also “genomic islands of divergence/speciation” (Wang et al., 2020) and have now been identified in a range of species (Fontsere et al., 2019).

There are multiple reasons why genetic differences that arose after the divergence of populations may not be well tolerated (Wolf and Akey, 2018): there could be negative selection on “archaic” variants (deleterious changes on the “archaic” lineage), or positive selection on human-specific variants (adaptive changes on the human lineage), or it may be due to drift. It is reasonable to expect, and indeed has been shown, that the X chromosome constitutes such a desertic region [not only in our species (Kuhlwilm et al., 2019; Martin and Jiggins, 2017)]. This could be due to repeated selective sweeps on this chromosome: genes involved in reproduction on this chromosome might act as strong reproductive barriers between populations (Fontsere et al., 2019).

In the case of modern humans, other genomic regions are devoid of Neanderthal and Denisovan introgression, for reasons that are perhaps less obvious, and therefore worth investigating further. A recent study (Chen et al., 2020) identifies four large deserts depleted of Neanderthal introgression, partially overlapping with a previous independent study (Vernot et al., 2016). As pointed out in (Kuhlwilm, 2018; Wolf and Akey, 2018), since it is likely that there were several different pulses of gene flow between us and our closest relatives (Iasi et al., 2021), the depletion observed in these four regions must have been reinforced repeatedly, and given the size of the deserts, it is reasonable that the “archaic” haplotype was purged within a short time after the gene flow event, as predicted by mathematical modeling on whole-genome simulations (Veller et al., 2021), and as evidenced in the analysis of genome-wide data from the earliest Late Pleistocene modern humans known to have been recovered in Europe (Hajdinjak et al., 2021).

The presence of FOXP2, a gene known for its role in language (Lai et al., 2001; Fisher, 2019), in one of these large deserts has attracted attention (Kuhlwilm, 2018), as it raises the possibility that the incompatibility between Homo sapiens and other hominin in such persistent introgression deserts may point to (subtle, but real) cognitive/behavioral differences. Indeed, the presence in such deserts of not only FOXP2 but also other genes like ROBO1, ROBO2, and EPHA3, all independently associated with language traits (St Pourcain et al., 2014; Wang et al., 2015; Eising et al., 2021; Mekki et al., 2022), together with an earlier observation in Vernot et al. (2016) that genes within large deserts are significantly enriched in the developing cerebral cortex and in the adult striatum, suggest a possible point of entry into some of the most distinctive aspects of the human condition (Pääbo, 2014). Such considerations, combined with independent evidence that introgressed Neanderthal alleles show significant downregulation in brain regions (McCoy et al., 2017), motivated us to focus on the brain in this study.

Specifically, we focused on the four largest genomic regions that resisted “archaic” introgression reported in (Chen et al., 2020), jointly with the most comprehensive catalog to date of signals of positive selection in our lineage (Peyrégne et al., 2017) (see Table 1), a combination that, to our knowledge, has not been previously studied in detail. Here, we tested if the genes that fulfill these two conditions (falling within large deserts of introgression and being under positive selection) follow particular (brain-region) expression trajectories that significantly deviate from that of other subsets of genes with evolutionary relevance or from control genomic regions. We characterized the gene expression dynamics (including genes falling within either deserts of introgression or positively-selection regions alone) by analyzing transcriptomic data from several brain regions encompassing multiple developmental stages from prenatal to adulthood. This dataset allows for greater resolution than the Allen Brain Atlas data used in (Vernot et al., 2016), especially at early stages of development (see Figure 1). Three of the brain regions under study showed marked transcriptomic divergence (i.e., a statistically significant difference when compared to all other regions, based on the Principal Component Analysis-derived Euclidean distances): the cerebellum, the striatum and the thalamus. Among the genes at the intersection of regions under positive selection and large deserts of introgression, we found CADPS2, ROBO2, or SYT6, involved in neurotrophin release, axon guidance and neuronal proliferation, and known to be expressed in the brain regions our analysis highlights.

There are two main findings to take away from our study: the importance of structures beyond the cerebral neocortex in the attempt to characterize some of the most derived features of our species’ brain, and the fact that some of the strongest effects in these regions takes place at early stages of development. In this way our work provides complementary evidence for the perinatal globularization phase as a species-specific ontogenic innovation (Gunz et al., 2010), and also provides new evidence for the claim that brain regions outside the neocortex (cerebellum, thalamus, striatum) significantly contribute to this phenotype (Boeckx and Benítez-Burraco, 2014; McCoy et al., 2017; Neubauer et al., 2018; Gunz et al., 2019; Weiss et al., 2021).

To our knowledge this is the first study to reveal the effect of the cerebellum in the context of large introgression deserts. For the striatum, previous studies have already highlighted the relevance of this structure: genes carrying Neanderthal-derived changes and expressed in the striatum during adolescence exhibit a higher McDonald-Kreitman ratio (Mafessoni et al., 2020). In addition, using a different range of introgressed regions and gene expression data from the Allen Brain Atlas (with lower temporal resolution than the database used in this study), it had already been noted (Vernot et al., 2016) that genes within large deserts are significantly enriched in the striatum at adolescence and adult stages, which converges with the life stages highlighted from our analysis using the most recent report of genomic regions depleted of archaic variants (Chen et al., 2020).

Naturally, the functional effects of these divergent developmental profiles for the cerebellum, the prenatal thalamus or the striatum remain to be understood, particularly in the context of the possible differences among Homo-species concerning regulation of the genes highlighted in this study. This is especially relevant in light of emerging evidence that selection against DNA introgression is stronger in regulatory regions (Vilgalys et al., 2021), which in addition have been found to be over-represented in putative positively-selected regions in Homo sapiens (Peyrégne et al., 2017; Petr et al., 2019). The fact that early developmental stages are critical holds the promise of using brain organoid technology to probe the nature of these differences, since such in vitro techniques best track these earliest developmental windows (Muchnik et al., 2019; Mostajo-Radji et al., 2020; Kyrousi and Cappello, 2020). Our level of analysis (mRNA-seq data, informed by paleogenetic studies) can be complemented with other omics data to finely resolve cell-type specificities of the genes considered here across brain areas, as with the use of single-cell RNA-seq data, or to infer gene regulatory networks (from differentially accessible and methylated regions and chromatin immunoprecipitation data) that underlie the divergent gene expression trajectories observed.

The fact that FOXP2 expression is known to be particularly high in the brain regions highlighted here (Lai et al., 2003) may help shed light on why FOXP2 is found in one of the large introgression deserts in modern human genomes. As pointed out in (Kuhlwilm, 2018), this portion of chromosome 7 is not a desert for introgression in other great apes, nor did it act as a barrier for gene flow from Sapiens into Neanderthals. As such, it may indeed capture something genuinely specific about our species.

Analyses were performed using R (R Core Team, 2019). Putative positively-selected regions were retrieved from the extended set of sweep regions in Peyrégne et al. (2017), built from two independent recombination maps using a Hidden Markov-based model applied to African and Neaderthal/Denisovan genomes. Coordinates for (large) deserts of introgression were retrieved from Chen et al. (2020), and genes within these two sets of regions were obtained using the BioMart R package version 2.42.1 (Durinck et al., 2009), using the respective genomic region coordinates as input and filtering by protein-coding genes.

mRNA-seq analysis. Publicly available transcriptomic data of the human brain at different developmental stages was retrieved from (Li et al., 2018) and analyzed using R (full code can be found at https://github.com/jjaa-mp/desertsHomo). Reads per kilo base per million mapped reads (RPKM) normalized counts were log-transformed and then subsetted to select genes either in large deserts of introgression or in both deserts and putative positively-selected regions. The complete log-transformed, RPKM normalized count matrix was subsetted to select genes with median expression value > 2 , as in (Li et al., 2018), while no median filtering was employed for the subsets of genes within deserts and positively-selected regions, due to the potential relevance of the outliers in these specific regions for the purposes of our study. To assess transcriptomic variability between brain regions accounted for by genes either in large deserts or in deserts and positively-selected regions, we performed principal component analysis and calculated the pairwise Euclidean distances between brain regions for each dataset [following (Li et al., 2018)]. We then statistically evaluated such differences at each developmental stage using pairwise Wilcoxon tests with Bonferroni correction. Significant differences were considered if p < 0.01. Our analysis based on statistically significant principal components did not make it possible for us to use the Allen Brain Atlas data for comparisons with the psychENCODE project dataset used in this study, due to the more limited resolution, especially at prenatal stages, offered by the former.

To evaluate the expression profile of genes from our regions of interest in comparison to other regions of the human genome, we generated sets of random regions of the same length and gene density (that do not overlap with the genomic coordinates of deserts on introgression). These served as control regions for comparisons of mean expression values using two-way repeated measures ANOVA, implemented in R. ANOVA tests were performed taking mean expression values as dependent value, with structure names as subject identifiers and the different regions of interest (datasource) as between-subjects factor variable. Posthocs tests were performed similarly but with the mean expression data grouped by the datasource, obtaining an ANOVA table for each structure, with a Bonferroni correction to account for repeated measures. The stage-version of the ANOVA grouped subject identifiers by stage. Two Kruskal-Wallis tests were used, one designed to detect whether non positively-selected genes in deserts of introgression have different mean expression levels than genes that are both in deserts and in positively-selected windows; and the second to determine whether any particular brain structure has a particularly different expression mean than the rest, regardless of selection. We also used two two-level linear mixed-effects regression models, to compare non-positively selected genes and positively selected genes within introgression deserts. These models consist of repeated measures of expression on different brain structures in three different groups: control, deserts of introgression, and deserts with selection signals. The same model applies when stages are taken into account, replacing structure identifiers. Tukey’s test was then used to fit the model.

Gene-specific expression trajectories. The R package Trendy version 1.8.2 (Bacher et al., 2018) was used to perform segmented regression analysis and characterize the expression trajectories of genes within both deserts of introgression and putative positively-selected regions (12 genes). The normalized RPKM values [from (Li et al., 2018)] in the form of a gene-by-time samples matrix was used to fit each gene expression trajectory to an optimal segmented regression model. Genes were considered if their adjusted R ² was > 0.5 . In addition, a maximum number of breakpoints (significant changes in gene expression trajectory) was set at 3, minimum number of samples in each segment at 2, and minimum mean expression, 2.

The permutation tests using gene expression data from (Li et al., 2018) were done using the regioneR package version 1.26.1 (Gel et al., 2016) at n = 1,000.