text:
Differential regulation of abundance and deadenylation of maternal transcripts during bovine oocyte maturation in vitro and in vivo

Background In bovine maturing oocytes and cleavage stage embryos, gene expression is mostly controlled at the post-transcriptional level, through degradation and deadenylation/polyadenylation. We have investigated how post transcriptional control of maternal transcripts was affected during in vitro and in vivo maturation, as a model of differential developmental competence. Results Using real time PCR, we have analyzed variation of maternal transcripts, in terms of abundance and polyadenylation, during in vitro or in vivo oocyte maturation and in vitro embryo development. Four genes are characterized here for the first time in bovine: ring finger protein 18 (RNF18) and breast cancer anti-estrogen resistance 4 (BCAR4), whose oocyte preferential expression was not previously reported in any species, as well as Maternal embryonic leucine zipper kinase (MELK) and STELLA. We included three known oocyte marker genes (Maternal antigen that embryos require (MATER), Zygote arrest 1 (ZAR1), NACHT, leucine rich repeat and PYD containing 9 (NALP9)). In addition, we selected transcripts previously identified as differentially regulated during maturation, peroxiredoxin 1 and 2 (PRDX1, PRDX2), inhibitor of DNA binding 2 and 3 (ID2, ID3), cyclin B1 (CCNB1), cell division cycle 2 (CDC2), as well as Aurora A (AURKA). Most transcripts underwent a moderate degradation during maturation. But they displayed sharply contrasted deadenylation patterns that account for variations observed previously by DNA array and correlated with the presence of a putative cytoplasmic polyadenylation element in their 3' untranslated region. Similar variations in abundance and polyadenylation status were observed during in vitro maturation or in vivo maturation, except for PRDX1, that appears as a marker of in vivo maturation. Throughout in vitro development, oocyte restricted transcripts were progressively degraded until the morula stage, except for MELK ; and the corresponding genes remained silent after major embryonic genome activation. Conclusion Altogether, our data emphasize the extent of post-transcriptional regulation during oocyte maturation. They do not evidence a general alteration of this phenomenon after in vitro maturation as compared to in vivo maturation, but indicate that some individual messenger RNA can be affected.

