Identifying Environmental Endocrine Disruptors Associated With the Age at Menarche by Integrating a Transcriptome-Wide Association Study With Chemical-Gene-Interaction Analysis.

Menarche is the first occurrence of menstrual bleeding and one of the most important events of female puberty. Alarmingly, over the last several decades, the mean age at menarche (AAM) has decreased. Environmental endocrine disruptors (EEDs) are chemicals that may interfere with the endocrine system, resulting in adverse developmental, immunological, neurological, and reproductive effects in humans. Thus, the effects of EEDs on fertility and reproduction are growing concerns in modern societies. In this study, we aimed to determine the influence of genetic and environmental factors on AAM. We used data from an AAM genome-wide association study of 329,345 women to conduct a transcriptome-wide association study (TWAS) with FUSION software. As references, we determined the gene-expression levels in the hypothalamus, pituitary gland, ovaries, uterus, and whole blood. We performed Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses using the significantly dysregulated genes identified by the TWAS. Using the STRING database, we also generated a protein-protein-interaction network to analyze common AAM-specific genes identified by the TWAS with different tissues. We performed chemical-related gene set enrichment analysis (CGSEA) and identified significant TWAS genes to uncover relationships between different chemicals and AAM. The TWAS identified 9,848 genes; among these, 1580 genes were significant (P < 0.05), and 11 genes were significant among the hypothalamus, pituitary, ovary, uterus, and whole blood. CGSEA identified 1,634 chemicals, including 120 chemicals significantly correlated with AAM. In summary, we performed a TWAS (for genetic factors) and CGSEA (for environmental factors) focusing on AAM and identified several AAM-associated genes and EEDs. The results of this study expand our understanding of genetic and environmental factors related to the onset of female puberty.

Introduction

Puberty is a complex process occurring between childhood and adulthood, producing internal and external physical changes that promote the development of primary and secondary sexual characteristics important for sexual reproduction (1). Sex hormones are responsible for the physical manifestations of female puberty, including thelarche, pubarche, and menarche (2). Menarche is the first menstrual bleeding and one of the most important events of female puberty. Over the last several decades, the mean age at menarche (AAM) has declined, which is concerning (3). Epidemiological evidence suggests that the onset of puberty is advancing in humans through undetermined mechanisms (4).

Environmental endocrine disruptors (EEDs) are chemicals that interfere with the endocrine system and evoke adverse developmental, immunological, neurological, and reproductive effects in humans (5). EEDs are common in human living environments and include pesticides, plasticizers, industrial by-products, drugs, and some naturally occurring phytochemicals (6). These exogenous chemicals can interfere with the complex endocrine system, causing adverse health effects, such as reproductive disorders, metabolic diseases, and various cancers (7). Chronic exposure to EEDs may play a role in accelerating or delaying the onset of menarche, and extensive research has shown that pesticides, phenols, polycyclic aromatic hydrocarbons, phthalates, and some heavy metals are responsible for hormone metabolism disorders occurring during puberty (8). Consequently, the effects of EEDs on fertility and reproduction are a growing concern in modern societies (9).

Recent genome-wide association studies (GWASs) have identified thousands of genetic variants associated with complex phenotypes and have provided insights into their genetic architectures. GWASs are also extremely well-suited for identifying common single-nucleotide polymorphism (SNP)-based variants (10). GWAS analysis has been conducted to research early puberty and identify the genetic characteristics of idiopathic central precocious puberty and validate the polygenic risk for early puberty (11). Felix R Day et al. have identified 389 independent, genome-wide significant signals for AAM, which explained ~7.4% of the population variance in AAM and corresponded to ~25% of the estimated heritability (11).

Genetic loci cause trait variations, ranging from growth and fitness in simple organisms to disease in humans. Determining the genetics of gene-expression differences has emerged as a key approach for linking DNA-sequence variations to phenotypes (12). Transcriptome-wide association study (TWAS) analysis has been used to identify significant expression-trait associations by integrating genotypes, gene-expression levels, and phenotypes in order to gain insights into the genetic basis of complex traits (13). In a recent study, a TWAS was performed to discover transcriptome differences that affect the age of natural menopause (ANM), and 34 ANM-associated genes were reported (14). The hypothalamic–pituitary–gonadal axis controls puberty and reproduction (15). Therefore, investigating the associated gene-regulation relationships may help identify important genes that are co-expressed in all tissue types during AAM.

