Genomic organization and evolution of the ULBP genes in cattle

Background The cattle UL16-binding protein 1 (ULBP1) and ULBP2 genes encode members of the MHC Class I superfamily that have homology to the human ULBP genes. Human ULBP1 and ULBP2 interact with the NKG2D receptor to activate effector cells in the immune system. The human cytomegalovirus UL16 protein is known to disrupt the ULBP-NKG2D interaction, thereby subverting natural killer cell-mediated responses. Previous Southern blotting experiments identified evidence of increased ULBP copy number within the genomes of ruminant artiodactyls. On the basis of these observations we hypothesized that the cattle ULBPs evolved by duplication and sequence divergence to produce a sufficient number and diversity of ULBP molecules to deliver an immune activation signal in the presence of immunogenic peptides. Given the importance of the ULBPs in antiviral immunity in other species, our goal was to determine the copy number and genomic organization of the ULBP genes in the cattle genome. Results Sequencing of cattle bacterial artificial chromosome genomic inserts resulted in the identification of 30 cattle ULBP loci existing in two gene clusters. Evidence of extensive segmental duplication and approximately 14 Kbp of novel repetitive sequences were identified within the major cluster. Ten ULBPs are predicted to be expressed at the cell surface. Substitution analysis revealed 11 outwardly directed residues in the predicted extracellular domains that show evidence of positive Darwinian selection. These positively selected residues have only one residue that overlaps with those proposed to interact with NKG2D, thus suggesting the interaction with molecules other than NKG2D. Conclusion The ULBP loci in the cattle genome apparently arose by gene duplication and subsequent sequence divergence. Substitution analysis of the ULBP proteins provided convincing evidence for positive selection on extracellular residues that may interact with peptide ligands. These results support our hypothesis that the cattle ULBPs evolved under adaptive diversifying selection to avoid interaction with a UL16-like molecule whilst preserving the NKG2D binding site. The large number of ULBPs in cattle, their extensive diversification, and the high prevalence of bovine herpesvirus infections make this gene family a compelling target for studies of antiviral immunity.


Background
The cattle Major Histocompatibility Complex Class I-like Gene Family A (MHCLA) was initially discovered in a cattle spleen cDNA library during a search for highly divergent mammalian genes [1]. Two transcripts, MHCLA1 and MHCLA2, were found to be members of the MHC Class I superfamily, encoding cell-surface transmembrane proteins containing α1and α2-like domains, but no α3-like domain. These molecules have peptide sequence similarity to their homologues in other mammalian species, including the ULBP and RAET1 molecules in humans [2,3] and the H60, RAE1 and MULT1 molecules in mice [4][5][6][7]. To establish consistency with the human nomenclature, the cattle MHCLA1 and MHCLA2 genes are renamed ULBP1 and ULBP2, respectively, in this study. The function of cattle ULBP molecules is not known, but the human and mouse homologues have been demonstrated to interact with the NKG2D receptor, leading to activation of natural killer (NK) cells and T cell subsets in antitumour and infectious disease immunity [8]. In vitro studies have demonstrated that the soluble human cytomegalovirus (hCMV) protein UL16 interferes with the ability of ULBP1 and ULBP2 to interact with NKG2D, and coexpression of UL16 with ULBP1 or ULBP2 results in cytoplasmic retention of the ULBP molecules [2,9,10].
Southern blot analysis revealed the existence of a high copy number of ULBP genes in the cattle genome and seven other ruminant genomes. It was thus hypothesized that the cattle ULBP genes evolved rapidly by duplication and sequence divergence in response to selective pressure exerted by a viral pathogen(s). Extensive duplication of the cattle ULBP genes may serve to increase the repertoire of ULBP molecules able to bind NKG2D to initiate an immune response even in the presence of a UL16-like molecule [1].
The purpose of the present study was to identify the number of ULBP genes in cattle and describe their genomic organization. Six cattle bacterial artificial chromosome (BAC) clones were sequenced, resulting in the identification of 30 ULBP loci organized in two gene clus-ters on BTA9. Sequence analysis of the paralogues revealed that extensive gene duplication led to the present organization of the ULBP gene clusters. Bioinformatics tools were employed to characterize domains and sequence motifs in ten ULBP genes predicted to encode cell surface molecules, the majority of which are predicted glycoproteins. Substitution analysis identified specific codons in these genes that appear to be under positive Darwinian selection, and these selected sites were interpreted in a structural context using homology modelling.

