extfieldsofstudy: Biology Medicine

s2fieldsofstudy: Biology Medicine

text:
Natural HLA-B*2705 Protein Ligands with Glutamine as Anchor Motif

Background: HLA-B27 is strongly associated with ankylosing spondylitis (AS). Results: A change in the position (P) 2 anchor motif was detected in the 3% of HLA-B27 ligands identified, which showed significant homology to pathogenic bacterial sequences. Conclusion: Unusual HLA-B*2705 ligands bind with different conformations to both AS-associated and non-AS-associated HLA-B27 subtypes. Significance: This could be sufficient to initiate autoimmune damage in patients with AS-associated subtypes. The presentation of short viral peptide antigens by human leukocyte antigen (HLA) class I molecules on cell surfaces is a key step in the activation of cytotoxic T lymphocytes, which mediate the killing of pathogen-infected cells or initiate autoimmune tissue damage. HLA-B27 is a well known class I molecule that is used to study both facets of the cellular immune response. Using mass spectrometry analysis of complex HLA-bound peptide pools isolated from large amounts of HLA-B*2705+ cells, we identified 200 naturally processed HLA-B*2705 ligands. Our analyses revealed that a change in the position (P) 2 anchor motif was detected in the 3% of HLA-B*2705 ligands identified. B*2705 class I molecules were able to bind these six GlnP2 peptides, which showed significant homology to pathogenic bacterial sequences, with a broad range of affinities. One of these ligands was able to bind with distinct conformations to HLA-B27 subtypes differentially associated with ankylosing spondylitis. These conformational differences could be sufficient to initiate autoimmune damage in patients with ankylosing spondylitis-associated subtypes. Therefore, these kinds of peptides (short, with GlnP2, and similar low affinity to all HLA-B27 subtypes tested but with unlike conformations in differentially ankylosing spondylitis-associated subtypes) must not be excluded from future researches involving potential arthritogenic peptides.

The presentation of short viral peptide antigens by human leukocyte antigen (HLA) class I molecules on cell surfaces is a key step in the activation of cytotoxic T lymphocytes, which mediate the killing of pathogen-infected cells or initiate autoimmune tissue damage. HLA-B27 is a well known class I molecule that is used to study both facets of the cellular immune response.
Using mass spectrometry analysis of complex HLA-bound peptide pools isolated from large amounts of HLA-B*2705 ؉ cells, we identified 200 naturally processed HLA-B*2705 ligands. Our analyses revealed that a change in the position (P) 2 anchor motif was detected in the 3% of HLA-B*2705 ligands identified. B*2705 class I molecules were able to bind these six GlnP2 peptides, which showed significant homology to pathogenic bacterial sequences, with a broad range of affinities. One of these ligands was able to bind with distinct conformations to HLA-B27 subtypes differentially associated with ankylosing spondylitis. These conformational differences could be sufficient to initiate autoimmune damage in patients with ankylosing spondylitis-associated subtypes. Therefore, these kinds of peptides (short, with GlnP2, and similar low affinity to all HLA-B27 subtypes tested but with unlike conformations in differentially ankylosing spondylitis-associated subtypes) must not be excluded from future researches involving potential arthritogenic peptides.
Proteolytic degradation of self-proteins and pathogenic proteins in the cytosol by the combined actions of the proteasomes and degradative peptidases generates peptides, mostly of 8 -10 residues, that are translocated to the endoplasmic reticulum lumen by the transporter associated with antigen processing (TAP) 2 molecules. These short peptides assemble with human leukocyte antigen (HLA) class I heavy chain and ␤ 2 -microglobulin molecules (1). Typically, this interaction is made possible by two major anchor residues at position 2 (P2) and the C terminus (C⍀) of the antigenic peptide (2,3) that are deeply inserted into specific pockets of the antigen recognition site of the HLA class I molecule (4,5). Finally, the stable trimolecular peptide-HLA-␤ 2 -microglobulin complexes are transported to the cell membrane and presented for cytotoxic T lymphocyte recognition (6). This recognition of pathogenic or self-peptide ligands can lead to the beneficial killing of pathogen-infected cells or initiate a pathologic autoimmune damage respectively.
One of the most interesting class I HLA alleles is the HLA-B27, which is strongly associated with ankylosing spondylitis (AS), a chronic inflammatory spondyloarthropathy (7). Furthermore, although most of the HLA-B27 subtypes are strongly associated with AS, HLA-B*2706 and -B*2709 are either not associated or perhaps only weakly associated with this disease (reviewed in Ref. 8), suggesting that this polymorphism subtype modulates disease susceptibility. Previous studies show an unambiguous functional distinction between the closely related AS-associated B*2704 and non-AS-associated B*2706 subtypes (9). However, after four decades of research, the basis for this association remains an enigma. Several hypotheses, each based on a particular feature of the HLA-B27 gene, have been proposed to elucidate this intriguing association, but none have yet satisfactorily explained the mechanism and the differential association of HLA-B27 subtypes with disease. The arthritogenic peptide hypothesis (10) assumes that HLA-B27 can present a microbial epitope, eliciting a normal cytotoxic T lymphocyte response against the pathogen. Unfortunately, some of these cytotoxic T lymphocytes would cross-react with an autol-* This work was supported by grants from the Programa Ramó n y Cajal and ogous self-ligand also presented by this molecule and show molecular mimicry with the primary pathogenic epitope. This cross-reaction promotes the autoimmune tissue injury and inflammation. In this hypothesis, the differential binding or immune recognition of either microbial or mimetic self-peptide(s) between the HLA subtypes explains the disparity in association with AS. The interest in this hypothesis has been highly renewed by recent studies on the role of the interaction between HLA-B27 and the endoplasmic reticulum aminopeptidase involved in the trimming of peptides for HLA class I antigen presentation, ERAP1 (11,12).
The major outstanding feature of HLA-B*2705 specificity is its almost mandatory requirement for Arg at P2 (SYFPEITHI database (3)). In contrast, two ligands with GlnP2 were identified (13,14). To broaden this study, by means of an immunoproteomics analysis of peptide pools isolated from HLA-B27 ϩ cells the current study identifies, in addition to several hundred HLA-B27 ArgP2 self-ligands, a small fraction of peptides contained Gln at the P2 anchor motif. These peptides bind to HLA-B27 subtypes with differential association with AS and show significant homology to arthritogenic bacterial sequences.
Synthesis of Peptides-Peptides were synthesized in a peptide synthesizer (model 433A; Applied Biosystems, Foster City, CA) and purified by reverse-phase HPLC. The molecular mass of peptides was established with MALDI-TOF MS, and composition was determined by LC-MS/MS.
Electrospray-Ion Trap Mass Spectrometry Analysis-Peptide mixtures recovered after the ultrafiltration step were concentrated with Micro-Tip reverse-phase columns (C 18 , 200 l, Harvard Apparatus, Holliston, MA). Each C 18 tip was equilibrated with 80% acetonitrile in 0.1% TFA, washed with 0.1% TFA, and then loaded with the peptide mixture. The tip was then washed with an additional volume of 0.1% TFA, and the peptides were eluted with 80% acetonitrile in 0.1% TFA. Peptide samples were then concentrated to ϳ18 l using vacuum centrifugation.
Recovered HLA class I peptides were analyzed in three HPLC procedures by LC-MS/MS using an Orbitrap XL mass spectrometer (Thermo Electron) fitted with a capillary HPLC column (Eksigent, Dublin, CA). The peptides were resolved on a homemade ReproSil C18 capillary column (75-m inner diameter) with a 7-40% acetonitrile gradient for 2 h in the presence of 0.1% formic acid, as in Ref. 22. The seven most intense masses that exhibited single-, double-, and triple-charge states were selected for fragmentation by collision-induced dissociation from each full mass spectrum.
Database Searches-Pep-Miner (23) was used for peak list generation of the LC-MS/MS data. The HLA peptides were identified using the following search engines: Pep-Miner (23); Proteome Discoverer 1.0 SP1 (Thermo Electron), combining the results of Sequest 3.31 and Bioworks Browser 3.3.1 SP1 (Thermo Electron) (24); and Mascot (server 2.2, Matrix Science) (25), using the human and the viral sections of the National Center for Biotechnology Information (NCBI) database, which includes 448,769 proteins. The search was not limited by enzymatic specificity, the peptide tolerance was set to 0.005 Da, and the fragment ion tolerance was set to 0.5 Da. Identified peptides were selected if the following criteria were met: Pep-Miner score above 80; Mascot score above 20; Sequest Xcorr Ͼ1.7 for singly charged, Ͼ2.2 for doubly charged, and Ͼ2.9 for triply charged peptides; P pep less than 1 ϫ 10 Ϫ4 with Bioworks Browser; Proteome Discovered score higher than 20; and mass accuracy of 0.005 Da (22). When the MS/MS spectra fitted more than one peptide, only the highest scoring peptide was analyzed. The false positive rate for peptide identification was set to 2% based on the search of a reversed database. In addition, the corresponding synthetic peptides were made, and their manually identified MS/MS spectra were used to confirm the assigned sequence of HLA-B27 ligands.
MHC/Peptide Stability Assays-The following synthetic peptides were used as controls in complex stability assays: Flu NP (SRYWAIRTR, HLA-B27-restricted) (26) and C4CON (QYDDAVYLK, HLA-Cw4-restricted) (27). RMA-S B*2705 transfectant cells, a cell line deficient in TAP that expresses low levels of cell surface MHC class I, were incubated at 26°C for 16 h in RPMI 1640 medium supplemented with 10% heat-inactivated FBS. This allows the expression of empty MHC class I molecules (without antigenic peptide) at the cellular membrane that are stable at 26°C but not at 37°C. The cells were washed and incubated for 2 h at 26°C with various concentrations of peptide in the same medium. The cells were maintained at 37°C for an additional 4 h and then collected for flow cytometry. This method allows empty MHC class I molecules to become internalized and can thus discriminate between bound and unbound peptides. MHC expression was measured using 100 l of hybridoma culture supernatant containing ME1 (anti-HLA-B27) mAb as described previously (28). Data were acquired on a FACSCalibur flow cytometer (BD Biosciences) and analyzed using CellQuest Pro 2.0 software (BD Biosciences). Cells incubated without peptides exhibited peak fluorescence intensities close to the background staining observed with secondary antibody alone. The fluorescence index was cal-culated for each time point as the ratio of peak channel fluorescence of the sample to that of the control incubated without peptide. Binding of peptides was also expressed as EC 50 , which is the molar concentration of the peptide at 50% of the maximum fluorescence obtained in a concentration range of 100 -0.001 M.
Statistical Analysis-To analyze statistical significance, an unpaired Student's t test was used. p values Ͻ 0.001 were considered to be significant.