Here, we aimed to determine the influences of genetic and environmental factors on AAM by performing a large-scale TWAS for AAM based on a GWAS data set. We investigated gene-expression levels in the hypothalamus, pituitary gland, ovaries, uterus, and whole blood. We also reevaluated the expression of TWAS-identified genes, functionally explored the genes, and identified AAM-associated EEDs.

Methods

Summary of the AAM GWAS Data Used in This Study

We used published GWAS summary data for female AAM (16). Briefly, Day et al. performed a meta-analysis of pooled GWAS data from multiple studies of 329,345 women of European ancestry, including 40 studies from the ReproGen consortium (N = 179,117), 23andMe (N = 76,831), and the United Kingdom Biobank (N = 73,397). The results were grouped into the relatively sparse HapMap 2 reference panel or attributed to the 1000 Genomes Project reference panel using gene-centric arrays. In each study, the associations of SNPs with AAM were based on a two-tailed additive linear-regression model and several factors, including the age at the study visit and other study-specific covariates. Day et al. performed an expanded genomic analysis of AAM in women that was nearly three times larger and used denser genomic data than previous studies. Information on the subjects, genotypes, responsibilities, and quality control were detailed in the published study (16). In addition, ethical approval was not applicable for this study as publicly available data were used for the analysis.

TWAS Analysis