Identification of the minor and major ULBP gene clusters
Four minimally overlapping ULBP-containing BACs were identified by hybridization-based screening with a fulllength cattle ULBP1 clone and then sequenced: RP42-147E22 [GenBank: AC092858], RP42-152A4 [GenBank: AC096629], RP42-146C17 [GenBank: AC098686] and RP42-194O5 [GenBank: AC098687]. Sequence alignment revealed that the former three BACs were overlapping, and the latter BAC was a singleton. Using BAC-end sequence data, two additional minimally overlapping BAC clones were identified: RP42-522F4 [GenBank: DQ405274] and CHORI240-21B24 [GenBank: DQ405273]. The overlapping clones were used to reconstruct two gene clusters, termed the "minor" ULBP cluster [GenBank: DP000082], spanning 331,973 bp, and the "major" ULBP cluster [Gen-Bank: DP000081], spanning 464,586 bp ( Table 1). The minor and major cluster sequences could not be further extended or joined by querying publicly available cattle genome sequence data [NCBI Build 2.0]. The ULBP1 locus [1] was not identified in this study, and therefore the major ULBP cluster sequence may be incomplete upstream of ULBP7.
Four ULBP loci were identified within the minor cluster, and 26 ULBP loci were identified in the major cluster. Nine loci represent coding sequences, and 21 loci are probable pseudogenes. Exons were identified by alignment and manual inspection ( Table 2, 3). Loci were designated as genes if they contained uninterrupted coding sequence in the signal peptide, α1 and α2 domains. Loci either lacking an exon corresponding to the signal peptide, α1 or α2 domains or containing a stop codon in the coding sequence of one of these three domains were designated as pseudogenes. Many of the pseudogenes contain exons with intact coding sequence ( Table 2, 3). It may be speculated that these pseudogenes serve as a repository for generating novel ULBP paralogues through gene conversion.
The nine ULBP genes identified in this study have a canonical five exon structure. An exception is ULBP21, which has six exons; the first two exons encode the signal peptide. All nine ULBP genes contain GU/AG exon splicing motifs. Because of the high degree of interlocus sequence identity among ULBP genes (e.g., ULBP9 and ULBP27 have 99.8% nucleotide identity over 1252 bp), the assignment of expressed sequence tags (ESTs) to any particular locus was problematic. Thus, EST data could not be used to definitively support ULBP gene annotation.

Comparative genome organization
Both the cattle minor and major ULBP clusters were localized to BTA9 using radiation hybrid mapping methods (data not shown). Comparative analysis showed that STXBP5 and SAMDC1 share a conserved orientation on HSA6q and BTA9 ( Figure 1). The cattle ULBP3, ULBP4, ULBP5 and ULBP6 loci located between STXBP5 and SAMDC1 likely originated by duplication and insertion of genes from the major ULBP cluster (see below).
The cattle NFYB gene is orthologous to human NFYB on HSA12 ( Figure 1). Comparison of nucleotide alignments in the cattle minor ULBP cluster and HSA12 genomic sequence demonstrates that sequence similarity is limited to the NFYB gene. The absence of genomic sequence similarity flanking this gene in humans and the lack of intronic sequence in cattle NFYB suggests that the cattle NFYB locus represents a retrotransposed gene. Although unlikely, a chimeric cattle BAC clone or sequence assembly error in the human genome cannot be ruled out as an explanation for these findings.
The discovery of at least 30 distinct cattle ULBP paralogues makes cattle the species with the largest number of ULBPlike genes identified to date ( Figure 1). Our findings confirm and extend previous Southern blot analysis indicating a large number of ULBP paralogues in cattle and seven other ruminant artiodactyl genomes [1]. In contrast, the more distantly related artiodactyls, swine and alpaca, appear to have relatively few ULBP genes [1,11].