HLA-B*2705
Peptidome-HLA-bound peptide pools were isolated from large amounts of B27-C1R cells. This peptide mixture was subsequently separated by reverse-phase HPLC and analyzed by mass spectrometry. Using several software technologies (see "Experimental Procedures"), 198 fragmentation spectra were resolved with high confidence parameters as peptidic sequences of different human cellular proteins (data not shown). As a control, a reverse database search of the same HLA-B27-bound peptide pool showed a 2% false positive rate.
Arg or Gln as the HLA-B27 P2 Anchor Motif-The classical anchor motif for HLA-B*2705 binding, Arg at the P2 residue (SYFPEITHI database (3)), was present in 193 (97%) of the detected ligands, and no differences with the HLA-B27 ligands previously described were found (data not shown). In contrast, this motif was absent in the other six ligands (Table 1). Although virtually all significant fragments of all MS/MS spectra were assigned as daughter ions of the tentative peptidic sequence (supplemental Figs. 1-6, upper panels), the bioinformatics assignment of each sequence was confirmed by identity with the MS/MS spectrum of each corresponding synthetic peptide (supplemental Figs. 1-6, lower panels). These confirmations of synthetic peptides indicate the correct database assignment of the HLA-B27 peptidome.
Two of the ligands with GlnP2 had been previously detected in an earlier study of the HLA-B*2705 peptidome (14) ( Table 1). One peptide was derived from the cartilage-related protein ADAM10, two nested peptides were derived from the monomorphic ␤2-microglobulin of the HLA class I molecule, and the last three were derived from a tumor antigen, a hypothetical protein, and the SF3B2 protein, respectively (Table 1).
In summary, these results indicate that a total of six ligands with absent ArgP2 anchor motifs were endogenously processed and presented in the HLA-B*2705 ϩ cell line. In addition, these data confirm that GlnP2 is an anchor motif derived from peptides bound to the HLA-B*2705 class I molecule.
GlnP2 Ligands Bind to the B*2705 Molecule-Following a similar strategy to the one used in this study, recent studies have identified several hundred HLA-B27 ligands by immunoprecipitation with the W6/32 mAb of HLA-B27-peptide complexes using the B27-C1R cell line (summarized in the SYFPEITHI database (3)). As these six ligands with GlnP2 do not possess the major HLA-B27 anchor motif ArgP2, one possibility is that they could have a defective interaction with the HLA-B27 class I molecule. To test this hypothesis, MHC/ peptide complex stability assays were carried out using TAPdeficient RMA-S cells transfected with HLA-B*2705 and three of the six noncanonical ligands. These included the larger (IQRTPKIQVY) and shorter (RQPQVSI) ligands as well as a nonamer (GQYGNPLNK). The GQYGNPLNK synthetic peptide induced similar numbers of HLA-peptide surface complexes to a well known HLA-B*2705 epitope from the influenza virus, whereas the other two peptides induced fewer HLA-B*2705-peptide complexes (Fig. 1A). In addition, the relative MHC class I affinity was determined for all peptides. The GQYGNPLNK peptide bound to HLA-B*2705 class I molecules with EC 50 values in the range commonly observed among natural ligands (Fig. 1B). In contrast, the HLA affinity was substantially less for the IQRTPKIQVY (EC 50 , 57 Ϯ 14) and RQPQVSI (EC 50 , 200 Ϯ 21) peptides, and thus, these peptides must be considered as medium and low affinity ligands, respectively (Fig. 1B). These data indicate that all ligands detected in the B27-C1R cell line were endogenously presented in association with the B*2705 molecule and that this molecule could bind peptides with a broad range of affinities.
Binding of the RQPQVSI Peptide to HLA-B*2705 Molecules-Peptides without an N-terminal interaction with the HLA molecule have been identified as endogenously bound to HLA-B39 (29) and the murine MHC class I molecule H-2D d (30), indicating that canonical MHC-peptide interactions in the P1 pocket are not always necessary for endogenous peptide presentation. Usually, the HLA-B*2705 molecule binds peptides of 8 -13 residues (3,31), but in this study, a 7-mer was found bound to this HLA class I molecule. One possibility for this is that the short RQPQVSI peptide binds to HLA-B27 molecules in the same manner that other peptides lacking the N-terminal binding residue do. Thus, new HLA-B*2705/peptide complex stability assays were carried out exchanging either the ArgP1 or the GlnP2 residues of the RQPQVSI peptide for alanine. Fig. 2 shows similar binding affinities for both the HLA-B*2705 natural ligands and the single Ala-substituted peptides. The exchange of both ArgP1 and GlnP2 residues by Ala abrogated binding to HLA-B*2705 molecules. Collectively, these data indicate that ArgP1 as well as GlnP2 residues of the RQPQVSI  (13,14).
peptide are sufficient to support interactions with class I molecules.
GlnP2 Ligands Bind to HLA-B27 Subtypes with Differential Association with Ankylosing Spondylitis-In contrast to HLA-B*2705, some HLA-B27 subtypes are not associated with AS (reviewed in Ref. 8). Thus, we selected a pair of HLA-B27 subtypes with differential association with this autoimmune disease. Although HLA-B*2706 is not associated with AS, its most related HLA-B27 subtype, B*2704, is an AS-associated subtype. Thus, a comparison of the specific binding of these HLA class I molecules to peptides with noncanonical GlnP2 residues was carried out. The binding to HLA-B*2706 was significantly less than that of B*2704 (and the B*2705) with both the GQYGNPLNK and the IQRTPKIQVY peptides; thus, these peptides must be considered low affinity ligands in the subtype associated with AS (Fig. 3). No differences in binding were found with the RQPQVSI peptide in both subtypes (Fig. 3). In summary, the three GlnP2 peptides studied bound to HLA-B*2706 with relatively low affinity. In addition, the differential binding between the B*2704 and B*2706 subtypes detected with the GQYGNPLNK and IQRTPKIQVY peptides, but not with RQPQVSI peptide, correlated with the previously reported repertoire of ArgP2 peptide specificity of these subtypes, where B*2704 but not B*2706 favors the binding of Tyr or Arg at residue carboxy-terminal of a peptide (P⍀) (32).
Differential Binding of the RQPQVSI Peptide to HLA-B27 Subtypes with Disparate Association with Ankylosing Spondylitis-The RQPQVSI peptide bound to the HLA-B*2704 and -B*2706 subtypes with similarly low affinities (Fig.  3). To study the requirements of the binding of these peptides to the related HLA-B27 subtypes, the exchange of alanine for either the ArgP1 or the GlnP2 residues of the RQPQVSI peptide was tested in complex stability assays. The relative binding to the B*2704 subtype of both natural ligands and Ala-substituted peptides was identical (Fig. 4, left panel). Thus, the RQPQVSI peptide bound in similar ways to two different AS-associated subtypes, B*2705 and B*2704. On the other hand, the ArgP1 residue of the RQPQVSI peptide is critical for binding to the B*2706 subtype (Fig. 4, right panel) but not to the B*2704 sub-  type (Fig. 4, left panel). Thus, in the interactions of this subtype with the RQPQVSI peptide, the ArgP1 residue is the anchor motif and mediates binding to this subtype, as do other peptides lacking the N-terminal binding residue (29,30). In summary, the RQPQVSI peptide bound with a comparably low affinity, but with different conformation, to HLA-B27 subtypes associated (B*2705 and B*2704) or not associated (B*2706) with ankylosing spondylitis ( Table 2).