We used FUSION (13) software (http://gusevlab.org/projects/fusion/) to analyze the GWAS summary data for the previous meta-analysis of AAM. The most popular TWAS methods, such as PrediXcan, TWAS-Fusion, and SMR, test causal relationships between gene-expression levels and complex traits (17), among which, the TWAS-Fusion method is used more often. Briefly, Bayesian sparse linear-mixed models (18) were used to calculate SNP expression weights for specific genes at the 1-Mb cis position and estimate the association of predicted expression levels with AAM using the following formula: Ztwas = w + Z/(w × [Lw]1/2) (13), where w denotes the weight, Z denotes the Z-score, and L denotes the SNP correlation matrix (definition, LD). We used the gene-expression weights for the hypothalamus (N samples = 108; N features = 2,315), pituitary (N samples = 157; N features = 4,402), ovary (N samples = 122; N features = 2,809), uterus (N samples = 101; N features = 2,135), and whole blood (N samples = 1,264; N features = 4,701) as references, which are available for download from the FUSION website (http://gusevlab.org/projects/fusion/). We estimated the transcriptome-wide significance as P = 5.08 × 10−6 (0.05/9,848) using the Bonferroni correction (14). Manhattan plot was made by “CMplot”(v. 3.6.2) in R package.

Functional Exploration of Genes

We performed Kyoto Encyclopedia of Genes and Genomes (KEGG) (19) and Gene Ontology (GO) (20) enrichment analyses to identify and confirm related biological processes. KEGG and GO enrichment were performed using the R packages “org.Hs.eg.db” and “clusterProfiler” (R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/).

Interaction-Network Analysis

We generated a protein–protein-interaction (PPI) network using the STRING database, v11.5 (STRING, https://string-db.org), requiring a confidence score of 0.15 and “active interaction sources” based on a previous study (21). Cytoscape (22) was used to visualize all interaction networks, and the Molecular Complex Detection (MCODE) plugin (23) was used for module analysis.

Chemical Gene Expression Annotation Data Set

The chemical-related gene-expression annotation data set used in this study was downloaded from the Comparative Toxicology Genomics Database (CTD) (http://ctdbase.org/downloads/). The CTD mainly provides four data sets, including a chemical–gene-interaction function, a chemical–disease association, a genetic disease association, and a chemical element-phenotypic association. The CTD integrates the four data sets to automatically construct a hypothetical chemical–gene phenotype disease network to illustrate the molecular mechanisms underlying diseases that affect the environment (24). Cheng et al. downloaded and used 1,788,149 chemical–gene pair annotation terms for humans and mice, generating 11,190 chemical substance-related gene sets (25). We also used that data set to perform our chemical-related gene set enrichment analysis (CGSEA).

CGSEA

CGSEA is a flexible tool for assessing associations between chemicals and complex diseases. Briefly, the software uses genome-wide summary data (e.g., a summary of the TWAS data and messenger RNA [mRNA]-expression profiles) to explore functional relationships among chemical substances and diseases from a genomics perspective, for many complex diseases and characteristics. We used the CTD Chemical Gene Interaction Network and TWAS Expression Association to test the AAM statistics and weighted Kolmogorov–Smirnov running sum statistics to explore the relationships between chemicals and AAM, as described previously in greater detail (26). Specifically, in this study, we performed 10,000 permutations to obtain the empirical distribution of the gene set enrichment analysis (GSEA) statistical data for each chemical and then calculated the P-value of each chemical based on the empirical distribution of the CGSEA statistical data. Based on previous findings, we excluded gene sets containing <10 or >200 genes to limit the influence of gene set sizes on the results (27). To avoid deviations from expression correlations between genes, a package of lme4qtl software (28) was used to adapt the mixed-model regression of the TWAS Z-score based on the number of members in each gene set, in order to consider the correlation of the Z-score between genes caused by LD (28). A detailed description of the analytical method used was provided previously (25).

Results

TWAS Analysis of AAM

TWAS analysis identified 9,848 genes from the GWAS summary data and of those 1,580 genes expression was associated with AAM (P <0.05) while 64 showed a significant association (P <5.08 × 10−6); 2,289 genes tested in hypothalamus and of 322 genes expression was associated with AAM (P <0.05) while 8 showed a significant association (P <5.08 × 10-6; ), 4,362 genes tested in pituitary and of 568 genes expression was associated with AAM (P <0.05) while 26 showed a significant association (P <5.08 × 10−6; ), 2,768 genes tested in ovary and of 366 genes expression was associated with AAM (P <0.05) while 14 showed a significant association (P <5.08 × 10−6; ), 2,104 genes tested in uterus and of 253 genes expression was associated with AAM (P <0.05) while 7 showed a significant association (P <5.08 × 10−6; ), and 4,671 genes tested in whole blood and of 604 genes expression was associated with AAM (P <0.05) while 9 showed a significant association (P <5.08 × 10−6; ), respectively ().

Figure 1

Manhattan plots of the association results from the AAM TWAS. The dashed horizontal lines represent P = 5.00 × 10−2. The solid horizontal lines represent P = 5.00 × 10−6 (Bonferroni correction). Each dot represents the genetically predicted expression of one specific gene in the hypothalamus, pituitary gland, ovary, uterus, and whole blood tissues. The X axis represents the chromosome (Chr) encoding the corresponding gene, and the Y axis represents the negative logarithm of the association P value. (A) Gene-expression weights for the hypothalamus. (B) Gene-expression weights for the pituitary gland. (C) Gene-expression weights for the ovaries. (D) Gene-expression weights for the uterus. (E) Gene-expression weights for whole blood.

Functional Exploration of TWAS-Identified Genes Associated With AAM

Tissues have unique gene-expression profiles. Thus, we performed an overlap analysis of the significant genes in different tissues to identify the most representative genes. shows the resulting Venn diagram, which indicates the number of genes expressed in one or more tissues. Overall, 163 TWAS-identified significant AAM-specific genes were associated with the hypothalamus: 38 significant genes were associated with the hypothalamus and pituitary gland; 7 significant genes were associated with the hypothalamus, pituitary, and ovary; and 11 significant genes were associated with the hypothalamus, pituitary, ovary, uterus, and whole blood. The 11 novel TWAS-significant AAM-susceptible genes identified in all five tissues were RBM6 (RNA-binding motif protein 6; chromosome 3), PILRB (paired immunoglobin-like type 2 receptor beta; chromosome 7), CPSF1 (cleavage and polyadenylation-specific factor 1; chromosome 8), INPP5E (inositol polyphosphate-5-phosphatase E; chromosome 9), MRPL43 (mitochondrial ribosomal protein L43; chromosome 10), HSD17B12 (hydroxysteroid-(17-β)-dehydrogenase 12; chromosome 11), TIPIN (TIMELESS-interacting protein; chromosome 15), FLYWCH1 (FLYWCH-type zinc finger 1; chromosome 16), EXOSC6 (exosome component 6; chromosome 16), ADORA2B (adenosine A2b receptor; chromosome 17), and SPATA20 (spermatogenesis-associated 20; chromosome 17). presents detailed information regarding these 11 genes, including the rsIDs of the most significant (i.e., best) GWAS SNPs in the locus (i.e., BEST.GWAS.ID) and the TWAS P-value (i.e., P ).

Figure 2

Functional exploration of the TWAS-identified genes associated with AAM. (A) Venn diagram revealing the overlap of TWAS-significant genes in different tissues. Red, hypothalamus; green, pituitary gland; purple, ovaries orange, uterus; blue, whole blood. (B) Bar plot of enriched GO terms for the overlapping genes.

Table 1

Significant TWAS-identified genes associated with AAM in all five tissues studied.

Gene	BEST.GWAS.ID	P TWAS
		Hypothalamus	Pituitary gland	Ovary	Uterus	Whole blood
RBM6	rs3905330	0.00271	0.00009	0.00041	0.00041	0.00303
PILRB	rs2950520	0.02887	0.03400	0.02900	0.03800	0.03621
CPSF1	rs35253356	0.00528	0.01188	0.01056	0.02224	0.01915
INPP5E	rs10448340	0.03503	0.02700	0.01299	0.04100	0.01770
MRPL43	rs11190901	0.00014	0.00017	0.00030	0.00038	0.00070
HSD17B12	rs6485443	0.00272	0.00061	0.00316	0.00089	0.00368
TIPIN	rs2113688	0.00420	0.00902	0.00147	0.00420	0.02590
FLYWCH1	rs1834026	0.01500	0.01500	0.00772	0.01210	0.02630
EXOSC6	rs7196842	0.00019	0.00075	0.00003	0.00001	0.00002
ADORA2B	rs178837	0.01740	0.00027	0.00027	0.00027	0.00110
SPATA20	rs989128	0.02210	0.04920	0.03560	0.00639	0.00593

We subjected the TWAS-identified genes to GO analysis (). Five enriched GO terms belonged to the biological process (BP) category, including cell cycle DNA replication, long-chain fatty-acyl-CoA metabolic process, DNA replication checkpoints, estrogen biosynthetic processes, and positive regulation of immune effector processes. Three significantly enriched GO terms belonged to the cellular component (CC) category, including the mRNA cleavage factor complex, mRNA cleavage and polyadenylation specificity factor complex, and nuclear exosome (RNase complex). In terms of the molecular function (MF) category, the enriched GO terms primarily involved estrogen metabolism (such as steroid dehydrogenase activity or steroid dehydrogenase activity acting on the CH-OH group of donors), nicotinamide adenine dinucleotide or nicotinamide adenine dinucleotide phosphate as an acceptor, and estradiol 17-beta-dehydrogenase activity.

PPI Network of the TWAS-Identified Genes

We used 1,580 TWAS-significant AAM-associated genes for PPI analysis and successfully transformed 1,056 protein-coding genes (). To effectively find densely connected regions of the PPI network, we formed seven MCODE clusters with the PPI network genes (). The hub genes identified using the MCODE plugin were further analyzed for functional exploration. MCODE cluster 1 (MCODE1) was related to ribosome biogenesis, female pregnancy, and blastocyst development. MCODE2 was characterized by genes in the human leukocyte antigen (HLA) family associated with the immune process. MCODE3, MCODE4, and MCODE5 were related to tRNA, energy metabolism processes, and biosynthetic processes. MCODE6 was associated with lipid oxidation and maternal processes involved in female pregnancy. MCODE7 contained some important significantly enriched terms, including fatty acid metabolic, gonadotropin-releasing hormone (GnRH), estrogen, oxytocin, and PI3K-Akt signaling pathways and breast cancer.

Figure 3

PPI network and significant modules. Red and blue circles indicate upregulated and downregulated TWAS-identified genes. (A) PPI network of the TWAS-identified genes. (B-H) Significant modules 1 to 7 of the PPI network and their functional exploration. GeneRatio: the ratio of the number of genes associated with a term linked to an enriched gene to the total number of enriched genes.

CGSEA of the TWAS-Identified Genes

We performed CGSEA to investigate environmental factors that influence the onset of puberty, which identified 1,634 chemicals, including 120 chemicals that correlated significantly with AAM. These significant chemicals included some drugs (e.g., fluoxetine), pesticides (e.g., ametryne), plant extracts (e.g., isoflavones), nutrients (e.g., cholesterol), phenols (e.g., cannabidiol), phthalates (e.g., monobutyl phthalate), heavy metals (e.g., uranium), and atmospheric pollutants (e.g., phosgene). illustrates our constructed network of EEDs and their target genes based on the TWAS-identified genes. Chemical with absolute normalize enrichment score (NES) values of >1 were considered significantly enriched according to GSEA. We identified 77 significantly enriched chemicals with an |NES| value of >1 and a P-value of <0.05 ().

Figure 4

CGSEA results. Network of EEDs and their target genes, based on the TWAS-identified genes.

Table 2

Significant chemicals identified by CGSEA of the significant TWAS-identified AAM-associated genes.

Chemical name	NES	P	Chemical name	NES	P
Chlorine	12.35459	0.00050	Acetaldehyde	3.24524	0.01399
Fluoxetine	8.68168	0.00050	Heliotrine	2.50818	0.01449
Proton pump inhibitors	63.99450	0.00050	Pseudocumene	3.38592	0.01499
Fexofenadine	18.12249	0.00100	Chlordecone	2.98902	0.01499
Torcetrapib	7.28264	0.00150	Gefitinib	1.74158	0.01549
Methylmethacrylate	6.88844	0.00250	Cholesterol	2.32618	0.01599
Chromous chloride	8.16988	0.00350	N-nitroso-tris-chloroethylurea	2.92314	0.01649
Thalidomide	4.83989	0.00350	Methylmercury cysteine	1.12554	0.01799
Pentosan sulfuric polyester	5.33529	0.00400	Pyrene	2.66545	0.01799
Anacardic acid	4.53182	0.00450	Dehydroxymethylepoxyquinomicin	1.58757	0.01799
Calcimycin	5.92628	0.00450	Casticin	2.55138	0.01899
Palbociclib	2.05821	0.00550	Squalestatin 1	2.60788	0.01949
Temozolomide	4.66128	0.00550	Periodate-oxidized adenosine	2.80821	0.01999
Uranyl acetate	3.30552	0.00600	Isoflavones	2.19951	0.02049
Catechol	4.40109	0.00650	Melphalan	2.67401	0.02099
Cannabidiol	3.31796	0.00650	Clomipramine	1.57621	0.02149
Amitriptyline	4.15778	0.00700	Uranium	2.49394	0.02199
Ibuprofen	3.81067	0.00700	Monobutyl phthalate	2.52550	0.02249
Uric acid	2.27066	0.00700	Rimonabant	2.71940	0.02349
Methylparaben	5.83928	0.00750	Toxaphene	1.21501	0.02399
Oleoyl-estrone	1.50555	0.00750	Ursodeoxycholic Acid	2.44936	0.02399
2-(4-nitrophenyl)-4-(4-fluorophenyl)-5-(4-pyridinyl)-1H-imidazole	1.72254	0.00750	Letrozole	1.83199	0.02449
Benzbromarone	3.93596	0.00750	4-toluidine	2.47743	0.02499
Ethionine	4.06789	0.00750	Carcinogens	2.27375	0.02649
SK&F 83959	3.56873	0.00800	Okadaic acid	2.35105	0.02649
Potassium perchlorate	3.55962	0.00850	Sulfasalazine	1.98058	0.02749
Phenanthrene	3.09849	0.00900	Pyrazolanthrone	2.25782	0.02899
SNX 2112	1.69629	0.00900	Pravastatin	1.71224	0.03098
Aripiprazole	1.82475	0.00950	Piperonyl butoxide	2.13714	0.03298
Asbestos	3.88868	0.00950	Potassium chloride	1.23405	0.03298
S-Adenosylmethionine	2.34165	0.00950	Bucladesine	1.80893	0.03498
Reserpine	2.09526	0.01049	3-(2-hydroxy-4-(2-methylnonan-2-yl) phenyl) cyclohexan-1-ol	2.16625	0.03598
Ciprofibrate	3.59698	0.01099	Caffeic acid	1.43915	0.04248
Probenecid	2.54975	0.01099	Methimazole	1.97181	0.04248
Pyrogallol	3.56596	0.01099	Methyl cellosolve	1.86438	0.04298
Homocysteine	3.54662	0.01149	Estrogens	1.47629	0.04448
Phosgene	3.01541	0.01199	Ciglitazone	1.26701	0.04598
Tebuconazole	2.51505	0.01349	Hyaluronic acid	1.84193	0.04698
Paraoxon	2.58233	0.01349

Discussion

Puberty is a developmental period affecting the body and behavior of children developing secondary sexual characteristics (29). Pubertal hormones play important roles in the adrenal, gonadal, and growth axes (30). Pubertal timing is determined by the activation of pulsatile hypothalamic GnRH secretion when pituitary gonadotropin secretion and downstream gonadal maturation are also initiated (31). Thus, various tissue types are affected by AAM, including the hypothalamus, pituitary gland, ovaries, uterus, and whole blood.

Up to now, previous studies have identified several genes associated with AAM, these results are similar to our studies (). Here, we performed a comprehensive TWAS to evaluate the relationship between AAM and predicted genes found in the hypothalamus, pituitary gland, ovaries, uterus, and whole blood. We identified 11 genes whose genetically predicted expression was associated with AAM (P ≤ 0.05), including 10 novel genes (RBM6, PILRB, CPSF1, INPP5E, MRPL43, TIPIN, FLYWCH1, EXOSC6, ADORA2B, and SPATA20) and one gene (HSD17B12) linked to AAM in a previous GWAS. The 11 protein-coding genes identified in our study have been implicated in estrogen metabolism. We also identified enriched fatty acid-metabolism pathways, similar to the findings of a recent study (36), indicating that different types of fatty acids may influence puberty timing. Additionally, individual fatty acids might have different physiological and metabolic effects, including the progression of pubertal development. Therefore, our study provides new information that improves our understanding of the genetics and etiology of AAM.

Table 3

Previously reported significant AAM-related genes.

Authors	SNP	Gene	Pa	PTWASb
John RB Perry, et al. (32)	rs10148448	MEG3	2.10E-04	2.52E-07
Diana L. Cousminer, et al. (33)	rs12917823	MAPK3	P < 5.00E-02	1.22E-03
	rs17046434	ADCY3	P < 5.00E-02	2.42E-02
	rs7759938	LIN28B	P < 5.00E-02	6.37E-49
Cathy E. Elks, et al. (34).	rs7759938	LIN28B	1.60E-58	6.37E-49
	rs10148448	BEGAIN	1.70E-10	2.55E-06
	rs4929923	TRIM66	2.40E-08	1.38E-09
	rs3905330	RBM6	1.40E-09	8.60E-05
	rs633715	SEC16B	1.50E-09	4.07E-03
	rs4955420	KLHDC8B	1.80E-09	6.13E-06
	rs16938437	PHF21A	1.40E-09	9.98E-09
	rs2687729	EEFSEC	1.00E-08	2.51E-11
	rs8104651	OLFM2	4.60E-10	1.10E-07
	rs939317	ECE2	2.30E-09	4.89E-03
	rs12641981	GNPDA2	8.70E-08	3.85E-02
	rs12056794	MSRA	2.40E-02	1.00E-02
Nicholas Mancuso, et al. (35).	rs6580698	CCDC65	P < 5.00E-08	7.30E-05
	rs7330016	COG6	P < 5.00E-08	4.40E-02
	rs12917823	INO80E	P < 5.00E-08	2.47E-05
	rs3761919	NUCKS1	P < 5.00E-08	4.86E-08
	rs4717903	PMS2P5	P < 5.00E-08	6.58E-04
	rs3761919	RAB7L1	P < 5.00E-08	9.47E-07
	rs4717903	STAG3L2	P < 5.00E-08	4.37E-08
	rs2274351	TMEM180	P < 5.00E-08	2.24E-04

P values of previously reported AAM-related genes.

Significant TWAS-identified genes associated with AAM in present research.

HSD17B12 encodes 17β-hydroxysteroid dehydrogenase (17β-HSD), which is crucial for converting estrone into estradiol and fatty acid elongation. Our results support those of Kemiläinen et al. (37), who found that HSD17B12 played a vital role in female fertility through arachidonic acid metabolism. Specifically, an investigation of the 17β-HSD enzyme expression in human and mouse ovaries revealed that female HSD17B+/- mice more often had frequent dysfunctional oogenesis and ovulation, leading to less frequent births. Under extreme circumstances, haploinsufficiency of the HSD17B12 gene in female mice resulted in subfertility. HSD17B12 can affect fatty acid elongation and ceramide accumulation in the serum (38). Previous data showed that central ceramides participate in the timing of female puberty (39), suggesting that HSD17B12 could be a candidate gene for female puberty. Notably, HSD17B12 was also related to endometrial (40), ovarian (41), and breast (42) cancers.

Previous results have shown that INPP5E is concentrated at the cilia base, where it helps control phosphoinositide metabolism (43). Primary cilia are present on virtually all cell types. Cilia appear to exert a crucial modulatory role in appropriate axonal wiring due to INPP5E-dependent activation of the PI3K–AKT signaling pathway, triggering an axonal Ca2+ wave (44). GnRH neurons in adults are multi-ciliated, and the percentage of GnRH neurons possessing multiple Kiss1r-positive cilia increases during puberty, correlating with sexual maturation (45). ADORA2B encodes an adenosine receptor and is involved in axon elongation. A previous report showed that ADORA2B transcripts were significantly downregulated in GnRH neurons during proestrus (46), while it has not been studied in pituitary. Therefore, these newly defined genes may influence the electrophysiology of puberty in hypothalamus and pituitary gland.

TIPIN is part of the replisome complex and binds the replication fork-protection complex TIMELESS, which is involved in circadian rhythm regulation (47). Circadian rhythms are well known to play key roles in animal reproduction (48). The existing body of research suggests that during the chronotype and the circadian timekeeping system change during the puberty, including endogenous rhythm period and sensitivity to environmental time cues (49). The results of several studies have established that melatonin promotes follicle-stimulating hormone in the pituitary gland and increases serum estrogen levels, thereby accelerating the onset of puberty (50). Those studies indirectly support a potential role of TIPIN in AAM.

CPSF recognizes the AAUAAA signal in pre-mRNA and interacts with other factors to facilitate RNA cleavage and the poly(A) addition. CPSF1 is the largest subunit of the CPSF complex. Previous data established that CPSF1 may promote ovarian cancer (51), cell proliferation, and triple-negative breast cancer (52). However, women have a higher risk for developing autoimmune diseases than men, which is attributed to sweeping endocrinological changes during puberty that considerably affect the immune system (53). CPSF1 was previously found to have a high affinity for HLA molecules (54). Thus, CPSF1 may be linked to the immune response in adolescent girls.

We extended the well-established GSEA approach to detect associations between environmental chemicals and AAM using published GWAS summary data sets and identified 120 chemicals, including drugs, organic compounds, inorganic compounds, plant extracts, nutrients, phenols, plasticizers, pesticides, herbicides, pollutants, and heavy metals. Aromatase is a rate-limiting enzyme in the conversion of androgens to estrogens, and letrozole is a selective aromatase inhibitor (55). Letrozole has been used to treat McCune–Albright syndrome, which has been associated with precocious puberty in girls (56) and delayed growth and puberty in boys (57). Further, fluoxetine is drug that is commonly used to treat adults with depressive disorders (58). The effect of fluoxetine on puberty has been a controversial issue. The results of one study demonstrated that fluoxetine exposure in utero delays puberty onset in female rats (59), whereas other data suggested that fluoxetine exposure during gestation did not alter plasma estrogen concentrations in peripubertal offspring (60).

Some plant extracts act as EEDs, as they are among a group of secondary metabolites with chemical structures similar to those of endogenous hormones. Isoflavones are phytoestrogens that are mainly produced in soybeans and can promote advancement of the vaginal opening in female rats after exogenous supplementation (61). Galangin is a naturally occurring flavonoid that inhibits the effects of flavonoids on human cytochrome P450 (62). Galangin was found to inhibit the aryl hydrocarbon receptor and is considered a potential drug for treating breast cancer (63, 64). Notable health benefits have also been reported for galangin (65), and previous data suggest that galangin is useful for treating precocious puberty. Zearalenone (ZEN) is a non-steroidal mycoestrogen that can exert adverse endocrine effects in mammals (66). To date, there is little agreement regarding the effect of ZEN on puberty. One group found that ZEN might trigger central precocious puberty (CPP) development in girls (67), whereas another group found an association between ZEN and normal pubertal development in adolescent girls (68).

A growing body of literature has associated being overweight or obese with early puberty (69); overweight and obese girls undergo menarche earlier than normal-weight girls (70). Several cross-sectional data have suggested that a substantial proportion of girls with CPP have hypertriglyceridemia (71). Cholesterol is an endogenous ligand of estrogen-related receptor alpha (72) that plays a role in activating estrogen receptors (73). Based on these data, we can infer that a high-cholesterol diet may also be a risk factor for precocious puberty. Plasticizers are additives used to produce or promote plasticity and flexibility in plastics, which are commonly used in everyday life. Further, parabens have been added to personal care products as antimicrobial preservatives (74), and methylparaben exposure increases glandular tissue sizes during critical developmental windows (75). A large survey showed that elevated phthalate metabolite levels pose potentially high health risks to children (76). In addition, monobutyl phthalate exposure might be associated with a risk for sexual precocity in girls (77). Animal experiments have demonstrated that two alkylphenols (4-nonylphenol and 4-tert-octylphenol) modify hormone biosynthesis and delay the onset of puberty (78).

This study has some limitations. First, the pooled GWAS data were obtained from the United Kingdom Biobank, and the study subjects were predominantly from European populations; these aspects affect the extrapolation of the results. Therefore, our results should be used cautiously when studying AAM in other populations. Second, although some of the genes we screened have been confirmed in other studies, our enrichment analysis used genes that did not pass correction for multiple testing (P <0.05), which may introduce some false positives so it requires caution when using our results. And some genes related to AAM susceptibility identified here have not been verified via molecular biology experiments, which should be performed in future studies. Further, some EEDs identified in this study were previously demonstrated to play a role in AAM, whereas others are not yet validated, which will require more clinical observations and cohort studies. However, to the best of our knowledge, this is the first large study using CGSEA to identify candidate EEDs related to AAM. Our TWAS analysis detected AAM-associated genes at the DNA level, and our CGSEA extended the classic GSEA approach to detect associations between environmental chemicals and AAM.

Conclusion

In this study, we aimed to determine the influences of genetic and environmental factors on AAM. Therefore, we performed a TWAS and CGSEA related to AAM and identified multiple AAM-associated genes and EEDs. The results of this study expand our understanding of the genetic and environmental factors affecting the timing of female puberty.

Data Availability Statement

The original contributions presented in the study are included in the article/. Further inquiries can be directed to the corresponding author.

Author Contributions

Author ML and RF collected and processed the data, as well as wrote this article. YQ and HD provided language help and writing assistance. BL proofreaded the article. CY and YX designed the study. All authors contributed to the article and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.