{"paper_id":"b8f5a995-45d6-4f58-b547-5735f0b93695","body_text":"Age at menarche (AAM), the onset of the first menstrual period in girls, is considered as a landmark of female pubertal development. Menarche generally occurs after a series of complex neuroendocrine events leading to full activation of the hypothalamic-pituitary-gonadal axis  [1] . Menarche is associated with physical, emotional, and social development  [2] . In addition, AAM was shown to be associated with the risk of various diseases. Early AAM is reported to be one of the significant risk factors for depression  [3] , eating disorders  [4] , obesity  [5] , diabetes  [6] , breast cancer  [7] , and coronary heart disease  [8] . On the other hand, late AAM has been associated with osteoporosis  [9]  and taller adult stature  [10] . Therefore, the identification of loci contributing to variation in AAM could lead to a better understanding of a wide range of phenotypes.\nAAM is known to be a complex trait determined by an array of genetic and environmental variables  [11] – [13] . Twin and family studies suggest a significant genetic contribution to AAM with a heritability of more than 50%  [12] ,  [14] . Several genetic variations within candidate genes such as the estrogen receptor genes ( ESR1  and  ESR2 )  [15] ,  [16] ,  CYP19A1 \n [17] , and the  SHBG  gene  [18]  were shown to be associated with AAM. To date, a number of genome-wide linkage analyses  [14] ,  [19] ,  [20]  and genome-wide association studies (GWAS) [21] – [25]  for genes underlying variation in AAM have been performed. In 2009, the association of genetic variations in  LIN28B  with AAM was identified by four independent groups. Currently, more than 30 loci have been shown to be significantly associated with AAM. However, most of these studies were conducted using women of European ancestry. Here we performed a large scale meta-analysis of GWAS using more than 15,000 Japanese female samples.\n\nA total of 15,495 Japanese female subjects from four GWAS using different SNP genotyping systems were enrolled in this analysis. Characteristics of samples and genotyping methods are summarized in  \n Table 1 \n . All the subjects were of Japanese origin and obtained from the Biobank Japan Project  [26] . Samples consist of patients that were classified into 33 disease groups. The average and S.D. of AAM in each disease cohort is shown in  Table S1 . Some diseases such as breast cancer and osteoporosis are likely to be associated with early or late AAM, as reported previously  [7] ,  [9] . We also found that AAM was negatively associated with birth year (p<0.0001). Thus we used disease status and birth year as covariates in this study. Genotyping was performed with over 500,000 SNP markers using Illumina HumanHap 550 Genotyping BeadChip, Illumina610-Quad Genotyping BeadChip, or Illumina Omni Express (Illumina, CA, USA). We applied stringent quality control criteria as mentioned in the methods section. We also conducted principal component analysis  [27]  to evaluate potential population stratification. To extend the genomic coverage and conduct meta-analysis, we subsequently performed a whole-genome imputation of the SNPs, using HapMap Phase II genotype data  [28] . After the imputation, we performed SNP quality control (minor allele frequency≥0.01 and an imputation score ( Rsq  value by MACH software  [29] )≥0.7) and found that more than two million autosomal SNPs satisfied these criteria.\nThe associations of these imputed SNPs with AAM were evaluated using a linear regression model  [30]  and meta-analysis. Quantile-Quantile plots of P-values indicated the Inflation factors to be as low as 1.039 ( \n Figure 1 \n ), suggesting no substantial population stratification existed in our study population. Although we could not identify significant association that satisfied the genome-wide significant threshold (P<5×10−8), SNP rs364663 which is located within intron 2 of the  LIN28B  gene at 6q21 indicated the strongest association with a P-value of 5.49×10 −7  ( \n Table 2 \n ,  \n Figure 2 \n ). SNP rs364663 exhibited the association with AAM in all four cohorts without significant heterogeneity in both effect sizes and directions (P-value for heterogeneity = 0.29;  \n Table 3 \n ), suggesting that the observed associations at  LIN28B  are not the result of false-positives due to study-specific bias. Regional p-value plots indicated that all of the AMM-associated SNPs were clustered around the  LIN28B  locus ( \n Figure 2 \n ). Previously reported SNPs rs314263 and rs7759938 near the  LIN28B  locus  21 , 24  also associated with AAM (P = 1.03×10−6 and 1.12×10−6, respectively) in the same direction. In addition to the 6q21 locus, 17 SNPs in seven genomic regions exhibited suggestive association with a p-value of<1×10−5 ( \n Table 2 \n ). To further investigate the physiological role of these loci, the associations of these variations with body mass index (BMI) and height were evaluated using previously published results in 26,620 Japanese subjects  [31] . SNP rs9404590 near the  LIN28B  locus associated with height with p-value of 0.0003 ( Table S2 ), but we did not find significant association (P<0.01) between these loci and BMI ( Table S2 ).\nA total of 15,495 female samples were analyzed in this study.\nUpper panel;  P  values of genotyped SNPs are plotted (as −log 10  values) against their physical location on chromosome 6 (NCBI Build 36). Estimated recombination rates from HapMap JPT shows the local LD structure. Inset; Colors of other SNPs indicate LD with rs2596542 according to a scale from  r \n 2  = 0 to  r \n 2  = 1 based on pair-wise  r \n 2  values from HapMap JPT. Lower panel; Gene annotations from the University of California Santa Cruz genome browser.\nGenotyping result of 15,495 Japanese subjects were anlayzed in this study. Imputed SNPs with R2 of less than 0.7 were excluded from this analysis. A1 frequency of JPT was those from release 24 Hapmap JPT.\nEffect size and SE of allele1 on age at menarche (year per allele) and P-values were obtained by inverse-variance method.\nEffect size and S.E. of allele1 on age at menarche (year per allele) and P-values were obtained by inverse-variance method.\nSince AAM is associated with various disease risks, we conducted separate analyses for each disease and then performed meta-analysis for the top 42 loci. As a result, some SNPs showed slightly stronger association, but none cleared the genome wide significant threshold ( Table S3 ). Similar to the current study, our group had previously conducted QTL analyses using disease status as a covariate and successfully identified many QTL loci  [31] – [35] . Therefore, different background due to disease status was unlikely to significantly affect the result of association analysis.\nAdditionally, we examined the loci already known to show significant association with AAM in women of European ancestry  [22] – [25] ,  [36] . We selected 37 SNPs for this candidate analysis and successfully obtained the genotyping results of 33 SNPs ( \n Table 4 \n ), and the risk allele was consistent with previous reports for 31 SNPs. In addition, eight SNPs in or near  RXRG ,  CCDC85A ,  LOC100421670 ,  TMEM38B ,  ZNF483 ,  ARNTL , and  CA10  indicated possible associations with AAM (P<0.05). Among them, rs4452860 and rs7028916 at the  TMEM38B  locus exhibited significant association even after Bonferroni's correction (P<0.0015 = 0.05/33). Taken together, these variations as well as  LIN28B  are likely to be common AAM loci, although their effect sizes are different between women of European ancestry and those of Japanese.\nGenotyping result of 15,495 Japanese subjects were anlayzed in this study. Imputed SNPs with R2 of less than 0.7 were excluded from this analysis. A1 frequency of JPT were those from release 24 Hapmap JPT. N.D.; no data. References: 1 Elks et al Nat Genet 2010, 2 He et al Nature Genet 2009, 3 Liu et al Plos Genet 2009. P-value of 0.0015 (0.05/33) was set at the significant threshold for this candidate analysis.\n:Effect size and S.E. of allele1 on age at menarche (year per allele) and P-values were obtained by inverse-variance method.\nConcordance of association direction between this study and the previous report.\n\nAAM is a complex trait that is influenced by both genetic and environmental factors. In recent years, genome wide association analyses have become a standard method to identify genetic factors related with various diseases and phenotypes. In 2002 our group performed the first GWAS for myocardial infarction and successfully identified  LTA  as a disease susceptibility gene. Using this method, we have identified a number of loci associated with various phenotypes and common diseases  [37] – [42] .\nThrough the analysis of more than 15,000 Japanese female subjects from the Biobank Japan Project, we here report the association of  LIN28B  with AAM. In our analysis, SNP rs364663 within intron 2 of the  LIN28B  gene exhibited the strongest association. Previously reported SNPs rs314263, rs7759938, rs314280, and rs314276 near or in the  LIN28B  gene also showed equivalent association with p-values of 1.03−1.57×10−6. Although the directions of associations were consistent between women of European ancestry and those of Japanese, the effect sizes among Japanese (0.085–0.087) were less than those among European ancestry (0.09–0.14). Thus  LIN28B  is a common AAM loci, but its effects are relatively small among the Japanese population.\nThe lin-28 gene was originally identified through C. elegans mutants showing abnormality in developmental timing  [43] . Deleterious mutations in lin-28 produce an abnormal rapid tempo of development through larval stages to adult cuticle formation  [43] . Two mammalian homologs, Lin28a and Lin28b, possess similar biochemical activities  [44] ,  [45] . Transgenic mice expressing Lin28a exhibited increased body size, crown-rump length and delayed onset of puberty due to increased glucose metabolism and insulin sensitivity  [46] . In humans, both  LIN28B  and  LIN28A  control preprocessing of the let7 microRNA family  [45] .  Lin28B  is highly expressed in the majority of human hepatocellular carcinomas and embryonic stem cells  [47] , and regulate cell pluripotency  [48]  as well as cancer growth  [47] . However the role of genetic variations in the  LIN28B  locus should be investigated in the future study.\nAmong the seven suggestive loci indentified in our GWAS,  ZFHX4  and  NKX2.1  loci were previously shown to be associated with height in the Caucasian population  [49] ,  [50] . These loci also exhibited association with height among the Japanese population (p = 0.053 and 0.023, respectively). The  NKX2.1  gene encodes a thyroid-specific transcription factor that was shown to be associated with hypothyroidism  [51] . Since delayed pubertal development is a common manifestation of hypothyroidism, variation at  NKX2.1  locus might affect AAM through the regulation of serum thyroid hormone level. The  ZFHX4  gene encodes a homeodomain-zinc finger protein.  ZFHX4  mRNA is highly expressed in brain, liver, and muscle, however the molecular mechanism whereby variations within this genetic locus affect AAM is not clear. In addition,  KIRREL3  locus was shown to be associated with breast size in women of European ancestry  [52] . The  KIRREL3  gene encode a member of the nephrin-like protein family.  KIRREL3  was highly expressed in brain and kidney, and shown to be associated with mental retardation autosomal dominant type 4  [53]  and neurocognitive delay associated with Jacobsen Syndrome  [54] . Taken together, these three loci seem to be common development traits in women of both European ancestry and Japanese.\nIn our candidate analysis, the association of SNPs rs4452860 and rs7028916 on 9q31.2 with AAM was also validated. This locus was repeatedly shown to be associated with AAM  [25] ,  [55] as well as serum ALT (alanine transaminase) level [56] . SNPs rs4452860 and rs7028916 were located more than 300 kb away from the  TMEM38B  gene. In mice,  TMEM38B  is strongly expressed in brain, and  TMEM38B  deficient mice is neonatal lethal [57] . In humans, a homozygous mutation of  TMEM38B  was associated with autosomal recessive osteosgenesis imperfecta [58] . These findings also suggest an important role of TMEM38B in developmental processes.\nIn this study, we observed early AAM among breast cancer patients compared with overall samples (12.99 vs 13.44,  Table S1 ). On the other hand, patients with osteoporosis exhibited late AAM (14.33 y.o.). These findings are concordant with previous reports. Since we do not have age-matched controls, the association of these diseases and AAM should be examined using more appropriate sample sets.\nAlthough more than 15,000 samples were employed in this study, no SNP cleared the genome wide significance threshold (p<5×10 −8 ). This could be explained by some limitations in our study. Firstly, the information about AAM was self-reported and might not be very accurate. Although this distinct event is often well recalled many years later  [59] , possible error due to age-associated memory impairment would increase the risk of false negative association. Secondly, we did not have a replication sample set, and the total sample size may be inadequate to detect variations with modest effects. At present, we unfortunately do not have sufficient samples for replication, since more than 2,900 additional samples are necessary to obtain a genome wide significant association for the top SNP rs364663. Previous studies comprising up to 17,510 women detected only one or two genome-wide significant signals  [22] – [25] . However, 30 significant loci were identified by using 87,802 subjects in the screening stage  [21] . Although no SNPs reached the genome wide significance in our study, analyzing an increased number of subjects in a relatively younger generation would improve statistical power as well as data accuracy, and subsequently enable us to identify other genetic factors with modest effects.\n\nThe subjects enrolled in the GWAS meta-analysis for age at menarche (AAM) (n = 15,495) consisted of female patients that were classified into 33 disease groups; cancer (lung, esophagus, stomach, colon, liver, cholangiocarcinoma, pancreas, breast, uterine cervix, uterine body, ovary, hematopoietic organ), diabetes, myocardial infarction, brain infarction, arteriosclerosis, arrhythmia, drug eruption, liver cirrhosis, amyotrophic lateral sclerosis, osteoporosis, fibroid, and drug response, rheumatoid arthritis. chronic hepatitis B, pulmonary tuberculosis, keloid, heat cramp, brain aneurysm, chronic obstractive lung disease, glaucoma, and endometriosis ( \n Tables 1 \n \n  and  \n \n 2 \n ). All subjects were collected under the support of the BioBank Japan Project  [26] , in which the individuals with any one of the 47 common diseases were enrolled between 2003 and 2008. Subjects under 17 years of age were not included in this study. Various life style information such as AAM was obtained through a face-to-face interview by trained medical coordinators using questionnaire sheets at each hospital. All participants provided written informed consent as approved by the ethical committees of each institute. For participants between the age of 17 and 20 years old, we obtained written informed consent from both the participant and her parents. AAM was collected by self-report on the questionnaire. The subjects with AAM between the ages of 10 and 17 were enrolled for the study. This project was approved by the ethical committees of the BioBank Japan Project  [26]  and the University of Tokyo.\nIn the four GWAS enrolled in the meta-analysis of AAM, all samples were genotyped at more than 500,000 loci using one of the following platforms: Illumina HumanHap550v3 Genotyping BeadChip, Illumina HumanHap610-Quad Genotyping BeadChip, or Illumina Omni express Genotyping BeadChip ( \n Table 1 \n ). Then we applied the following quality control for each GWAS separately: exclusion criteria; subjects with call rates<0.98, SNPs with call rates<0.99 or with ambiguous clustering of the intensity plots, or non-autosomal SNPs. We then excluded subjects whose ancestries were estimated to be distinct from East-Asian populations using principle component analysis performed by EIGENSTRAT version 2.0  [27] . Subsequently, the SNPs with MAF<0.01 or P-value of the Hardy-Weinberg equilibrium test<1.0×10 −7  were excluded.\nAfter the quality control criteria mentioned above were applied, genotype imputation was performed by MACH 1.0.16  [29]  using the genotype data of Phase II HapMap JPT and CHB individuals (release 24) [28]  as references, in a two-step procedure as described elsewhere  [60] . In the first step of the imputation, recombination and error rate maps were estimated using 500 subjects randomly selected from the GWAS data. In the second step, imputation of the genotypes of all subjects was performed using the rate maps estimated in the first step. Quality control filters of MAF≥0.01 and  Rsq  values≥0.7 were applied for the imputed SNPs.\nAssociations of the SNPs with AAM were assessed by linear regression assuming the additive effects of the allele dosages on AAM, using mach2qtl software  [29] . Years of birth and affection statuses of the diseases were used as covariates. Meta-analysis of all four GWAS was performed using an inverse-variance method from the summary statistics of beta and standard error (SE), using the Java source code implemented by the authors  [61] . Genomic control correction was applied for each GWAS separately, and applied again for the results of GWAS meta-analysis. We set the P - values of 5.0×10 −8  and 0.0015 ( = 0.05/33, Bonferroni's correction based on the numbers of the evaluated loci) as significance thresholds in GWAS and candidate gene analysis, respectively. Heterogeneity of the effect sizes among the studies was evaluated using Cochran's Q statistics.\nThe associations of 42 SNPs with BMI and height were evaluated using previously published results, in which a total of 26,620 subjects with 32 diseases from Biobank Japan were enrolled  [31] . Of the 26,620 subjects, 12,350 were also included in this AAM study. Genotyping was performed using the Illumina HumanHap610-Quad Genotyping BeadChip. Genotype imputation was performed using MACH 1.0. BMI was calculated based on self-reported body weight and height data. A rank-based inverse-normal transformation was applied to the BMI values of the subjects. Associations of the SNPs with transformed values of BMI were assessed by linear regression assuming additive effects of allele dosages (bound between 0.0 and 2.0) using mach2qtl software using gender, age, smoking history, the affection statuses of the diseases and the demographic classifications of the medical institutes in Japan where the subjects were enrolled were used as covariates.\nThe URLs for data presented herein are as follows. BioBank Japan Project,  http://biobankjp.org\nMACH and mach2qtl software,  http://www.sph.umich.edu/csg/abecasis/MACH/index.html  International HapMap Project,  http://www.hapmap.org  PLINK software,  http://pngu.mgh.harvard.edu/~purcell/plink/index.shtml  EIGENSTRAT software,  http://genepath.med.harvard.edu/~reich/Software.htm \n R  statistical software,  http://cran.r-project.org .\n\nAge at menarche in each disease cohort.\n(DOCX)\nClick here for additional data file.\nThe result of association analysis of TOP 42 SNPs from Japanese AAM study with body mass index and height from previous GWAS by Okada et al.\n(DOCX)\nClick here for additional data file.\nThe result of top 42 SNPs by separated analysis.\n(DOCX)\nClick here for additional data file.","source_license":"CC-BY-4.0","license_restricted":false}