{'fraction_non_alphanumeric': 0.036173097200404555, 'fraction_numerical': 0.0234554883690904, 'mean_word_length': 4.961005643919959, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Species-Specific Changes in a Primate Transcription Factor Network Provide Insights into the Molecular Evolution of the Primate Prefrontal Cortex

Abstract The human prefrontal cortex (PFC) differs from that of other primates with respect to size, histology, and functional abilities. Here, we analyzed genome-wide expression data of humans, chimpanzees, and rhesus macaques to discover evolutionary changes in transcription factor (TF) networks that may underlie these phenotypic differences. We determined the co-expression networks of all TFs with species-specific expression including their potential target genes and interaction partners in the PFC of all three species. Integrating these networks allowed us inferring an ancestral network for all three species. This ancestral network as well as the networks for each species is enriched for genes involved in forebrain development, axonogenesis, and synaptic transmission. Our analysis allows us to directly compare the networks of each species to determine which links have been gained or lost during evolution. Interestingly, we detected that most links were gained on the human lineage, indicating increase TF cooperativity in humans. By comparing network changes between different tissues, we discovered that in brain tissues, but not in the other tissues, the human networks always had the highest connectivity. To pinpoint molecular changes underlying species-specific phenotypes, we analyzed the sub-networks of TFs derived only from genes with species-specific expression changes in the PFC. These sub-networks differed significantly in structure and function between the human and chimpanzee. For example, the human-specific sub-network is enriched for TFs implicated in cognitive disorders and for genes involved in synaptic plasticity and cognitive functions. Our results suggest evolutionary changes in TF networks that might have shaped morphological and functional differences between primate brains, in particular in the human PFC.


Introduction
Understanding why humans have unique cognitive abilities requires the identification of morphological and molecular aspects that are unique to the human brain. Unique morphological features of the human brain include its larger size (Povinelli and Preuss 1995;Koechlin et al. 2003;Schoenemann et al. 2005;Enard 2015;Smaers et al. 2017;Donahue et al. 2018) its cell-type compositions (Sherwood et al. 2006;Oberheim et al. 2009;Spocter et al. 2012), and specific cortical architectural structures (Buxhoeveden et al. 2006;Smaers et al. 2011). At the molecular level, there are several genes with brain functions that have been shown to evolve under positive selection on the human lineage, making them prime candidates for having contributed to the evolution of human-specific features. For example, ASPM (Zhang 2003;Mekel-Bobrov et al. 2005;Montgomery and Mundy 2012) and MCPH1 (Ponting and Jackson 2005;Voight et al. 2006;Pulvers et al. 2015), which determine brain size, and FOXP2, which when mutated causes severe cognitive and speech deficits (Enard et al. 2002b;Fisher and Scharff 2009;Konopka et al. 2009). Moreover, evolutionary young KRAB zinc-fingers (ZNFs) genes have been shown to be preferentially expressed in the human developing prefrontal cortex (PFC) (Nowick et al. 2011) and to evolve rapidly in sequence and expression in primates (Nowick et al. 2009;Nowick et al. 2011), suggesting that this gene family has played an important role during the evolution of the human brain. In line with these findings, several studies identified expression differences in the human compared with the chimpanzee brain that might be linked to human-specific traits (Enard et al. 2002a;C aceres et al. 2003;Somel et al. 2009;Babbitt et al. 2010;Liu et al. 2012). Since expression changes are often controlled by transcription factors (TFs) that are activating or repressing the expression of target genes it seems likely that TFs are responsible for driving some of the expression pattern differences and hence morphological differences between humans and other primates.
Despite the importance of TFs, only a limited number of studies so far have focused on evolutionary changes in TFs or TF networks in primates (Nowick et al. 2009;Schmidt et al. 2010;Schwalie et al. 2013;Ballester et al. 2014). These studies were limited in that the network analysis was based on ubiquitously expressed genes (i.e. not being able to reveal brain-specific differences) and only included human and chimpanzee samples (i.e. not being able to distinguish between changes on the human or chimpanzee lineage). Recent studies have investigated co-expression differences between humans and chimpanzees using an outgroup (Konopka et al. 2012;Sousa et al. 2017;Xu et al. 2018). This work has revealed human-specific co-expression modules especially in the frontal lobe, which are enriched for genes involved in neuronal processes and psychiatric diseases. However, due to the methods used in these studies (WGCNA; Langfelder and Horvath 2008), it is restricted to an analysis at the level of modules, thus not providing the possibility to pinpoint particular network links that have changed during evolution. Moreover, while progress in uncovering the biological cascades that take place during mammalian brain development has been made, how the morphological and functional differences of the human brain are determined is still not well understood.
To gain more insights into the gene regulatory processes that might underlie human-specific brain evolution, we investigate here how a TF co-expression network evolves in the primate PFC. To do so, we analyzed genome-wide expression data from PFC samples of humans, chimpanzees, and rhesus macaques to first determine the genes that are specifically changed in each species. In total, we identified 645 genes coding for TFs that show lineage-specific expression, among them 134 known to be involved in brain development, functions, and/or diseases. We then derived weighted topological overlap (wTO) networks from the changed TFs and their correlated genes and compare these networks between the three species to infer the ancestral network and evolutionary network changes in the human and chimpanzee lineages. To further evaluate which evolutionary changes might be specific to the brain, we used genome-wide expression data from multiple tissues. We found increased cooperation of TFs in the human brain compared to the chimpanzee brain, but not in the other tissues. We further showed that the network of the human PFC is enriched for TFs implicated in crucial brain functions and regulates genes involved in brain development, neuronal functions, synaptic plasticity, cognition, and others. In addition to these functionally characterized TFs, we suggest that at least some of the other TFs with prominent positions in the human PFC network or drastically changed connections are very good candidates for being also involved in human-specific functions of the PFC.

Data Sets
Raw RNA-Seq, Tag-Seq, and microarray data were downloaded from Gene Omibus Database (http://www.ncbi.nlm. nih.gov/geo/). For the differential expression profiling, we used Tag-Seq data of the PFC of 5 adult human, chimpanzee, and rhesus macaque individuals (GSE50782). For the correlation and network analyses, we used a microarray data set of PFC samples from which we selected 12 individuals for each species with different ages (GSE22521) (Somel et al. 2011) (supplementary table S8, Supplementary Material online). For a comparable age collection, we implemented a linear model using the specific life traits of each species such as sexual maturity, first reproduction, age at gestation, litter per year, weaning, and maximum life expectancy (Somel et al. 2011). We used a second data set to confirm network patterns: an RNA-Seq data set comprising multiple tissues of six adult human, chimpanzee, and rhesus macaque individuals (GSE49379) (Bozek et al. 2014).
Expression Profiling RNA-Seq and microarrays were analyzed using the R programming language and Bioconductor packages. RNA-Seq reads were aligned to primate genomes (hg19, panTro3, rheMac3) using segemehl (Hoffmann et al. 2009). Unmapped and multi-mapped reads were further removed. Counts and Reads Per Kilobase of transcript, per Million mapped reads were calculated using GenomicRanges and biomaRt implemented on R (Durinck et al. 2005;Lawrence et al. 2013). LiftOver tool where implemented to translate the non-human primates coordinates into hg19 coordinates (Hinrichs et al. 2016). We retained expressed genes with RPKM > 0.5 in at least one species. Only orthologous genes were used and human gene names were selected for further analysis. Differential expression was calculated using the DESeq package (Anders and Huber 2012). Genes were defined as differentially expressed between species if having jlog2FCj > 0.3 and False Discovery Rate < 0.05. For the microarray data set, we first analyzed the RNA-degradation patterns with the Affy package in R (Gautier et al. 2004). We identified no significant differences between the human, chimpanzee, and rhesus macaque samples. We next performed a computational mask procedure using the maskBAD package (http://bioinf.eva.mpg.de/masking/) (Dannemann et al. 2009). This removed probes with binding affinity differences between species. For further analysis, we only considered the probe sets with more than four probes left after masking. We determined gene expression levels (Robust MultiArray Average values) and MAS5 detection P value from the remaining probes using the "affy" package (Gautier et al. 2004). We considered only the probesets significantly detected in at least one individual (P < 0.05). Furthermore, for genes represented by more than one expressed probeset, we calculated the mean of the expression values of all its probesets.

Gene Sets
The list of all TFs was taken from TFcheckpoint (Chawla et al. 2013) in which they selected and manually curated genes coding for TFs in the human genome. For our analysis, we only included TFs with orthologs expressed in all three species. The "Brain TF" gene set was manually curated from different and independent sources, databases, and studies (Polymeropoulos 2000;Inlow and Restifo 2004;Greydanus and Pratt 2005;Hamosh et al. 2005;Buxhoeveden et al. 2006;Ropers 2008;Bertram 2009;Banerjee-Basu and Packer 2010;Jia et al. 2010;Kaufman et al. 2010;Darnell et al. 2011;van Bokhoven 2011;Voineagu et al. 2011;Iossifov et al. 2012;Lill et al. 2012;Parikshak et al. 2013; Schizophrenia Working Group of the Psychiatric Genomics Consortium 2014) and contains candidates of genes coding for TFs implicated with brain development, neurogenesis, and brain disorders (Berto et al. 2016). We required a positive signal from more than 2 GWAS studies for Alzheimer's disease, Parkinson's disease, and Schizophrenia (Bertram et al. 2007;Allen et al. 2008;Lill et al. 2012). The Autism gene set was downloaded from SFARI (Banerjee-Basu and Packer 2010).

Enrichment of cis-Regulatory Elements
Human PFC H3K27ac peaksets where downloaded from GSE67978. Genome coordinates were translated from hg38 to hg19 using liftOver (Hinrichs et al. 2016). A consensus peakset was calculated using bedtools (Quinlan and Hall 2010), retaining only peaks present in all the three replicates. Peak annotation was performed using HOMER (Heinz et al. 2010), resulting in 5,967 putative CREs located near genes. TF annotation and enrichment was performed using Enrich (Chen et al. 2013).

Correlation Analysis
We performed Spearman rank correlations between the expression values of each expression-changed TF and all expressed genes. To derive the wTO networks incorporating all significantly (P < 0.05) correlated genes, we calculated the wTO values as previously described (Nowick et al. 2009;Berto et al. 2016). Briefly, we calculated a wTO matrix starting from the adjacency matrix A ¼ [a ij ], with a ij ¼ Corr ðijÞ [À1, 1] or a ij ¼ 0 if i¼ j, where i and j represent the differentially expressed TFs. In contrast with method previously described (Zhang and Horvath 2005), our method incorporates the correlations of two TFs associated gene sets denoted as u. Our approach further considers positive and negative correlations as following: a ij [À1, 1] when a ij ! 0 ! a iu a uj ! 0 for all u and a ij [À1, 1] when a ij 0 ! a iu a uj 0 for all u. Inserting the weighted connectivity of a node i as ja ij j, the wTO (x) is calculated as:

Network Robustness Tests
To test the robustness of the networks, we performed two different methods: Firstly, we performed a permutation test shuffling 1,000 times the expression values of all expressed genes for each individual. We then calculated the wTO values with these randomized expression values. The randomized networks showed fewer links and high structural differences compared with the empirical networks of all species for all tested cutoffs jwTOj¼ [0.2.0.6], resulting in a P value of 0.001. None of the randomized networks displayed edges with jwTOj > 0.4. Therefore, we chose jwTOj > 0.4 as cutoff for the network comparisons. Robustness of the networks was examined through permutation and leave-one-out methods (jwTOj > 0.4; permutation test, P value < 0.001; supplementary fig. S3, Methods, Supplementary Material online). For the EC-sub-network, we also performed permutation tests by randomizing the expression values of all expressed genes 1,000 times. Also here, the structures of the randomized networks were considerably different from the empirical network for all tested cutoffs jwTOj¼ [0.2.0.6]. Since none of the 1,000 shuffled networks presented wTO values higher than 0.3 (P < 0.001), we chose jwTOj > 0.3 as cutoff for the EC-sub-networks. To test if the high number of species-specific links could be an artifact of particular individuals, we recalculated the networks using the "leave-one-out" method. This resulted in 12 networks per species constructed from 11 individuals each. All these networks clustered according to species, demonstrating that the strong divergence in network links between species is robust (supplementary figs. S7 and S8, Supplementary Material online).

Additional TF Enrichment
We performed the TF motif enrichment using the Jaspar and TRANSFAC databases (Matys et al. 2003;Mathelier et al. 2014). We compared the 5 kb upstream promoter regions to three different background data sets: 5 and 2 kb promoter regions of all human genes and human CpG islands. To perform the motif enrichment, we used the MEME suite (Bailey et al. 2009).

Other Statistics
To test for enrichment of "Brain-TFs" in the EC-sub-networks, P values were calculated with one-sided Fisher's exact tests (confidence level ¼ 0.99, simulated P value with 1,000 replicates). A one-sided Wilcoxon ranked test was implemented to evaluate the difference of the network connectivity between human compared with non-human primates (alternative ¼ "g," confidence level ¼ 0.99, paired ¼ FALSE). P values for the overlaps of differentially expressed TFs between data sets (Somel et al. 2011;Liu et al. 2012) and for differentially expressed TFs with "Brain-TFs" were calculated with the hypergeometric test. In each case we used an independent background for population size based on genes expressed in the human, BrainSpan data set (15,585 genes; Parikshak et al. 2013). Two-way permutation tests of 10,000 were employed to validate the overlaps. First we randomized the external gene sets (e.g. human-specifically changed genes) by randomly selecting the same number of genes from an independent brain expressed genes list (BrainSpan gene set ¼ 15,585 genes) and subsequently calculating the overlap P values with the TF gene set. The second approach randomized the internal gene sets (e.g. TF gene set) by randomly selecting the same number as TFs detected from the expressed genes and subsequently calculating the overlap P values. Moreover, we used a permutation test to evaluate the detected differentially expressed genes, randomizing 1,000 times the RNAseq data and recalculating the differentially expressed genes detecting that none of the permuted data showed the same number of differentially expressed genes.

Interactive Networks and Input Tables
Using D3.js, we developed a method to visualize the wTO networks and associated conservation or specificities of the links per each lineage. The interactive networks are associated with the manuscript as additional files and can be downloaded here: http://www.nowick-lab.info/? page_id¼470. Input tables containing TFs and their correlated genes are available upon request.
We ranked all genes based on the number of TFs the genes are correlated with and used the Wilcoxon ranked test implemented in FUNC for testing for enrichment of GO groups. We report GO groups with enrichment P values < 0.05 before and after refinement.

Results
TFs with Known Functions in the Brain are Enriched among Human-Specifically Changed TFs in the PFC To identify species-specific expression patterns we analyzed genome-wide expression data, derived from PFC samples of 5 adult human, chimpanzee, and rhesus macaque individuals (see Materials and Methods section). Genes were defined as species-specifically changed if their difference in expression was significant (DESeq, FDR < 0.05, jlog2 fold changej > 0.3) in one species compared to the other two species, but not significant between the other two species ( fig. 1; supplementary table S1, Supplementary Material online). Among the genes with species-specific expression changes, we found 645 genes coding for TFs, consisting of 103 humanspecifically changed TFs, 80 chimpanzee-specifically changed TFs, and 462 rhesus macaque-specifically changed TFs. This represents a significant enrichment of TFs among differentially expressed genes (8%, chi-square test, P ¼ 0.02).
Due to its distant evolutionary relationship with great apes, we found the highest number of specifically changed genes in rhesus macaques. However, when we normalized the number of species-specific expression changes for divergence time, we found about equal numbers of changes in all three lineages (e.g. number of TF/divergence time in million years between species and common ancestor ¼ H: 103/6 ¼ 17.2; C: 80/6 ¼ 13.5, R: 462/25 ¼ 18.5), suggesting that overall gene expression changes are similar between lineages.
We validated the species-specific expression changes with data of an independent primate comparative study of the frontal lobe (Konopka et al. 2012), and observed significant overlap of species-specific differential expression (human-specific: P ¼ 0.04; chimpanzee-specific: P ¼ 0.02; rhesus macaque-specific: P ¼ 0.001; hypergeometric test followed by permutation test, P < 0.0001). In addition, humanspecifically changed genes are also over-represented among the genes with human-specific PFC expression detected in a study that used multiple different brain regions (Sousa et al. 2017) (P ¼ 2.40eÀ06; hypergeometric test followed by permutation test, P < 0.0001), supporting our findings.
To conjecture potential impacts of the species-specific TF expression changes on species differences in brain functions, we first asked how many of the changed TFs are known to have a role in the brain. Our literature review discovered 134 changed TFs that are described to have a function during brain development or are implicated in a brain disease (Berto et al. 2016) (we will refer to them here in short as "Brain-TFs"; see Materials and Methods section). "Brain-TFs" are specifically overrepresented among human-specifically changed TFs (27 out of 103; Fisher's exact test, P ¼ 0.028) but not in chimpanzee (14 out of 80; Fisher's exact test, P ¼ 0.62) and rhesus macaque (93 out of 462; Fisher's exact test, P ¼ 0.15). Among these human-specifically changed "Brain-TFs" (fig. 2), are for example CLOCK, a circadian regulator involved in cognitive disorders (Vitaterna et al. 1994;Gekakis et al. 1998;Menet and Rosbash 2011;Fontenot et al. 2017) and a hub in a human-specific PFC co-expression module (Konopka et al. 2012); CC2D1A, which is implicated in non-syndromic mental retardation (Basel-Vanagaite et al. 2006;Rogaeva et al. 2007) with humanspecific expression (Konopka et al. 2012); NR1D1, involved each species has been used to identify the species-specifically differentially expressed genes and TFs. Differentially expressed TFs and Genes were used for calculating two types of networks, (B) for inferring network evolution and (C) species-specifically changed EC-sub-networks, using a stringent criteria for correlation and nominal P value cutoffs. (B) Network evolution: We calculated Spearman rank correlations for each of the TFs with species-specific expression with all expressed genes. Correlated genes were filtered according to the criteria shown in red in each box, whereby pval stands for the P value of the correlation and rho for the correlation strength, which needed to have the same sign (positive or negative) in the species to inferred that a link was present in the networks of the ancestors of these species. We then calculated a wTO network from all genes that passed the respective filtering criteria for humans, chimpanzees, rhesus macaques, the HC-, and the HCR macaque-ancestor. A comparison of these five networks allowed us to investigate the evolution of network links. (C) Species-specific EC-sub-networks: For the species-specific EC-sub-networks we only considered TFs that were changed in expression in the respective species. Their correlated genes (Spearman rank correlation, P<0.05) were filtered for also being species-specifically expressed in the same species and for displaying an expression change that is in the direction that is in agreement with the direction of the expression change of the TF and the sign of the correlation to that TF (see text). In blue, species-specifically upregulated TFs and correlated genes; in orange, species-specifically downregulated TFs and correlated genes. The wTO of the species-specific sub-networks were calculated from the genes that passed this filter. in brain development and autism (Goto et al. 2017); and EGR1, a gene implicated in brain plasticity, cognitive disorders, and social behavior (Robinson et al. 2008;Duclot and Kabbaj 2017).

Species-Specific Changes in TF Networks
Only about 200 TFs are functionally characterized (Matys et al. 2003;Wang et al. 2013;Mathelier et al. 2014). Accordingly, the functions of many of the TFs with species-specific expression change are currently unknown (Encode Project Consortium 2004;Wang et al. 2013;Mathelier et al. 2014). Thus, we analyzed co-expression patterns of TFs to gain more insight into the functions of the species-specifically changed TFs and into the potential phenotypic impact of their expression changes. Further, since TFs with similar sets of coexpressed genes are likely functionally related, we aimed at capturing the co-expression patterns of the changed TFs and their similarities using a network approach. To not bias this analysis by the expression patterns that exist in the data sets we used for discovering differentially expressed genes, we utilized another independently derived data set for our network analyses (Somel et al. 2011). From this data set, we selected 12 age-matched individuals per species (see Materials and Methods section). We compared these TF networks between the three species aiming to answer two main questions: 1) How did the TF network evolve? 2) What underlies species-specific expression changes?

Connectivity of the TF Network Increased in the Human Brain
To answer the first question, we first identified for each of the species-specifically expressed TFs the genes with correlated expression patterns across the individuals of a species ( fig. 1). Since TFs can activate or repress the expression of genes, we calculated positive and negative correlations. To analyze the overlap in the correlated gene sets between the TFs, we calculated the wTO using a method we developed previously that considers both, positive and negative correlations (Nowick et al. 2009). This allowed us constructing a wTO network for each species in which the nodes represent the TFs and the links the correlations between the TFs including the commonality of the TFs in their sets of correlated genes. From a biological perspective, TFs that are linked in the wTO network might cooperatively regulate a significant set of potential target genes.
We inferred the human-chimpanzee (HC) and the humanchimpanzee-rhesus (HCR) macaque ancestral networks based on correlations that are present in humans and chimpanzees or in all three species, respectively ( fig. 1B). Using the rhesus macaque as outgroup, we further determined the network links that are likely specific to either the human or chimpanzee network ( fig. 3). Only 531 links are shared between the three species (HCR ancestor) and 239 links between humans and chimpanzees (HC ancestor). In contrast, the human network contains 2,238, the chimpanzee network 1,113, and the rhesus macaque network 389 specific links (supplementary table S2, Supplementary Material online). Robustness and correlation Quality Control were tested by investigating different wTO cutoffs and using a "leave-one-out" test (supplementary methods, Supplementary Material online). In addition, to rule out the possibility that the differences in the number of links is driven by a general species difference in the number of correlations across all genes (e.g. caused by biological or technical differences), we compared the correlation distributions of all genes after P value cutoff (supplementary fig. S1, Supplementary Material online). All three species showed similar patterns of that distribution, suggesting that the species difference in connectivity is due to the number of genes correlated with TFs rather than an overall difference in correlations.
The human network has significantly higher connectivity (number of links per nodes; c ¼ 13.2; jwTOj > 0.4) than the chimpanzee network (c ¼ 8.1, jwTOj > 0.4; Wilcoxon test, P value ¼ 2.47 Â 10 À10 ) and the rhesus macaque network (c ¼ 3.9, jwTOj > 0.4; Wilcoxon test, P value ¼ 2.2 Â 10 À16 ). Higher connectivity points to an increased cooperation or FIG. 2.-Human-specific differential expression. Expression patterns of human-specifically changed "Brain TFs," that is, TFs that are known to be involved in brain functions and disorders are displayed. Shown are Z-scores. Red ¼ high expression, blue ¼ low expression in human PFC. coordination of TFs in the human compared to the chimpanzee and rhesus macaque brain. Taken together, this indicates that the network complexity increased on the human lineage.
Our data further allow us to follow how the network architecture has been changed during evolution. While, for example, BBX ( fig. 4A), FOXG1, RAB37, ZMAT3, and ZNF436 are hubs in almost all networks, STAT6, is only a hub in the human, chimpanzee, and HC ancestor network, suggesting that it derived this status on the lineage to great apes. Several hubs are human-specific, that is, possess a relatively large number of human-specific links, such as CC2D1A ( fig. 4B), GLIS3, KLF5, MEF2D, and ZNF286A. We invite the reader to further explore our interactive network visualization FIG. 3.-Ancestral and species-specific links in the TF wTO network of the PFC. In light blue are links common to the human, chimpanzee, and rhesus macaque network; in purple are links common to the human and chimpanzee network; in green are rhesus macaque-specific links; in red are chimpanzeespecific links; and in blue are human-specific links. Highlighted are the hubs of each network. Note that we cannot predict links that have been lost during evolution. (supplementary data 1, Supplementary Material online) to discover more details about changes in this network during evolution.
Having found an increased connectivity between TFs in the human network, we asked whether this increase can also be found in other brain areas or in other tissues (Bozek et al. 2014). Networks clustered well by tissue ( fig. 5A). We found a higher extent of connectivity changes in human brain regions compared with the non-human primates ( fig. 5B). Human networks of all examined brain areas (PFC, cerebellum, and visual cortex) were characterized by higher connectivity compared to chimpanzee networks. However, we did not observe an increase in connectivity in the examined other tissues, muscle, and kidney ( fig. 5C; supplementary

Human-Specific Expression Changes are Associated with Higher TF Cooperation
While all TFs in the presented networks have changed in expression, it does not mean that the genes correlated with these TFs have also changed in expression. In fact, many of the lineage-specific network changes could have evolved to compensate for other mutations to keep the expression of the associated genes conserved. Since the genes with speciesspecific expression changes are most likely the ones that drive phenotypic differences between the three species, we determined next which TF correlated genes have speciesspecifically changed in expression.
To this end, we filtered all TF correlated genes requiring that their expression change in the PFC is consistent with the expression change of the TF in the PFC (supplementary methods, Supplementary Material online; fig. 1C). We then constructed another TF wTO network for each species. In contrast to the networks above, these TF wTO networks contain only the TFs specifically changed in expression in that species, which implies that the node sets do not overlap between the networks of the three species. Further, the wTO is only measured from correlated genes with species-specific expression changes in the respective species. Despite the conceptual difference between the two types of networks, the node set of the second type is a subset of the first type. The link set of the second network type is not an exact subset of the first network type, but wTO values correlate between the networks (Pearson correlation, r ¼ 0.13, P < 2.2 Â 10 À16 ). We thus call the second type of networks species-specific expression-changed networks (EC-sub-networks), albeit they are not strict sub-networks of the first ones.
The EC-sub-networks are densely connected with high clustering coefficients (human C ¼ 0.  . 6). The human EC-sub-network has significantly higher connectivity than the chimpanzee EC-sub-network (human c ¼ 7.8, chimpanzee c ¼ 3.4, Wilcoxon test, P ¼ 4.85 Â 10 À05 ). Consistently with this result, we also found higher connectivity in the EC-sub networks derived from the PFC samples of Bozek et al. (2014) (human c ¼ 15.7, chimpanzee c ¼ 6.8, Wilcoxon test, P ¼ 3.3 Â 10 À09 ) and significant overlap with the links determined with the Bozek et al. (2014) data set (human: hypergeometric test, P ¼ 4.78eÀ29, chimpanzee: hypergeometric test, P ¼ 2.00eÀ05; supplementary figs. S3 and S4, Supplementary Material online). Taken together, this indicates that human-specific expression changes in the PFC are, at least in part, driven by a particular set of TFs with human-specific expression levels in the PFC and higher cooperativity among them.

Binding Sites and Enhancers
Several of the links between TFs, we discovered in the human PFC network had already been discovered experimentally, such as the interactions between MEF2C and HIRA (Yang et al. 2011), MEF2C andHDAC9 (Haberland et al. 2007;Potthoff and Olson 2007) and MEF2D and SP1 (Park et al. 2002), supporting the validity of our network inferences.
To provide further support for our network inferences, we tested for enrichment of TF binding sites in the promoters of genes that are correlated with the respective TF and have changed in expression (see Materials and Methods section). For instance, for the genes contributing to the human EC-sub-network we found enrichment for binding sites of several human-specifically changed TFs (e.g. CLOCK, EGR1, HNF4A, LMO2, PRDM14, and SMAD2), lending support for the biological relevance of the inferred EC-sub-networks.
In addition, we evaluated the relationship between genes contributing to the EC-sub-networks and cis-regulatory elements (CREs) regulated by H3K27ac in the human brain (Vermunt et al. 2016) (see Materials and Methods section). We detected a significant overlap between such CREs and genes of the human EC-sub-network (161 of 470, P ¼ 0.02, hypergeometric test), and confirmed this significant overlap also with the EC-sub-network derived from the Bozek et al. (2014) data set (195 of 529, P ¼ 3 Â 10 À4 , hypergeometric test). For 15 TFs of the human EC-sub-network exists direct evidence of binding sites (Lambert et al. 2018). To refine the CRE enrichment, we next tested whether binding sites for these TFs are within the CREs near genes of the human EC-sub-networks. For EGR1, CLOCK, KLF8, and SP4 we detected significant enrichment of their binding sites (see Materials and Methods section). These results reflect that genes in the human ECsub-networks are associated with active enhancers that are potentially bound and regulated by TFs of the humanspecific EC-sub-network.

Human-Specific Network Changes Seem to be Related to Cognitive Functions
To gain insights into the TF network functions, we tested for enrichment of GO groups among the genes correlated with the TFs (see Materials and Methods section). The wTO networks of all three species, as well as the putative ancestral networks of humans and chimpanzees and of all three species are enriched for axonogenesis, synaptic transmission, learning and memory, and other brain functions (supplementary table  S4 and fig. S5, Supplementary Material online). This suggests that, although species-specific links exist, overall the functions and pathways regulated by the TFs in the PFC are conserved. Still, the strongest enrichment for brain function related GO groups was for networks of the human PFC, while that enrichment was lower in the chimpanzee and rhesus macaque PFC and in the other brain areas and tissues, indicating that human-specific links are involved in regulating functions important for brain development and cognition.
To further test this possibility, we asked which functional groups are over-represented among the genes of the EC-subnetworks. The human EC-sub-networks are enriched for genes involved in axon guidance, myelination, and cell differentiation. Such functions are not over-represented in the chimpanzee EC-sub-networks (supplementary table S5, Supplementary Material online). This is remarkable given that very similar GO groups have been enriched in the chimpanzee and ancestral networks built from all correlated genes. While the overall function of the PFC network seems to be  Bozek et al. (2014); in blue, the cerebellum (CBC) samples of Bozek et al. (2014); in black, the PFC; in brown, the kidney (KD) samples of Bozek et al. (2014); in pink, the muscle (MSC) samples of Bozek et al. (2014). The networks of the rhesus macaque PFC and CBC are the most different ones. (B) Change in connectivity per million years of all TF wTO networks. lsPFC refers to the data set from Somel et al. (2011), while adPFC represents the Bozek et al. (2014) data set. Human networks have a higher number of changes compared with the other primates in brain regions. (C) One-sided Wilcoxon signed rank test comparing connectivity between human and non-human primates. Humans showed a greater connectivity compared with chimpanzee and rhesus macaques in all brain regions whereas not in kidney and muscle. conserved since the HCR ancestor, TF genes with humanspecific expression changes seem to particularly affect the expression of genes involved in certain brain processes such as axon guidance, myelination, and cell differentiation.

Species-Specific Changes in a Primate Transcription Factor
Indeed, "Brain TFs" had more links than other TFs in the human PFC (Wilcoxon test, P ¼ 0.035; Bozek et al. 2014: Wilcox test, P ¼ 0.006) compared with chimpanzee PFC (Wilcoxon test, P ¼ 0.65; Bozek et al. 2014: Wilcox test, P ¼ 0.53). Additionally, brain-TFs have more links in the human PFC compared with the chimpanzee PFC (Wilcoxon test, P ¼ 0.08; Bozek et al. 2014: Wilcoxon test, P ¼ 1.04 Â 10 À05 ), suggesting a more central role for those "Brain-TFs" in the human PFC network (supplementary fig. S6, Supplementary Material online). Examples of hubs in the human EC-subnetworks are the aforementioned CC2D1A, and ZNF24 and ZNF536, two zinc finger genes implicated in maintenance of neural progenitor cells (Khalfallah et al. 2009) and neuronal differentiation (Qin et al. 2009 To further evaluate the association of EC-sub-networks with brain functions, we used cell-type markers from singlecell RNA-Seq (sc-RNA) (Zhang et al. 2014). Genes of the human EC-sub-networks are significantly enriched for genes expressed in myelinating oligodendrocytes (Fisher's exact test, P ¼ 0.0008;Bozek et al (2014):, Fisher's exact test, P ¼ 0.003; Benjamini-Hochberg adjusted). In contrast, genes of the chimpanzee EC-sub-networks did not show any significant enrichment in any cell-type. We further confirmed the enrichment for cortical oligodendrocytes expressed genes in the human EC-sub networks with CSEA (Xu et al. 2014). Taken together, these data support the notion that the network with genes with human-specific expression changes is involved in myelination.
Because the EC-sub-networks are bi-modular, we also tested for GO enrichment among the genes of each module (see Materials and Methods section). While the chimpanzee modules did not show any significant enrichment, one human modules was enriched for genes involved in cellular differentiation and morphogenesis and the other one for axon guidance, synaptic plasticity, learning and memory, cognition, and brain development (supplementary table S7

Discussion
While a few studies have already analyzed co-expression networks in primate brains (Oldham et al. 2006;Konopka et al. 2012;Sousa et al. 2017), we particularly focused here on the evolution of TF species-specific networks in PFC. We found that TFs with a known function in the brain are enriched among the TFs with human-specific differential expression in the PFC. Using our system-level network approach, we determined potential interactions and target genes of all TFs with species-specific expression changes. We represent them in two types of networks. The first network type allowed us to directly compare links between all species and to inferred the ancestral HC and HCR networks as well as species-specific network features. The second type of networks are built only from TFs and other genes with species-specific expression changes. These networks enabled us to investigate TF interactions that might be responsible for driving lineage-specific changes in expression patterns. Higher connectivity in humans was observed with both types of networks, suggesting that cooperativity between TFs is stronger in the human brain. We further showed that this increase in connectivity is specific to the brain, while it was not observed in the investigated nonbrain tissues. Similar GO groups have been found to be enriched among genes of all PFC networks; however, functions related to brain development and cognition were most pronounced in the human networks, and "Brain-TFs" were only enriched in the human EC-sub-networks. Our findings provide context to previous observation that gene co-expression might have undergone strong evolutionary remodeling in primates and point out not only TFs that had already previously been linked with human-specific networks but also novel candidates for human brain evolution and cognitive uniqueness.
Our findings provide context to previous observations of increased network connectivity in the human compared to non-human primates (Konopka et al. 2012). In our results, higher connectivity between TFs means that the TFs overlap more strongly in their potential target genes and interaction partners. Such higher connectivity can be interpreted as an increase in the complexity of gene regulatory mechanisms in the human brain. The fact that these two different methods applied to different data sets both found higher connectivity in the human PFC compared to other primates strongly supports the idea that the transcriptional landscape in human has undergone marked rewiring during evolution. This connectivity increase might not be restricted to the PFC, as we have observed it also in other brain areas, albeit not in other tissues. However, at date, accessibility to non-human primates tissues and sample size are limiting factors. Future studies that assess TF networks in different human brain regions compared with non-human primates should address these possibilities.
Our method allowed us to directly compare the network links between species, inferring ancestral networks, and the species-specific gain/loss of TF connections. For example, CC2D1A, a gene implicated in non-syndromic mental retardation (Basel-Vanagaite et al. 2006;Rogaeva et al. 2007), has gained 91 human-specific links. Loss of CC2D1A reduces dendritic complexity and regulates NF-kB signaling pathways implicated in neuronal morphology (Manzini et al. 2014). Therefore, the increased connectivity of CC2D1A might underline a more complex transcriptional regulation in the human lineage compared with non-human primates. We therefore hypothesized that CC2D1A might have had an important role during the evolution of the human PFC transcriptional landscape.
Despite the evolutionary differences, the ancestral as well as the species-specific TF networks are predicted to regulate the expression of similar sets of genes involved in axonogenesis, synaptic transmission, learning, and memory. This indicates that overall the functional output of the network is conserved. Our finding agrees with the hypothesis that many new network links can arise non-adaptively and can be added to or deleted from the network without changing its function (Sorrells and Johnson 2015).
An alternative explanation is that some of the rewiring we observed might compensate for other mutations during primate evolution to overall keep expression of genes in the brain conserved. This implies that in part some molecular pathways are regulated differentially across different primate species.
This prompted us to specifically investigate which TF changes might be related to expression changes of their potential target genes, reasoning that such changes should be more relevant for changing phenotypes. With this assumption in mind we calculated EC-sub-networks consisting of only the TFs and their potential targets that show species-specific expression patterns. These EC-sub-networks show an organization in two modules, with the human networks having higher connectivity than the chimpanzee networks. This data emphasize the complexity of the transcriptional regulation, often represented by modules of TFs that bind cooperatively or competitively promoter regions of target genes (Berman et al. 2002;Nowick et al. 2009). Interestingly, many hubs in the human EC-sub-network are "Brain-TFs." Genes potentially regulated by TFs in the human EC-sub-network are enriched for axon guidance, myelination, and cell differentiation. Furthermore, genes that are specifically correlated in human but not in chimpanzee brains are significantly enriched for genes expressed in myelinating oligodendrocytes. This result is intriguing since the evolutionary dynamic of the human brain and cognitive enhancement are subjected to prolonged myelination when compared with close relatives such as chimpanzees (Roth and Dicke 2005;Miller et al. 2012). One of the modules of the human EC-sub-network, which includes for instance BBX, CLOCK, and ZNF24 as hubs, is linked to functions such as cellular proliferation and migration, implicating TFs in proliferative function in the human brain. The other module of the human EC-sub-network, with for instance CC2D1A, HOXD1, and KCNH3 as hubs, was strongly and specifically enriched for GO groups associated with synaptic plasticity, learning and memory, and cognition making that module an excellent candidate for setting the stage for the evolution of human-specific cognitive abilities.
Despite the fact that the samples we analyzed are bulk RNA-Seq from brain tissue, this allowed us to highlight genes that represent potential target genes or interaction partners of these TFs by co-expression. In fact, our results suggest that a network of TFs that regulates genes involved in brain development and cognitive processes has changed during primate evolution. Our work not only highlights the complexity of transcriptional networks in human brain regions with a focus on the PFC but also adds to previous findings on humanspecific morphological changes in the PFC (Rilling et al. 2008;Semendeferi et al. 2011) and human-specific gene expression changes in the PFC (Somel et al. 2009;Babbitt et al. 2010;Konopka et al. 2012;Liu et al. 2012;Sousa et al. 2017;Xu et al. 2018), by suggesting candidate TFs and interactions that might drive these human-specific changes. In this respect our study points to TFs that should be further investigated in the future for understanding the evolution of the human brain.
Moreover, our observation of enrichment for and higher connectivity of "Brain TFs" specifically in the human PFC network, suggests that human-specific changes in expression and network integration of some "Brain-TFs" might in part be associated with the evolution of human-specific cognitive abilities. Interestingly, among the most highly connected TFs in the human EC-sub-network are TFs that have been implicated in cognitive disorders. It is possible that these TFs became risk genes for brain disorders, because they moved into such central position in the human network.
We yet have to better understand the complexity of gene regulatory networks and their phenotypic consequences, but the TF network changes we identified here might have changed the expression of genes that are involved in determining human-specific traits, such as bigger brain size, particular cognitive abilities, behavior, and brain disorders.

Supplementary Material
Supplementary data are available at Genome Biology and Evolution online.