Identification of genes involved in the evolution of human intelligence through combination of inter-species and intra-species genetic variations

Preparation of bait genes

To exploit the genes related to human intelligence, six major works were collected from GWAS or GWAS meta-analyses on human intelligence with a large sample size from the last five years (various intelligence-related phenotypes including general cognitive, reaction time and verbal-numerical reasoning). A total of 271 loci associated with intelligence in the human genome (Table S1) were identified in these six studies.

The six GWAS studies included: (i) meta-analyses of general cognitive function (n = 53949, Davies et al., 2015); (ii) GWAS of cognitive function and educational attainment (n = 112151, Davies et al., 2016); (iii) meta-analyses of calculated Spearman’s g or a primary measure of fluid intelligence (n = 78308, Sniekers et al., 2017); (iv) meta-analysis and gene-based analysis of human cognition using 24 cohorts (n = 35298, Trampush et al., 2017); (v) GWAS using human populations with extremely high intelligence (n = 1238, Zabaneh et al., 2018); (vi) a recent meta-analysis of 14 independent epidemiological cohorts with intelligence assessed (n = 269867, Savage et al., 2018).

All the independent, significantly associated SNPs (IndSigSNPs) nearest genes (based on ANNOVAR annotations) were integrated with redundancies (the same gene identified in more than one study) removed. We took these 549 genes as “bait genes” (Table S2), which were candidates for human intelligence-related genes.

It should be noted that some GWAS works underlying human intelligence that have been published very recently may be not included in this study. This would not affect our analyses because this study aimed to provide some candidate genes for human intelligence evolution and could not identify all related genes at one time. In addition, because genotype imputation was not performed for the X chromosome in some cohorts (e.g., the UKB cohort including 195,653 samples with the assessed phenotype verbal and mathematical reasoning), the potential genes related to intelligence in the X chromosome were not included in our “bait genes”.

Preparation of variation pond

The new genome assemblies of chimpanzees, gorillas and orangutans have become fully available recently through single-molecule, real-time (SMRT) long-read sequencing, which improved the resolution of large and complex genomic regions. We incorporated the human-specific structural variations (from intermediate size to large size) into a “variation pond”, including exon gain/loss, short tandem repeats (STRs), insertion/deletion (indels) of more than 50 bp, and inversions. Moreover, considering the important role of the human-specific segmental duplications (HSDs, > 1 kb sequence with > 90% similarity, indicating the large recent duplication events, Bailey et al., 2001) in new gene function and human evolution (Dennis & Eichler, 2016), the HSDs identified recently from the comparative genomic information of both macaque and mouse (Dennis et al., 2017) were added into the “pond”. The “baits” and “ponds” were then integrated to detect whether there were human-specific variations hit by “bait genes”, thus generating the “prey genes” (Tables S3–S5), which refer to the intelligence-related genes with human-specific variations. The human-specific variation that was located within each of the 549 genes was left for further analysis using a window of 1 Kb (the locations of the variations and genes on the chromosomes were divided by 1000 for the computation efficiency).

The genome sequences of three great apes were downloaded from the NCBI database. The reference genome sequences from Pan troglodytes (chimpanzee), Pongo abelii (Orangutan) and Gorilla (Western Lowland Gorilla) were downloaded from: https://www.ncbi.nlm.nih.gov/assembly/GCF_002880755.1/, https://www.ncbi.nlm.nih.gov/assembly/GCF_002880775.1/ and https://www.ncbi.nlm.nih.gov/assembly/GCF_000151905.2/, respectively.

Localized alignments of the target gene sequences were performed to filter out the false positive human-specific variations from previous reports (McLean et al., 2011; Kronenberg et al., 2018; Dennis et al., 2017). The local sequences of the human genome were retrieved from chromosomes found in the GRCh38 version, and were then aligned with the great ape genomes using BLASTn (ncbi-blast-2.2.28+ version) with the parameters “-evalue 1E-50 -dust no”. The human-specific variation that was undetectable with a local BLAST was then removed for subsequent analyses.