The cattle ULBP loci evolved through extensive gene duplication
The cattle minor and major ULBP clusters were analyzed for internal nucleotide sequence similarity ( Figure 2 and Figure 3, respectively) in order to identify duplicated segments. The largest was a duplication of ULBP28, ULBP2, ULBP29, ULBP30, and ULBP31 to form ULBP16, ULBP17, ULBP18, ULBP19, and ULBP20 ( Figure 3). The directionality of this duplication event was determined from the expansion of a novel cattle-specific repeat (see below) in the first intron of cattle ULBP17 as compared to the smaller corresponding repeat region in the first intron of ULBP2. There appear to be four tandem duplications involving blocks containing ULBP9 and ULBP10, ULBP11 and ULBP12, ULBP13 and ULBP14, and ULBP15 ( Figure  3). However, similarity to ULBP27 was observed for ULBP9, ULBP11, ULBP13, and ULBP15, thus providing evidence that ULBP27 was likely also part of the large duplication involving ULBP28 through ULBP31 described above. In addition to the duplication events described above, there are two other segments that contain duplicated genes: i) ULBP7, ULBP8, and ULBP9 are related to ULBP22, ULBP23, and ULBP24, and ii) ULBP21 and ULBP22 are related to ULBP25 and ULBP26 ( Figure 3).
Known repetitive elements were identified in the minor and major ULBP clusters (Table 4). An additional novel genomic repeat was identified within the first introns of ULBP17 and ULBP2. The novel repeat spans 11,938 bp in the first intron of ULBP17 and 2,100 bp in the first intron of ULBP2 ( Figure 3) [GenBank: DP000081]. These repeats are specific to the cattle major ULBP cluster and are not found elsewhere in the cattle genome. The large size of the ULBP17 repeat region relative to the corresponding repeat region in ULBP2 suggests active repeat expansion. A full understanding of the means by which these repeats contributed to the evolution of the ULBP gene family awaits complete genomic sequencing of this region and sequencing of additional haplotypes.

Structure and evolution of ULBP proteins
Conceptual translations of the nine cattle ULBP genes identified in this study and the previously identified cattle ULBP1 [1] are shown in Figure 4. Each molecule contains a 24 to 42 amino acid (aa) signal peptide sequence, an 88 aa α1 domain, an 84 aa α2 domain and a 25 to 30 aa connecting peptide region followed by a hydrophobic segment. Peptide sequence identity was determined within the α1 and α2 domains for each cattle ULBP and the por-  cine ULBP ( Figure 5) [11]. ULBP9, ULBP21, and ULBP27 have glycophosphatidylinositol (GPI) anchor sites (P < 0.01, P < 0.05, and P < 0.001, respectively). The other seven ULBPs have predicted transmembrane domains of 23 to 25 aa followed by cytoplasmic tails ranging from 27 to 73 aa in length. The signal sequences and transmembrane domain or GPI anchor motifs indicate that all 10 of the expressed ULBPs are localized extracellularly. Each protein has predicted N-glycosylation motifs, with the exceptions of ULBP2 and ULBP17, suggesting that at least eight cattle ULBPs are glycoproteins.
Pairwise substitution analyses of ten ULBP genes showed an average global nonsynonymous to synonymous substitution ratio (ω t ) of 0.934 ( Figure 6). Values of ω > 1.0 are regarded as indicating that positive selection has operated on the sequences analyzed [12]; however, global substitution analysis is stringent and may mask evidence of positive selection in molecular subregions [13]. Heterogeneity in selection intensity was investigated within the ULBP α1 and α2 domain regions (Table 5). In model 2 (M2), a continuous positive selection model with an additional (third) ratio of nonsynonymous to synonymous substitutions (ω 2 ) estimated from the data, ω 2 is 3.17, but represented only a small proportion (p 2 = 0.08 out of 1.0) of codon sites. The log likelihood test of M2 vs. M1, the continuous neutrality model, was not statistically significant. In M3, the unconstrained discrete positive selection model, ω 2 is 1.90 with p 2 = 0.28. The log likelihood test of M3 vs M1 was significant, providing evidence of heterogeneity in ω ratios among codon sites. Model 8, a beta distribution with an added ω class estimated from the data, was compared to M7, a beta distribution that did not allow for positively selected sites. The log likelihood test of M8 vs M7 was significant, allowing the detection of positively selected codon sites (Table 5). Thirteen codon sites were determined to be under a high degree of positive selection (> 90% probability, Figure 4, Table 5). Twelve of the positively selected sites in the cattle ULBPs were mapped onto the three dimensional structure of human ULBP3 [PDB: 1KCG, chain C] (Figure 7). Eleven of the twelve positively selected residues were located at outwardly directed positions, indicating that positive selection acted at the level of interaction between the ULBPs and another molecule. On the basis of the structural data, fourteen human ULBP3 sites interact with Internal sequence identity plot of the minor ULBP cluster  Table 1.
Internal sequence identity plot of the major ULBP cluster Figure 3 Internal sequence identity plot of the major ULBP cluster. Sequence numbers displayed on the X-and Y-axes indicate alignment orientation, originating in the lower left corner. The center line represents identity; other lines indicate regions of internal sequence identity. Genes are annotated above the figure with arrows indicating orientation. Nucleotide sequences used to construct the major ULBP cluster are listed in Table 1. The duplication of ULBP28-ULBP31 to form ULBP16-ULBP20 is within the red shaded area. The cattle specific repeat regions within ULBP17 and ULBP2 are indicated by red asterisks. The duplication of ULBP7-ULBP9 corresponding to ULBP22-ULBP24 is within the blue shaded area. The duplication of ULBP21-ULBP22 corresponding to ULBP25-ULBP26 is within the green shaded area. The tandem duplications of ULBP9 and ULBP10, ULBP11 and ULBP12, ULBP13 and ULBP14, and ULBP15 are within the yellow shaded area, and lines showing their similarity to the ULBP27 region are within the violet shaded area.
NKG2D [14], and only one of these binding residues was found to overlap with the cattle ULBP sites under positive selection ( Figure 7). Therefore, the positively selected cattle ULBP sites, located outside of the predicted NKG2Dbinding residues on the basis of the homology modelling data, appear to interact with molecules other than NKG2D.
Several members of the Herpesviridae, which is the taxonomic family to which hCMV belongs, infect cattle, including bovine herpesviruses-1 through -5 and bovine lymphotrophic herpesvirus. The sequenced genomes of bovine herpesvirus-1, -4, and -5 do not encode molecules with detectable sequence similarity to UL16 of hCMV (HHV5); however it is conceivable that peptides encoded by the bovine herpesviruses or other viral pathogens may disrupt ULBP cell surface expression or the molecular interactions mediated by ULBP molecules. Thus, the rapid expansion of the ULBP gene family and the maintenance of such a large gene cluster are likely adaptive, serving to provide cattle with at least ten ULBP molecules through which an immune activation signal can be transmitted, even in the presence of an inhibitory pathogen-derived peptide.

