Accumulation of CTCF-binding sites drives expression divergence between tandemly duplicated genes in humans

Background During eukaryotic genome evolution, tandem gene duplication is the most frequent event giving rise to clustered gene families. However, how expression divergence between tandemly duplicated genes has emerged and maintained remain unclear. In particular, it is unknown if epigenetic regulators have been involved in the process. Results We demonstrate that CCCTC-binding factor (CTCF), the master epigenetic regulator and the only known insulator protein in humans, has played a predominant role in generating divergence in both expression profiles and expression levels between adjacent paralogs in the human genome. This phenomenon was not observed for non-paralogous adjacent genes. After tandem duplication events, CTCF-binding sites gradually accumulate between paralogs. This trend was more prominent for genes involved in particular functions. Conclusions The accumulation of CTCF-binding sites drives expression divergence of tandemly duplicated genes. This process is likely targeted by natural selection. Our study reveals the importance of CTCF to the evolution of animal diversity and complexity. Electronic supplementary material The online version of this article (doi:10.1186/1471-2164-15-S1-S8) contains supplementary material, which is available to authorized users.


Background
Gene duplication is a major driver for the emergence of organismal complexity and evolutionary innovations [1,2]. Tandem duplication is the most common route to the formation of clustered paralogous genes [3,4]. A newly duplicated gene must diverge from its progenitor gene in coding sequence or expression, or it will degenerate due to redundancy [1,5,6]. To acquire novel transcription patterns, tandemly duplicated genes need to interrupt expression similarity due to shared upstream cis-elements upon origin [3,7] and transcriptional interference contributed by physical proximity [8]. Despite the challenges, functionally important gene clusters consisting of paralogs with distinct expression patterns are found in a wide range of species, including human [9,10]. Therefore, it is important to understand the origin and maintenance of expression divergence between tandem paralogs.
CCCTC-binding factor (CTCF), the only known human insulator protein, plays a master role in determining the transcriptional landscape of genomes. When bound to insulator sequences (CTCF-binding sites), CTCF can prevent repressive heterochromatin from spreading into neighboring regions [11,12]. In addition, CTCF interferes with enhancer-promoter communication [13] and guides long-range chromatin interactions [14]. Because changes in epigenetic marks play roles in regulatory divergence of duplicated genes [15][16][17], we hypothesize that CTCF also plays a role. Using human RNA-seq and ChIP-seq data (see Methods), we examined how CTCF drives regulatory divergence of duplicated genes, especially tandemly arrayed paralogs.

