Selection of study sample

The sample of this study consists of 100 women from the BARCOS [¹⁰] cohort. This cohort is composed of about 1,500 Spanish postmenopausal women monitored at the Hospital del Mar de Barcelona. Women diagnosed with osteomalacia, Paget's disease, some metabolic or endocrine disorder, or those undergoing hormone replacement therapy or drug treatment that could affect bone mass were excluded from the cohort. Women with early menopause (before age 40) were also excluded. The data collected for each sample were BMD, age, age of first menstruation, age of menopause, years since menopause, weight and height. Blood samples and informed, written consents were obtained from each patient according to the regulations of the Marital Health Park Clinical Research Ethics Committee. BMD (g/cm2) was measured at the femoral neck and lumbar spine. A dual energy X-ray densitometer was used to carry out the measurements.

Two groups of 50 samples with extreme BMD values were selected according to the Z-score value. Specifically, the groups consisted of the 50 samples with the highest Z-score (range: 2.98 to 0.73) and the 50 samples with the lowest Z-score (range: -2.41 to -4.26) of the BOATS cohort.

Preparation of genomic samples

Each woman's DNA was extracted from peripheral blood samples. The concentration and quality of the DNA samples (260/280 and 260/230 ratios) were measured by spectrophotometry on a NanoDrop ND-1000 (NanoDrop Products) instrument. To determine DNA integrity, 5 µl of each sample was analyzed by 1% agarose gel electrophoresis. Finally, the samples were normalized to a concentration of 100 ng/µl.

Long-Range PCR (LR-PCR)

A 28 kb genomic region (containing the FLJ42280 gene [22 kb] along with 3.8 kb of flanking region at 5 'and 2 kb of flanking region at 3') was divided into 7 overlapping fragments (Figure 2). The sizes and coordinates of these 7 fragments and the primer pairs used to amplify them are shown in Table 1. The fragments, 2 to 5 kb, were amplified by LR-PCR. Each LR-PCR reaction included: 100 ng of genomic DNA, 5 µl of "Magnesium +" Ex Taq buffer (20 mM Mg++; Takara) x10, mixture of dNTPs (2.5 mM each), Ex Taq polymerase 5 U/µl) and primers (20 µM), in a final volume of 50 µl. The reactions were carried out in a GeneAmp® PCR System 2700 (Applied Biosystems) thermocycler. Each fragment required distinct elongation time and hybridization temperature conditions. The total number of amplicons was 700 (100 samples x 7 fragments). The quantity and quality of all amplicons were checked by 1% w/v agarose gel electrophoresis in TBE x1 buffer.

Table 1 Amplicons used to sequence the FLJ42280 region 

Figure 2 In the genomic region of FLJ42280 previous GWA studies have identified 4 SNPs showing association with BMD, details of which are shown in the table included in the figureTo better understand the variability of this locus, the region was subdivided into 7 overlapping fragments (yellow arrows) for resequencing in women with very high or very low BMD 

Purification and quantification of samples

To remove the residues from the PCR reagents, the PCR products were purified using 96-well filter plates with a suitable pore size (Pall Corporation). The vacuum (Vacuum Manifold, Merck Millipore) was applied and the DNA retained in the filter was resuspended in 35 µl of milliQ water. The PCR products were then quantified using the Quant-iT PicoGreen dsDNA Reagent Kit (Life Technologies) following manufacturer's instructions. Briefly, a standard curve of concentrations was constructed by fluorescence emission measurements at 520 nm after the DNA was excited at 480 nm. The curve was then used to calculate the samples' DNA concentration.

Table 2 Number and location of single nucleotide variants found in this study 

^tfn4*includes flanking regions, 5'UTR, 3'UTR and introns; LFV: low frequency variant (MAF <5%).

Standardization of sample concentrations and pooling

The samples on the plates were normalized to a concentration of 5 ng/µl and then 5 µl of each sample of one plate (one tube per plate) was mixed in a single tube using the epMotion® 5075 Liquid Handling Workstation (Eppendorf). Thus, 14 tubes with 250 µl each were obtained, two for each PCR fragment (high BMD and low BMD). The 14 pools were pooled up to 5 times using the Genevac EZ-2 evaporator (Genevac SP Scientific) and each tube was quantified using Qubit® 2.0 Fluorometer (Life Technologies). Finally, the PCR fragments were mixed equimolarly in two tubes, one for high BMD and one for low BMD.

Massively parallel sequencing

The mass sequencing of the samples was carried out at the Genomics Service of the Scientific and Technological Centers of the University of Barcelona, using the GS 454 Junior System (Roche). Briefly, the DNA was fragmented by nebulization, two labeled libraries were prepared with adapters (sequences of 10 nucleotides), one for each group, which were mixed in a single tube. The mixture was then amplified by emulsion PCR and the final library was loaded onto a picotiter plate for pyrosequencing. Four sequencing runs were carried out, corresponding to 140 Mb of final data (35 Mb/stroke). This volume of data provides a theoretical coverage of 40x for each initial sample.

Processing of sequencing data and variant selection

The readings obtained from the sequencing were preprocessed based on their quality and aligned against the reference genome (GRCh37) using the GS Mapper program (Roche). The readings were indexed and filtered using SAMtools. The variants present in the two groups were detected by GATK using standard filtering parameters [¹¹]. The variants found were prioritized according to the following criteria: variants were selected with a coverage of at least 1,000 readings, present in 1% of readings and with a low bias strand. The number of readings of the variants that passed the filters was normalized by the coverage and the variants were classified between common (with a frequency greater than 5%) and rare or low frequency (with frequency less than 5%).