Homology of GlnP2 Ligands with Pathogenic Bacterial
Sequences-A comparison of the sequences of six natural HLA-B27 GlnP2 ligands identified in the current study with 12 species of bacteria associated with HLA-B27-dependent reactive arthritis (33) was carried out. Table 3 shows diverse bacterial sequences that differ from endogenous HLA-B27 peptides by 2-3 residues. Most of these changes are in surface-exposed regions of proteins; thus, they could modulate differential interactions with the T-cell receptor as compared with endogenous cellular HLA-B27 peptides as shown in the HLA-B*2705 crystal structure (34) (PDB: 1HSA). In addition, four proteins with HLA-B27 GlnP2 ligands, ADAM10, ␤ 2 -microglobulin, cutaneous T-cell lymphoma antigen, and SF3B2, were included in a network of connective tissue disorders and inflammatory disease following analysis with the Ingenuity Pathway Analysis software.

DISCUSSION
In this study, we have investigated several issues concerning the processing of natural HLA-B*2705 ligands, including the analysis of a refined binding motif and the binding of GlnP2 ligands to HLA-B27 subtypes differentially associated with AS. Previously, some HLA-B*2705 ligands with GlnP2 had been described (13,14). In the present study, this change in the P2 anchor motif was detected in the 3% of HLA-B*2705 ligands identified. Thus, these data confirm the refined HLA-B*2705 binding motif. In addition, five of six ligands identified by mass spectrometry analysis show the preferential residue usage both P1 and/or P⍀ positions previously described from canonical HLA-B*2705 ligands (14,35). Thus, the GlnP2 ligands could interact with the pockets of the HLA class I-presenting molecule in the same way those ArgP2 canonical ligands.
Five of the six ligands with GlnP2 were derived from four proteins included in a network of connective tissue disorders. ADAM10 is a chondrocyte-derived metalloproteinase involved in joint pathology (36,37). The ␤ 2 -microglobulin is associated with HLA class I heavy chains and reduces HLA-B27 misfolding promoting arthritis and spondylitis in HLA-B27-transgenic rats (38). HD-CL-01 is an antigen of cutaneous T-cell lymphoma, and polyarthritis in the presence of this lymphoma is a phenomenon previously described (39 -41). Lastly, SF3B2 is a component of the spliceosome complex (42) that has been linked to rheumatoid arthritis and several connective tissue diseases (43,44).
All these human HLA-B27 GlnP2 ligands have homology with distinct arthritogenic bacterial sequences that differ from these endogenous HLA-B27 peptides by 2-3 residues. Remarkably, several arthritogenic bacterial sequences possess changes in the ArgP1 of the RQPQVSI sequence that could alter the binding conformations between the different subtypes of HLA-B27 associated or not associated with AS.
The observation that HLA-B*2704 and -B*2706 subtypes share a very broad fraction of ϳ90% of their peptide repertoires (32) suggests that putative arthritogenic peptides could be confined in a first approximation to a relatively small portion of ligands shared by AS-associated subtypes such as B*2705 and B*2704 that are missing from the B*2706 peptide repertoire.  The self-derived ligands with GlnP2, analyzed in the current study, bind to AS-associated class I subtypes (HLA-B*2705 and -B*2704) with a broad range of affinities, whereas they have shown a low affinity in their interaction with the non-AS-associated HLA-B*2706 allele. Thus, the differences of affinity in self-derived ligands could trigger autoimmune tissue injury in AS-associated versus non-AS-associated HLA-B27 class I positive individuals. In addition, our study indicates that some unusual HLA-B*2705 ligands, such as the RQPQVSI peptide, bind with similar affinity but different conformations to both AS-associated and non-AS-associated HLA-B27 subtypes. The HLA-B27 subtype-dependent conformation of canonical ArgP2 self-peptides was previously described (45,46). Thus, these conformational differences, which could be detected by the fine specificity of autoreactive cytotoxic T lymphocytes, could be sufficient to initiate autoimmune damage in patients with AS-associated subtypes. Therefore, these kinds of peptides (short, with GlnP2, and with similar low affinity to all HLA-B27 subtypes tested but with dissimilar conformations in differentially AS-associated subtypes) must not be excluded from future research involving potential arthritogenic peptides.
In summary, these results inform us about the interaction between unusual GlnP2 self-derived ligands and different HLA-B27 subtypes in the pathogenic role of these class I molecules in the triggering of AS autoimmune disease.

TABLE 3
Alignment of natural HLA-B27 peptides with sequences of bacteria associated with B27-dependent reactive arthritis


text:
Genome-Wide Detection of Gene Coexpression Domains Showing Linkage to Regions Enriched with Polymorphic Retrotransposons in Recombinant Inbred Mouse Strains

Although gene coexpression domains have been reported in most eukaryotic organisms, data available to date suggest that coexpression rarely concerns more than doublets or triplets of adjacent genes in mammals. Using expression data from hearts of mice from the panel of AxB/BxA recombinant inbred mice, we detected (according to window sizes) 42−53 loci linked to the expression levels of clusters of three or more neighboring genes. These loci thus formed “cis-expression quantitative trait loci (eQTL) clusters” because their position matched that of the genes whose expression was linked to the loci. Compared with matching control regions, genes contained within cis-eQTL clusters showed much greater levels of coexpression. Corresponding regions showed: (1) a greater abundance of polymorphic elements (mostly short interspersed element retrotransposons), and (2) significant enrichment for the motifs of binding sites for various transcription factors, with binding sites for the chromatin-organizing CCCTC-binding factor showing the greatest levels of enrichment in polymorphic short interspersed elements. Similar cis-eQTL clusters also were detected when we used data obtained with several tissues from BxD recombinant inbred mice. In addition to strengthening the evidence for gene expression domains in mammalian genomes, our data suggest a possible mechanism whereby noncoding polymorphisms could affect the coordinate expression of several neighboring genes.

genetics of gene expression quantitative trait loci clustering of coexpressed genes transposable elements structural variants A significant component of gene expression variability in populations is due to variations in the DNA sequence (Schadt et al. 2003). Accordingly, genetic mapping studies have led to the identification of quantitative trait loci (QTL) linked to the expression levels of particular genes within cells and/or tissues (Göring et al. 2007). When the expression of a given gene associates with a genetic polymorphism that maps close to that gene's locus, the corresponding "expression QTL" (eQTL) is identified as a proximal eQTL (also called "cis-eQTL"). In other cases, when the expression level of a gene associates with a locus clearly distinct from that of the gene itself, it is defined as a distal eQTL (also called "trans-eQTL"). In the case of trans-eQTLs, it has been observed that a single genetic locus can show linkage to the abundance of the mRNA transcript of several genes across the genome, forming so-called trans-eQTL "hotspots" or "bands" (Grieve et al. 2008;Wu et al. 2008). One example would, for instance, be that of a cis-eQTL regulating the expression of a transcription factor, which in turns regulates the expression of many other genes belonging to a hotspot of trans-eQTLs.
Contrary to trans-eQTLs, investigators typically associate cis-eQTLs with the expression level of just one gene. The premise is that the cisacting polymorphism that is located in close proximity to the gene is likely to affect the regulatory machinery of that same gene and that machinery is unlikely to be shared by genes other than the ones that are immediately adjacent. Nonetheless, if the genome contained features that could influence the expression of several neighboring (and not necessary immediately adjacent) genes, one consequence from a genetic standpoint would be the clustering of several cis-eQTLs within a narrow genetic interval.
Although clustering of cis-eQTLs has not been reported yet, there is evidence that gene coexpression domains exist in several eukaryotic organisms (Michalak 2008;Elizondo et al. 2009). For instance, approximately 20% of the genes in Drosophila are arranged into clusters of similarly expressed genes, with the clusters spanning intervals from 20 to 200 kb and containing 10 to 30 genes each (Spellman and Rubin 2002). In mammals, coexpressed genes have been reported to cluster both at either short-range (1 Mb) or long-range (.10 Mb) levels (Woo et al. 2010). One particular case of short-range coexpression clusters concerns that of conserved clusters of paralogous genes arising from tandem duplication (such as, for instance, Hox, globin, and major histocompatibility complex genes). Beyond these paralogous clusters, it was reported in humans that the overall level of coexpression of genes was greater than expected by chance when the genes are located within distances smaller than 1 Mb, although the level of expression did not exceed that of more distant genes by a very large margin (Lercher et al. 2002). Likewise, in other studies reporting on the clustering of co-expressed genes in mammals, it was found that coexpression rarely concerned more than doublets or triplets of immediately adjacent genes (Sémon and Duret 2006;Purmann et al. 2007). In such cases, coexpression is generally believed to derive from the sharing of one regulatory element by adjacent genes. Nonetheless, clusters containing on average two to six coordinately regulated genes within 1-Mb intervals have been observed under special circumstances, such as in fibroblasts during replicative senescence (Zhang et al. 2003).
One limitation of many studies to date is that they have been not been performed within the framework of panels of individuals (or strains) with well-characterized genetic backgrounds. If "short-range" clustering of coexpressed genes could derive from physical elements in the genome, the impact of the latter would be easier to detect in situations in which they are polymorphic, such as in animals from genetic crosses. Moreover, because clusters of cis-eQTLs would all map to precise genomic regions, further analysis of these regions might reveal the nature of the polymorphisms associated with coordinate changes in gene expression.
To complement previous genetic studies reporting on QTL linked to cardiac morphologic characteristics, we used Illumina microarrays to obtain the profiles of cardiac gene expression in a panel of 24 mouse recombinant inbred strains (RIS). When performing linkage analysis to detect eQTLs for all detected genes, we observed several instances in which three or more cis-eQTLs clustered within small genomic intervals. Because such clustering of cis-eQTLs had not been reported previously, we used our dataset to analyze the characteristics of corresponding regions in a more systematic fashion. Likewise, to test to which extent this observation could be generalized to other tissues and/or crosses, we compared our findings with those obtained in other tissues from either the same or other mouse RIS panels. We found that clusters of cis-eQTLs could be detected in all the tissues from all genetic mouse crosses we tested and that coexpression of cis-eQTLs within these clusters reached very high levels. Further analysis of these regions revealed that they showed enrichment for particular types of structural polymorphisms.