Results and discussion
In the genome, adjacent genes are coexpressed due to a common origin [18], cofunctionality [19], or deleterious  [8]. Using RNA-seq data from six human adult tissues [20], we measured gene expression dissimilarity between adjacent genes in terms of expression profile or expression level using ExpD 1-r or ExpD Euc (see Methods), respectively. ExpD 1-r focused on changes in the shape across the tissue dimension, while ExpD Euc focused on summed changes in abundance. We studied divergences in both profiles and abundances because previous studies have shown that the observations from the two aspects did not necessarily produce consistent results [21], possibly due to different underlying mechanisms controlling the properties [22].
Mammalian genes close to each other have similar expression profiles [8,19]. CTCF-binding sites can prevent undesirable crosstalk between active and inactive genomic regions [23]. DNA methylation upstream of a gene can inhibit its transcriptional initiation [24]. To determine if the effect of intergenic distance, the number of CTCF-binding sites, and DNA methylation is related to shared evolutionary origin, we measured ExpD 1-r and ExpD Euc in non-paralogous adjacent genes. Non-paralogous gene pairs with a longer intergenic distance (d, calculated as the distance in nucleotides between the transcriptional start sites, TSS), more overlapping CTCF-binding sites (#CTCF, see Methods), or a larger difference in upstream DNA methylation (ΔCpG O/E , see Methods) had greater ExpD 1-r or ExpD Euc (see ρ in Table 1 Table 1). Although CTCF binding can vary among cell types [25], when we define #CTCF using joint CTCF ChIP-seq peaks (see Methods) instead of overlapping CTCF peaks, the results did not change (Table S1 in Additional file 1).
Mechanisms generating coding sequence divergence, such as change in protein structure or splicing, between paralogous genes have been intensively investigated [2,[26][27][28]. By contrast, mechanisms generating expression divergence have garnered less attention. Although #CTCF was not the strongest determinant of expression divergence in non-paralogous adjacent gene pairs (Table 1), for adjacent paralogs, #CTCF had the greatest direct influence on ExpD 1-r and ExpD Euc , followed by Table 2). Repeating analysis based on an independently generated and unpublished RNA-seq dataset including 16 human tissues (Illumina BodyMap 2.0 project, see Methods) produced a result consistent with Table 1 and 2 ( Table S2 in Additional file 1), suggesting the robustness of the pattern found.
We classified adjacent paralogs into three groups according the orientations: head-to-head, head-to-tail and tail-to-tail [29]. Although reduced sample sizes resulted in reduced statistical significance, the pattern of the greatest impact #CTCF on expression divergences held regardless of the orientation of paralogs (Table 3). Using microarray expression data, a previous study found that intergenic distance was related to expression divergence for tandemly arrayed paralogs [30]. However, microarray data is known to have cross-hybridization related biases [17], and we found no direct association between intergenic distance and expression divergence for tandemly arrayed paralogs after controlling for #CTCF and ΔCpG O/E (see ρ p of d vs. ExpD 1-r or ExpD Euc , Table 2). Taken together, these results suggest that CTCF-binding sites play a very significant, if not primary, role in driving expression divergence of tandemly duplicated genes.
There are two hypotheses to explain the influence of CTCF-binding sites in driving expression divergence of tandem paralogs. First, tandem paralogs that arose in genomic regions with high densities of CTCF-binding sites nearby are more likely to be preserved due to immediate independence of gene regulation. Second, CTCF-binding sites accumulated over time to enhance independent gene regulation of tandem paralogs, especially those have been functionally diverged. If the first Table 1 Rank correlations (r) and partial rank correlations (r p ) of the genomic properties with expression dissimilarities (measured by ExpD 1-r or ExpD Euc ) of the non-paralogous adjacent genes  Table 4). By contrast, the insignificant partial correlation of d with respect to d S (or T phy ) indicated that the increase d over time can be explained by the increase in the number of CTCFbinding sites and the associated changes in DNA methylation.
To determine whether tandem paralogs with specific functions tend to have a greater number of intervening CTCF-binding sites, we performed enrichment analyses in Gene Ontology (GO) terms. To eliminate the potential effect of duplicability in GO analysis [4], paralogous and non-paralogous gene pairs were analyzed separately.
To control for the potential effect of gene density [31], #CTCF/d was used instead of #CTCF (although using #CTCF produced a consistent result, which is not shown). Gene pairs in the top quartile of CTCF-binding site density (#CTCF/d) were compared against the bottom three quartiles. Enriched/depleted GO terms for the paralog group were substantially different from those in the non-paralog group (Table 5). Only two GO terms (GO:0010033, response to organic substance; GO:0031012, extracellular matrix) had the same enrichment status ("enriched") in both groups. Tandemly duplicated genes with a higher density of intervening CTCF-binding sites tended to specifically encode pro-  Table 5). This result implied that the densities of CTCF Table 2 Rank correlations (r) and partial rank correlations (r p ) of the examined genomic properties with expression dissimilarities (measured by ExpD 1-r or ExpD Euc ) of the paralogous adjacent genes  Table 3 Rank correlations (r) and partial rank correlations (r p ) of the examined genomic properties with expression dissimilarities (measured by ExpD 1-r or ExpD Euc ) of the adjacent paralogous genes of different orientations.
Genomic properties a ExpD 1-r ExpD Euc binding sites between tandem paralogs were not contributed by the genomic background. Genomic regions of high densities of CTCF binding sites can emerge through stochastic evolutionary processes. To examine if high CTCF-binding site density between tandem duplicates is the outcome of gradual CTCF-binding site accumulation by natural selection, we focused on the subset of 278 paralogous pairs in which at least one paralog had one of the abovementioned enriched GO categories. In this subset of tandem paralogs, ρ p of #CTCF vs. d S (or T phy ) after controlling for d and ΔCpG O/E was stronger (#CTCF vs. d S : ρ p =0.276, P<10 -5 ; #CTCF vs. T phy : ρ p =0.243, P<10 -4 ) than that observed in the full set of adjacent paralogs (#CTCF vs. d S : ρ p =0.189, P<10 -9 ; #CTCF vs. T phy : ρ p =0.178, P<10 -9 ) ( Table 4). When #CTCF was defined by joint CTCF ChIP-seq peaks rather than overlapping CTCF peaks, the results were similar (Table S3 in Additional file 1). Therefore, the trend to accumulate CTCFbinding sites following a tandem duplication event resulted in more CTCF-binding sites (both absolute number and density) between paralogs, especially for those with enriched GO-categories shown in Table 5.

Conclusions
Combining human genomic and transcriptomic data, this study demonstrates that CTCF and its binding sites play a major role driving the expression evolution of tandemly duplicated genes. Following tandem duplication events, CTCF-binding sites gradually accumulate between the paralogs to increase their divergences in expression profile and their divergences in expression level. The role of CTCF-binding sites is not limited to the insulation of DNA methylation domains, because the effects of #CTCF on ExpD 1-r (or ExpD Euc ) were still significant even after controlling for ΔCpG O/E (Table 2). Thus, CTCF, a conserved regulatory protein [32], affects the expression evolution of adjacent genes in genomes from flies [33] to humans and is important for the evolution of organismal complexity in animals.

Methods
Annotation of the human genome (Ensembl v72), including gene coordinates, TSS, and paralog divergence, was retrieved through BioMart (http://www.biomart.org/) [34]. Orientation for each paralogous pair (head-to-head, head-to-tail, or tail-to-tail) was determined using strand information. For each paralogous pair, the rate of synonymous changes (d S ) was calculated using PAML [35], and the phylogenetic age (T phy ) was assigned according to Table S4 in Additional file 1 based on Ensembl's annotation of the most recent common ancestor [36]. A smaller d S or T phy indicated a more recent divergence time. The RNA-seq-based gene expression signals in the human brain, cerebellum, heart, kidney, liver and testis (GSE30352) [20] were obtained from NCBI Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/ ). The raw reads of RNA-seq data in 16 human tissues (adipose, adrenal, brain, breast, colon, heart, kidney, liver, lung, lymph node, ovary, prostate, skeletal muscle, testes, thyroid, white blood cells) by Illumina BodyMap 2.0 project were downloaded from GSE30611 of GEO and were processed following our previous studies [17,37] to obtain expression signals. Upstream DNA methylation for a gene was calculated based on the 500 nucleotides upstream of the TSS by CpG O/E =P CpG /(P C × P G ), where P CpG , P C , and P G are the frequencies of CpG dinucleotides, C nucleotides, and G nucleotides, respectively [17,38]. The difference in upstream DNA methylation, ΔCpG O/E , was calculated by the absolute value of the difference in CpG O/E between the two adjacent genes compared. In total, we obtained 30,164 non-paralogous and 1,256 paralogous gene pairs based on 32,164 human genes with detectable expression and estimable CpG O/E .  ExpD 1-r , which represents the dissimilarity in expression profile, was calculated by 1-Pearson's correlation coefficient of the expression signals of the six tissues. ExpD Euc , which represents the summed difference in expression levels, was calculated by the Euclidean distance of the expression signals of the six tissues. CTCF-binding sites identified by ChIP-seq experiments [25] in 13 non-cancerous human cells (Table S5 in Additional file 1) were obtained from broadPeak data deposited as GSE30263 at NCBI GEO. We generated two sets of CTCF-binding sites. The number of overlapping CTCF-binding sites, which are CTCF-binding regions present in all examined cell types (Table S5 in Additional file 1), was determined using BEDTools [39]. The number of joint CTCF-binding sites, the union of CTCF-binding sites from all cell types, was determined using a custom Perl script. Enrichment analyses on GO terms were performed using FatiGO [40]. Partial correlation analyses were performed using modules of the "ppcor" package (v.1.0) [41] for R (http:// www.r-project.org/).