Investigation of HIEGs

Most “prey” genes contained the human-specific variation in introns, far away from the exon-intron junction site. These variants are less likely to affect the gene functions, so only genes that contained variations in the coding regions were considered, named human intelligence evolution genes (HIEGs). These HIEGs included all the genes containing exon-gain/loss (2 genes), hCONDEL (28 genes), or HSD (1 gene), and genes with exon-located indels (8 genes) or STRs (4 genes). 40 non-redundant genes were finally identified, which were associated with human intelligence and which carried significant human-specific variations.

The gene transcript information was obtained from Ensembl Release 95 (http://asia.ensembl.org/index.html). For exon alignment, all transcript isoforms of one gene were compared both intra-species (within human) and inter-species (between human and the great apes) in order to confirm the specificity of the new transcript in humans. The human principal isoform, the human variant isoform, the chimpanzee isoform, the gorilla isoform, and the orangutan isoform of the PCCB gene are ENST00000469217 (NM_001178014), ENST00000466072, ENSPTRT00000028811, ENSGGOT00000005484 and ENSPPYT00000016431, respectively. The human principal isoform, the human variant isoform, the chimpanzee isoform, the gorilla isoform, and the orangutan isoform of the STAU1 gene are ENST00000371856, ENST00000340954, ENSPTRT00000050938, ENSGGOT00000048652 and ENSPPYT00000012903, respectively.

The analysis for protein domain was performed using the Simple Modular Architecture Research Tool (SMART) in the normal mode (http://smart.embl-heidelberg.de). The protein accession numbers in humans, chimpanzees, gorillas, and orangutans for KMT2D are NP_003473, XP_016778992, XP_018894141, XP_024112209, respectively. The protein accession numbers in human, chimpanzee, gorilla, and orangutan for TRIOBP are NP_001034230, XP_016794633, XP_004063488, and H2P4B5 (UniProt), respectively. We performed multiple alignments using the software Constraint-based Multiple Alignment Tool (COBALT, https://www.ncbi.nlm.nih.gov/tools/cobalt/re_cobalt.cgi).

Expression analysis for the candidates for intelligence evolution

We assessed the expression patterns of the HIEGs using transcriptome data for humans and their closest relatives. The HPA RNA-seq data were downloaded from the Human Protein Atlas (http://www.proteinatlas.org), including 102 samples of 37 tissues, in which TPM (transcripts per million) were used for the evaluation of expression level. In order to compare the expression levels between humans and great apes, the RNA-Seq data (NCBI ID: (100796) from 107 samples in 8 brain regions of humans, chimpanzees, gorillas, and gibbons were used (Xu et al., 2018). The brain tissues included five neocortical areas and three other brain tissues. The five neocortical areas are the anterior cingulate cortex (ACC), the dorsolateral prefrontal cortex (DPFC), the ventrolateral prefrontal cortex (VPFC), the premotor cortex (PMC), the primary visual cortex (V1C). The three other brain tissues are the hippocampus (HIP), the striatum, and the cerebellum (CB). Hierarchical cluster analysis was applied using the RPKM (reads per kilobase per million) of 39 genes in 8 brain tissues of humans and the other three primates, and was displayed using the software MeV4.2 (http://mev.tm4.org/).

Identification of candidate genes for human intelligence evolution

The strategy that was used for detecting the candidate genes in the evolution of human intelligence is briefly described (Fig. 1). Hundreds of genetic loci for human intelligence and its related traits have been identified by GWAS using large population data. We integrated six high-quality GWAS works and large-scale GWAS meta-analyses from the last five years. This enabled the identification of a total of 271 loci in the human genome underlying intelligence-related phenotypes (Table S1). A total of 549 human genes, distributed across the whole human genome, were found to be located within the 271 intelligence-associated loci (Fig. 2, Table S2). We take the 549 genes as “bait genes”, which constituted intelligence-related genes from GWAS, underlying the intra-species intelligence-related variations related to intelligence in human populations.

In the meantime, the recently released long-read sequence assembly of chimpanzee, orangutan, and gorilla provided a large number of high-quality sequence differences between the human genome and the great ape genomes. We incorporated the human-specific structural variations, including exon gain and loss, STR, indels, hCONDEL, HSDs, and large structural variations, into a “variation pond” which are related to human evolution.

The “bait genes” and “variation pond” were then integrated to detect whether there were “bait genes” hit by the human-specific variation. After putting the “bait genes” into the “variation pond”, we found 406 sequence variations physically located within 213 genes related to intelligence, which is considered to be “prey genes” (Fig. 2). The 213 “prey genes” linked the inter-species and intra-species variations and may be related to the evolution of human intelligence. Additionally, the potential effects of the human-specific variation through local sequence comparisons and gene structure analyses were carefully checked. Following in-depth analyses, there were 40 strong candidate genes that were identified as containing human-specific variations, probably changing either the coding or the expression of intelligence-related genes, which were named as human intelligence evolution genes (HIEGs, Table 1).

Exon gain and loss on the intelligence-related genes

There were only two genes classified as the “prey genes” with exon gain or loss events, PCCB and STAU1 (Table 1). The PCCB gene encoding the propionyl-CoA carboxylase subunit beta, was located in the locus on chromosome 3 identified by GWAS for intelligence (P = 1.956 ×10⁻⁹, Savage et al., 2018). Human PCCB contained 15 protein coding isoforms according to the Ensembl release 95. The principal isoforms generated proteins of 539 aa, 559 aa and 570 aa, respectively, all of which have identical isoforms in the chimpanzee. However, a new transcript variant was human-specific. Compared with the principal isoform, the variant lost a 60-bp exon (exon 4 of principal isoform) but gained another new 60-bp exon (exon 11 of the variational isoform). The gained 60-bp exon did not appear in any transcripts of the PCCB of the three great apes (Fig. 3A). This human-specific transcript variant can be detected in the human cerebral cortex at a lower level than that of the principal isoform (transcriptome data were obtained from the Human Protein Atlas, HPA). However, the relative cerebral cortex expression level of the PCCB variational transcript was higher than that of the principal transcript (relative to the average expression level of the transcript in 37 tissues) (Figs. 3B and 3C). These results suggested the generally low expression level and the relative enrichment of the PCCB variant in the cerebral cortex, implying a potential role of the human-specific PCCB transcript variant in the cerebral cortex. Although there were no reports about the mechanism of PCCB function in neuron system development, mutations in PCCB are one of the major causes of the genetic disease propionic academia (PA). Neurological complications, such as intellectual disability, brain structural abnormalities, optic neuropathy, and cranial nerve abnormalities are significant symptoms of PA (Schreiber et al., 2012). Moreover, there were reports that patients carrying PCCB mutations exhibited intellectual disabilities (Witters et al., 2016).

STAU1 was also located in a locus identified by GWAS for intelligence (Savage et al., 2018). STAU1 encoding the double-stranded RNA-binding protein which regulates RNA metabolism. There was a total of 10 protein-coding isoforms of human STAU1. Compared with the longest principal isoform, a human-specific insertion resulted in a gain of 123-bp exon (exon 2, located within the 5′UTR of the gene) in one transcript variant. The isoform, with the addition of a 123-bp exon, was not detected in any transcripts of chimpanzee, gorilla or orangutan (Fig. 3D). The human-specific isoform was expressed in many human tissues. In the human cerebral cortex, the expression of the new isoform was equivalent to ∼21% of that of the principal isoform of the gene (Figs. 3E and 3F). Previous functional studies found that STAU1 plays a role on mRNA transport in neuronal dendrites (Broadus, Fuerstenberg & Doe, 1998).

STR variations on the intelligence-related genes

In the variation pond, there were a total of 1,465 human-specific STR contractions and 4,921 human-specific STR expansions. These STRs hit 100 “bait genes” (26 genes containing the STR contractions and 74 genes containing the STR expansions, Table S3). Most of the human-specific STR variations were located within the intron regions or intergenic regions. Only four human-specific STRs were located within the exonic regions, highlighted as HIEGs (Table 1). Among them, three STR expansions were located within the exon of non-protein coding isoforms (processed transcripts) of three genes (ARIH2, STAB1 and TSNARE1), and one STR contraction was located within the 39th exon of the KMT2D gene.

KMT2D, also known as MLL2, encoded an H3K4 histone methyltransferase made up of 5537 amino acids. GWAS for intelligence detected one associated locus (P = 2.518 × 10⁻¹⁴, Savage et al., 2018) on chromosome 12 containing the candidate gene KMT2D. Compared with the KMT2D sequence in the great apes, the STR region contracted for 60-bp in the c.11838 of human KMT2D. This region was also polymorphic among the human population. The contraction led to a 16-aa discontinuous deletion for the human KMT2D protein from the p.3958 position (Fig. 3G). The peptide coding by this STR region was in the coiled coil region for KMT2D, which could affect the protein structure through wrapping the hydrophobic residues and forming an amphipathic surface (Mason & Arndt, 2004). Mutations in KMT2D are the main cause of the genetic disease Kabuki syndrome, and several mutations in the 39th exon have been found in Kabuki patients. Kabuki syndrome affected mental capabilities and most of these patients showed various levels of intellectual disability (Lehman et al., 2017).

Indels on the intelligence-related genes

Among 5,894 human-specific deletions and 11,899 human-specific insertions in the “variation pond”, we found 94 “prey genes” with 148 deletions and 144 genes with 298 insertions (Table S4). Furthermore, it was found that there were 8 insertions and 1 deletion for the exonic regions, as highlighted by the HIEGs (Table 1). The insertion in ARIH2 was also identified as STR expansion which has been described above. Insertions in PDE4D, NRXN1, EXOC4, FUT8, and ZNF584, and the deletion in SLC27A5 affect the lengths of the non-protein coding isoforms (processed transcripts) of the six genes. Moreover, GNB5 and TRIOBP carried an insertion in the 3′UTR region and an insertion in the exon region, respectively.

The insertion in TRIOBP resulted in a gain of 675-bp coding region (c.887-1560) in the 5th exon of the longest isoform when compared with that in the chimpanzee. This variation resulted in a 238-aa discontinuous insertion in the region p.296–811 of the human TRIOBP protein (Fig. 3H). The GWAS for intelligence (P = 3.582 × 10⁻⁸, Savage et al., 2018) and the GWAS for underlying brain ventricular volume also identified the gene TRIOBP as a strong candidate (Vojinovic et al., 2018). Biochemistry experiments have shown that TRIOBP could physically interact with TRIO, which is an important protein involved in neural tissue development (Seipel et al., 2001). Mutations in TRIOBP could cause autosomal recessive deafness-28 (DFNB28), and surprisingly, several causal mutations were located within the human-specific insertion regions (e.g., R347X and Q297X) (Shahin et al., 2006).

We also searched for human conserved deletions (hCONDELs) near the “bait genes” that had been identified (McLean et al., 2011; Kronenberg et al., 2018). These sequences are lost in the human genome but are highly conserved among other species (including the great apes, the macaque, and the mouse). In total, 28 genes were identified as containing the hCONDELs. Although none of the 28 hCONDELs were located in the coding regions, they were all taken into HIEGs (Table 1), as the high lineage specificity of hCONDELs serving as important signs of intelligence evolution. The list included NRXN1, GRID2 and GRIA4, which were all involved in neurotransmission and the formation of synaptic contacts.

The GNB5 gene in the intelligence-associated locus on chromosome 15 (P = 2.47 × 10⁻¹¹, (Savage et al., 2018)), which encoded the beta subunit of the heterotrimeric GTP-binding proteins (G proteins). Aligned with the GNB5 sequences of the humans and the great apes, we found that there was a 292-bp human-specific insertion in the 3′UTR of the gene and a 1,472-bp hCONDEL in the intron (2.7 kb distance to the third exon) of the gene (Table 1). GNB5 was expressed in the brain tissues and participated in neurotransmitter signaling (e.g., through the dopamine D2 receptor) (Xie et al., 2012). In human populations, mutations in GNB5 have been reported to cause several diseases affecting intelligence, including language delay, ADHD/cognitive impairment with or without cardiac arrhythmia, and intellectual developmental disorder with cardiac arrhythmia (Lodder et al., 2016). The human-specific variation in GNB5 might participate in the evolution of human intelligence as well.

Large structural variations on the intelligence-related genes

The inversion variation is the rearrangement in which the genomic segment is reversed. Based on the previous report, there were a total of 625 inversions in our “variation pond” ranging from 9 kb to 8.4 Mb in size. Among these, 31 of them hit the “bait genes” (Table S5). However, none of these genes were located in the breakpoint of the human-specific inversions.

There were 218 human-specific duplications (HSD) the length of which were more than 5-kb (Prado-Martinez et al., 2013). Among them, one 24.6-kb HSD was detected to be overlapped with the AFF3 gene. The AFF3 also contained an hCONDEL around the intron-exon junction regions (31-bp distance to the exon). The region around AFF3 has been identified to be an intelligence-associated locus in found in two distinct meta-analyses conducted on intelligence (P = 1.56 × 10⁻⁸ , Sniekers et al., 2017 or P = 3.41 ×10⁻¹⁰, Savage et al., 2018), but no functional reports on neuron development are available.

Expression profiling analysis of the HIEGs

The transcriptome data of 37 human tissues from the Human Protein Atlas database (Uhlen et al., 2015) was used to investigate the tissue expression patterns of the 40 highlighted HIEGs (Fig. 4A). Of these there were 23 genes with higher expression levels in the cerebral cortex than their average expression levels in all 37 tissues. Furthermore, there were 12 genes with more than twofold expression levels in the cerebral cortex versus their average expression levels (Fig. 4A). These were regarded as genes that were potentially involved in the development of the cerebral cortex. Except FUT8, the remaining 11 genes are the hCONDEL-containing genes. The expression patterns of the 11 genes exhibited three types. Four genes, including GRIA4, NRXN1, CADM2, and CALN, showed the most significant high expressions in the cerebral cortex and low expressions in the other tissues. SGCZ, DCC, and GRID2 showed only relatively high expressions in the cerebral cortex compared with the low expressions in all other tissues. NCAM1, FBXL17, FUT8 and GNB5 were generally expressed in all tissues (Fig. 4A).

A transcriptome dataset sampling eight brain regions (five neocortical areas, hippocampus, striatum, and cerebellum) of both humans and four primate species (chimpanzees, gorillas, gibbon, and macaques) was reported previously (Xu et al., 2018), which enabled a comprehensive interspecies comparison. This dataset was used to examine whether the genes containing human-specific variation showed expression changes in the cortex between humans and the great apes. The expression data for 39 of the 40 HIEGs, with the exception of STAU1, could be found in the transcriptome dataset. Among them, the expression levels of the genes AFF3, SKAP1, REEP3, DCC, and SGCZ in the human neocortical areas were much lower than those in the neocortical areas of the great apes (fold change < 0.5), while the expressions of STAB1 in the human neocortex were much higher than those in the great apes (fold change = 4.8). Hierarchical clustering was also performed for the 39 gene expressions in all samples. We found that 5 neocortex tissues in the same species could be always clustered (that is, one clade for one species), while CB, STR, and HIP were generally clustered based on their tissue types (Fig4 B). This result suggested that the expression profiles of the HIEGs in the neocortex tissues displayed a strong species specificity, which was in contrast to the profiles in these non-neocortex tissues (e.g., CB). Taken together, the intelligence-related genes, containing the human-specific variations, expressed differently in the neocortex tissues between the humans and the great apes.