BACKGROUND
gene duplication is a major driver for the emergence of organismal complexity and evolutionary innovations  <cit> .
tandem duplication is the most common route to the formation of clustered paralogous genes  <cit> .
a newly duplicated gene must diverge from its progenitor gene in coding sequence or expression, or it will degenerate due to redundancy  <cit> .
to acquire novel transcription patterns, tandemly duplicated genes need to interrupt expression similarity due to shared upstream cis-elements upon origin  <cit>  and transcriptional interference contributed by physical proximity  <cit> .
despite the challenges, functionally important gene clusters consisting of paralogs with distinct expression patterns are found in a wide range of species, including human  <cit> .
therefore, it is important to understand the origin and maintenance of expression divergence between tandem paralogs.
ccctc-binding factor , the only known human insulator protein, plays a master role in determining the transcriptional landscape of genomes.
when bound to insulator sequences , ctcf can prevent repressive heterochromatin from spreading into neighboring regions  <cit> .
in addition, ctcf interferes with enhancer-promoter communication  <cit>  and guides long-range chromatin interactions  <cit> .
because changes in epigenetic marks play roles in regulatory divergence of duplicated genes , we hypothesize that ctcf also plays a role.
using human rna-seq and chip-seq data , we examined how ctcf drives regulatory divergence of duplicated genes, especially tandemly arrayed paralogs.
RESULTS
in the genome, adjacent genes are coexpressed due to a common origin  <cit> , cofunctionality  <cit> , or deleterious transcriptional interference  <cit> .
using rna-seq data from six human adult tissues  <cit> , we measured gene expression dissimilarity between adjacent genes in terms of expression profile or expression level using expd1-r or expdeuc , respectively.
expd1-r focused on changes in the shape across the tissue dimension, while expdeuc 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  <cit> , possibly due to different underlying mechanisms controlling the properties  <cit> .
mammalian genes close to each other have similar expression profiles  <cit> .
ctcf-binding sites can prevent undesirable crosstalk between active and inactive genomic regions  <cit> .
dna methylation upstream of a gene can inhibit its transcriptional initiation  <cit> .
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 expd1-r and expdeuc in non-paralogous adjacent genes.
non-paralogous gene pairs with a longer intergenic distance , more overlapping ctcf-binding sites , or a larger difference in upstream dna methylation  had greater expd1-r or expdeuc .
although d, #ctcf, and Δcpgo/e were interrelated , partial correlation analyses  suggested that the effects of d, #ctcf, and Δcpgo/e on generating expression dissimilarity between adjacent genes was independent.
d had the largest direct influence on expd1-r, and Δcpgo/e had the largest direct influence on expdeuc.
#ctcf had a weak but significant effect on expd1-r and an intermediate effect on expdeuc .
although ctcf binding can vary among cell types  <cit> , when we define #ctcf using joint ctcf chip-seq peaks  instead of overlapping ctcf peaks, the results did not change .table  <dig> rank correlations  and partial rank correlations  of the genomic properties with expression dissimilarities  of the non-paralogous adjacent genes
genomic properties a
expd
1-r
expd
euc
ρ
b
ρ
p
b,c
ρ
b
ρ
p
b,c
d
#ctcf
Δcpg
o/e
a "d", intergenic distance; "#ctcf", number of overlapping ctcf-binding sites; "Δcpg
o/e", difference in upstream dna methylation.
b
p values show the probabilities of the observations under the hypothesis of no correlation.
c spearman's partial correlation coefficient ρp is computed by controlling for the other two genomic properties listed in a
mechanisms generating coding sequence divergence, such as change in protein structure or splicing, between paralogous genes have been intensively investigated .
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 , for adjacent paralogs, #ctcf had the greatest direct influence on expd1-r and expdeuc, followed by Δcpgo/e .
repeating analysis based on an independently generated and unpublished rna-seq dataset including  <dig> human tissues  produced a result consistent with table  <dig> and  <dig> , suggesting the robustness of the pattern found.table  <dig> rank correlations  and partial rank correlations  of the examined genomic properties with expression dissimilarities  of the paralogous adjacent genes
genomic properties a
expd
1-r
expd
euc
ρ
b
ρ
p
b,c
ρ
b
ρ
p
b,c
d
#ctcf
Δcpg
o/e
a "d", intergenic distance; "#ctcf", number of overlapping ctcf-binding sites; "Δcpg
o/e", difference in upstream dna methylation.
b
p values show the probabilities of the observations under the hypothesis of no correlation.
c spearman's partial correlation coefficient ρp is computed by controlling for the other two genomic properties listed in a
we classified adjacent paralogs into three groups according the orientations: head-to-head, head-to-tail and tail-to-tail  <cit> .
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 .
using microarray expression data, a previous study found that intergenic distance was related to expression divergence for tandemly arrayed paralogs  <cit> .
however, microarray data is known to have cross-hybridization related biases  <cit> , and we found no direct association between intergenic distance and expression divergence for tandemly arrayed paralogs after controlling for #ctcf and Δcpgo/e .
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.table  <dig> rank correlations  and partial rank correlations  of the examined genomic properties with expression dissimilarities  of the adjacent paralogous genes of different orientations.
genomic properties a
expd
1-r
expd
euc
ρ
b
ρ
p
b,c
ρ
b
ρ
p
b,c
head-to-head
d
d
#ctcf
Δcpg
o/e
head-to-tail
d
d
#ctcf
Δcpg
o/e
tail-to-tail
d
d
#ctcf
Δcpg
o/e
a "d", intergenic distance; "#ctcf", number of overlapping ctcf-binding sites; "Δcpg
o/e", difference in upstream dna methylation.
b
p values show the probabilities of the observations under the hypothesis of no correlation.
c spearman's partial correlation coefficient ρp is computed by controlling for the other two genomic properties listed in a
d there were  <dig>   <dig> and  <dig> head-to-head, head-to-tail and tail-to-tail adjacent paralogs, respectively.
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 hypothesis is correct, there would be no correlation between the divergence time of paralogs  and #ctcf.
however, we observed a positive correlation between ds or tphy and #ctcf , suggesting that the second hypothesis is correct.
#ctcf had stronger rank correlation , which is positive, to ds  than d and Δcpgo/e .
partial correlation analyses further suggested that the increase of #ctcf with ds  is not caused by the change in d or Δcpgo/e or their combined effect over time .
by contrast, the insignificant partial correlation of d with respect to ds  indicated that the increase d over time can be explained by the increase in the number of ctcf-binding sites and the associated changes in dna methylation.table  <dig> rank correlations  and partial rank correlations  of the examined genomic properties with the divergence time of the paralogous adjacent genes.
genomic properties a
d
s
t
phy
ρ
b
ρ
p
b,c
ρ
b
ρ
p
b,c
all adjacent paralogs
d
#ctcf
Δcpg
o/e
adjacent paralogs associated with go terms in which high
#ctcf/d
genes were specifically enriched
d
#ctcf
Δcpg
o/e
a "d", intergenic distance; "#ctcf", number of overlapping ctcf-binding sites; "Δcpg
o/e", difference in upstream dna methylation.
b
p values show the probabilities of the observations under the hypothesis of no correlation.
c spearman's partial correlation coefficient ρp is computed by controlling for the other two genomic properties listed in a
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  terms.
to eliminate the potential effect of duplicability in go analysis  <cit> , paralogous and non-paralogous gene pairs were analyzed separately.
to control for the potential effect of gene density  <cit> , #ctcf/d was used instead of #ctcf .
gene pairs in the top quartile of ctcf-binding site density  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 .
only two go terms  had the same enrichment status  in both groups.
tandemly duplicated genes with a higher density of intervening ctcf-binding sites tended to specifically encode proteins involved in gene expression , metabolic processes , or cellular processes  through dna binding , smad binding , growth factor binding , or kinase interaction  in receptor complexes  or intracellular regions  .
this result implied that the densities of ctcf binding sites between tandem paralogs were not contributed by the genomic background.table  <dig> enriched/depleted go terms at level  <dig> for genes with high #ctcf/d flanking regions
go terms
a
go terms
molecular functions
dna binding
transcription factor binding
smad binding
growth factor binding
kinase inhibitor activity
peptidase inhibitor activity
protein kinase regulator activity
cellular components
intracellular
biological processes
regulation of metabolic process
regulation of catalytic activity
regulation of cellular process
nucleobase, nucleoside, nucleotide and nucleic acid metabolic process
gene expression
cellular biosynthetic process
programmed cell death
a specifically enriched go terms are highlighted in bold fonts
b the symbol ↑or ↓ represents the enrichment status of "enriched" or "depleted" of the go term, respectively.
c
p values were bonferroni-corrected for multiple tests
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  <dig> 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. ds  after controlling for d and Δcpgo/e was stronger  than that observed in the full set of adjacent paralogs  .
when #ctcf was defined by joint ctcf chip-seq peaks rather than overlapping ctcf peaks, the results were similar .
therefore, the trend to accumulate ctcf-binding sites following a tandem duplication event resulted in more ctcf-binding sites  between paralogs, especially for those with enriched go-categories shown in table  <dig>
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 expd1-r  were still significant even after controlling for Δcpgo/e .
thus, ctcf, a conserved regulatory protein  <cit> , affects the expression evolution of adjacent genes in genomes from flies  <cit>  to humans and is important for the evolution of organismal complexity in animals.

 gene duplication is a major driver for the emergence of organismal complexity and evolutionary innovations. tandem duplication is the most common route to the formation of clustered paralogous genes. a newly duplicated gene must diverge from its progenitor gene in coding sequence or expression, or it will degenerate due to redundancy. to acquire novel transcription patterns, tandemly duplicated genes need to interrupt expression similarity due to shared upstream cis-elements upon origin and transcriptional interference contributed by physical proximity. despite the challenges, functionally important gene clusters consisting of paralogs with distinct expression patterns are found in a wide range of species, including human. therefore, it is important to understand the origin and maintenance of expression divergence between tandem paralogs. ccctc-binding factor, the only known human insulator protein, plays a master role in determining the transcriptional landscape of genomes. when bound to insulator sequences, ctcf can prevent repressive heterochromatin from spreading into neighboring regions. in addition, ctcf interferes with enhancer-promoter communication and guides long-range chromatin interactions. because changes in epigenetic marks play roles in regulatory divergence of duplicated genes, we hypothesize that ctcf also plays a role. using human rna-seq and chip-seq data, we examined how ctcf drives regulatory divergence of duplicated genes, especially tandemly arrayed paralogs.  in the genome, adjacent genes are coexpressed due to a common origin, cofunctionality, or deleterious transcriptional interference. using rna-seq data from six human adult tissues, we measured gene expression dissimilarity between adjacent genes in terms of expression profile or expression level using expd1-r or expdeuc, respectively. expd1-r focused on changes in the shape across the tissue dimension, while expdeuc 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, possibly due to different underlying mechanisms controlling the properties. mammalian genes close to each other have similar expression profiles. ctcf-binding sites can prevent undesirable crosstalk between active and inactive genomic regions. dna methylation upstream of a gene can inhibit its transcriptional initiation. 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 expd1-r and expdeuc in non-paralogous adjacent genes. non-paralogous gene pairs with a longer intergenic distance, more overlapping ctcf-binding sites, or a larger difference in upstream dna methylation had greater expd1-r or expdeuc. although d, #ctcf, and Δcpgo/e were interrelated, partial correlation analyses suggested that the effects of d, #ctcf, and Δcpgo/e on generating expression dissimilarity between adjacent genes was independent. d had the largest direct influence on expd1-r, and Δcpgo/e had the largest direct influence on expdeuc. #ctcf had a weak but significant effect on expd1-r and an intermediate effect on expdeuc. although ctcf binding can vary among cell types, when we define #ctcf using joint ctcf chip-seq peaks instead of overlapping ctcf peaks, the results did not change.table rank correlations and partial rank correlations of the genomic properties with expression dissimilarities of the non-paralogous adjacent genes genomic properties a expd 1-r expd euc ρ b ρ p b,c ρ b ρ p b,c d #ctcf Δcpg o/e a "d", intergenic distance; "#ctcf", number of overlapping ctcf-binding sites; "Δcpg o/e", difference in upstream dna methylation. b p values show the probabilities of the observations under the hypothesis of no correlation. c spearman's partial correlation coefficient ρp is computed by controlling for the other two genomic properties listed in a mechanisms generating coding sequence divergence, such as change in protein structure or splicing, between paralogous genes have been intensively investigated. 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, for adjacent paralogs, #ctcf had the greatest direct influence on expd1-r and expdeuc, followed by Δcpgo/e. repeating analysis based on an independently generated and unpublished rna-seq dataset including human tissues produced a result consistent with table and, suggesting the robustness of the pattern found.table rank correlations and partial rank correlations of the examined genomic properties with expression dissimilarities of the paralogous adjacent genes genomic properties a expd 1-r expd euc ρ b ρ p b,c ρ b ρ p b,c d #ctcf Δcpg o/e a "d", intergenic distance; "#ctcf", number of overlapping ctcf-binding sites; "Δcpg o/e", difference in upstream dna methylation. b p values show the probabilities of the observations under the hypothesis of no correlation. c spearman's partial correlation coefficient ρp is computed by controlling for the other two genomic properties listed in a we classified adjacent paralogs into three groups according the orientations: head-to-head, head-to-tail and tail-to-tail. 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. using microarray expression data, a previous study found that intergenic distance was related to expression divergence for tandemly arrayed paralogs. however, microarray data is known to have cross-hybridization related biases, and we found no direct association between intergenic distance and expression divergence for tandemly arrayed paralogs after controlling for #ctcf and Δcpgo/e. 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.table rank correlations and partial rank correlations of the examined genomic properties with expression dissimilarities of the adjacent paralogous genes of different orientations. genomic properties a expd 1-r expd euc ρ b ρ p b,c ρ b ρ p b,c head-to-head d     d     #ctcf    Δcpg o/e head-to-tail d     d     #ctcf    Δcpg o/e tail-to-tail d     d     #ctcf    Δcpg o/e a "d", intergenic distance; "#ctcf", number of overlapping ctcf-binding sites; "Δcpg o/e", difference in upstream dna methylation. b p values show the probabilities of the observations under the hypothesis of no correlation. c spearman's partial correlation coefficient ρp is computed by controlling for the other two genomic properties listed in a d there were and head-to-head, head-to-tail and tail-to-tail adjacent paralogs, respectively. 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 hypothesis is correct, there would be no correlation between the divergence time of paralogs and #ctcf. however, we observed a positive correlation between ds or tphy and #ctcf, suggesting that the second hypothesis is correct. #ctcf had stronger rank correlation, which is positive, to ds than d and Δcpgo/e. partial correlation analyses further suggested that the increase of #ctcf with ds is not caused by the change in d or Δcpgo/e or their combined effect over time. by contrast, the insignificant partial correlation of d with respect to ds indicated that the increase d over time can be explained by the increase in the number of ctcf-binding sites and the associated changes in dna methylation.table rank correlations and partial rank correlations of the examined genomic properties with the divergence time of the paralogous adjacent genes. genomic properties a d s t phy ρ b ρ p b,c ρ b ρ p b,c all adjacent paralogs     d     #ctcf    Δcpg o/e adjacent paralogs associated with go terms in which high #ctcf/d genes were specifically enriched     d     #ctcf    Δcpg o/e a "d", intergenic distance; "#ctcf", number of overlapping ctcf-binding sites; "Δcpg o/e", difference in upstream dna methylation. b p values show the probabilities of the observations under the hypothesis of no correlation. c spearman's partial correlation coefficient ρp is computed by controlling for the other two genomic properties listed in a 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 terms. to eliminate the potential effect of duplicability in go analysis, paralogous and non-paralogous gene pairs were analyzed separately. to control for the potential effect of gene density, #ctcf/d was used instead of #ctcf. gene pairs in the top quartile of ctcf-binding site density 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. only two go terms had the same enrichment status in both groups. tandemly duplicated genes with a higher density of intervening ctcf-binding sites tended to specifically encode proteins involved in gene expression, metabolic processes, or cellular processes through dna binding, smad binding, growth factor binding, or kinase interaction in receptor complexes or intracellular regions. this result implied that the densities of ctcf binding sites between tandem paralogs were not contributed by the genomic background.table enriched/depleted go terms at level for genes with high #ctcf/d flanking regions go terms a go terms molecular functions dna binding transcription factor binding smad binding growth factor binding kinase inhibitor activity peptidase inhibitor activity protein kinase regulator activity cellular components intracellular biological processes regulation of metabolic process regulation of catalytic activity regulation of cellular process nucleobase, nucleoside, nucleotide and nucleic acid metabolic process gene expression cellular biosynthetic process programmed cell death a specifically enriched go terms are highlighted in bold fonts b the symbol ↑or ↓ represents the enrichment status of "enriched" or "depleted" of the go term, respectively. c p values were bonferroni-corrected for multiple tests 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 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. ds after controlling for d and Δcpgo/e was stronger than that observed in the full set of adjacent paralogs. when #ctcf was defined by joint ctcf chip-seq peaks rather than overlapping ctcf peaks, the results were similar. therefore, the trend to accumulate ctcf-binding sites following a tandem duplication event resulted in more ctcf-binding sites between paralogs, especially for those with enriched go-categories shown in table  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 expd1-r were still significant even after controlling for Δcpgo/e. thus, ctcf, a conserved regulatory protein, affects the expression evolution of adjacent genes in genomes from flies to humans and is important for the evolution of organismal complexity in animals. 
however, how expression divergence between tandemly duplicated genes has emerged and maintained remain unclear. the accumulation of ctcf-binding sites drives expression divergence of tandemly duplicated genes. after tandem duplication events, ctcf-binding sites gradually accumulate between paralogs. this phenomenon was not observed for non-paralogous adjacent genes. this process is likely targeted by natural selection. we demonstrate that ccctc-binding factor, 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.
 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.  we demonstrate that ccctc-binding factor, 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.  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 contains supplementary material, which is available to authorized users. keywords epigeneticsdna methylationgenome architecturechip-seqrna-seqissue-copyright-statement© biomed central ltd 2014 