Functional and statistical analysis of variants

The frequencies of each variant were compared between the two groups using an exact Fisher's test, applying the Bonferroni correction for multiple comparisons. The functional analysis of variants consisted in looking at whether they were described in databases such as dbSNP and 1000 Genomes and, if so, searching for their MAF in the European and Iberian population. In addition, for the exonic variants it was observed what amino acid change they assumed and their severity predicted by SIFT, PolyPhen and Provean. For the intronic variants, we analyzed the region containing the variant: sites of hypersensitivity to DNAse, binding of transcription factors, DNA methylation, histone modifications and regulatory regions. All these data were obtained from databases and repositories such as Ensembl, UCSC Genome Browser, ENCODE, BioMart, MatInspector. HaploReg was also used to search for regulatory annotations. Finally, all the variants found were analyzed with the Ensembl Variant Effect Predictor and UCSC and with the SNP function prediction from the US National Institute of Environmental Health Sciences.

Analysis of linkage disequilibrium

To calculate the linkage disequilibrium between all variants of the genomic region of FLJ42280 we used the genotypes of the SNPs present in the region and the haplotypes of the individuals of HapMap phase 3. To calculate such an imbalance and generate a graph the software was used HaploView.

EQTL Analysis

The SNPs that gave significant differences in the GWAs and the SNP rs4613908 were assessed as possible eQTLs using two approaches: using the GTEx project portal and using the genotypes of those SNPs in HapMap individuals and the levels of cis gene expression in the same individuals. Specifically, the SNPs genotypes were obtained from 210 unrelated HapMap phase 1 and 2 individuals and the levels of expression of the SHFM1, SLC25A13 and DLX5 genes from a lymphoblastoid cell line from the same individuals obtained.

Variants found and clues about its function

he genomic region of FLJ42280 (28 kb) was massively resequenced in two DNA pools corr

The genomic region of FLJ42280 (28 kb) was massively resequenced in two DNA pools corresponding to the 50 women with the highest BMD and the 50 women with the lowest BMD of the BARCOS cohort (see details in Material and Methods above) at a high depth (about 3,600x per group). We compared the number and frequency of the variants found in each group. A total of 110 variants were identified, of which 18 were new and 59 were rare or low frequency variants (Table 2). It was observed that the number of low frequency variants between the two extreme groups was balanced. Likewise, it was observed that the frequency differences of all variants were below the statistical power of the design, although 9 showed a tendency.

For each variant, its overlap was analyzed with functional elements annotated in the genome by the ENCODE project. Four of the variants overlapped with possible transcription enhancing sequences (or enhancers) of osteoblasts and one of them [SNP rs4613908; MAF(CEU)=0.39] overlapped with an active enhancer in osteoblasts (Figure 3).

Figure 3 A) Schematic of the genomic region of FLJ42280 and location of some promising variants. Common variants are indicated in red and low frequency variants (LFV) or rare in gray.Blue squares represent enhancers of osteoblasts; The darker represents an active enhancer of the osteoblasts. B) Detail of active enhancer in osteoblasts. The sequence conservation profile in vertebrates, transcription factor binding sites and histone modifications are shown. These data were extracted from UCSC Genome Browser -GRCh37- and from ENCODE (data referred to osteoblasts) 

Linkage disequilibrium analysis

The linkage disequilibrium among all variants common in this region was also studied. We plotted linkage disequilibrium (LD) using HaploView and haplotype information from the HapMap project (Figure 4) and noted that there is a large LD block that includes almost the entire gene (with the exception of the 3 'UTR region) And that by the upstream part of the gene extends 5 kb beyond the resequenced region. It was also found that the SNPs rs4613908 and rs4727338 (GWAs meta-analysis of Estrada et al. [⁷]) present a large linkage disequilibrium between them.

Figure 4 Graph of linkage disequilibrium present in the genomic region of FLJ42280.A large block (24 kb) is seen that includes much of the gene except its 3'UTR region. The 4 SNPs that have been associated with BMD in the various GWAs are marked with green ovals. The SNP overlapping with the active enhancer in osteoblasts is marked with a purple oval. The blue horizontal line indicates the resequenced genomic region in this study 

eQTL Analysis

To complete the functional analysis, an eQTL analysis was carried out. Having the genotypes of the four SNPs associated with BMD, and the SNP rs4613908 of 210 individuals from the HapMap project and the levels of gene expression of a genomic array in lymphoblastoid lines from these same individuals, we determined whether the different alleles or genotypes Of the SNPs correlated with the levels of gene expression of the genes located in the genomic region of FLJ42280. None of the SNPs showed influence on the expression levels of the SHFM1, SLC25A13 or DLX5 genes (in the array there is no information on expression levels of FLJ42280). We also accessed the GTEx database to collect eQTL information for the same SNPs and the result was negative for all of them. Finally, we searched for regulatory annotations in HaploReg. This latter analysis confirmed that the sequence surrounding SNP rs4613908 is highly conserved among mammals and that in several cell types, including primary osteoblasts, contains chromatin markers typical of enhancer sequences (H3K4me1, H3K27ac). On the other hand, HaploReg highlighted alteration of regulatory motifs of this SNP and rs10429035, but showed no effect of these SNPs on gene expression.