insulator proteins confer activity to insulators sequences , a class of functional elements with important roles in the regulation of chromosomal organization . in this study , we investigated the origin of the seven known proteins associated with insulator activity in d. melanogaster . to find potential orthologs of these proteins in other organisms , we followed two independent strategies . on the one hand , we performed blast searches and analyzed a large number of candidate sequences in phylogenetic experiments ( table 1 ; figs .
s1s4 ) . on the other hand , we combined hmm searches and mcl clustering to examine which of the candidates share orthology with a known drosophila insulator ( table s4 ) . with both methods we found robust evidence that the known drosophila insulators ( except ctcf ) are restricted to arthropods and have been acquired successively during arthropod evolution ( fig . 2 ;
, we did not find evidence for the existence of known drosophila insulator proteins outside the arthropods by analyzing seven genomes from four protostome phyla ( annelida , mollusca , platyhelminthes , nematoda ) .
it is unlikely that a lack of sensitivity is responsible for this result as we find ctcf orthologs in the annelid capitella spi , in the mollusc lottia gigantea , and in the nematode trichinella spiralis , faithfully replicating the results of previous work with an independent method ( heger et al .
the detection of ctcf in all 19 analyzed arthropod genomes , including so far unknown orthologs ( e.g. , from the myriapod strigamia maritima and the strepsipteran m. moldrzyki ) , further confirms the specificity and sensitivity of our approach .
second , the zf proteins ctcf and su(hw ) are consistently present in arthropods from which genome sequences are available plus in some additional orders with a significant number of ests . in both cases , the presence of orthologous sequences in the three subgroups of arthropods indicates that the common ancestor of arthropods already had these proteins .
such an assumption is well supported for the ctcf protein that has been found in all bilaterians ( heger et al .
the origin of the su(hw ) protein is less clear . as it is absent in nematodes ( table s4 ; heger et al .
2009 ) , it could have evolved in the ancestor of arthropods or in the common ancestor of arthropods and a closely related ecdysozoan sister group , for example , tardigrades or onychophorans .
third , in contrast to ctcf and su(hw ) , the proteins zw5 and beaf-32 are restricted to a remarkably limited subset of arthropods .
we obtained with both methods only few beaf-32 candidates outside the drosophila genus ( 64 and 95 , respectively ) , and none of them could be placed into the drosophila beaf-32 orthology group .
although the domain composition of zw5 issued a much higher number of candidates ( 28,431 ) across the 26 genomes ( table s4 ) , we could only recover zw5 orthologs in drosophila and g. morsitans , another brachyceran fly .
a zw5 ortholog could not be found in nematoceran dipteres and other arthropods despite the existence of several sequenced genomes and large amounts of ests .
these results are based on the most comprehensive study undertaken so far and provide consistent evidence that zw5 and beaf-32 likely emerged in or close to the common ancestor of the genus drosophila .
finally , our results with respect to the btb domain containing insulators suggest that at least some of them have evolved at intermediate stages when compared with the  ubiquitous  proteins ctcf / su(hw ) and the  restricted 
the mapping of cp190 , for example , shows that it is present in all insects with a sequenced genome and in three crustacean orders .
the inability to detect this protein in a myriapod and two mite genomes and in a large amount of ests from chelicerates suggests that cp190 evolved in the ancestor of pancrustacea .
gaga factor , on the other hand , formed a highly supported cluster containing exclusively sequences from diptera and hymenoptera in phylogenetic experiments ( fig . s2 ) .
the hmm - based approach confirmed this result , but assigned in addition sequences from the hemipteran rhodnius prolixus to the gaga factor orthology group ( table s4 ) . in all
further insect , crustacean , and chelicerate genomes , orthologs were not detectable , indicating that this protein likely emerged in the common ancestor of hemiptera , hymenoptera , and diptera .
importantly , these findings imply that gaga factors do not exist outside the arthropod phylum . to elucidate the conflict of these findings with the proposed existence of vertebrate gaga factors ( matharu et al . 2010 ; kumar 2011 ) , we performed phylogenetic analysis with two sets of candidates , acquired by blast and hmm searches .
however , neither the four proposed vertebrate gaga factors nor our additional candidates were placed to the insect orthology group ( figs .
although our study provides evidence for a consecutive gain of insulator proteins in arthropods , it also suggests that insulator evolution is dynamic in terms of losses .
although it is difficult to prove the absence of a gene in potentially inaccurate genome assemblies , our results indicate that inference of at least some secondary losses is reasonable .
the two best examples are the repetitive loss of gaga factor in coleoptera / strepsiptera and in lepidoptera and the loss of su(hw ) in lepidoptera , each supported by the analysis of two genomes with both methods ( fig .
we were not able to discover an expected ortholog in a single genome ( table s4 ) .
this may be a consequence of incomplete genome assemblies , but could alternatively reflect a dynamic nature of insulator evolution in arthropods .
we did not observe a loss of the more ancient ctcf insulator in any arthropod genome .
this is compatible with the idea that ctcf function is needed for fundamental processes in the bilateria ( heger et al .
2012 ) and might have been supplemented and modified by additional components in the arthropod lineage .
although our results have been established by two fairly independent methods , we can not formally prove the absence of orthologous proteins from certain clades . with growing databases and dependent on the sensitivity of homology detection tools , the exact placement of the origin of some proteins may still change .
however , irrespective of uncertainties in the exact time of gain and loss , our observations reveal a consistent model for the evolution of the known d. melanogaster chromatin insulators through a series of successive acquisitions .
it has been confirmed repeatedly that insulator proteins colocalize and interact with each other ( gerasimova et al .
2007 ; bartkuhn et al . 2009 ; negre et al . 2010 ; van bortle et al .
2012 ) , thereby creating a network of dependencies that is thought to contribute to cell - specific differences in nuclear organization and gene expression ( yang and corces 2011 ) .
moreover , there is recent evidence that each d. melanogaster insulator subclass , determined by the presence of su(hw ) , ctcf , beaf , or gaf , shares the cp190 protein and possibly also mod(mdg4 ) ( van bortle et al .
our findings point out that the mechanisms , interactions , and components of insulator complexes must be considerably different from drosophila not only in the great majority of arthropods ( that share two of the seven factors ) , but also in other bilaterian phyla that only have ctcf in common . according to our findings ,
the presence of different insulator systems and the complex interactions between their components seen in d. melanogaster today are a result of ongoing evolution and diversification .
these processes started more than  600 million years ago in the ancestor of bilaterians ( www.timetree.org ) with ctcf , the oldest known chromatin insulator of multicellular animals ( heger et al . 2012 ) . at or close to the root of arthropods ,
these two basal systems experienced subsequent additions and modifications , with the gain of beaf-32 in the genus drosophila ,  60 million years ago , being the most recent acquisition . to what extent the successive acquisition of new insulator genes is involved in adaptive processes , shall be the topic of future investigations . although we can describe in detail an expansion of insulator proteins only in the drosophila history ,
millions of arthropod species trace back to the same common ancestor that was equipped with the ctcf and su(hw ) insulators .
it is therefore possible that other arthropods simultaneously increased their initial repertoire of insulator proteins during the past 600 million years , giving rise to a plethora of unexplored territories .
thus , the expansion of insulator mechanisms that happened in d. melanogaster history might in fact be a general characteristic of arthropods and other animals .
the description of non - ctcf insulators in echinoderms exemplifies that such an expansion could indeed apply to a greater variety of animals ( yajima et al .
if so , this characteristic could provide a mechanism to modulate an ancient insulator system and fine - tune gene expression in a lineage - specific way .
with the known drosophila insulator proteins as query ( dctcf , gi:21356747 ; su(hw ) , gi:33860216 ; cp190 , gi:23171337 ; gaf , gi:83287912 ; mod(mdg4 ) , gi:158030328 ; zw5 , gi:45549097 ; beaf , gi:17647187 ) , standard blastx , blastp , or tblastn ( altschul et al . 1997 ) searches were conducted in publicly available sequence databases at ncbi ( http://blast.ncbi.nlm.nih.gov/ ) . to minimize the chance of missing an ortholog , we performed parallel searches in different databases ( nucleotide , est , and protein ) and subdivided a search into smaller entities if a given taxonomic range reported > 500 hits below threshold (= download restriction at the ncbi blast web interface ) .
we set the threshold for btb domain protein candidates to 10 because in preliminary experiments this value effectively incorporated btb domain proteins distinct from btb domain containing insulators , for example , tramtrack or broad complex .
collected nucleotide sequences were translated to the appropriate reading frame using emboss ( rice et al .
we obtained from the initial dataset a collection of unique sequences using the emboss tool 
all other candidates from the initial dataset were passed to the silix clustering algorithm ( miele et al .
2011 ) as  partial  sequences ( command line parameter  -p  ) .
, we took all reference sequences plus candidates that clustered to one of the references .
for the zw5 , su(hw ) , and ctcf proteins , the described drosophila orthologs served as positive control for cluster identification .
similarly , we included the known drosophila orthologs of the btb - domain proteins gaga factor , cp190 , and mod(mdg4 ) to specify these clusters .
outgroup for the analysis of zf insulators was the widely distributed zf transcription factor yy1 .
for btb domain insulators , we used as outgroup a set of lola - like orthologs from diverse arthropods .
phylogenetic trees resulting from the alignments were computed under the maximum likelihood criterion using parallel raxml version 7.2.6 ( stamatakis 2006 ) with 100 bootstrap resamplings .
as optimal models of sequence evolution we used the wag+ or dcmut++f model ( zf domain proteins ) , the jtt++f model ( btb domain proteins ) , and the wag++f model ( vertebrate gafs ) as selected by prottest3 ( darriba et al .
likelihood trees were visualized and arranged with figtree ( http://tree.bio.ed.ac.uk/software/figtree/ ) and then graphically edited with adobe illustrator software . to verify our results with an independent method , we combined a search on the basis of insulator - specific hidden markov models ( http://hmmer.org/ ) with the markov cluster algorithm ( van dongen 2000 ) .
we downloaded from various sites ( table s3 ) genome and , if available , proteome and unplaced reads data of 19 different arthropod and seven outgroup species .
the 26 selected species represent maximal diversity while avoiding over - representation of well - sampled groups such as dipterans . as subsequent comparisons relied on protein sequences
, we translated the 26 genomes and the respective unplaced reads data into all six reading frames ( > 90 nucleotides ) and combined the resulting open reading frames of each species with the corresponding protein set as offered by the sequencing center . to obtain specific hmmer profiles
, we selected  with attention to maximal diversity and length  for each insulator protein 815 previously verified orthologs as representatives of a cluster .
we then calculated and manually refined multiple alignments of these sequences using the mafft  einsi 
algorithm ( katoh et al . 2005 ) and derived a representative full length hmmer profile for each insulator protein .
all 26 orf sets were scanned with the six custom made hmmer profiles and all sequences below the default inclusion threshold of hmmsearch ( e - value < 0.01 ; 39,573 unique sequences ) were fed to a dedicated orthomcl pipeline ( li et al .
we used the recommended inflation parameter  1.5  for the mcl step . before clustering , we removed duplicates and supplemented each of the 26 sequence collections with reference insulator sequences of that species ( verified by phylogenetic experiments ) to facilitate cluster recognition .
zw5 , and beaf , we detected a single orthology group whereas su(hw ) ( 2 ) , cp190 ( 4 ) , and gaf ( 3 ) orthologs split to more than one group .
like described above , we constructed two multiple sequence alignments from representative members of the gaga cluster ( full length and btb domain only ) and uploaded these alignments to hhpred ( http://toolkit.tuebingen.mpg.de/hhpred ) for highly sensitive profile  profile searches .
we searched with default parameters in the proteomes of h. sapiens and mus musculus , the only deuterostome datasets available . after removing duplicates , we analyzed the remaining 45 candidates with phylogenetic methods ( see above ) .
with the known drosophila insulator proteins as query ( dctcf , gi:21356747 ; su(hw ) , gi:33860216 ; cp190 , gi:23171337 ; gaf , gi:83287912 ; mod(mdg4 ) , gi:158030328 ; zw5 , gi:45549097 ; beaf , gi:17647187 ) , standard blastx , blastp , or tblastn ( altschul et al . 1997 ) searches were conducted in publicly available sequence databases at ncbi ( http://blast.ncbi.nlm.nih.gov/ ) . to minimize the chance of missing an ortholog , we performed parallel searches in different databases ( nucleotide , est , and protein ) and subdivided a search into smaller entities if a given taxonomic range reported > 500 hits below threshold (= download restriction at the ncbi blast web interface ) .
we set the threshold for btb domain protein candidates to 10 because in preliminary experiments this value effectively incorporated btb domain proteins distinct from btb domain containing insulators , for example , tramtrack or broad complex .
collected nucleotide sequences were translated to the appropriate reading frame using emboss ( rice et al .
we obtained from the initial dataset a collection of unique sequences using the emboss tool 
all other candidates from the initial dataset were passed to the silix clustering algorithm ( miele et al .
2011 ) as  partial  sequences ( command line parameter  -p  ) .
, we took all reference sequences plus candidates that clustered to one of the references .
for the zw5 , su(hw ) , and ctcf proteins , the described drosophila orthologs served as positive control for cluster identification .
similarly , we included the known drosophila orthologs of the btb - domain proteins gaga factor , cp190 , and mod(mdg4 ) to specify these clusters .
outgroup for the analysis of zf insulators was the widely distributed zf transcription factor yy1 . for btb domain insulators , we used as outgroup a set of lola - like orthologs from diverse arthropods .
phylogenetic trees resulting from the alignments were computed under the maximum likelihood criterion using parallel raxml version 7.2.6 ( stamatakis 2006 ) with 100 bootstrap resamplings . as optimal models of sequence evolution we used the wag+ or dcmut++f model ( zf domain proteins ) , the jtt++f model ( btb domain proteins ) , and the wag++f model ( vertebrate gafs ) as selected by prottest3 ( darriba et al .
likelihood trees were visualized and arranged with figtree ( http://tree.bio.ed.ac.uk/software/figtree/ ) and then graphically edited with adobe illustrator software .
to verify our results with an independent method , we combined a search on the basis of insulator - specific hidden markov models ( http://hmmer.org/ ) with the markov cluster algorithm ( van dongen 2000 ) .
we downloaded from various sites ( table s3 ) genome and , if available , proteome and unplaced reads data of 19 different arthropod and seven outgroup species .
the 26 selected species represent maximal diversity while avoiding over - representation of well - sampled groups such as dipterans . as subsequent comparisons relied on protein sequences , we translated the 26 genomes and the respective unplaced reads data into all six reading frames ( > 90 nucleotides ) and combined the resulting open reading frames of each species with the corresponding protein set as offered by the sequencing center . to obtain specific hmmer profiles , we selected  with attention to maximal diversity and length  for each insulator protein 815 previously verified orthologs as representatives of a cluster .
all 26 orf sets were scanned with the six custom made hmmer profiles and all sequences below the default inclusion threshold of hmmsearch ( e - value < 0.01 ; 39,573 unique sequences ) were fed to a dedicated orthomcl pipeline ( li et al .
we used the recommended inflation parameter  1.5  for the mcl step . before clustering , we removed duplicates and supplemented each of the 26 sequence collections with reference insulator sequences of that species ( verified by phylogenetic experiments ) to facilitate cluster recognition .
zw5 , and beaf , we detected a single orthology group whereas su(hw ) ( 2 ) , cp190 ( 4 ) , and gaf ( 3 ) orthologs split to more than one group .
like described above , we constructed two multiple sequence alignments from representative members of the gaga cluster ( full length and btb domain only ) and uploaded these alignments to hhpred ( http://toolkit.tuebingen.mpg.de/hhpred ) for highly sensitive profile  profile searches .
we searched with default parameters in the proteomes of h. sapiens and mus musculus , the only deuterostome datasets available . after removing duplicates , we analyzed the remaining 45 candidates with phylogenetic methods ( see above ) .
additional supporting information may be found in the online version of this article at the publisher s website : figure s1 .
download location of 26 genomes and proteomes used for a hmm / orthomcl - based search for insulator proteins .
 alteration of regulatory dna elements or their binding proteins may have drastic consequences for morphological evolution . 
 chromatin insulators are one example of such proteins and play a fundamental role in organizing gene expression . 
 while a single insulator protein , ctcf ( ccctc - binding factor ) , is known in vertebrates , drosophila melanogaster utilizes six additional factors . 
 we studied the evolution of these proteins and show here that  
 in contrast to the bilaterian - wide distribution of ctcf  all other d. melanogaster insulators are restricted to arthropods . 
 the full set is present exclusively in the genus drosophila whereas only two insulators , su(hw ) and ctcf , existed at the base of the arthropod clade and all additional factors have been acquired successively at later stages . 
 secondary loss of factors in some lineages further led to the presence of different insulator subsets in arthropods . 
 thus , the evolution of insulator proteins within arthropods is an ongoing and dynamic process that reshapes and supplements the ancient ctcf - based system common to bilaterians . 
 expansion of insulator systems may therefore be a general strategy to increase an organism s gene regulatory repertoire and its potential for morphological plasticity . 
Discussion
Methods Summary
Blast-Based Search for Candidate Insulator Proteins
Clustering and Multiple Sequence Alignment
Phylogenetic Analysis
HMM-Based Candidate Search and MCL-Clustering of Orthologous Groups
HHpred-Based Search for GAGA Factor Orthologs
Supporting Information