Conclusion
This study provides insights into the genomic organization and evolution of the cattle ULBP genes, a recently  Sequence identity in the extracellular domains ofcattle and swine ULBPs Figure 5 Sequence identity in the extracellular domains ofcattle and swine ULBPs. Percent pair-wise sequence identity between the α1 and α2 domains of ten cattle ULBP proteins and porcine ULBP [GenBank: AAP81932]. Sequence alignments were made using BLASTP and edited manually.  Figure 4 Protein alignment of the cattle ULBP family. Alignment gaps are represented by periods, and aligned ULBP sites sharing identity with cattle ULBP1 are represented by dashes. The signal peptide region (sp) is represented by dashes within brackets above the alignment. The α1 and α2 domains are designated by an underscore above the alignment. Universally conserved cysteine, proline, and tryptophan residues are annotated with Cs, Ps and Ws, respectively, beneath the alignment. Predicted Nglycosylation motifs are represented by underlined text in the alignment. Transmembrane domains are represented by black shaded background. Predicted GPI anchor motifs are represented by gray shaded background, and the associated downstream hydrophobic regions are indicated by italicized text with gray shaded background. Positively selected sites are indicated by asterisks above the alignment, and posterior probabilities associated with each positively selected site (see Table 2.5) are represented within the alignment by normal text (probability > 0.90), bold text (probability > 0.95), and italicized bold text (probability > 0.99).  , and centrifuged at 10,000 × g for 10 min. An 860 µl aliquot of cleared lysate was combined with 600 µl isopropanol, placed on ice for 15 min, and centrifuged at 10,000 × g for 5 min. The supernatant was decanted, and the pellet was washed with 500 µl 70% ethanol before Global substitution analysis of ten ULBP paralogues Figure 6 Global substitution analysis of ten ULBP paralogues. The bold numerical values correspond to ω t ratios; raw d N /d S values are listed in parentheses. Accession numbers associated with each sequence analyzed are listed in the Methods. centrifugation at 10,000 × g for 5 min. The dried pellet was suspended in 40 µl of a solution containing 10 mM Tris and 1 mM EDTA.