MATERIALS AND METHODS
Detection of eQTLs in hearts from AxB/BxA mouse RIS The AxB/BxA mouse RIS originate from reciprocal crosses between the two parental C57BL/6 and A/J inbred strains and were derived from 20 generations of inbreeding of the F2 progeny of these two strains (Marshall et al. 1992). We had previously used a set of 24 RIS from that panel to detect QTLs linked to cardiac left ventricular mass (Llamas et al. 2007). Using the same 24 RIS, we extracted total RNA from cardiac left ventricles from four male mice for each strain and used Illumina MouseRef-8 v2.0 BeadChips to obtain the profile of gene expression in the tissues, as described previously (Llamas et al. 2009). The raw data were obtained by using the BeadStudio software (Illumina) and imported into the R programming environment. The data were processed and normalized using the Limma software (Smyth 2004). After filtering out genes not detected across the chips (by retaining only genes detected in more than 50% of the biologic replicates for at least one strain), we selected a set of 8725 genes for further analysis. Processed data have been submitted for public access to GeneNetwork (www.genenetwork.org; accession number GN421).
For genomic mapping, genomic DNA was extracted from spleens of all corresponding 24 RIS using the DNeasy tissue kit (QIAGEN, Mississauga, ON). All samples were hybridized at The Jackson Laboratory on the Affymetrix Mouse Diversity Array, which contains 623,124 single-nucleotide polymorphic (SNP) and 916,269 invariant genomic probes (Yang et al. 2009). Signal intensities were extracted from CEL files using the MouseDivGeno package (Didion et al. 2012), and genotyping was performed by comparing the intensity and contrast of signals in a given line with that in the parental strains (Simecek et al. 2011). In total, we detected 977 informative SNPs (meaning that they were polymorphic for at least one strain among all 24 strains from the panel) defining intervals averaging 2.59 6 2.95 Mb. The average value of the r 2 coefficient (calculated as a descriptor of linkage disequilibrium for all pairs of adjacent informative SNPs) was r 2 = 0.8 (where 0 is the value for perfect equilibrium and 1 is the value obtained when two markers have identical information).
To detect and map eQTLs, the data were analyzed with the "R-QTL" tool, and all other statistical analyzes were performed with the statistical language R. We used a detection threshold corresponding to a "logarithm-of-the-odds" score of 3.3, as suggested previously (Lander and Kruglyak 1995). For each eQTL, we then determined whether the transcription was regulated in cis or in trans by defining cis-eQTLs as those whose peak eQTL was within 1 Mbp of the physical location of the corresponding gene start. Of note, artifactual detection of eQTLs might theoretically occur when polymorphic SNPs occur within sequences corresponding to the probes used by the microarray. Although the vast majority of SNPs within probes have been shown to have no significant effect on hybridization efficiency for Illumina microarrays (Schurmann et al. 2012), we nonetheless used data from the Sanger website to detect all high-quality (score . 100) polymorphic SNPs, compared their positions with those of all probes in the microarray (as annotated in GeneNetwork), and verified that the polymorphisms had no impact the eQTL analysis.