Background
The fully grown oocyte (FGO) ability to sustain genome reprogramming and embryo development is associated with a specific and complex pattern of gene expression. This has prompted high-throughput studies of the mammalian oocyte transcriptome in mouse [1,2], in human [3] or bovine [4]. In most tissues, the transcriptome reflects the contemporary pattern of RNA synthesis and protein synthesis. In this respect, the fully grown and maturing oocytes are atypical models, since at these stages the oocyte is essentially transcriptionally quiescent. Rather, the FGO transcriptome results from active RNA synthesis during the oocyte growth phase, followed by post-transcriptional regulation. A subset of maternal transcripts undergoes post-transcriptional regulation through shortening of the polyA tail, and is stored protected from degradation and from the translation machinery. At the time when the protein is required, during maturation, fertilization or the early stages of embryo development, such dormant transcripts can undergo cytoplasmic polyadenylation and support translation, then are rapidly degraded. The underlying molecular mechanisms were shown to involve cis regulatory elements within the 3' untranslated region, including cytoplasmic polyadenylation elements (CPE) (for reviews, see [5][6][7]). However, during maturation most maternal transcripts are expected to follow a default deadenylation and degradation pathways [8][9][10].
Besides a complex post-transcriptional regulation pathway, the mammalian oocyte transcriptome is characterized by specific transcripts not found in somatic tissues, as first evidenced in mouse [11][12][13]. The hypothesis that socalled oocyte-specific genes would be important for fertility was confirmed by functional genomics experiments in mouse (reviewed in [12]). A first class of genes is required for normal folliculogenesis and fertilization. Other genes are maternal effect genes, i.e. they are dispensable until fertilization but they are required for proper embryo development: most embryos from mothers null for Zygote arrest 1 (Zar1), Maternal antigen that embryos require (Mater, official symbol Nalp5 as a member of the NACHT, leucine rich repeat and PYD containing, (NALP) family), or Stella (official name Developmental pluripotency associated 3 (Dppa3)) are blocked before the blastocyst stage [14][15][16][17]. Functional importance of other oocyteenriched genes remains to be assessed, such as Oocyte secreted protein 1 (Oosp1), Maternal embryonic leucine zipper kinase (Melk), the Oogenesin genes and several Nalp family members [13,[18][19][20][21] (Studies in bovine have identified orthologous genes, as well as potentially novel oocyte-specific genes in this species [22][23][24][25][26][27]. The main objective of this study was to investigate posttranscriptional regulation of maternal transcripts during bovine oocyte maturation in relation with its quality, i.e. its ability to undergo meiosis and sustain embryo development. Bovine is definitely a pertinent scientific model to analyze the relationship between oocyte developmental competence and the regulation of maternal transcripts. First, embryonic genome major activation is delayed up to the 8/16 cell stage about 3.5 days after fertilization, whereas in mouse it is already completed 24 hrs after fertilization at the end of the 2-cell stage. This chronology highlights the importance of timely control of maternal RNA. Most importantly, several physiological and experimental models of oocyte quality have been defined in this species (reviewed in [28]). In particular, in vivo matured oocytes present the highest developmental competence (as estimated by the blastocyst rate), while in vitro maturation remains a suboptimal alternative [29][30][31][32]. As nuclear maturation is not affected, the difference is usually attributed to an inadequate reproduction of the cytoplasmic features of oocyte maturation in vitro.
Here, using real time PCR, we have analyzed the variation of maternal transcripts during bovine oocyte maturation and embryo development, in terms of abundance and polyadenylation. We selected genes likely involved in oocyte quality, based on their restricted expression in the germ cell or their previous isolation as markers of maturation [33]. Three oocyte specific genes (MATER, NALP9 and ZAR1) were previously identified. Four additional genes are characterized here for the first time: the bovine MELK and STELLA genes, as well as two genes whose oocyte preferential expression was not previously reported in any species, ring finger protein 18 (RNF18) and breast cancer anti-estrogen resistance 4 (BCAR4). We also examined the cell-cycle regulator Aurora A (AURKA) whose transcript is abundant in meiotically arrested oocytes. In addition, we followed 6 transcripts previously identified as down-or up-regulated during oocyte maturation: cyclinB1 (CCNB1) and cell division cycle 2 (CDC2) encoding components of the M-phase promoting factor, as well as peroxiredoxin (PRDX) 1 and 2, and inhibitor of DNA binding (ID) 2 and 3 [33]. Overall, our results show a moderate degradation of most transcripts during maturation, and evidence contrasted deadenylation patterns that mostly account for variations observed previously by DNA array. As a differential model of developmental competence, we compared regulation of maternal transcripts during in vitro and in vivo maturation. During embryonic development, distinct profiles were observed, and most oocyte marker genes remained silent after embryonic genome activation.

Characterization of four bovine genes preferentially expressed in the oocyte
We characterized transcripts represented by five distinct EST in our library of EST preferentially represented in bovine oocyte as compared to somatic tissues [23]. One EST presented a high sequence homology with the human MELK gene, while four others did not return a match. Virtual Northern blots were run and hybridized using the considered EST as probes. Then, full length transcript sequence was determined by Rapid Amplification of cDNA Ends (RACE). Gene structure and chromosomal localization was deduced by aligning the full length cDNA and the bovine genome sequences. We were able to identify the three remaining genes as STELLA, RNF18 and BCAR4.

MELK
One EST was 88% homologous to human MELK. As observed by RT-PCR, MELK is not really oocyte specific; rather, it is expressed at low level in most tissues, but strongly expressed in both the male and female gonads and in the oocyte (Fig 1A). A virtual Northern blot revealed a single 2.5 kbp cDNA (not shown). The full length 2445 bp sequence deduced from RACE (genbank accession hoEF446902) is assembled from 18 exons, including a first non coding exon that is lacking in the bovine MELK transcript that has since been predicted within the genome onto BTA8. The full-length sequence contains a 1953 nucleotide long open reading frame (ORF), encoding a predicted 650 aa/74 kD protein wse sequence is 90% and 80% identical to its human and mouse counterparts, respectively; it presents similar serine/threonine kinase domain and kinase associated domain at the N-terminus and C-terminus respectively.

STELLA
Two independent EST initially could not be identified based on their sequence, including previously described EST5B4 [23]. By virtual Northern blot we observed a doublet: the corresponding 1590 and 1347 bp long cDNA, as deduced from RACE, are generated by alternative splicing at their 5' end ( Fig 1B) (genbank accession EF446904, EF446905). They present the same 492 nt ORF encoding a 163 aa putative protein. We identified a bovine orthologue of the germ cell marker Stella, despite limited sequence homology with the human and mouse counterparts (30% and 26% respectively). By RT-PCR, both variants were undetected in somatic tissues. Variant 1 was detected in testis, ovary and oocyte, while variant 2 was detected only in oocytes ( Fig 1A).

RNF18
Initially our EST could not be identified by sequence homology with sequences in public databases. By RT-PCR, the transcript was specifically detected in testis, ovary and oocyte (Fig. 1A). A 1723 bp cDNA, observed by virtual Northern blot, was cloned by RACE (genbank accession EF446903). The 1353 nt long ORF encodes a 450 aa/52 kD putative protein with 3 major regions: a Nterminal ring finger, a B-box, and a C-terminal SPRY domain also found in ryanodine receptor. This structure identifies the protein as a member of the tripartite motif (TRIM) family. The corresponding locus LOC525113 on BTA29 has since been predicted from the bovine genome with the annotation "similar to ring finger protein 18 (Testis-specific ring-finger protein)", whose human orthologue was found primarily expressed in testis [34] and is also named TRIM49. Identity is confirmed by analogous gene structure, with six coding exons of similar length in human and bovine. As compared to the predicted bovine transcript, our cDNA contains three additional exons,

BCAR4
Two transcripts (genbank accession EF446906, EF446907) observed by virtual Northern blot were cloned by RACE starting from our unidentifiable EST. They are 638 and 1317 bp long respectively, and are generated through the alternative use of two polyadenylation signals, downstream a putative cytoplasmic polyadenylation element (Fig. 1C). The 363 nt long ORF encodes a putative 120 aa/14 kD protein. The bovine gene is located onto BTA25, between two loci annotated "similar to zinc finger CCCH-type domain containing 7" and "ribosomal L1 domain containing 1"respectively. By conserved synteny we identified the human ortholog breast cancer antiestrogen resistance 4 (BCAR4). Both the bovine and human ORF are included within a single exon, and encode proteins with two putative transmembrane domains; yet they do not display significant sequence identity. In human, ESTs have been isolated mainly from placenta, and to our knowledge expression of BCAR4 has not been reported in the gonads or the oocyte. By RT-PCR the transcript was detected specifically in the ovary at trace level and in the oocyte (Fig. 1A).

Variation of maternal transcripts during in vitro maturation
We followed the variation of selected maternal transcripts during in vitro oocyte maturation, in terms of abundance and polyadenylation. Real-time PCR was run onto random hexamer-primed RT products, in order to estimate their abundance, and oligo(dT) 15 -primed RT products, which are affected by both abundance and polyadenylation. For each transcript, the ratio to the mean value in immature oocytes is represented (Fig 2). For most transcripts, including 18S ribosomal RNA, the abundance did not vary significant during maturation, although an average 20% decrease was observed. PRDX1, ID2 and MATER transcripts appeared the most affected (over 33% degraded). As expected due to transcription repression during maturation, no transcript displayed a significant increase. The patterns of polyadenylated transcripts revealed sharp contrasts. CDC2, PRDX1, PRDX2, ID2, MATER, ZAR1 transcripts were strongly affected, as detection level dropped below one third of their initial level. For all these genes, variations of the total and polyadenylated forms were significantly different during maturation ( Fig. 2, *). For NALP9 and BCAR4 variant 1, the decrease (41 and 49% respectively) was still more important than the decline observed for the total form (27 and 23%), but the difference was not statistically significant.
On the other hand, seven genes appeared hardly affected; the polyadenylated transcript followed a parallel evolution as the total form, and decreased by less than 30%.

Comparison of in vitro maturation and in vivo maturation
In order to compare the variation of maternal transcripts during in vitro and in vivo maturation, a subset of genes was followed in oocytes collected by ovum pick up (OPU) (Fig. 2). BCAR4 variant 1, STELLA variant 2, ID2 and ID3 transcripts were not quantified in these minute samples due to requiring more substrate. Again, for most transcripts the total form displayed a limited and statistically non significant variation during maturation. An exception was PRDX1, whose detection level dropped by 62% during in vivo maturation: this was statistically different from both the level in immature oocytes and the 33% drop during in vitro maturation ( Fig. 2, ≠). As observed during in vitro maturation, polyadenylated transcripts displayed more contrasted patterns. ACTB, PRDX1, PRDX2, CDC2, MATER, NALP9, ZAR1 decreased significantly (Fig. 2, *). Indeed, for polyadenylated transcripts, the mean values were often remarkably close after in vitro and in vivo maturation. Interestingly, PRDX1 polyadenylated form displayed a similar variation in both cases.

Profile of maternal transcripts during preimplantation development
Then we analyzed the variation of oocyte markers during the early stages of embryo development (Fig. 3). As examples of housekeeping genes, we quantified 18S and ACTB transcripts. Their levels started to increase between the 5/ 8-cell and morula stages, as expected since major genome activation occurs in bovine at the 8/16-cell stage, and rose further during blastocyst formation and hatching. By contrast, detection of all oocyte markers decreased between the 5/8-cell and morula stages, and all but MELK decreased further in blastocysts, often below detection level. MELK was indeed activated in blastocysts.
Finally, we did not observe significant differences between the total and polyadenylated forms of any transcript during embryo development from zygote to hatched blastocysts.

Regulatory elements in the 3' untranslated region
We searched the 3' untranslated region of the transcripts for some of the regulatory sequences that are known to be involved in the post-transcriptional regulation of maternal RNA in model species. There is no unique consensus sequence for the CPE, but rather variants from the basic U 4-5 A 1-2 U sequence were identified in xenopus and mouse within 100 nucleotides upstream the nuclear polyadenylation signal [5,9,35] Based on this definition, potential CPE-like sequences were found in ID2, ID3, CCNB1, MELK, RNF18, STELLA and AURKA (Table 1). Except for ID2, their total and polyadenylated forms were not differentially regulated during in vitro or in vivo maturation (Fig. 2, no *). Reciprocally, no such CPE could be recognized within the 3' untranslated region of CDC2, PRDX1, PRDX2, MATER, NALP9 and ZAR1 transcripts, or much farther from the nuclear polyadenylation signal for  * * * * ≠ ACTB; their polyadenylated forms all displayed a significant decline as compared to their total form during in vivo maturation (Fig. 2, *), as well as in vitro maturation although below significance threshold for NALP9 and ACTB. CPE-like sequences were present in BCAR4 transcripts, either close to the nuclear polyadenylation signal for variant 2, or 600 nucleotides upstream for variant 1 ( Table 1). BCAR4 variant 1 exibited a similar profile as NALP9 during in vitro maturation. When both variants were analysed together, neither the total nor the polyadenylated transcripts appeared much affected (Fig. 2). Of note, STELLA and ACTB contained U 10 CU 8 and U 12 , U 15 and U 11 sequences respectively. These might be embryonic CPE, i.e. elements that target transcripts for cytoplasmic polyadenylation after fertilization [6]. Yet no such polyadenylation was evidenced by real-time PCR.  We also searched the 3' untranslated region for so called A/U rich element (ARE) that can target corresponding transcripts for active deadenylation, as opposed to the default deadenylation pathway (Table 1). Although there is no strict consensus sequence, AREs have tentatively been grouped into three major classes [36,37] Class I and II ARE contain one or several dispersed or adjacent copies of AUUUA motifs in a U-rich context; such motifs were found in bovine ACTB, CDC2, ID2, ID3, MELK, RNF18, STELLA and CCNB1 (Table 1). Class III ARE are described as over 50 nucleotide long U-rich regions; several transcripts may fall into this category. Finally, the presence of multiple UGU motif in a U/G-rich context has also been reported to target mRNA for deadenylation [38]; several bovine transcripts indeed contained at least 4 copies of the triplet within 45 nucleotides. But altogether, no obvious correlation appeared between the presence of putative ARE motifs or UGU-rich sequences and the transcripts profile during maturation or development.

Discussion
In this study, we have characterized four novel bovine genes with a preferential oocyte expression. Then, using real-time PCR, we have analysed these and other maternal transcripts for changes in their abundance and polyadenylation during in vitro or in vivo maturation as well as in vitro development up until the hatched blastocyst stage.

Influence of polyA tail onto transcripts detection by real time PCR and microarrays
A few studies have commented on the different apparent transcript levels deduced from real time PCR onto cDNA synthesized from oligo(dT) or random primers in bovine oocytes and embryos [39][40][41] Here we extend this obser-vation to a number of oocyte-specific or non specific transcripts. The discrepancy is attributed to variations in the polyA tail, which is of major interest in this model due to the predominantly post-transcriptional regulation of gene expression in FGO and pre-MET embryos. We may also suspect a similar bias when polyadenylated RNA are purified [42][43][44][45] Indeed, the influence of polyA tail has consequences beyond the analysis of real time PCR data, into the field of embryogenomics. Due to limited biological samples, amplified probes are usually generated from mammalian oocytes/embryos for hybridizing onto microarray, through protocols that involve reverse transcription from oligo(dT). Therefore, a similar bias associated with the polyA tail is expected in these transcriptomic analyses. Although this had been suggested before, experimental validation analyzing both the total and polyadenylated transcripts was lacking. In this study we have included six genes previously identified by DNA array as markers of oocyte maturation: CDC2, PRDX1 and PRDX2 transcripts were shown to be negatively regulated, while CCNB1, ID2 and ID3 were positively regulated [33]. Here, using real time PCR, we observed only a minor degradation of five transcripts (all but ID2) during maturation. But we did observe a dramatic decrease in polyadenylated transcripts level for CDC2, PRDX1 and PRDX2, while ID3 and CCNB1 polyadenylated transcripts appeared hardly affected. These results indeed confirm that the array experiment was influenced by the polyadenylation status of these transcripts. Similarly, down regulation of MATER during oocyte maturation, as observed on a microarray [46], must reflect the significant decrease of the polyadenylated transcript evidenced here. In our experiment, ID2 did not follow the expected pattern: while among up-regulated transcripts revealed by DNA array screening, the none none polyadenylated transcript displayed a significant decline when assessed by real time PCR. We have currently no satisfying explanation for this discrepancy. Hybridization of a non-specific probe due to cross-species hybridization appears unlikely due to the high sequence conservation (94%) between the bovine and human coding sequences. Alternatively, a different transcript variant might have hybridized onto the array than the one studied by real time PCR.
Altogether, our results should prompt to reconsider interpretation of DNA array and real-time PCR data in models with potential regulation of the polyA tail. In mouse and bovine oocytes and embryos, array validation experiments were based on reverse-transcription initiated either from oligo(dT) [44,45,47] or from a mix of both oligo(dT) and random hexamers [48]. In order to minimize the impact of variable polyA tails, a modified protocol for probe amplification was recently developed [49].
On the other hand, detecting transcripts polyadenylation status on an array may be viewed as an asset, since it provides a footprint of mRNA that can be translated at a given stage, circumventing the need for purifying polyadenylated transcripts [50].

Selective deadenylation of maternal transcripts during in vitro maturation
Maturation-associated degradation of stored maternal transcripts has been known for a long time. In the mouse, earlier studies reported a 19% decrease in RNA content, as measured by ethidium bromide fluorescence; most of it was accounted for by the degradation of rRNA [51,52] Here, bovine 18S ribosomal RNA was found to decrease during in vitro maturation by 22% (not statistically significant), a moderate degradation that could not be observed previously by Northern blot [39]. Our value appears fairly close to the value reported in the mouse. Yet, based on different developmental kinetics in both species, one might have expected a lesser extent of degradation in the bovine, as maternal rRNA will be required to sustain translation in embryo for several days versus only a few hours in mouse. A higher absolute RNA content in bovine oocytes as compared to mouse oocytes [53] might compensate for this similar global degradation rate.
Non ribosomal RNA is also affected, but the destruction of transcripts during maturation in mice appears a selective rather than promiscuous process: this was observed on individual genes in earlier studies [54] and was recently confirmed using a microarray [49]. Here, most transcripts displayed a limited degradation during maturation. Within our set of genes, MATER, PRDX1 and ID2 transcripts were most affected, over 33% degraded. Degradation of mouse Mater transcript in ovulated oocytes was also recently reported [49]. Contrasted profiles were observed for the polyadenylated forms. In this respect, oocyte specific transcripts did not exhibit a uniform or a particular pattern; a similar heterogeneity was observed for housekeeping genes [40]. For a few genes, protein data are available. During in vitro maturation, the level of CDC2 and MATER proteins was not much altered as observed by Western blot [55,56], while an isoform of PRDX2 was shown to be synthesized through a proteomics approach [57]. Here, all three polyadenylated transcripts displayed a major decline. This apparent discrepancy may only reflect complex post-transcriptional and translational control, and illustrates the uncoupling of transcriptome and proteome in FGO and mature oocytes.

Comparison of in vitro and in vivo maturation
Whether considering variation of total or polyadenylated transcripts, for most genes we did not evidence a different pattern after in vitro vs in vivo maturation. Indeed, for polyadenylated transcripts mean values were often remarkably close. Overall, our data support the notion that degradation and deadenylation are not severely affected by in vitro maturation as compared to in vivo maturation. Lower developmental competence of in vitro matured oocytes should not be attributed to a general alteration of post-transcriptional regulation of maternal transcripts. Yet the following points have to be kept in mind. First, we considered only relative variations between immature and mature oocytes; should a transcript be more abundant at both stages in one model (in vitro or in vivo maturation), the corresponding ratio might not be affected. Second, we collected oocytes before and after maturation, but not in the course of maturation, when transient differences might have been revealed. Third, it can not be excluded that superovulation treatment affects oocyte transcriptome, so that OPU oocytes may not reflect normal physiology.
An intriguing exception was PRDX1, whose total form declined more severely during in vivo as compared to in vitro maturation. This transcript was reported to increase in oocytes from growing follicles over 8 mm in diameter [58]; however, direct comparison with our own results is difficult since in this paper, there is no mention of whether transcript total form or only polyadenylated form was quantified, and data were normalized to an internal rather than an exogenous reference. In our study, immature oocytes were meant to be collected both in vivo and post-mortem from follicles under 8 mm, and thus the follicle size should not impact our data. Altogether, these studies may indicate that PRDX1 level in oocytes is particularly sensitive to environmental changes during late folliculogenesis and maturation. Analysis in other physiological models will reveal whether it is a general marker of bovine oocyte developmental competence.
In a previous study, a correlation between shorter polyA tail and low developmental competence had been reported; the difference in polyA tail could be observed both before and after in vitro maturation [59]. The authors considered a different model of oocyte quality, based on the morphology of the ovary: their competent oocytes, fairly similar to our post mortem collected oocytes, were compared with oocytes with really poor developmental (5% blastocyst rate) [60].

Profile of maternal transcripts during embryo development
Profiles during embryo development reflected the degradation of maternal transcripts between the immature oocyte and the 5/8-cell stage. In morulae and blastocysts, detection of 18S and ACTB increased, reflecting embryonic genome activation. Strikingly, oocyte markers continued to decline in morulae, often below detection level. Thus they were not activated at the time of major genome activation, as previously reported for MATER, OOSP1, BMP15, GDF9, ZAR1 and NALP9 [22,24,27,56] and for the vast majority of mouse oocyte specific genes [2]. All but MELK remained silent in expanded and hatched blastocysts. Interestingly, Melk is also activated in mouse 8cell embryos [2,19] Most transcripts displayed a rapid decline during the earlier stages of embryo production, i.e. oocyte maturation, fertilization and zygote formation. BCAR4 was fairly constant during maturation, and then dropped concomitantly with zygote formation. Whether this decline follows translation of a protein required for early embryogenesis remains to be investigated. By contrast, MELK, RNF18, STELLA variant 1 and AURKA exhibited a progressive decrease from immature oocytes to 5/8-cell embryos. Again, it is tempting to speculate that these transcripts support translation at all stages to ensure continuous production of a short half-lived protein necessary for development; but they may just as likely be degraded without supporting translation.
For all transcripts, the total and polyadenylated forms displayed parallel profiles throughout embryo development from zygote to hatched blastocyst. We did not evidence polyadenylation of a transcript at a specific stage when it would be recruited for translation. However, this may be too transient a phenomenon to be observed in a population of several embryos. Rather, individual analysis of numerous embryos might be required to reveal the few embryos containing a transcript with an elongated polyA tail.

CPE and the control of maternal transcripts during maturation
Globally, we observed a correlation between the presence of a CPE-like sequence and the regulation of polyade-nylated transcripts during maturation. The profiles of transcripts without a CPE were coherent with them following the expected default deadenylation pathway. Reciprocally, the polyadenylated forms of transcripts with a putative CPE appeared stable, protected from further deadenylation, with the exception of ID2. The divergent profile of ID2 may indicate that a CPE is necessary, but not sufficient for stabilizing the polyA tail. Alternatively, it may have been generated by a distinct transcript variant with a longer 3' untranslated region (and thus a far upstream CPE, see BCAR4).
Altogether, our results are in agreement with conclusion from a recent microarray study in Xenopus tropicalis: CPE were under-represented in transcripts that were deadenylated during maturation, and over-represented in transcripts that were polyadenylated during maturation or substrate for post-fertilization deadenylation [61]. Here we did not observe a significant increase that would result from cytoplasmic polyadenylation. This is not unexpected as the polyadenylated form may be transient and no longer present in mature oocytes; a time-course analysis might be necessary, especially during the first 6 to 10 hours when cytoplasmic polyadenylation is believed to initiate [62,63] Alternatively, the elongation of the polyA tail may not be sufficient to be detected by this methodology; in fact, in a previous study polyadenylation of one CCNB1 variant was observed using a PCR based polyA test, but not by real-time PCR [64].

Conclusion
In this study, we have analyzed the regulation of selected maternal transcripts during bovine oocyte maturation, and its correlation with the presence of a putative CPE. As a model of differential developmental competence, we have compared maturation in vitro and maturation in vivo. The control of degradation and deadenylation did not appear globally compromised by in vitro manipulation, as only PRDX1 exhibited strikingly different patterns in the two conditions.

Post mortem oocyte collection and embryo production
Ovaries from adult cows were collected at a slaughterhouse and the oocyte-cumulus complexes were aspirated from 3-8 mm follicles, selected based on morphological criteria and washed in saline solution (for details see [65]). Denuded germinal vesicle stage immature oocytes and in vitro matured metaphase II oocytes (polar body extrusion was verified) were obtained as previously described [33]. Briefly, oocyte-cumulus complexes were denuded by mechanical treatment either before or after in vitro maturation in TCM199 (Sigma, Saint Quentin Fallavier, France) supplemented with 10 ng/ml epidermal growth factor and 10% fetal calf serum for 22 h at 39°C in water-saturated air with 5% carbon dioxide. Embryos were produced as previously described [56]. Groups of 10 embryos were collected over preimplantation development: zygotes and 2-cell embryos (day 2), 4-cell, 5 to 8cell embryos (day 3), morulae (day 5), expanded blastocysts (day 7) and hatched blastocysts (day 8). All oocyte and embryo samples were stored frozen in RNA later (Ambion, UK).

In vivo oocyte collection
All procedures were approved by the Agricultural and Scientific Research Government Committees in accordance with the guidelines for Care and Use of Agricultural Animals in Agricultural Research and Teaching (approval PB380). Oocytes were collected in vivo from eight superovulated Montbeliard cows by OPU at 12 or 60 h after prostaglandin injection to collect immature and in vivomatured oocytes, respectively. Oocytes were denuded as described above, and stored frozen in RNAlater.

Characterization of full-length cDNA
Gonads were collected from animals slaughtered in a commercial abattoir. Other organs were collected from an animal bred and killed in an experimental setting; all procedures were approved by the Agricultural and Scientific Research Government Committees in accordance with the guidelines for Care and Use of Agricultural Animals in Agricultural Research and Teaching (approval A37801). RNA was isolated from biopsies using Trizol (Invitrogen, Cergy Pontoise, France) and DNAse-treated. RACE-PCR and virtual Northern blot (i.e. Southern blot of full length cDNA) were carried out as previously described [22]. Sequences are available from Genbank and from the Ovogenae programme database [66]. Primers sequences are given in Table 2.

Reverse-transcription for real-time PCR
After adding 1 pg of luciferase mRNA (Promega, Charbonnières-les-bains, France) per oocyte/embryo as an exogenous standard, total DNAse-treated RNA was purified from 4 independent pools of 10 oocytes or embryos at each stage (or 5 to 16 OPU oocytes) using the Picopure RNA isolation kit (Alphelys, Plaisir, France). Reverse transcription was performed at 37°C for 50 min using oligo(dT) 15 primers (Promega) or random hexamers (Promega) by mouse Moloney leukaemia virus reverse transcriptase (Invitrogen) onto RNA equivalent to 2.5 in vivo collected oocytes and 5 post-mortem collected oocytes or in vitro produced embryos.

Real time PCR
Target cDNA were quantified by real-time PCR using iQ SYBR green supermix (Bio-Rad, Marnes la Coquette, France) with a iCycler or MyiQCycler system (Bio-Rad) using specific primers (sequences are presented in table   2). Primers for CCNB1 did not discriminate variants present in bovine oocytes [64]. A 3-step protocol (95°C for 30 sec, 60°C for 30 sec, 72°C for 20 sec) was repeated for 40 to 50 cycles, followed by acquisition of the melting curve. Amplicons were sequenced after cloning into pCRII vector using TA cloning kit (Invitrogen). The standard curve was deduced from serial dilutions of a plasmid including the target sequence, adjusted so that samples would fall within the considered range. Depending on the target gene, a cDNA amount equivalent to 0.0005 to 0.2 oocyte or embryo was used in triplicate PCR reactions, and the median value was considered (0 when not detected). For each sample, the data were normalized to the median value for exogenous luciferase. The gene/luciferase value in mature oocytes was then compared to this same ratio in immature oocytes (mean of four oocyte or embryo collections); data are presented as mean ± SEM. For each gene, Mann-Whitney comparison test was used to compare a) detection in mature vs immature oocytes, b) data obtained with random hexamers primed vs oligo(dT) 15 primed RT products, c) regulation in in vitro matured vs in vivo matured oocytes. For embryonic development from mature oocytes onwards, after analysis of variance, differences were evaluated using Tukey comparison test. In all instances, differences were considered statistically significant at the 95% confidence level (P < 0.05).