BAC selection, isolation and sequencing
The BAC DNA was digested using the HindIII restriction enzyme, separated by electrophoresis on 1X TAE agarose, stained using SYBR Green (Invitrogen), and visualized using image analysis (Typhoon Visual Imaging System, Molecular Dynamics) according to an established protocol [15]. Gel images showing restriction fragments were analyzed semiautomatically using IMAGE v3.10b [16], and band migration information was analyzed using FPC v6.0 [17,18] to determine clone overlap for contig assembly as previously described [15].
The first round of sequencing was performed for four minimally overlapping ULBP-containing BACs [GenBank: AC098687, AC092858, AC096629 and AC098686]. BAC DNA isolation, shotgun cloning, sequencing, quality analysis and sequence assembly were performed using established protocols [19,20]. The sequenced BACs were aligned using BLASTN v2.2.13 [21] to generate one contiguous genomic DNA sequence containing three BACs and one singleton. Repetitive elements in the contiguous sequences were masked using REPEATMASKER v3. 1.3 [22] before the sequences were used to query publicly available cattle genome trace sequences [NCBI Build 2.0] to identify additional minimally overlapping cattle BACs. Two additional ULBP-containing BAC clones were identified [GenBank: DQ405273 and DQ405274], and a second round of sequencing was performed. Shotgun cloning and sequencing was performed using the Topo Shotgun Subcloning Kit (Invitrogen) according to the manufacturer's instructions. PHRED [23,24] was used to remove low quality sequence (PHRED score < 20). CROSSMATCH and PHRAP [25] were used to remove vector sequence and assemble BAC subclone sequences, respectively. For the two BACs sequenced in the second round, sequence gaps were closed by primer walking. Contiguous sequences [GenBank: DP000081 and DP000082] were constructed from overlapping full-clone BAC sequences using BLASTN. REPEATFINDER [26] was used to identify genomic sequence repeats not identified by REPEATMASKER. [GenBank: AY160681] sequences were aligned to the genomic sequences using BLASTN to assist in the annotation of the ULBP genes. Exon/intron boundaries were verified by manual inspection and editing. For each locus identified, all exons were joined and conceptually translated using SIXFRAME [27] to identify open reading frames. Homologous positions in the human genome [NCBI Build 35] were identified using the UCSC Genome Browser [28].

Substitution analysis
Ratios of nonsynonymous to synonymous substitutions (d N /d S or ω) were determined for the ULBP genes using the PAML software package [37]  showing the spatial arrangement of homologous cattle ULBP residues under positive selection (>90% probability) as well as the human ULBP3 residues that interact with the NKG2D molecule. The human ULBP3 backbone appears in blue and green. Twelve of 13 cattle ULBP sites under positive selection are mapped onto the structure, and eleven appear as red space-filling residues. Thirteen of the fourteen human ULBP3 sites that interact with NKG2D appear as yellow space-filling residues. One site corresponding to both a selected cattle ULBP site and a human ULBP3 site that interacts with NKG2D appears in orange.
ULBP9, ULBP11, ULBP13, ULBP15, ULBP17, ULBP21, ULBP27 and ULBP2 [annotated in GenBank: DP000081]. Only the extracellular α1 and α2 domain regions were analyzed. The YN00 program in PAML was used to estimate ω t for each group of aligned sequences using the method of Yang and Nielsen [38]. The CODEML program in PAML was used to identify variation in selection intensity. The data were modelled using maximum likelihood methods [39], and the results were compared to obtain a test statistic. Three comparisons were performed. Model M1, a neutrality model that constrained ω to be either 0 or 1, was compared to both M2, a selection model that added an additional ω ratio class estimated from the data, and M3, a selection model that used an unconstrained discrete distribution to model classes of ω ratios. This analysis used three discrete classes for M3. In addition, M7, a continuous distribution neutrality model that estimates ω using a beta function limited to the interval from 0 to 1, was compared to M8, a continuous distribution selection model that adds an additional class of sites with ω estimated from the data and not constrained to the interval between 0 and 1. A test statistic of twice the negative value of the difference between the log likelihood values generated under each model was compared to a χ 2 distribution with degrees of freedom calculated from the difference in the number of model parameters (M2 vs M1, df = 2; M3 vs M1, df = 4; M8 vs M7, df = 2). Posterior probabilities for ULBP sites under positive selection were generated under M8.