Origin of datasets
Lists of transposable elements (TEs) that are polymorphic between the C57BL/6J and the A/J mouse strains were obtained from two different sources. The first one corresponded to a supplementary file from the recent publication of Nellåker et al. (2012), in which the authors developed a comprehensive catalog of TE variants across 18 mouse strains (Nellåker et al. 2012). The second source corresponded to the MouseIndelDB database (http://variation.osu.edu/mouse_indel/ index.html), which reports on structural variants that show polymorphisms between four inbred strains (Akagi et al. 2008(Akagi et al. , 2010. From both databases, we extracted the locations of elements that showed either "insertion" (i.e., present in C57BL/6J but absent in A/J) or "deletion" (i.e., present in A/J but absent in C57BL/6J) vs. the mm9 reference sequence of the whole genome from C57BL/6J. For simplicity, the aforementioned publications used the convention of referring to these two types of structural variants as polymorphic "insertion-deletions" (indels). Although the number of polymorphic TEs reported by MouseIndelDB is lower than that reported by Nellåker et al. 2012. (Supporting Information, Table S1), sequences in MouseIndelDB have been characterized in greater detail and contain useful annotations. Of note, the vast majority of TE sequences present in a given species are "fixed," so that in mice, they will be present in both C57BL/6J and A/J and thus nonpolymorphic between the two strains. To obtain data on the abundance of "fixed" TEs in mice, lists of TEs located within all protein-coding genes were downloaded from the TranspoGene database (http://transpogene.tau.ac.il). Recombination rates across the mouse genomes corresponded to the values calculated in a recent report (Brunschwig et al. 2012).
Data on genomic binding sites for several factors were obtained from several sources. A list of 33,172 regions associated with CCCTCbinding factor (CTCF) in chromatin from the hearts of adult C57BL/6J mice (as assessed by chromatin immunoprecipitation and massively parallel sequencing) was obtained from the ENCODE/LICR database of transcription factor binding sites using a custom track in the UCSC Genome browser. These data (Release 2, April 2012) correspond to the results of experiments performed by the laboratory of Bing Ren at the Ludwig Institute for Cancer Research. Regions corresponding to binding sites of transcription factors Gata4, Mef2A, Nkx2.5, Srf, and Tbx5 in cardiac chromatin (amounting to either 16,753, 1337, 20,573, 23,806, or 55,582 regions, respectively) corresponded to those published by Schlesinger et al. (2011). The 10,486 regions corresponding to the abundance of acetylated histone 3 sites and the 3596 binding sites for the p300 histone acetylase in cardiac chromatin corresponded to those published by Blow et al. (2010) and Schueler et al. (2012), respectively.
Lists of other eQTLs obtained in RIS were downloaded from the www.genenetwork.org web site using the Genograph tool. In short, the tool uses WebQTL to detect all eQTLs associated to the expression levels of genes within a given dataset (Chesler et al. 2005). Datasets used are listed in Table S2. In addition to data for whole eyes from AxB/BxA mice, we used data obtained with five tissues (eye, kidney, hippocampus, hypothalamus, and cerebrum) from BxD mouse RIS. The latter originate from crosses between the parental C57BL/6J (B6) and DBA/2J strains. After analysis, eQTLs were selected using a false-discovery rate threshold of 0.2. As for our own data, we defined cis-eQTLs as those whose peak eQTL was within 1 MB of the physical location of the corresponding gene start. Using the same parameters as for detection of cis-eQTL clusters (i.e., boxes containing at least three cis-eQTL separated by maximum interval of 500 kb), we calculated for each pair of analyzed tissues which proportion of cis-eQTL2containing regions overlapped between the two datasets.

Selection and comparative analysis of genomic regions
Clusters of cis-eQTLs were detected by defining regions in which cis-eQTLs were separated by maximum distances of either 250, 500, or 750 kb. Control clusters were defined using the same maximum intervals between genes detected by the Illumina array in mouse hearts and imposing a maximal limit on the overall size of control clusters to obtained clusters whose size was not significantly different from that of matching cis-eQTL clusters (Table S3). To further verify that both types of clusters had similar properties, we calculated the number of "Entrez" genes in each cluster using the biomaRt R package (version 2.10.0) (Durinck et al. 2009) interfaced to BioMart databases. Coexpression levels were quantified by calculating the absolute value of the Pearson correlation coefficient among expression levels of detected genes in the cis-eQTL clusters and compared with the coexpression levels observed in two other kinds of boxes: (1) control clusters (whose characteristics were similar to cis eQTL clusters in terms of size, number of genes detected in the heart by the Illumina microarray, total number of genes, and overall level of expression of detected genes); and (2) random regions (corresponding to boxes of similar size chosen randomly within the genome).
For comparisons of the abundance of structural variants and/or binding sites of regulatory factors, the regions analyzed were slightly larger than the clusters themselves and were selected by adding flanking regions of either 250, 500 or 1000 kb to four types of boxes: (1) the same cis-eQTL and control clusters defined above (using maximum intervals between cis-eQTLS or detected genes of either 250 or 500 kb; and (2) regions with the same size as the previous ones but either centered around single cis-eQTLs or selected randomly throughout the genome. Data calculated represented the number of features per Mb in each different region (cis-eQTL cluster, control cluster, and random region). For easy comparison across different types of features, all data were normalized by dividing them by the mean number of features (i.e., structural variants or binding sites) found in the random group. Accordingly, the mean normalized number of features in random groups was 1 (6 SD), and the values in other regions corresponded to "fold difference" compared to random regions.
Motif searches were performed in the sequences of polymorphic short interspersed elements (SINE) and long-terminal repeat (LTR)-TEs, using the HOMER bioinformatic package (http://biowhat.ucsd. edu/homer/motif). Sequences analyzed corresponded to those that were present in C57BL/6 but absent in A/J (referred to as C57 (+)/A/J(2) in Table 1), as these were the only ones where full sequence information was available. To test whether the regions containing polymorphic SINEs had specific characteristics, we use the "annotate_peaks" function provided by the HOMER package.

Statistics
Comparisons between groups were performed by Student's t-test (in the case of 2 groups) or by one-way analysis of variance followed be Tukey's honestly significant difference post-hoc multiple comparison tests (in the cases of comparisons between more than two groups). Differences in the relative abundance of total and polymorphic TEs in different genomic regions were tested by two-way analysis of variance followed by Tukey's honestly significant difference post-hoc tests.

Detection of cis-eQTL clusters
Using the Illumina MouseRef-8 microarray, we detected a total of 8725 genes expressed in hearts from the AXB/BXA mice. Further genomic mapping revealed that 777 of these genes were linked to cis-eQTLs, and several of them formed clusters of three or more cis-eQTLs within genomic intervals of a few hundred kilobases. We thus tested several window sizes to best define cis-eQTL clusters and control clusters (Table S3). By using maximum intervals of 250 kb, we detected a total of 42 cis-eQTL clusters (containing in average 4.23 6 1.9 detected genes, ranging from 3 to 11, within intervals averaging 221.9 6 130 kb), and 188 control clusters (containing in average 4.75 6 1.85 detected genes, ranging from 3 to 13, within intervals averaging 248 6 77.6 kb). By using maximum intervals of 500 kb, we detected a total of 53 cis-eQTL clusters (containing in average 4.9 6 3.5 cis-eQTL genes, ranging from 3 to 19, within intervals averaging 467 6 486 kb), and 59 control clusters (containing in average 5.2 6 2.5 detected genes, ranging from 3 to 17, within intervals averaging 456 6 119 kb). There were no significant differences between cis-eQTL and control clusters for any of the aforementioned values (Table S3). Because genes detected by the Illumina array in heart extracts do not correspond to all genes present in the genome, we also verified the density of all Entrez-annotated genes in the same regions to estimate total gene density, and found no significant difference between cis-eQTL and control clusters. When using maximum intervals of 750 kb, we detected a total of 61 cis-eQTL clusters but detected only 21 matching control clusters (Table S3). Further analyses were thus restricted to the "250 kb" and "500 kb" clusters. None of the cis-eQTL clusters corresponded to clusters of paralogous genes known to arise from tandem duplication. The coordinates of all 42 "250 kb" cis-eQTL clusters are listed in Table S9, along with the symbols of corresponding cis-eQTL genes. All control clusters and random regions are listed in Table S10 and Table S11.
Although SNPs within probes are not likely to affect the hybridization efficiency of Illumina microarray probes (Schurmann et al. 2012), we nonetheless used data from the Sanger website to detect all high-quality (score . 100) SPs that are polymorphic between A/J and C57Bl/6J and verified that no cis-eQTL within the clusters could represent an artifact resulting from a SNP polymorphism. Accordingly, we found a total of 91 SNPs falling within probe sequences. Among them, eight corresponded to cis-eQTLs, but none of them corresponded to those found in the cis-eQTL clusters.
Coexpression levels were quantified by calculating the absolute values of the Pearson correlation coefficients between each pair of genes within clusters. Within each cis-eQTL cluster, overall coexpression levels were calculated by the averaging all pairwise coexpression values and then compared with values obtained in either corresponding control clusters or in 500 "random groups." The latter comprised either 3 to 11 genes (for the "250 kb" clusters) or 3 to 19 genes (for the "500 kb" clusters), all genes being chosen randomly throughout the genome. For cis-eQTL genes within "250 kb" clusters, coexpression level was 0.755 6 0.07 (mean 6 SD), this value being significantly greater (P , 10e-16) than that obtained for detected genes in control clusters (0.23 6 0.09). Although mean coexpression level in control clusters was about 17% greater than in random groups of genes (0.196 6 0.04, P = 7.4e-07), it was also more than three times lower than that found in cis-eQTL clusters ( Figure 1A). Very similar results were obtained in terms of coexpression levels and intergroup differences when using the "500" kb clusters ( Figure 1B). The increased coexpression of genes within the cis-eQTL clusters was not due to an overall greater level of expression of genes in the cluster: the mean log2 value of expression level of genes within the cis-eQTL clusters was 9.023, this value being not significantly different (P = 0.9) from the level of expression of genes in control clusters (i.e., 9.02).
Given the average size of cis-eQTL clusters (i.e., 248 to 456 kb), the average size of intervals between polymorphic SNPs in the RIS panel (2.59 6 2.95 Mb) and the high level of linkage disequilibrium between adjacent informative SNPs (r 2 = 0.8), the vast majority of neighboring and coexpressed cis-eQTLs within cis-eQTL clusters are likely to have the same allelic origin. In contrast, cis-eQTL genes within cis-eQTL clusters did not show homogeneity either in terms of regulation (because consistent up-or down-regulation of cis-eQTL genes was found in only 26% of the cis-eQTL clusters) or in terms of genomic strand origin (because both strands of genomic DNA contributed to the sequences of neighboring and co-expressed cis-eQTL in 77% of the cis-eQTL clusters). There was little evidence to suggest that the cis-eQTL clusters could correspond to either recombinant blocks or to regions with different recombination rates. When defining minimal haplotype blocks as regions flanked by polymorphic markers, we found a total of 930 blocks whose average size (2.59 6 2.95 Mb) was considerably larger than that of cis-eQTL clusters (248 to 456 kb). Moreover, the average coexpression value of detected genes within these minimal haplotype blocks was 0.25 6 0.11 and thus much lower than that of cis-eQTLs within cis-eQTL clusters (0.755 6 0.07). Thus, high coexpression levels were found only for cis-eQTL genes within cis-eQTL clusters, and not for all detected genes throughout the haplotype blocks. Finally, the distribution of recombinant rate values in cis-eQTL clusters did not appear to be different from that of detected genes in control clusters ( Figure S1).

Structural characteristics of regions containing cis-eQTL clusters
The detection of clusters of cis-eQTLs suggested that genetic polymorphisms in some regions could associate with changes in the expression levels of several genes in the same region. Considering that it would be unlikely that such coordinate changes in the regulation of several neighboring genes could result from SNPs each affecting the expression of corresponding cis-eQTL genes in a proportionate manner, we mined databases to question whether cis-eQTL cluster regions could show enrichment in structural variants (with potential of affecting expression of all cis-eQTLs in the region). Because the majority (i.e., 98%) of mouse structural variants have been reported to correspond to TE variants (Akagi et al. 2008;Yalcin et al. 2011), we used data from the most recent report that established a catalog of TE variants across mouse strains (Nellåker et al. 2012) to test whether cis-eQTL clusters and their surrounding regions would contain more TE variants than either control clusters or regions of similar size centered around single cis-eQTLs.
The respective abundances of TEs that were reported as polymorphic between A/J and C57BL/6J are listed in Table S1. Considering that regulatory regions can be located either upstream or downstream of the genes under consideration, we defined the regions to be analyzed by adding flanking sequences with lengths of either 250, 500, or 1000 kb to the regions corresponding to both types of clusters, thus corresponding to regions of six different sizes (Table S4 and Table S5). For each of the six sizes of regions, comparisons were made between the three types of "defined" regions (cis-eQTL clusters, control clusters, and regions centered on single cis-eQTLs) and the fourth type of region, consisting of random regions of matching size. Overall, the same interregion differences were found regardless of how the regions were defined. For simplicity, the regions corresponding to the "250 kb" clusters augmented by flanking regions of 250 kb were chosen as representative data for presentation ( Figure 2). Long interspersed elements (LINEs) and LINE fragments showed some minor (usually nonsignificant) differences in abundance between the four types of regions. Both polymorphic and fixed SINEs and LTRs were more abundant in all three defined regions than in random regions. For SINEs, this finding is in keeping with the fact that these elements are generally more abundant in gene-rich than in gene-poor regions (Jurka et al. 2005). However, for polymorphic TEs, only SINEs were significantly greater (P , e-5) in cis-eQTL clusters compared with both control clusters and regions centered on single eQTLs. Moreover, the foldenrichment and the significance of these differences were of greater magnitude for polymorphic SINEs than for fixed SINEs, which indicated that the greater abundance of polymorphic SINEs in cis-eQTL clusters was not a mere consequence of total abundance of TEs but rather a consequence of the genetic differences between C57BL/6J and A/J mice. These differences were not due to differences in total gene density because the abundance of total genes in cis-eQTL regions (10.1 6 17) was not significantly different than that in control regions (10.18 6 3.5). Representative examples that compare two cis-eQTL clusters and control clusters of matching sizes are shown in Figure 3. Altogether, polymorphic SINEs appeared to be a signature characteristic of cis-eQTL regions. We also found that the density of polymorphic SINEs in cis-eQTL cluster regions (taken along with their 250 kb flanking regions) was calculated to be 10.3 6 10.7 polymorphic SINEs/MB; this value was significantly greater (P , 0.002) than that found in other regions of corresponding haplotype blocks outside of the cis-eQTL clusters (5.1 6 8.1). This finding provided additional evidence that the cis-eQTL clusters had features that differentiated them from the haplotype blocks that contained them.
To test the possible functional impact of polymorphic SINEs, motif enrichment analyses were performed for TE sequences that are present in C57BL/6 and deleted in A/J mice (because the full sequences of TEs deleted in C57BL/6 and present in A/J are not available yet). The list of most significantly enriched binding sites is shown in Table S6. For polymorphic SINEs, the binding site that was most statistically enriched (P = 1e-1283) was that previously reported to be bound by BORIS, a CTCF paralogue that binds a CTCF-like binding site (Pugacheva et al. 2010). This site matched the composition of the M1 moiety of the full CTCF binding site recently described by Schmidt et al. (2012). Most other significantly enriched sites corresponded to binding sites for transcription factors from several families, in addition to another CTCF binding site that in fact matches the composition of the M2 moiety of the full CTCF binding site (Schmidt et al. 2012).
According to MouseIndelDB annotation, 19% of SINEs that are polymorphic between C56BL/6J and A/J correspond to B1 SINEs, with the remainder corresponding to B2 SINEs. We also investigated the nature of regions harboring polymorphic TEs (Table 1). Notwithstanding a few exceptions, the great majority of polymorphic TEs fell into either intronic or intergenic regions. Not surprisingly, C57(+)/ A/J(2) LTR-TEs and SINEs fell into C57 regions previously annotated as having these characteristics. Interestingly, C57(2)/A/J(+) LTR-TEs sometimes fell into regions annotated as containing LINEs or SINEs in C57BL/6J, whereas C57(2)/A/J(+) SINEs sometimes fell into regions already containing LTR-TEs or LINEs. Beyond polymorphic TEs, we obtained from previous publications (reporting the results of Chip-Seq experiments performed on mouse heart chromatin) lists of all regions corresponding to either binding sites for transcription factors or chromatin modifications. We then tested whether all corresponding genomic features showed differential abundance between cis-eQTL, control, and random regions by using the six different types of window sizes described previously (Table S7 and Table S8). Overall, the same interregion differences were found regardless of how the regions were defined, but differences tended to be more pronounced for regions surrounding and comprising the "250 kb" clusters. Accordingly, the regions corresponding to the "250 kb" clusters augmented by flanking regions of 250 kb were chosen as representative examples for presentation ( Figure 4). Overall, the abundance of all binding sites was lower in random regions than in all types of defined regions. For some regulatory factors (CTCF, H3Ac, SRF, and Tbx5), their abundance was much greater in cis-eQTL regions than in the other two types of defined regions control regions and more abundant in the latter than in random regions (Figure 4).

Comparisons with other panels of RIS
To test to which extent cis-eQTL clusters would be conserved across tissues, we questioned whether regions containing cis-eQTL clusters for cardiac genes would overlap with regions containing cis-eQTL clusters for genes expressed in another tissues (with clusters of cis-eQTLs being defined on the basis of maximum intervals of 250 Mb between each cis-eQTL). We first analyzed gene expression data obtained in AxB/BxA eyes with Illumina microarrays (Table S2), which allowed us to detect a total of 35 cis-eQTL clusters (while verifying, as explained previously, that none of the cis-eQTL genes could represent an artifact due to the presence of a SNP polymorphism within the sequence of the probes). Despite the fact that the latter data were obtained by other investigators, 12 of the 42 regions containing cis-eQTLs clusters for genes expressed in AxB/BxA hearts also contained cis-eQTL clusters for genes expressed in eyes from the same RIS.We also analyzed gene expression data for other tissues from the BxD RIS mouse panel, where gene expression was analyzed (Table  S2) with either the Affymetrix MoGene 1.0 ST microarray (for hypothalamus) or with the Affymetrix Mouse Genome 430 microarray (for eye, kidney, hippocampus, and cerebellum). Of note, consultation of the list of SNP polymorphisms between C57BL/6J and DBA/2 mice revealed a total of 0.005% and 0.0006% of probes used by the Affymetrix MoGene 1.0 ST and Affymetrix Mouse Genome 430 microarrays, respectively, are affected by such polymorphisms. Corresponding genes were excluded from the cis-eQTL analysis. Analysis of the gene expression data from BxD RIS tissues allowed us to detect 52, 279, 260, 77, and 36 cis-eQTL clusters for kidneys, eyes, hippocampus, Figure 3 Representative examples illustrating two cis-eQTL clusters and control clusters of matching size. For each cluster, the first line (on top) represents respective genomic scales (in Mb); the second line represents the genomic positions of either cis-eQTLs or detected genes; and the third and fourth lines represent the genomic positions of polymorphic LTR-TE and SINEs. On the second line, cis-eQTL genes are represented at their respective position by arrows; the arrow for the first gene on the left points upwards; other genes are represented by either upwards or downward arrows, according to whether they correlated positively or negatively, respectively, with the change in expression of the first gene. On the third and fourth lines, plain vertical lines correspond to TEs that are C57(+)/A/J(2); dashed lines correspond to TEs that are C57(2)/A/J(+).
hypothalamus, and cerebrum, respectively. Again, high proportions of cis-eQTL2containing regions were shared across tissues. For instance, of the 52 regions detected for genes expressed in kidneys, the number of regions also detected for genes from other tissues amounted to 44 for eyes, 44 for hippocampus, 15 for hypothalamus, and 15 for cerebrum. Finally, to test whether regions containing cis-eQTL clusters in one tissue would overlap with regions containing cis-eQTL clusters in one same tissue from two RIS panels, we compared data obtained for eye genes from AxB/BxA and BxD RIS. Of the 35 regions containing cis-eQTL clusters for genes from AxB/BxA eyes, 16 corresponded to regions containing cis-eQTL clusters for genes from BxD eyes.

DISCUSSION
Despite evidence suggesting the existence of gene coexpression domains in several eukaryotic organisms (Michalak 2008;Elizondo et al. 2009), data available to date suggested that coexpression in mammals rarely concerns more than doublets or triplets of immediately adjacent genes (Sémon and Duret 2006;Purmann et al. 2007). Our data show that within panels of mouse RIS, where gene expression variability is due in part to genetic polymorphisms, it is indeed possible to detect short-range clustering of more than three neighboring (but not necessarily adjacent) coexpressed genes. Moreover, all coexpressed genes within these domains showed linkage to loci having the same position as that of the domains, thus showing that genetic polymorphisms can associate with the expression levels of several neighboring genes within these domains. Depending on the sizes of windows used to define them, cis-eQTL clusters detected for AxB/BxA hearts contained in average 4.2 to 4.9 cis-eQTL genes within intervals averaging 248-456 kb. The number of cis-eQTL clusters detected for hearts from AxB/BxA panel (i.e., 53) was within the same range as the numbers of cis-eQTL clusters detected with three other tissues from the distinct BxD RIS panel (i.e., 60 for cerebrum, 74 for kidneys, and 96 for hypothalamus). Moreover, a great proportion of the regions containing cis-eQTL for AxB/BxA hearts were identical to those containing cis-eQTL clusters for eyes from either AxB/BxA or BxD RIS, and many of the regions containing cis-eQTL clusters for BxD kidneys overlapped with those detected for four other tissues from the same panel. Altogether, the data indicate that cis-eQTL clusters occur in a fairly robust and consistent manner across different tissues and/or genetic backgrounds in mouse RIS. Of note, all tissues do not necessarily express the same genes, which may explain in part why some cis-eQTL clusters are distinct between tissues. With some tissues (eyes and hippocampus from BxD RIS), up to 300 cis-eQTL clusters were detected. A comprehensive list of all genomic regions having potential to contain cis-eQTLs and a full understanding of the extent to which they overlap would require data from a more exhaustive list of tissues.
One distinct feature of cis-eQTL clusters was that cis-eQTL genes in these regions showed a level of coexpression that greatly exceeded that found for detected genes in control regions with similar gene density. For instance, in regions selected as controls for the "500 kb" cis-eQTL clusters, only 2 of 59 regions displayed coexpression levels greater than 0.56 (i.e., the value corresponding to the lowest value for gene coexpression in cis-eQTL clusters). The cis-eQTL clusters thus corresponded to "gene coexpression domain" QTLs. The fact that all coexpressed genes within domains showed linkage to a common locus suggested that the genome contains polymorphisms that can alter coexpression levels. Reasoning that identification of the nature of such polymorphisms could reveal insights as to what drives coexpression of genes within coexpression domains, we mined databases to test whether the cis-eQTL regions harbored any particular types of polymorphic structural variants. Accordingly, we found that polymorphic SINEs [either C57(+)/A/J(2) or C57(2)/A/J(+)] were significantly more abundant in cis-eQTL clusters than in any other type of control regions. Given the possibility that polymorphic SINEs could in fact drive these high levels of coexpression in corresponding domains, we tested whether they showed enrichment for particular motifs. In addition to binding sites for various transcription factors belonging to different families, polymorphic SINEs showed enrichment for two sites corresponding to CTCF-binding regions. This finding is an agreement with the previous report in which the authors showed that CTCF-binding regions in the mouse genome were preferentially embedded in B2 SINE elements (Bourque et al. 2008). B2 SINEs constitute types of SINEs that are specific to rodent genomes, where they have undergone waves of amplification (Kass et al. 1997). Accordingly, 81% of the SINEs that are polymorphic between C57BL/6J and A/J were in fact B2 SINEs. Moreover, we compared the relative abundance of either chromatin marks or binding sites for either CTCF or cardiac transcription factors in cis-eQTL clusters vs. all three other types of regions. Accordingly, we found that binding sites for the transcription factors SRF and TBX5, the chromatin-organizing factors CTCF and p300, and the H3Ac chromatin were all significantly enriched in cis-eQTL clusters vs. all three other regions.
SINEs and CTCF binding sites are of particular interest. The structural and regulatory organization of the mammalian genome is fundamentally dependent on CTCF, which has been dubbed the "master weaver of the genome" (Phillips and Corces 2009). CTCF generally acts as an insulator preventing the spread of inactive heterochromatin, and is often associated with open chromatin (The ENCODE Project Figure 4 Relative abundance of binding sites for CTCF, SRF, and Tbx5 and of the H3Ac chromatin marks previously reported in chromatin from mouse hearts cis-eQTL clusters, control clusters, single cis-eQTL regions, and random boxes. Consortium 2012). This may be particularly pertinent in the context of our current data further documenting the existence of genetically controlled gene coexpression domains. As a matter of fact, recent genome-wide studies on chromatin structure have revealed that mammalian genomes are organized into topological domains, the boundaries of which show enrichment for SINEs and CTCF binding, and where the spread of heterochromatin is constrained (Dixon et al. 2012). Likewise, it was recently shown that in several mammalian species, CTCF-binding events are associated with waves of retrotransposon expansion, thus revealing the mechanism by which there are born (Schmidt et al. 2012). Our data extend these previous reports in several ways: (1) in addition to interspecies differences, polymorphic SINEs can cause differential abundance of CTCF binding sites between strains of one same species; (2) the CTCF-dependent organization of genomic mammalian domains may not be static because polymorphic SINEs could reshape that organization; and (3) polymorphic SINEs and CTCF-binding sites may constitute a mechanism defining gene coexpression domains, for which there was so far little evidence in mammals beyond doublets of triplets of genes.
Of note, the allelic origin of the cis-eQTL cluster regions did not affect the expression of all corresponding cis-eQTL in a consistent manner (as illustrated in Figure 3). However, because cis-eQTL clusters contained an average of five polymorphic TEs, changes in corresponding chromatin domains may be more complex than simply corresponding to a chromatin structure that is entirely open or closed for the whole region. Moreover, TEs can affect single gene expression by a variety of different possible mechanisms, as for instance by providing alternative promoters or enhancers, serving as insulators or transcriptional silencers, disrupting the exon-intron structure and/or causing premature transcriptional termination (Gogvadze and Buzdin 2009). TE might also affect gene regulation by mechanisms other than just providing regulatory elements. For instance, we found that some C57(2)/A/J(+) polymorphisms fell themselves within regions containing other TEs, and might thus disrupt in A/J the organization of some TE elements that have regulatory effects within the C57BL/6J strain. In combination with complex changes in chromatin structure, all the above mechanisms might account for the nonuniform effects of the allelic origin of cis-eQTL clusters on gene expression.
Close to half of mammalian genomes is derived from ancient TEs (primarily retroelements) (van de Lagemaat et al. 2003). Given their great abundance in the genome and the emerging recognition of their role in gene regulation, TEs that show polymorphism between inbred strains are increasingly being recognized as potential players in the genetics of quantitative traits (Nellåker et al. 2012). Accordingly, some classes of polymorphic TEs have been reported to show small but significant enrichment in refined genomic intervals selected on the basis of previous detection of quantitative trait loci for mouse quantitative traits (Nellåker et al. 2012). In terms of gene regulation, recent reports on the effects of TEs have so far concerned mostly single genes (Chernova et al. 2008;Palmer and Dulawa 2010;Li et al. 2012). However, complex quantitative traits are usually considered to result from the combined regulatory effects of several genes rather than from the highly penetrant effect of single mutations (Manolio et al. 2009). Gene coexpression domains might thus be of particular interest, as coordinate dysregulation of the expression levels of several genes within cis-eQTL clusters could possibly have greater effects on the phenotypic expression of complex traits than dysregulation of single genes.
Regulatory sequences located outside of coding regions are also interesting in the light of evidence that: (1) up to one third of noncoding sequence variation contributes causally to the traits under investigation in genome-wide association studies (Visel et al. 2009); (2) a great proportion of regulatory variants of gene expression are found at a fairly great distance from transcription start sites (Teng et al. 2011);and (3) chromatin structure plays important roles in the organization and regulation of our genes (The ENCODE Project Consortium 2012). Genome-wide approaches allowing the discovery and functional characterization of such elements might improve our understanding of their role in human biology and disease susceptibility (Teng et al. 2011). However, if polymorphic TEs turn out to have important consequences on gene regulation, appropriate technologies to detect them genome-wide will need to be developed because technologies based on detection of SNP polymorphisms are not sufficient in this regard.


text:
Detection of Ascitic Feline Coronavirus RNA from Cats with Clinically Suspected Feline Infectious Peritonitis

ABSTRACT Ascitic feline coronavirus (FCoV) RNA was examined in 854 cats with suspected feline infectious peritonitis (FIP) by RT-PCR. The positivity was significantly higher in purebreds (62.2%) than in crossbreds (34.8%) (P<0.0001). Among purebreds, the positivities in the Norwegian forest cat (92.3%) and Scottish fold (77.6%) were significantly higher than the average of purebreds (P=0.0274 and 0.0251, respectively). The positivity was significantly higher in males (51.5%) than in females (35.7%) (P<0.0001), whereas no gender difference has generally been noted in FCoV antibody prevalence, indicating that FIP more frequently develops in males among FCoV-infected cats. Genotyping was performed for 377 gene-positive specimens. Type I (83.3%) was far more predominantly detected than type II (10.6%) (P<0.0001), similar to previous serological and genetic surveys.

Feline infectious peritonitis (FIP) is an immune-mediated progressive infectious disease of domestic cats and wild felids caused by infection with feline coronavirus (FCoV), a single-stranded RNA virus classified as Alphacoronavirus [1,28]. It is divided into effusive FIP, in which exudate is observed in the body cavity, and non-effusive FIP, in which multiple pyogranuloma lesions are formed, but differences in lesions are influenced by the immunity of individuals [3]. This disease is considered to be induced by a mutant of feline enteric coronavirus (FECV), FIP virus (FIPV) [8,26]. The incidence is generally about 1% in FCoV-infected cats, but it varies depending on the age, breed, environment and superinfection with other viruses [7,8,17,23]. There are 2 genotypes of FCoV, and FCoV type II is considered to arise by recombination between FCoV type I and canine coronavirus (CCoV) [9,27]. Based on serological and genetic investigations, FCoV type I is overwhelmingly dominant over FCoV type II, and mixed infection with both types is not rare [6,11,14,15]. Since FIPV and FECV cannot be distinguished serologically or genetically, it is generally difficult to diagnose FIP [1], but it has been stated that the demonstration of FCoV RNA in ascites is one of most reliable diagnostic indicators of FIP [3,10]. Several results concerning the state of FIP incidence have been reported overseas [6,15,16,18,29]. In Japan, there have been several reports on the prevalence of FCoV infection [11,13,20,24], but, to our knowledge, no studies on FIP prevalence have been reported. Thus, in this study, to investigate the state of FIP incidence in Japan, we examined FCoV RNA in ascites of domestic cats suspected of having FIP using RT-PCR and investigated the genotype of detected FCoV RNA.
Ascites was obtained from a total of 854 privately owned domestic cats that were clinically suspected of having FIP from veterinary clinics throughout Japan in the period of June 2008 to October 2012. Viral RNA was extracted using a QIAamp Viral RNA Mini Kit (Qiagen, Tokyo, Japan). Detection of FCoV RNA was performed using a QIAGEN On-eStep RT-PCR Kit (Qiagen) with the P205-P211 primer pair providing a 223 bp amplicon of the 3'-untranslated region (UTR) gene of FCoV [10], according to the manufacturer's instructions. FIPV strain 79-1146 (American Type Culture Collection, Rockville, MD, U.S.A.) was used as a positive reference virus. Furthermore, for samples in which the 3'-UTR gene was detected, genotyping was performed with the nested RT-PCR reported by Addie et al. [4], providing 360 and 218 bp amplicons of FCoV types I and II, respectively. FIPV strain UCD1 (kindly provided by Azabu University) and FIPV strain 79-1146 (same as above) were used as positive reference viruses for the types I and II, respectively. Gene positivity was analyzed by chi-square test (StatView, Adept Scientific, Herts, U.K.), and a P value of <0.05 was considered statistically significant.
The 3'-UTR gene was detected in 377 of 854 cats (44.1%). Grouped by breed, the gene positivity was significantly higher in purebreds (62.2%; 181 of 291) than in crossbreds *CorrespondenCe to: Soma, T., Veterinary Diagnostic Laboratory, Marupi Lifetech Co., Ltd., 103 Fushiocho, Ikeda, Osaka 563-0011, Japan. e-mail: takehisa-soma@ds-pharma.co.jp # Contributed equally to the work and therefore should be considered equivalent authors. ©2013 The Japanese Society of Veterinary Science (34.8%; 196 of 563) (P<0.0001, χ 2 =58.4). Thus, age-related changes in positivity were investigated in purebreds and crossbreds. As shown in Fig. 1, the highest positivity was noted in cats younger than 1 year of age in both breeds and decreased with aging after 1 year of age. The positivity was higher in the purebreds than in the crossbreds until 5 years of age, and a significant difference was noted between the purebreds and crossbreds at 1 year of age and younger (P=0.0001, χ 2 =14.6 and P=0.0119, χ 2 =6.3, respectively).
Then, the positivity was investigated in various purebreds. As shown in Table 1, the positivities in the Norwegian forest cat (92.3%) and Scottish fold (77.6%) were significantly higher than the average for purebreds (62.2%) (P=0.0274, χ 2 =4.9 and P=0.0251, χ 2 =5.0, respectively), whereas that in the American shorthair (30.0%) was significantly lower than the average for purebreds (P<0.0001, χ 2 =26.3); the positivity in the American shorthair was also lower than that in crossbreds, although the difference was not significant. The  The breeds that consisted of fewer than 10 cats were classified together as others. A significance test (χ 2 test) was performed between each breed and the average of purebreds (62.2%).
positivity was low in spring-summer and high in fall-winter, although no significant differences were observed (38.8, 40.9, 47.3 and 49.5% in spring, summer, fall and winter, respectively). Grouped by gender, the positivity was significantly higher in males (51.5%; 235 of 456) than in females (35.7%; 142 of 398) (P<0.0001, χ 2 =21.7). The positivity was investigated in 10 regions of Japan (Hokkaido, Tohoku, Shinetsu, Kanto, Hokuriku, Tokai, Kinki, Chugoku, Shikoku and Kyushu), but no marked regional difference was noted (data not shown). Furthermore, the genotype was investigated in 377 specimens positive for the 3'-UTR gene. As shown in Table 2, FCoV type I (83.3%) was more predominantly detected than FCoV type II (10.6%) (P<0.0001, χ 2 =399.8), and both types were detected in 3.4% of specimens. In 36 specimens, no amplicons were observed; therefore, their genotypes could not be identified. There were no significant differences in the positivity of either type between the purebreds and crossbreds older than 1 year of age. But, at younger than 1 year of age, the positivity of FCoV type I was significantly higher in the crossbreds (94.9%) than in the purebreds (75.0%) (P=0.0003, χ 2 =13.3), whereas that of FCoV type II was significantly higher in the purebreds (23.5%) than in the crossbreds (3.8%) (P=0.0002, χ 2 =13.9). Of 40 cats in which FCoV type II was detected, 31 (77.5%) were purebreds younger than 1 year of age.
Since demonstration of FCoV RNA in ascites of FIPsuspected cats using RT-PCR was a useful indicator to diagnose FIP [3,10], the results of this study likely represent the epidemiology of effusive FIP in domestic cats in Japan. However, whether or not cats negative for the 3'-UTR gene were infected with FCoV was unknown. This status needs to be investigated further, such as by antibody test and quantitative RT-PCR.
In previous surveys, the incidence of FIP was higher in purebreds than in crossbreds, and it was particularly high in specific breeds [16][17][18][19]29]. Similarly, in this study, the positivity of the 3'-UTR gene was significantly higher in purebreds than in crossbreds, particularly in the Norwegian forest cat and Scottish fold, whereas the positivity was lower in the American shorthair. However, in previous reports from the U.S.A. and Australia, the incidence of FIP was high in the Abyssinian, Bengal, Birman, Himalayan, Ragdoll and Rex and the Australian mist, British shorthair, Burmese and Rex, respectively, and low in the Burmese, Exotic shorthair, Manx, Persian, Russian blue and Siamese and the Persian, respectively [16,18]. These breeds are inconsistent with our findings. The breeds that had high and low positivity in this study are consistent with FCoV antibody prevalence in Japan [24]. These findings indicate that countries and regions vary in the incidence of FIP by breed, and the prevalence of FCoV in breeding colonies seems to be closely related to the incidence, rather than bloodlines.
The positivity was highest in cats younger than 1 year of age and markedly decreased with aging. These findings are consistent with the incidence of FIP in previous reports [2,6,16,17,19,29]. At younger than 1 year of age, the positivity was significantly higher in purebreds than in crossbreds, indicating high prevalence of FCoV in multi-cat households, such as pet shops and breeding catteries [8]. The positivity was also significantly higher in males than in females, as previously reported [6,16,17,29]. However, no gender difference has generally been noted in the antibody prevalence in previous studies [5,12,20,24], indicating that FIP more frequently develops in males among FCoV-infected cats, although no serological data were shown in this study. This may be due to the influences of the stress of fighting and infectious diseases, such as feline leukemia virus and feline immunodeficiency virus infections, since these viruses more frequently infect males than females [22]. On investigation by season, the positivity tended to increase in fall-winter. A similar tendency has been reported in the U.S.A. [19], suggesting the seasonality of FIP development. This may be related to weather and the estrus cycle.
FCoV type I was far more predominant than FCoV type II, and this is consistent with previous serological and genetic surveys [6,11,14,15]. However, most of the FCoV type II cases (77.5%) were detected in purebreds younger than 1 year of age. Most purebred cats were likely to have been closely maintained with dogs in pet shops. Considering the fact that the CCoV infection rate is very high in dogs maintained in multi-dog households [21,25], this finding epidemiologically may indicate the possibility that FCoV type II is a recombinant of FCoV type I and CCoV [9,27], although further detailed analyses are necessary for clarifying this relationship. The above findings may help to clarify the state of FIP incidence in pet cats in Japan. But, the results were from cats with effusive FIP and are not necessarily applicable to the non-effusive FIP. Moreover, genotype-associated differences in the pathogenicity could not be investigated, because the severity and course of symptoms of the cats from which the test specimens were collected were unclear. It may be necessary to perform epidemiological surveys based on further detailed information and to investigate specimens from cats with non-effusive FIP.