Authors' contributions
All authors were involved in biological material collection. AT participated in experimental design, data analysis, figure editing and carried out most experimental work: oocyte in vitro maturation, embryo production, RNA extraction, real-time PCR for oocyte marker genes, statistical analysis and characterization of oocyte markers MELK and RNF18. PP carried out real-time PCR experiments onto genes previously identified as oocyte maturation markers. CP carried out in vitro maturation and embryo production. SP characterized STELLA and BCAR4. SU and JMT supervised and participated in the study of AURKA. PH and CJ supervised and carried out animal ovarian stimulation and OPU. RDT designed and supervised the study, analysed data and wrote the manuscript. All authors have read and approved the final manuscript.


extfieldsofstudy:
Computer Science
Biology
Medicine

text:
EPGD: a comprehensive web resource for integrating and displaying eukaryotic paralog/paralogon information

Gene duplication is common in all three domains of life, especially in eukaryotic genomes. The duplicates provide new material for the action of evolutionary forces such as selection or genetic drift. Here we describe a sophisticated procedure to extract duplicated genes (paralogs) from 26 available eukaryotic genomes, to pre-calculate several evolutionary indexes (evolutionary rate, synonymous distance/clock, transition redundant exchange clock, etc.) based on the paralog family, and to identify block or segmental duplications (paralogons). We also constructed an internet-accessible Eukaryotic Paralog Group Database (EPGD; http://epgd.biosino.org/EPGD/). The database is gene-centered and organized by paralog family. It focuses on paralogs and evolutionary duplication events. The paralog families and paralogons can be searched by text or sequence, and are downloadable from the website as plain text files. The database will be very useful for both experimentalists and bioinformaticians interested in the study of duplication events or paralog families.

INTRODUCTION
The occurrences and consequences of gene and genome duplication events have been discussed for a long time (1,2). The duplication of genes and large genome regions (or entire genomes) is proposed to be an important mechanism for the evolution of phenotypic complexity, diversity and innovation, and as an origin of novel gene functions. To uncover the evolutionary trajectories of duplicated genes, previous studies have integrated transcriptomic, interactomic and other data (1). Such integrated approaches, focusing on gene duplications in genomes, have already contributed robust insights into important evolutionary questions, such as the complexity of genes (3), the evolution of genome architecture (4), growth of gene networks (5), the 2R hypothesis (6) and diversity of gene expression (7). Moreover, the duplicated genes can be used to investigate diverging gene functions, which, when allied with computational methods, may provide useful information for experimental approaches. An example is the analysis of the molecular basis of the adaptive evolution of the duplicated pancreatic ribonuclease gene in leaf-eating monkeys with both computational and experimental approaches (8).
As more genomes are examined, increasing evidences support the dominating role of gene duplication events in the expanding of genome content (2,9). A crucial step in the study of gene duplications is to identify duplicated genes (known as paralogs) in genome sequences and to distinguish these from genes that have similar sequence but arisen from convergent evolution or other mechanisms. Algorithm-based homology detection from primary sequences is the preferred approach to detect paralogs or paralogous regions (4).
In contrast to ortholog databases, there are only a few specific paralog databases available in the public domain. Even though several general homolog databases, such as Inparanoid (10), Ensembl Compara (11), NCBI homologene (12), include some paralog information, they did not comprehensively summarize and display the evolution information of paralogs. In order to construct a stable web resource that supports easy browsing and downloading of evolutionary information on paralogous genes, we created EPGD (Eukaryotic Paralog Group Database; http://epgd.biosino.org/EPGD/). Several steps used to identify the paralogs contained in the EPGD were used previously to detect the duplication events in the family of animal transmembrane genes (13). Using this work (13) as a basis, we developed a semi-automatic procedure for collecting the within-species paralog families from genomes and pre-calculating several evolutionary indexes of these families. We collected the paralogs only from eukaryotes, as they are known to have a higher rate of gene duplication than Prokaryotes (14) and are more widely studied in this field.
A pioneer in the construction of paralog database is paraDB (15). A highlight of paraDB is the display of paralogons, which have been thoroughly investigated in the human genome (16) and are reviewed by Van de Peer (4). EPGD inherits this feature and adopts the term 'paralogon', defined as homologous genomic segments created by partial or complete genome duplication. EPGD focuses on families of paralogs and integrates spatial and temporal data to diagnose gene duplication processes comprehensively (17). The ratio of dN (the rate of non-synonymous substitutions) to dS (the rate of synonymous substitutions) (18), synonymous distance/clock, transition redundant exchange (TREx) clock (19), paralogons and several other features were generated by computational methods and deposited in the database.
In the current EPGD version, 26 eukaryotic genomes were processed and 35 991 paralog families and 29 480 paralogons were identified and stored (Table 1). To our knowledge, it is one of the most extensive paralog databases in public domain. All data can be browsed, searched and downloaded directly from the website.

CONSTRUCTION AND CONTENT
EPGD is implemented through MySQL relational database (http://www.mysql.com) and JavaServer Pages technology (http://java.sun.com/products/jsp/). The raw datasets of 26 eukaryotic genomes (Table 1) in GeneBank flat file format (GBK) were downloaded from the NCBI FTP site (ftp://ftp.ncbi.nih.gov/genomes) in March 2007. Proteins, coding sequences (CDS) and gene location information were extracted from these GBK files with a PERL script.
Overview of the procedure A total of 531 715 coding sequences and corresponding proteins were obtained after preprocessing. Only the protein sequences were used to construct the paralog families. The procedure is briefly described below: (i) Pairwise alignments of the proteins using gapped BLAST (20), with filtering for low sequence complexity regions using SEG (21). The default parameters were used, except for the threshold E-value of 10 À5 . (ii) Definition of the homologous genes. Four criteria must be satisfied. (a) all high-scoring segment pairs (HSPs) in the target sequence have to be arranged in the same order as in the query protein sequence (22); (b) the remaining HSPs cover more than 80% of the protein length; (c) the similarity of each HSPs is more than 50% (two amino acids are considered similar if their BLOSUM62 similarity score is positive) (22) and (d) these conditions are symmetrical for both genes. (iii) Single linkage clustering of homologous genes (13).
Generation of the primary paralog families.  (25) and the Yang and Nielsen methods (26), which were carried out using yn00 from the PAML (Phylogenetic analysis by maximum likelihood) packages (27). The TREx distances were computed based on the definition (19): the fractional identity of silent sites in conserved 2-fold redundant codon sites, which was implemented by ourselves. (viii) Construction of the arithmetic average (UPGMA) trees for grouping the proteins in a paralog family. These trees were derived from the dS matrix, because the synonymous substitutions are thought to be approximative neutral molecular markers. (ix) Identification of the paralogons using the algorithm developed by McLysaght et al. (16). Paralogons are two genomic segments that share a set of paralogous genes (4,16). After tandem duplications were masked, a greedy search algorithm was used to identify all paralogons between all pairs of chromosomes, based only on gene content but not gene order (4). Two criteria must be satisfied for a pair of paralogons. (a) they should contain at least two pairs of paralogous genes; (b) the gap size between two neighboring paralogous points in either chromosome should be less than the average length of 30 genes (16).

Content in the database
Large datasets were obtained when the procedure was applied to 26 genomes. We housed the data in a MySQL relational database. The kernel tables in the schema of EPGD are the table of paralog families and the table of paralogons. The peripheral tables, i.e. evolutionary indexes and annotation information, surround these two core tables. A summary of the data in EPGD is shown in Table 1.

Web interface
The web interface was implemented using Java and JavaServer Pages technologies. The user can inspect the datasets in the EPGD and see a summary of the current version. The records of paralog families, paralogons and genes ( Figure 1) are randomly selected each time when 'Glance' page is visited (http://epgd.biosino.org/ EPGD/glance.jsp). As shown in Figure 1, if the gene record is obtained, the corresponding paralog family and paralogons can be linked from this page. The main content of the gene page ( Figure 1A) starts with basic information of this gene (NCBI gene ID, taxonomy, EPGD family ID, location in the chromosome and simple description), followed by EPGD paralogons, which include or cover this gene. We defined that a gene is 'included' in a paralogon if it has at least one corresponding paralog in this paralogon region (paralogon-defining gene), while a gene is 'covered' by a paralogon if it does not have any corresponding paralog in this paralogon region (paralog-intervening gene). The coding sequences of the gene are listed at the bottom of the page.
The outline of the family page is similar to that of gene page ( Figure 1B). Multi-aligned sequences in protein or codon level, pre-calculated evolution indexes [dN, dS, TREx (19), etc.] and UPGMA tree based on dS are displayed on this page. The multi-alignments can be viewed in plain text or be displayed with the Jalview alignment viewer (28) (Figure 1). In the page which is hyperlinked from 'Evolution indexes of Pairwise CDSs', a row with a dN/dS different from the neutral expectation of 1 (z score > 1.96 or z score < À1.96) is color coded orange (Figure 1). The z score is computed using equation (18) where z is the z score, dN is the rate of non-synonymous substitutions, dS is the rate of synonymous substitutions, SE dN and SE dS are the standard errors of dN and dS, and Cov(dN, dS) is the covariance of dN and dS. We assume that the non-synonymous substitutions and the synonymous substitutions are independent and set Cov(dN, dS) to zero (18). The main part of the paralogon page contains basic information (taxonomy, locations in the chromosomes, average block length, average block density, number of links) of the paralogon, followed by an image thumbnail displaying a graphic view of the paralogon. Here, 'the average block density' is the arithmetic mean of the ratio of paralogon-defining genes to all genes in both sides of the paralogon; 'number of links' is the number of unique paralog families linked in the paralogon region. When the mouse hovers over this thumbnail, an enlarged view of this image pops up. Gene names and their regions in the enlarged graphic view of this paralogon are hyperlinked to the gene records in database.
The user can access the records in the EPGD with customized queries (Figure 2). From the 'iSearch' webpage (Figure 2A), 'any text' and nucleic acid or protein sequences can be searched without setting any parameter. Advanced Search pages with numerous input options ( Figure 2B and C) can be accessed via the links ('Advanced Text Search' or 'Advanced Sequence Search') from 'iSearch' page. The sequence search is powered by NCBI Blast package (20). Each search returns a result list of records in the database, which provides the hyperlinks to detailed pages ( Figure 2D).

DATA AVAILABILITY
The EPGD is available for download through the 'DOWNLOADS' link in the website as a FASTA file containing all proteins, family members lists, evolutionary indexes and paralogon regions in plain text files.

RESULTS AND DISCUSSION
The properties of the paralog family spaces in EPGD Table 1 gives a summary of the content of the current EPGD version. The proportions of duplicated genes in eukaryotes collected by EPGD range from 9% (Plasmodium falciparum) to 52% (Strongylocentrotus purpuratus), and are smaller than previously reported (e.g. Homo sapiens, 38%; Arabidopsis thaliana, 65%; Drosophila melanogaster, 41%; Caenorhabditis elegans, 49%; Saccharomyces cerevisiae, 30%) (2). This is due to the rigorous criteria for paralog definition used to construct the EPGD and because many duplicated genes have eliminated characteristic signatures from their sequences during their evolution history (2). Since evolutionary indexes are highly unreliable for ancient gene duplications, rigorous criteria are essential for our database.
The size of the paralog families tends to be smaller than five genes. The distributions of paralog family size in all species of EPGD follow power law (data not shown) (29,30). As an example, Figure 4A displays the distribution of paralog family sizes in H. sapiens and the corresponding log-log diagram. The power law distribution indicates the robustness of our family detection method and the quality of gene prediction in the original data (29).
Consistent with previous studies on Bacteria and a small set of Eukarya (9,29,31), large genomes possess more paralog families and a higher proportion of genes belonging to paralog families than small genomes ( Figure 3A and C). We find, however, only a weak correlation between the average size of families and the genome sizes ( Figure 3B, r = 0.26, P = 0.19), in contrast to the finding in Bacteria that average family size increases with genome size (31). This result suggests that the higher percentage of paralogs in large eukaryotic genome stems mainly from the emergence of new paralogon families. An expansion of existing gene families is not evident in Eukarya ( Figure 3B).
The number of paralogons increases with the genome size ( Figure 3D, r = 0.86, P = 3.356 Â 10 À8 ), indicating the effect of duplication of large genome segments on the evolution of genome size. Furthermore, the distribution of the paralogon size is also a skewed distribution (e.g. Figure 4D). Most of the paralogons have less than five linked familes (98% of all human paralogons), because of the high level of gene loss after duplication, as well as recombination, chromosomal rearrangements and recombination. Still, the identification of putative paralogons provides many insights into evolutionary mechanisms (4).

The example of H. sapiens
Taking H. sapiens as an example (Figure 4), we plotted the distribution of paralog family size ( Figure 4A), a scatter diagram of TREx distance versus dS ( Figure 4B), a log-log graph of dN versus dS ( Figure 4C) and the distribution of paralogon size (the number of linked families) ( Figure 4D). Transition redundant exchange (TREx) processes at the position of conserved 2-fold codon sites are thought to offer an approximation for a neutral molecular clock (19). We calculated the TREx distances for each paralog family, which provide a more homogeneous molecular clock than that provided by the dS. If the time since two genes diverged is long relative to the reciprocal of the rate constant with which these silent sites suffer transition substitutions, the TREx distance approximates 0.5. As seen from Figure 4B, TREx distances are negatively correlated with dS ( Figure 4B, r = À0.89, P < 2.2 Â 10 À16 ). Therefore, the TREx distance can be used as an alternative of dS.
Similar to the work of Lynch et al. (32), dN was plotted as a function of dS ( Figure 4C). The accumulation of nonneutral points when dS increases ( Figure 4C) confirms the gradual increase of selective constraint on duplicates during evolutionary history (32). When dS is greater than 2, there are more points around the neutral expectation ( Figure 4C). This is an artifact, resulting from the saturation effects in the estimation of dN and dS (33).

PERSPECTIVES
We plan to update EPGD every six months. As new eukaryotic organisms are fully sequenced and annotated, they will be added to EPGD using our procedure. In the future, ortholog annotation information will also be included. However, the development of the utilities for EPGD will still focus on tools for the analysis of duplication events, such as statistical tests of the paralogons (unpublished data) and chromosome ideograms. Furthermore, we will thoroughly analyze the data in EPGD and present insights into the effect of duplication events on genome evolution. The procedure to build the EPGD is currently semi-automatic. We will make the procedure totally automatic and start an open source project in the future.


extfieldsofstudy: Medicine Biology

s2fieldsofstudy: Biology Medicine

text:
A case-control study of rheumatoid arthritis identifies an associated single nucleotide polymorphism in the NCF4 gene, supporting a role for the NADPH-oxidase complex in autoimmunity

Rheumatoid arthritis (RA) is a chronic inflammatory disease with a heritability of 60%. Genetic contributions to RA are made by multiple genes, but only a few gene associations have yet been confirmed. By studying animal models, reduced capacity of the NADPH-oxidase (NOX) complex, caused by a single nucleotide polymorphism (SNP) in one of its components (the NCF1 gene), has been found to increase severity of arthritis. To our knowledge, however, no studies investigating the potential role played by reduced reactive oxygen species production in human RA have yet been reported. In order to examine the role played by the NOX complex in RA, we investigated the association of 51 SNPs in five genes of the NOX complex (CYBB, CYBA, NCF4, NCF2, and RAC2) in a Swedish case-control cohort consisting of 1,842 RA cases and 1,038 control individuals. Several SNPs were found to be mildly associated in men in NCF4 (rs729749, P = 0.001), NCF2 (rs789181, P = 0.02) and RAC2 (rs1476002, P = 0.05). No associations were detected in CYBA or CYBB. By stratifying for autoantibody status, we identified a strong association for rs729749 (in NCF4) in autoantibody negative disease, with the strongest association detected in rheumatoid factor negative men (CT genotype versus CC genotype: odds ratio 0.34, 95% confidence interval 0.2 to 0.6; P = 0.0001). To our knowledge, this is the first genetic association identified between RA and the NOX complex, and it supports previous findings from animal models of the importance of reactive oxygen species production capacity to the development of arthritis.

Introduction
Rheumatoid arthritis (RA) is a chronic inflammatory disease that leads to erosion and deformation of the joints. The prevalence in the general population is 0.5% to 1%, and women are at two to three times greater risk for developing the disease. Twin studies show a concordance rate of 12% to 15% in monozygotic twins and 4% in dizygotic twins, and the genetic heritability is estimated at 60% [1]. Despite much effort to identify arthritis causing genes, only few genetic loci have been confirmed to be associated with RA, among which are the human leucocyte antigen (HLA) locus and the protein tyrosine phosphatase non-receptor 22 gene (PTPN22) [2][3][4].
RA is a heterogeneous disease with symptoms and disease progression that vary between patients. Classification of typical RA requires fulfilment of four out of seven criteria established by the American College of Rheumatology [5]. The heterogeneity of the symptoms in RA probably mirrors the CCP = cyclic citrullinated peptide; CEPH = Centre d'Etude du Polymorphisme Humain; CI = confidence interval; EIRA = Epidemiological Investigation of Rheumatoid Arthritis; HWE = Hardy-Weinberg equilibrium; NOX = NADPH oxidase; OR = odds ratio; RA = rheumatoid arthritis; RF = rheumatoid factor; ROS = reactive oxygen species; SE = shared epitope; SNP = single nucleotide polymorphism. underlying genetic contribution; hence, the various combinations of symptoms observed in patients are probably caused by different combinations of risk alleles.
In the search for markers that can predict development of the disease before its onset, several autoantibodies, including rheumatoid factors (RFs) and antibodies against cyclic citrullinated peptides (anti-CCP antibodies), have been detected [6][7][8]. RFs, autoantibodies that recognize the Fc part of immunoglobulins [9], are present in 60% to 70% of RA patients and have been found to be associated with more severe clinical manifestations [10,11]. Anti-CCP antibodies are present in 50% to 60% of RA patients and have also been shown to predict more severe disease [10,[12][13][14]. Interestingly, recent data suggest that the HLA-DRB1 locus, which has been shown to be associated with RA, is only associated with the presence of anti-CCP antibodies, and this association is independent of both RA development and the presence of RFs [15]. Together with recent data describing the interaction between environmental factors and genetic predisposition [16,17], these findings support the hypothesis that several independent disease mechanisms may lead to development of RA. They also emphasize the need to use relevant subgroups of RA patients in genetic association studies.
Various approaches have been used to identify genes that contribute to common diseases such as RA [18,19]. Because of increasing knowledge of gene functions and immunological pathways, a candidate gene approach can be more efficient in terms of both cost and time than traditional linkage analysis or genome-wide approaches. The selection of genes in a candidate gene study can be based on previous knowledge of gene functions as well as on disease or immunological mechanisms. It can also be based on findings in animal models.
The Pia4 locus in rats has been shown to reduce the severity of pristane-induced arthritis (an arthritis model) [20]. A single nucleotide polymorphism (SNP) in the Ncf1 gene was found to be responsible for the protective effect [21]. Ncf1 encodes the protein p47 phox , which is part of the NADPH-oxidase (NOX) complex that produces reactive oxygen species (ROS) in response to infectious stimuli. Rats carrying the risk allele have reduced capacity to produce ROS, which is linked to increased risk for development of arthritis [21,22]. Apart from NCF1, the NOX complex is composed of five other proteins encoded by the genes NCF2, NCF4, CYBB, CYBA and RAC2 [23,24]. In phagocytes massive ROS production, called the oxidative burst, takes place in the phagosome or endosome (intracellular) or in the plasma membrane (extracellular) after ingestion of invading pathogens or after stimulation of innate immune receptors [25][26][27]. However, the ability to generate ROS extends to cells other than classical phagocytes, such as dendritic cells, suggesting that the NOX complex has additional functions in the immune system [28]. Interestingly, recent findings from our group show that the redox balance of T-cell membranes has an important effect on the activation and proliferation of T cells [29].
Results obtain in animal models make NCF1 and the other genes of the NOX complex candidate genes for human RA. However, transferring animal data to the human setting is not straightforward in this case. In contrast to rats, the genomic organization of NCF1 in humans is very complex, which makes SNP-based association analysis difficult [30,31]. However, because of the complex interplay between the proteins of the NOX complex, it is likely that genetic changes in any of the genes could have the same effect on ROS production, as is seen for Ncf1 in rats.
We used a candidate gene approach to study the association with RA of 51 SNPs in five genes included in the NOX complex. We found a SNP, located in intron 4 in NCF4, to be associated with RA in a Swedish case-control cohort. This supports the hypotheses that the NOX complex is involved in the development of RA and that ROS could be an important regulator of the immune system.

Selection of single nucleotide polymorphism
The HapMap genome browser [32] was used as the primary source for selection of SNPs. Primarily, HapMap-validated SNPs were selected. However, in order to obtain evenly dispersed SNPs across the genes, SNPs were also selected from dbSNP [33]. The selection was based on several criteria; the minor allele frequency should be more than 5% in the Centre d'Etude du Polymorphisme Humain (CEPH) or other European cohort, the SNP should be located in a nonrepetitive sequence, and the validation status at dbSNP should be at least two-hit. We also consulted the Linkage disequilibrium (LD) maps of the HapMap CEPH samples in order to disperse evenly the SNPs with respect to the LD structure. Sixty-seven SNPs were selected (see Additional file 1): six in CYBB, eight in CYBA, 22 in NCF2, 21 in NCF4 and 10 in RAC2. The number of SNPs selected in each gene reflects both the size of the gene and the availability of validated SNPs.

Samples
Samples from the Epidemiological Investigation of Rheumatoid Arthritis (EIRA), a population based case-control study, were analyzed. The study base comprised the population, aged 18 to 70 years, in a geographically defined area in the middle and southern parts of Sweden during the period from May 1996 to 2005. A case was defined as a person in the study base who, during the study period and for the first time, received a diagnosis of RA based on the American College of Rheumatology criteria of 1987 [5]. All cases were examined and diagnosed by a rheumatologist at a participating centre. All public rheumatology units in the study area and almost all of the (very few) private units participated in the study. The present study includes data from 1,842 RA patients and (page number not for citation purposes) 1,038 controls, randomly selected from the study base and matched on age, sex and residential area. Of the cases and controls, 71% and 72%, respectively, are female; the mean age was 51 ± 13 years in cases and 54 ± 12 in controls. Information about RF and anti-CCP antibody status was available for 1,315 of the patient samples; 63% of these were RF positive and 61% were anti-CCP positive.

Genotyping
Assay design, validation and genotyping were performed by Kbiosciences (London, UK) using a fluorescence resonance energy transfer based competitive allele-specific polymerase chain reaction system (KASPar

Hardy-Weinberg analysis
The Haploview software calculates P values for deviations from Hardy-Weinberg equilibrium (HWE) for each marker on the complete uploaded dataset [34]. A significance threshold of P ≤ 0.001 was used. In order to obtain HWE P values for the case and control groups separately, each dataset was uploaded separately.

Single marker analysis
Contingency tables were created for each SNP using the JMP 5.0 software (SAS Institute, Cary, NC, USA), and P values for association with RA were calculated using χ 2 tests for genotype frequencies. The sample set was stratified for sex in the initial analysis, and for sex and RF or anti-CCP status in the stratified analysis. Because CYBB is located on the X chromosome, allele frequencies were used to estimate association in the male samples.
To obtain a corrected α level in the stratified analysis, the Bonferroni correction method was applied, as implemented on the simple interactive statistical analysis website [35].
We also performed logistic regression analysis for the rs729749 SNP, in order to adjust for age, sex and living area. Logistic regression analysis was conducted for a subset of the material, in which all information regarding genetic factors, antibodies and matching variables was available. This same sample set was used for the frequency analysis of rs729749. We used the SAS software for Windows (version 9.1; SAS Institute, Cary, NC, USA) to perform logistic regression analysis.

Haplotype analysis
Haplotype blocks were calculated using Haploview [34]. SNPs that departed from HWE in both cases and controls were excluded from the analysis. The haplotype predictions were based on the CI method proposed by Gabriel and coworkers [36]. However, the other available methods yielded the same result. Haplotypes were assessed using the WHAP software [37]. Sub-haplotypes were investigated for associations based on the block structure predicted in Haploview. Because the haplotype block analysis is based on all genotyped SNPs, we used only the genotype information from the subset of the EIRA cohort comprising 1,069 cases and 634 controls. The permuted P values for the haplotypes are based on 5,000 permutations.

Genome analysis
The region surrounding rs729749 was investigated for the presence of regulatory elements, transcription factor binding sites and conserved regions using the University of California, Santa Cruz (UCSC) Genome Browser [38]. The ESPERR Regulatory potential (seven species), Conservation, and TFBS Conservation options were used. P values and odds ratios (OR) with 95% CIs for the genetic models were calculated using the GraphPad Prism software (GraphPad Software, San Diego, CA, USA).

Fifty-one single nucleotide polymorphisms were successfully genotyped in a Swedish rheumatoid arthritis cohort
To investigate the genes in the NOX complex for association with RA, 67 SNPs evenly dispersed in the five candidate genes (Table 1) were selected from the HapMap genome browser and dbSNP (National Center for Biotechnology Information; Additional file 1). Fifty-one out of the 67 SNPs were successfully assayed and genotyped in the Swedish EIRA cohort [39,40], with a 98% call rate (Additional file 1). Deviations from HWE were estimated for each SNP in all samples combined and in cases and controls separately. The SNP rs3788524 (NCF4) significantly deviated from HWE (P < 0.001) in both the case and control groups, indicating a possible genotype failure.

Genotype analysis indicates male-specific associations with rheumatoid arthritis in NCF2, NCF4 and RAC2
To evaluate the SNPs for association with RA, contingency tables were created for genotype frequencies and P values were calculated using χ 2 tests. Because RA is more frequent in women than in men, the disease mechanisms might be sex dependent. The samples were therefore analyzed either all together or stratified by sex. The initial cut-off P value for significance was set at 0.05 to reduce the risk for missing subgroup-specific associations, which would be weak in the complete sample set. Three SNPs fulfilled this criterion for association with RA, all in men; rs789181 in NCF2 (P = 0.03), rs729749 in NCF4 (P = 0.001) and rs1476002 in RAC2 (P = 0.05; Table 2).

Stratification based on autoantibody status reveals a highly significant association of rs729749 with rheumatoid factor negative rheumatoid arthritis in men
Because RA is a heterogeneous disease, it is likely that the disease contributing genes differ in different subgroups of patients. Presence or absence of autoantibodies such as RF and anti-CCP antibodies might reflect different underlying disease mechanisms. In order to determine whether the detected associations are specific for a subclass of RA, we stratified the male cases by RF or anti-CCP antibody status and performed an analysis on genotype frequencies for the associated SNPs.
We found the SNPs rs789181 (NCF2) and rs1476002 (RAC2) to be significantly associated in RF-positive men (P = 0.03 and 0.02, respectively). The SNP rs1476002 was also significantly associated in anti-CCP antibody positive men (P = 0.05). However, the P values are only slightly improved com-pared with the initial analysis, suggesting that the associations are not specific for any of these subgroups. On the other hand, the association for rs729749 in NCF4 is highly significant in RF-negative men (P = 0.0002), whereas no association was detected in RF-positive men (P = 0.0713; Table 3). Furthermore, a comparison of the frequencies in RF-negative against RF-positive men also yielded a significant association with RFnegative disease (P = 0.01). We also detected a weaker association in anti-CCP antibody negative men (P = 0.003; Table  3), although the anti-CCP antibody negative versus anti-CCP antibody positive comparison is not significant (P = 0.08). Logistic regression analysis of the CT or TT genotype against the CC genotype, adjusted for age and residential area, also indicated a sex and subgroup specific effect for rs729749 (Table 4). With regard to RF-negative RA among men, the OR for CT against CC is 0.34 (95% CI = 0.20 to 0.60; P = 0.0001). Here we also see a strong association in anti-CCP antibody negative men (OR = 0.4, 95% CI = 0.2 to 0.7; P = 0.0004).
To determine the validity of the P value, we estimated a corrected α value using the Bonferroni correction method. The corrected α level for the number of tests performed (159) was estimated at 0.0003. Hence, only the association detected for rs729749 in RF-negative males remained significant using the corrected α level.

Haplotype analysis reveals an associated haplotype in NCF4 caused by the rs729749 single nucleotide polymorphism
We wished to estimate haplotypes and haplotype blocks to determine whether rs729749 is included in a conserved and possibly associated haplotype. The Haploview software was used to determine the haplotype block structures of the genes. The analysis showed that the NCF4 gene is divided into three blocks ( Figure 1). One block of 1 kilobase contains rs729749 as well as three other genotyped SNPs: rs760519, rs2284027 and rs17811059. To test the association of the haplotype blocks, we used the WHAP software. The block containing rs729749 is borderline significant in the nonstratified analysis (P = 0.055), but excluding rs2284027 and rs17811059 yields a significant association (empirical P = 0.04 in all and 0.03 in men; Table 5). In order to determine how  (11) CCP, cyclic citrullinated peptide; RF, rheumatoid factor. Values indicate the genotype frequencies and number (in brackets). important the impact of the rs729749 SNP is for the identified associated haplotype, we conducted a conditional analysis in which the association of the haplotype was estimated while controlling for the haplotype background of rs760519 and rs2284027. This test gave an empirical P value of 0.01, which is lower than for the associated haplotype, indicating that the rs729749 SNP itself or a SNP in strong linkage disequilibrium is causing the haplotype association.

Homozygosity could explain the genetic risk associated with rs729749
The rs726749 SNP is noncoding and located in the beginning of intron 4 in NCF4. Analysis of the region using the UCSC genome browser (March 2006 assembly) [38] revealed that it is not located at any known transcription factor binding site or regulatory element; however, there is a conserved regulatory region (as predicted by the ESPERR regulatory potential option) approximately 300 base pairs downstream of rs729749. Nonetheless, the genetic effect of this SNP is not obvious.
By looking at the genotype frequencies in RF-negative males for rs729749, it appears as though both homozygous genotypes are more common in the RA patients (Table 3). This led us to perform an analysis to determine which genetic model best fits the data. Not surprisingly, the homozygous model (CC + TT versus CT) gave the best fit in the subgroup of RFnegative men (OR = 2.67, 95% CI = 1.60 to 4.46; P = 0.0002; Table 6). However, because the allele frequency of the T-allele is quite low and few individuals are homozygous for this allele, there is a degree of uncertainty in this prediction. The CC genotype is clearly more common in the patients and homozygosity at this locus also gives a significant genetic model (Table 6).

Discussion
Here we report an association of the SNP rs729749, located in the gene NCF4, with RA in Swedish RF/anti-CCP antibody negative male patients. NCF4 is part of the NOX complex, and our results support recent findings that ROS play a role in the development of RA. However, because this is the first report of an association between a gene in the NOX complex and RA, genetic replication as well as functional studies will be required before it is possible to determine conclusively whether NCF4 is an RA susceptibility gene. Furthermore, the results from this study do not exclude the possibility that the other genes in the NOX complex could affect RA susceptibility. There could, for example, be epistatic effects between the different genes in the complex.
There is a lack of consensus in the field about statistical significance levels in association studies. In this study we initially used an α level of 0.05 to reduce the risk for overlooking associations that are specific for a subgroup of RA patients. There are strong medical and biological arguments in favour of subgrouping RA patients, and an analysis of nonstratified samples might conceal a true association. In the initial analysis of the genotype frequencies we found indications of association with three SNPs, and stratification of the samples indicates that the association for rs729749 is male specific and predominates in RF-negative disease. We also found a weaker association in anti-CCP antibody negative disease, but at this stage it is not possible to say whether this is due to the fact that the presence of RFs and anti-CCP antibodies is dependent [41] or whether the anti-CCP antibodies independently influence genetic risk.
The detected association of rs729749 is strongest in RF-negative males, and therefore it is possible that this particular SNP modulates the affected disease pathway so that it mainly influences the clinical disease in RF-negative RA, specifically in men. However, it is not impossible that this pathway play a role also in other subtypes of RA. The results from this study strengthen the view that different combinations of genes are involved in the disease progression in different RA subclasses defined, for instance, by sex or autoantibody status. The most striking aspect of the identified association of the rs729749 SNP is the clear sex specificity, indicating that the disease pathway affected by the rs729749 SNP is specific to men. The fact that RA is three times more common in women than in men suggests that men and women respond differently to factors that trigger the onset of RA and that some disease pathways could be sex specific. The RF-negative specificity of the association is not as striking as the sex specificity, but it might still point toward which RA subgroup is affected by the rs729749 SNP. Because both RFs and anti-CCP antibodies can be detected before the appearance of any disease symptoms [7,42], they might reflect or be part of initial disease mechanisms that are only present in anti-CCP antibody or RF positive patients. Furthermore, it has repeatedly been shown that patients who are negative for these autoantibodies have a less severe disease outcome than do anti-CCP antibody or RF positive RA patients [10][11][12]. This further supports the view that autoantibody negative RA is another variant of RA. Also, the few genes confirmed to be associated with RA have been shown mainly to affect subclasses of RA.
The well established association of the HLA-DRB1 alleles, which encode the so-called shared epitope (SE), was recently suggested to be restricted to anti-CCP antibody but not RF positive RA [15,16,41]. However, it has also been shown that SE alleles are associated with higher titres of autoantibodies against type II collagen [43]. Interestingly, the HLA-DR3 locus, which does not encode the SE, appears to be associated only Haplotype blocks in NCF4 Haplotype blocks in NCF4. Three haplotype blocks were identified in NCF4. Block 2 contains the associated single nucleotide polymorphism (SNP) rs729749. Colour scheme of the linakge disequilibrium LD map is based on the standard D'/LOD option in the Haploview software. The LD blocks are calculated based on the CI method.
with anti-CCP negative RA [41,44]. The recently discovered PTPN22 risk allele has been shown to be associated with both RF and anti-CCP antibody positive RA [43,45,46], but there are conflicting data regarding the autoantibody restriction of this association [47]. Furthermore, the reported association of the PD-1.2A allele of the PDCD1 gene is restricted to RF negative as well as SE negative RA [48]. Hence, even though little is known about the precise disease mechanisms in RA, it is obvious that the symptoms can be caused by several different distinct disease pathways.
Because of the lack of certain autoantibodies, it is tempting to speculate that RF or anti-CCP antibody negative RA patients have a lesser B-cell component, and reports indicating that Bcell depleting therapy is a less effective treatment for RF-negative RA supports such a view [49,50]. Several recent clinical studies, however, have yielded conflicting results on B-cell depletion efficiency [51], and more extensive studies are required before any conclusions can be drawn about the Bcell dependence in RF and anti-CCP antibody negative RA. In addition, it is fair to assume that the autoantibody negative forms of RA, in themselves, represent a heterogeneous collection of pathways related, for example, to osteoarthritis, arthropathies, or lupus.
The rodent arthritis models are heterogeneous as well but they may relate to more defined arthritis conditions than the human disease. Interestingly, the pristane-induced arthritis model used to identify the polymorphism in Ncf1 in rats is B-cell independent. Studies have shown that transfer of T-cells from genetically susceptible and immunized rats into genetically resistant and un-immunized rats is enough to induce arthritis in the resistant rat [52].
Both the HLA-DRB1 locus and PTPN22 are believed to affect T-cell activation. HLA-DRB1 is expressed by antigen-presenting cells and restricts the presentation of antigens to T cells. PTPN22 on the other hand has an intracellular effect, and the amino acid shift to the disease-associated tryptophan impairs the function of the Lyp protein, encoded by PTPN22, rendering it a less potent negative regulator of T-cell activation [53].
Studies in animal models as well as recent findings suggest that the NOX complex could also be involved in T-cell activation [29,54]. ROS production has been found in antigen-presenting cells, and it has also been shown that T-cell receptor signalling is affected by oxidation of the T-cell membrane as well as by intracellular levels of ROS [55,56]. Interestingly, hydrogen peroxide generated from • O 2produced by the NOX complex can readily cross the cell membrane and inactivate protein tyrosine phosphatases, including that encoded by PTPN22, through oxidation of a specific cysteine residue [54]. ROS production of antigen-presenting cells, such as dendritic cells, has been found to be crucial for the activation of T-cells [28] and could therefore determine the immune response to an antigen. Furthermore, recent work from our group highlights the importance of the redox balance on the surface of T cells for activation. T-cells originating from the Ncf1 mutated rat DA have an increased number of reduced thiol groups on the cell surface, which increases proliferation and activation [29]. Gelderman and coworkers [29] showed that T-cells originating from the nonmutated E3 rats could be made arthritogenic by increasing the number of reduced thiol groups on the cell surface. Because NOX production could not be detected in Tcells, it was concluded that the redox balance of T-cell surfaces is determined by other cells, such as macrophages or dendritic cells.
In the rat a clear difference in the capacity to produce oxidative burst is the result of a single SNP in Ncf1, resulting in the shift from threonine to the disease-promoting methionine at position 153 in the p47 phox protein [21]. The consequences of the amino acid shift have not yet been identified, although position 153 does not coincide with any critical binding sites. It is most likely that the shift affects the three-dimensional structure of the protein, impairing the binding capabilities to the other proteins in the complex. The activation of the NOX complex is initiated through phosphorylation of the three cytoplasmic proteins p47 phox (NCF1), p67 phox (NCF2) and p40 phox (NCF4) [23,26,57]. The phosphorylations lead to conformational changes of p47 phox , and subsequent translocation of the three proteins to the membrane, where it interacts with and activates the membrane bound complex flavocytochrome b 558 , which is composed of gp91 phox (CYBB) and p22 phox (CYBA) [24,58]. RAC2 also translocates to the membrane after dissociation from Rho-guanine nucleotide dissociation inhibitor and interacts with flavocytochrome b 558 [23,59]. Studies show that p67 phox and Rac2 are critical for the function of the complex [25,60], whereas p40 phox and p47 phox function as adaptor proteins and are responsible for the assembly of the complex [61]. The precise role of p40 phox is still debated, and there is evidence of both positive and negative regulation of the NOX complex [61][62][63]. Recent studies have suggested that p40 phox , as well as p47 phox , functions by 'carrying' p67 phox to the membrane through the interaction with phospholipids [62,64]. Studies of mutations causing human chronic granulomatous disease show that mutations in NCF1, CYBB, CYBA, or NCF2 lead to a complete lack of function of the NOX complex [65,66]. However, SNPs in less critical regions, like the 153 SNP in Ncf1 found in rat, apparently reduce but do not completely abolish ROS production, which leads to an increased susceptibility to arthritis [21]. These facts indicate that the SNPs of interest for examination in an RA association study of the NOX complex will probably be located at noncritical positions and will not completely abolish but only slightly modify binding to the other proteins.
The rs729749 SNP is located in the middle part of intron 4 in NCF4. By studying the genome maps using the UCSC genome browser, we could not find any evidence of it affecting a transcription factor binding site or a known regulatory element. However, there is a conserved regulatory region predicted approximately 300 base pairs downstream of rs729749. The haplotype analysis indicates that rs729749, or a SNP in very high LD, is responsible for the association detected for the haplotype. Also, because the predicted regulatory region is very close to rs729749, it could very well contain the causative SNP.
The genetic analysis of the association shows that the homozygous model gives the best significance, indicating that this could be the genetic effect. However, the low number of samples in the subgroup of RF-negative males together with the low frequency of the T allele make this assumption somewhat uncertain. Furthermore, it is not obvious what the functional consequences of a homozygous effect would be. We therefore chose to present all models tested in order to provide the full picture.
In the logistic regression analysis, weak associations were also detected for RF and anti-CCP antibody positive disease. One possible explanation for these findings is that some patients could have been 'misdiagnosed' regarding autoantibody status. The diagnosis is made by measuring the autoantibody titres, and certain threshold are used to determine positive versus negative status. However, sometimes autoantibody status changes as the disease progresses, and therefore there is a degree of uncertainty regarding autoantibody status.
A method or a strategy to apply multiple testing corrections accurately in a case-control candidate gene association study has not yet been established [67,68]. The conventional multiple correction methods, such as Bonferroni, are considered by some to be too stringent for large-scale studies [69,70]. The Bonferroni correction method is used to correct for the number of independent tests performed. However, in an association analysis neither the SNPs tested nor the different stratified analyses are truly independent. The allele frequencies of the SNPs are affected by the LD of the region, and the stratified analyses are based on the same SNPs and samples as the initial analysis and can therefore not be considered to be independent from each other. Nonetheless, in order to obtain an indication of the validity of the P values obtained in these analyses, we used the Bonferroni correction method to estimate an α level correcting for the number of tests (159) performed in the frequency analysis. The 159 tests reflect the sex-stratified analysis of all 51 SNPs (153) plus the stratified analysis of RF and anti-CCP antibody status performed for rs729749, rs789181 and rs1476002 in the male subset (6). Correcting for 159 tests gives a new α level of 0.0003, which means that only the association of rs729749 in RF-negative males passed the significance threshold using this stringent method.
The EIRA study included mainly individuals born in Sweden. Taking those born outside Sweden into consideration also, 97% were of Caucasian origin. In order to minimize potential bias resulting from population stratification, we performed logistic regression analyses in which ORs were adjusted for age and residential area.

Conclusion
We found evidence of an RA-associated SNP in the NCF4 gene of the NOX complex. The association is specific for male patients, with the strongest association found in RF-negative RA. This finding supports the notion of RA being a heterogeneous disease caused by a variety of disease pathways that are regulated by a variety of contributing risk genes. The detected association with a component of the NOX complex also strengthens previous evidence obtained in animal models that suggests that the NOX complex and ROS have a major impact on inflammation and autoimmunity. It is hoped that this finding will help to elucidate the complex genetics that underlie RA and autoimmunity.
contributed with intellectual property to the study. RH conceived of the study, participated in its design and helped to draft the manuscript.