Human chromatin remodeler cofactor, RNA interactor, eraser and writer sperm RNAs responding to obesity

ABSTRACT

In the United States almost 33% of adults are considered obese (BMI > 30 kg/m²). Both animal models and to a lesser extent human studies, have associated BMI, a measure of obesity, with alterations in sperm DNA methylation and RNAs. Sperm RNAs from the Assessment of Multiple Gestations from Ovarian Stimulation trial, were isolated and sequenced. A Generalized Linear Model identified 487 BMI associated human sperm RNA elements (short exon-sized sequences). They partitioned into four patterns; a continual increase with BMI, increase once obese (BMI>30 kg/m²); a steady decrease with BMI; and decrease once overweight (BMI 25 - 30 kg/m²). Gene Ontology revealed a unique relationship between BMI and transcripts associated with chromosome organization, adipogenesis, cellular stress and obesity-related inflammation. Coregulatory networks linked by Chromatin remodeler cofactors, RNA interactors, Erasers and Writers (CREWs) were uncovered to reveal a hierarchical epigenetic response pathway.

KEYWORDS:

• BMI
• sperm RNA
• obesity
• CREWs - chromatin remodeler cofactor RNA interactor eraser and writer RNAs

Introduction

Obesity, commonly measured as Body Mass Index (BMI), is associated with the pathogenesis of a variety of diseases including diabetes and cardiovascular disease [1–3]. More recently, obesity and male infertility have been associated [1,4,5]. Established hallmarks and associations of obesity including hormonal changes, metabolic and cellular stress, and inflammation have been suggested as endpoints [4]. Epididymal fat pads are a type of white adipose tissue (eWAT) surrounding the epididymis and testis [6]. As with other WAT depots, eWAT responds to changes in nutrition and is highly regulated [6,7]. Uniquely, eWAT, is important for spermatogenesis [8] and may serve as a source of obesity-related sperm signaling.

Infertility impacts approximately 15% of couples worldwide [5,9] with approximately one third of all cases attributed solely to the male partner [1]. Current sperm diagnostic tests rely on the assessment of standard semen parameters (such as sperm concentration and motility) and in some instances, DNA analysis (such as DNA fragmentation index and high DNA stainability) [10,11]. Associations between BMI and both sperm DNA fragmentation and common semen parameters, including semen volume, sperm concentration and total sperm count, have also been suggested [4,5,12–15]. Highlighting the importance of the paternal contribution to a child’s health, these have been extended to include paternal obesity associations with birth success, birth weight and a child’s phenotypic traits such as metabolic health and over-all longevity [14,16–21]. Animal models have identified alterations in sperm DNA methylation and RNAs in response to high fat and protein diets [21–24], some of which are inherited by the next generation resulting in metabolic disorders [23,25–27]. Studies by Ortega et al. [28] identified a number of circulating miRNAs altered in response to obesity, while Donkin et al. [29] identified altered DNA methylation and small RNAs in sperm when lean and obese men were compared.

We have previously identified sperm RNAs required for live birth [30]. Gene ontology along with developmental and disease associations of sperm RNAs and DNA methylation analyses suggest a potential impact on child health and development [29,31,32]. We sought to define the relationship between obesity as measured by BMI and human sperm RNAs to expand our understanding garnered from animal models of obesity [33,34] and other human studies [4,29,35]. Here we present sperm RNAs responding to BMI as 8 co-regulatory networks, each representative of obesity-related adipogenesis, cellular stress, and inflammatory pathways. Chromatin remodelers, RNA interactors, Erasers and Writers (CREWs) are likely regulators driving the network.

Results

Subset demographics

Paternal obesity is associated with live birth, birth weight, a child’s metabolic health, and other phenotypic traits [16,18–20] playing a significant role in a child’s life course. Mouse models [23,24,26,34], and to a lesser extent human studies [28,29], have suggested a correlation between obesity and DNA methylation and some small sperm RNAs. Longer RNAs, i.e., those above 100 nucleotides including mRNAs have remained an enigma. To clearly define the relationship of RNAs and obesity, the composition of idiopathic infertile male sperm RNAs from the Reproductive Medicine Network (RMN) Assessment of Multiple Gestations from Ovarian Stimulation (AMIGOS) trial [36] was examined as a function of BMI. A representative subset of 32 semen samples from subjects with an average BMI of 28.55 kg/m² (range 19.25–40.68 kg/m²) was selected (Table 1, Table S1). This approximated the AMIGOS study average BMI of 28.94 kg/m² (range 15.79–63.70 kg/m²). This group was segregated into three BMI categories: healthy (BMI < 25 kg/m², 8 samples), overweight (BMI 25–30 kg/m², 11 samples) and obese (BMI > 30 kg/m², 13 samples).

Table 1. Demographic characteristics of AMIGOS males. All values are reported as the mean (range) for A) the 32-sample subset and full AMIGOS study and B) the three BMI category groups; healthy (BMI < 25 kg/m²), overweight (BMI 25–30 kg/m²) and obese (BMI > 30 kg/m²).

A. AMIGOS subset compared to full study population demographics
	AMGIOS subset (n = 32)	AMIGOS study (n = 900)
BMI (kg/m^2)	28.548 (19.25–40.68)	28.94 (15.79–63.70)
Age (years)	34.22 (29–44)	34.49 (19–63)
Testosterone (ng/dL)	404.87 (34.10–841.85)	425.85 (9.20–1131.90)
Glucose (mg/dL)	85.02 (60.60–121.80)	89.77 (18.30–364.70)
Insulin (mIU/L)	14.84 (2.00–70.70)	19.82 (2.00–272.00)
SHBG (nmol/L)	29.14 (10.10–47.30)	29.43 (5.97–136.00)
FSH (mIU/mL)	4.03 (1.08–11.20)	4.50 (0.10–24.63)
LH (IU/L)	3.92 (0.80–7.80)	3.76 (0.10–48.85)
TSH (mIU/mL)	1.92 (0.52–7.39)	1.83 (0.26–17.90)
INHB (pg/mL)	210.65 (120.02–340.48)	191.6 (8.80–566.20)
B. AMIGOS subset demographics by BMI category
	Healthy BMI (n = 8)	Overweight BMI (n = 11)	Obese BMI (n = 13)
BMI (kg/m^2)	22.104 (19.25–24.37)	26.93 (25.10–19.16)	33.89 (30.30–40.68)
Age (years)	32.75 (29–39)	36.09 (30–44)	33.54 (30–39)
Testosterone (ng/dL)	498.40 (403.80–633.9)	353.69 (244.50–444.50)	380.87 (34.10–841.85)
Glucose (mg/dL)	86.49 (63.70–119.90)	83.92 (60.6–102.40)	84.86 (75.9–121.8)
Insulin (mIU/L)	10.64 (2.00–13.60)	19.18 (5.6–70.7)	14.38 (4.07–44.40)
SHBG (nmol/L)	34.48 (19.40 − 47.30)	21.86 (10.80–29.40)	31.04 (16.10–45.90)
FSH (mIU/mL)	3.93 (3.15–5.13)	3.73 (1.08–6.96)	4.32 (2.30–11.20)
LH (IU/L)	4.05 (1.99–6.42)	3.43 (0.80–6.98)	4.21 (1.29–7.80)
TSH (mIU/mL)	1.27 (0.52–2.89)	2.86 (1.30–7.39)	1.64 (0.91–2.80)
INHB (pg/mL)	209.78 (163.85–316.01)	201.10 (120.02–281.03)	218.39 (141.59–340.48)

SREs associated with BMI

Sperm RNA was isolated and sequenced as described [30]. Each sequence was then directly mapped to the human genome and a series of exon sized RNA elements (REs) were defined using the RNA Element Discovery algorithm, REDa [37]. REDa identifies novel transcribed orphan elements, elements within 10 kB of an exon, intronic elements as well as annotated exons. It does not solely rely on the prior knowledge of transcripts or splicing structure. This strategy simultaneously affords both hypothesis testing while providing an unbiased measure to evaluate the impact of BMI on the distribution of sperm REs. Using a Generalized Linear Model (GLM) [38], SRE abundance was assessed in response to increasing BMI. This resolved a set of 487 BMI associated SREs (BH adjusted p-value < 0.05), comprised of 12 novel orphan SREs, 20 SREs within 10 kB of an exon, 43 intronic SREs and 412 exonic SREs, which in total correspond to 421 genes. As a function of increasing BMI, 333 SREs (283 genes) were upregulated and 155 SREs (155 genes) were down regulated (Table S2). They encoded 184 proteins from the human sperm proteome [39], of which 132 were associated with upregulated transcripts and 57 with down regulated transcripts (Table S2). Genomatix Pathway Systems (GePS) Gene Ontology (GO) of the 421 SRE-associated gene name set highlighted cilium organization (GePS adjusted p-value < 0.001), cilium assembly (GePS adjusted p-value < 0.001), ciliopathies (GePS adjusted p-value < 0.001) and enrichment in chromosome organization (GePS adjusted p-value < 0.001). Signaling pathways related to the cellular response to stress (GePS adjusted p-value = 1.60e-2) were also identified.

To define the molecular pathways, the set of 421 BMI associated genes were evaluated across the Molecular Signature Database (MSigDB) that contains well-defined functional gene sets. This revealed a series of downstream adipogenesis responsive genes following transcription factor (TF) peroxisome proliferator activated receptor gamma (PPARG) and retinoid X receptor alpha (RXRA) activation [40] (MSigDB p-value = 1.93e-11, Table S3). Significant enrichment was also observed with respect to inflammatory processes within the MSigDB immunological gene sets. These included those genes downstream from the activation of effector CD4+ Th1 cells [GSE32901, 41] (MSigDB p-value = 2.52e-9, Table S3), toll-like receptor-2 (TLR2) mediated signaling [GSE25146, 42] (MSigDB p-value = 3.27e-9, Table S3) and interferon-alpha (IFN-α) compared to interferon-gamma (IFN-γ) signaling [GSE3920, 43] (MSigDB p-value = 2.12e-8, Table S3). Together, the above was indicative of a complex epididymal fat pad (eWAT) modifier regulatory network.

To identify membership within the co-regulated networks, complexity was reduced and modeled using the Mfuzz algorithm. RNA levels were considered as a function of BMI partitioned into three categories, healthy, overweight and obese respectively [2,5,29]. As shown in Figure 1 and summarized in Table S2, this resolved 8 unique expression patterns from 235 SREs, segregating into 4 general patterns. Three resolved as patterns derived from 88 SREs – 86 genes – 41 sperm proteome proteins that continually increased with increasing BMI (Figure 1(a)); 2 patterns derived from 69 SREs – 68 genes – 32 sperm proteome proteins that increased once obese Figures 1(b) and 2 patterns derived from 54 SREs – 51 genes – 20 sperm proteome proteins that steadily decreased with increasing BMI (Figure 1(c)). Pattern 1D, was characterized by 24 SREs from 26 genes for which 6 sperm proteome proteins were identified. This was the only group of SREs that decreased once overweight. However, insufficient membership precluded meaningful ontological analysis.

Figure 1. BMI responsive sperm RNA element (SRE) clusters. Response clusters identified in Mfuzz include 235 SREs that follow 8 unique cluster patterns. General responses are comprised of (a) 3 cluster patterns that respond with a general steady increase with body mass index (BMI) indicated by the red section of the pie chart in the upper left-hand corner. (b) Two cluster patterns that follow a general BMI response of increasing once obese. These clusters are indicated by the orange pie chart section in the upper left-hand corner. (c) Two clusters indicated by the blue pie chart section in the lower left-hand corner comprise a general steady decrease as BMI increases. (d) One cluster indicated with the green slice of the pie chart in the lower left-hand corner responded to BMI with a decrease once overweight. Clusters were determined from the significant 487 SREs identified by the generalized linear model (GLM). BMI category; healthy (< 25 kg/m²), overweight (25–30 kg/m²) and obese (30+ kg/m²), group means were assigned. The Mfuzz algorithm was set to cluster only those SREs with mean reads per kilobase per million (RPKM) > 10 in at least one BMI group and follow the cluster pattern with a membership of 0.4. Graphs are visualized across BMI groups based on expression fold-change of the cluster pattern.

Figure 2. Obesity related impact on sperm. Increased BMI induces the obesity-related hallmarks adipogenesis, stress and low-grade systemic inflammation within the epididymal white adipose tissue (eWAT) to alter the sperm environment. The altered sperm environment changes sperm RNAs regulated by Chromatin remodeler cofactors, RNA interactors, Erasers and Writers (CREWs) of methylation, acetylation and ubiquination to impact sperm characteristics and function.

SREs that steadily increased with BMI shown in Figure 1, cluster A_II (37 SREs/34 genes/18 sperm proteome) were attributed to ciliary processes related to cilium assembly (GePS adjusted p-value < 0.001) and cilium organization (GePS adjusted p-value < 0.001). The remaining A clusters were associated with a cellular response to stress (GePS adjusted p-value = 2.00e-3) and all A clusters associated with genes downstream of transforming growth factor – beta 1 (TGF-β1) as well as interleukin-6 (IL-6) (MSigDB p-value = 4.74e-6, Table S3) [GSE39820, 44]. Reiterating the enrichment in chromosome organization found within the 421 BMI associated genes, cluster B genes were enriched in the regulation of histone H3-K9 acetylation (GePS adjusted p-value < 0.001), chromosome organization (GePS adjusted p-value < 0.001) and regulation of chromosome organization (GePS adjusted p-value = 2.00e-3). Adipogenesis related gene set enrichment was evident in both B clusters that showed increased expression once obese in response to adipogenic hormones (MSigDB p-value = 2.04e-5, Table S3) [45] and C clusters that steadily decreased as BMI increased bound by the TFs PPARG and RXRA [40] (MSigDB p-value = 8.44e-6, Table S3), suggesting a possible TF coregulatory network.

The MSigDB TF binding motif gene set module was then used to identify enriched TF binding motifs in the promoter region of the 421 BMI associated genes and the 224 Mfuzz BMI responsive cluster genes (Table S2). In total, 70 TF binding motif associated gene sets were identified as enriched in both the 421 correlated and 224 responsive gene sets, with 55 assigned to a known TF [46,47]. Of the 55 binding motifs, only 8 with evidence of TF transcript presence were considered further (Table S4). These 8 TFs were evident in the sequencing data outside of the 421 BMI correlated genes; specificity protein 1 (SP1), lymphoid enhancer binding factor 1 (LEF1), GA binding protein (GABP, subunits alpha (GABPA), beta 1 (GABPB1) and beta 2 (GABPB2)), MYC associated zinc finger protein (MAZ), paired box 4 (PAX4), nuclear factor of activated T cells (NFAT, family members NFATC1, NFATC2, NFAT3, NFAT4, NFAT5 and NFAT2 interacting protein (NFATC2IP)), E2F transcription factor 1 (E2F1) and meis homeobox 1 (MEIS1), with MAZ also present within the sperm proteome. Of these lymphoid enhancer binding factor 1 antisense RNA 1 (LEF1-AS1), GABPB2, MAZ (short isoform only), NFAT3 and NFAT4 were full length (Table S5). Interestingly, LEF1-AS1, a non-coding RNA that acts to attenuate LEF1 expression [48], is an important mediator of Wnt signaling involved in early development [49]. Only the TFs MAZ, NFAT3 and NFAT4 retain functional capacity following fertilization as GABP is an obligate heterotetramer of the GABPA subunit and one GABPB, GABPB1 or GABPB2 subunit [50]. BMI response of the genes within these 8 TF:gene associated pairs was assessed. GABP is associated with 7 TF A cluster associated genes, MAZ with 11 and PAX4 with 8, that predominantly respond with a steady increase with BMI (Figure 1). Similarly, of those pairs primarily responding with an increase once obese (Figure 1), SP1 is associated with 15, LEF1 with 9, and MEIS1 with 5 B cluster responsive genes. In comparison, of the TF:gene associated pairs that steadily decreased as BMI increased (Figure 1), E2F1 is associated with 5 and NFAT with 8 C cluster responsive genes.

While the corresponding TF encoding SREs were present within the total sperm SRE profile (343,102 SREs), they were not within the 487 BMI associated SREs suggesting that, while necessary, they are likely not a master regulator of the BMI response. CRREWs (Chromatin remodeler cofactors, RNA interactors, Readers, Erasers and Writers) were then considered since they act prior to TFs, together regulating chromatin state as necessary for expression [51–54]. To determine whether CRREWs might function in the regulation of these pathways, the 421 BMI associated, and 224 BMI responsive genes were tested for the presence of CRREW function (Table 2). A list of 650 potential CRREWs [39,52,53,55–62] that are also reported in the human sperm proteome [39] were compiled. Of the 421 genes that associate with BMI, 24 are CREWs, Chromatin remodeler cofactors, RNA interactors, Erasers and Writers (p = 0.001, Table 2) as no Readers were identified, with 12 responding to the BMI responsive patterns (p = 0.0231, Table 2, Figure 1). Further, 21 of the 24 CREWs are reported within the sperm proteome with 9 exhibiting a unique BMI response (Figure 1, Table 2). Transcript integrity was evaluated (Table 2) and 16 of the 24 CREWs were full length, with 6 directly observed within the BMI responsive clusters as necessary for translation.

Table 2. Transcript integrity for Chromatin remodeler cofactors, RNA interactors, Erasers and Writers (CREWs) also found in the sperm proteome that are associated with BMI. MAZ and NFAT are indicated with an * as the NFAT members NFAT3 and NFAT4, along with MAZ are present as a full-length transcript. Slope refers to the change in SRE abundance (RPKM) as BMI is increased by a single integer while BMI Response informs whether the SRE is present within one of the four over-all responses; steady increase with BMI, increase once obese, steady decrease with BMI and stable decrease once overweight, or is only associated (positive, negative) with BMI from the generalized linear model.

Gene	Gene Name	Sperm Proteome	Function	Slope	BMI Response	Transcript Integrity	Presence in TF^a Binding Motif Gene Set
ARID4B	AT-rich interaction domain 4B	Yes	Chromatin remodeler cofactor	0.886816	Positive	Full	n/a
CHD2	chromodomain helicase DNA binding protein 2	Yes	Chromatin remodeler cofactor	−0.70262	Negative	Full	LEF1, MAZ*, NFAT*, MEIS1
CHD5	chromodomain helicase DNA binding protein 5	Yes	Chromatin remodeler cofactor	−1.22435	Negative	Full	n/a
CUL4B	Cullin 4B	Yes	Writer	2.730699	Increase once obese	Full	n/a
DDX5	DEAD (Asp-Glu-Ala-Asp)-box helicase 5	Yes	Chromatin remodeler cofactor	6.28644	Increase once obese	Full	SP1, LEF1, MAZ*
HUWE1	HECT, UBA And WWE domain containing 1, E3 ubiquitin protein ligase	Yes	Writer	1.41212	Steady increase with BMI	Full	n/a
KDM2B	lysine demethylase 2B	Yes	Eraser	3.963266	Positive	Full	n/a
KDM4B	lysine demethylase 4B	Yes	Eraser	−1.69199	Negative	Full	n/a
MYO1C	myosin 1C	Yes	Chromatin remodeler cofactor	0.824644	Positive	Full	SP1, PAX4, NFAT*
NELFE	negative elongation factor complex member E	Yes	RNA interactor	1.392059	Increase once obese	Full	n/a
SMARCC1	SWI/SNF related, matrix associated, actin dependent regulator of chromatin subfamily C member 1	Yes	Chromatin remodeler cofactor	−6.63593	Negative	Full	LEF1, MAZ*
SMCHD1	structural maintenance of chromosome flexible hinge domain containing 1	Yes	Chromatin remodeler cofactor	1.700137	Increase once obese	Full	n/a
TLE4	transducing like enhancer of split 4	Yes	Chromatin remodeler cofactor	0.816017	Positive	Full	SP1, MEIS1
USP15	ubiquitin specific peptidase 15	Yes	Eraser	0.686767	Positive	Full	SP1, MEIS1
USP7	ubiquitin specific peptidase 7	Yes	Chromatin remodeler cofactor/ Writer/Eraser	−1.42733	Negative	Full	n/a
KDM5A	lysine demethylase 5A	No	Eraser	0.681963	Steady increase with BMI	Full	n/a
HNRNPC	heterogeneous nuclear ribonucleoprotein C (C1/C2)	Yes	RNA interactor	1.681745	Increase once obese	5ʹ depleted	LEF1, PAX4, NFAT*
NIPBL	cohesion loading factor	Yes	Chromatin remodeler cofactor	1.610834	Steady increase with BMI	5ʹ depleted	SP1, GABP, MAZ*, MEIS1, E2F1
PHF12	PHD finger protein 12	Yes	Eraser	0.869446	Steady increase with BMI	5ʹ depleted	SP1, LEF1, GABP, MAZ*, PAX4, MEIS1, E2F1
FBL	fibrillarin	No	Writer	4.87741	Positive	5ʹ depleted	SP1, NFAT*
				4.374552	Steady increase with BMI
HDAC5	histone deacetylase 5	No	Eraser	1.397323	Increase once obese	5ʹ depleted	n/a
DNMT1	DNA methyltransferase 1	Yes	Writer	3.582606	Positive	3ʹ depleted	GABP, MAZ*, MEIS1, E2F1
ZNF451	zinc finger protein 451	Yes	Chromatin remodeler cofactor	4.449028	Increase once obese	3ʹ depleted	SP1
				2.234207	Increase once obese
				1.747638	Positive
DPY30	histone methyltransferase complex regulatory subunit	Yes	Writer	4.12142	Positive	Depleted	SP1

^a TF; Transcription Factor

As summarized in Table 2, of the 11 Chromatin remodeler cofactors, 9 were full-length transcripts and all 11 were within the human sperm proteome. Of these Chromatin remodeler cofactors, 7 were positively correlated with BMI and 4 were negatively correlated with BMI. Cluster membership of the Chromatin remodeler cofactors segregated within the B clusters (DEAD (Asp-Glu-Ala-Asp)-box helicase 5 (DDX5) SRE within cluster B_II, structural maintenance of chromosome flexible hinge domain containing 1 (SMCHD1) SRE within cluster B_I, and zinc finger protein 451 (ZNF451) SREs within cluster B_I and cluster B_II) and cohesion loading factor (NIPBL) SRE was within the A_II clusters (Figure 1, Table S2).

RNA interactors constituted of the smallest group containing the full-length negative elongation factor complex member E (NELFE) transcript and heterogeneous nuclear ribonucleoprotein C (C1/C2) (HNRNPC) (Table 2), which are within the sperm proteome. Both RNA interactors were positively correlated with BMI, and part of the Figure 1(b) clusters (NELFE SRE within cluster B_II and HNRNPC SRE within cluster B_I, Table S2). Seven methylation, acetylation or ubiquitin Erasers were identified, of which 5 were full-length and 5 were identified within the sperm proteome (Table 2). Based on cluster membership (Figure 1) histone deacetylase 5 (HDAC5) increased once obese (Figure 1 cluster B_II, Table S2) whereas lysine demethylase 5A (KDM5A) and PHD finger protein 12 (PHF12) increased steadily as obesity increased (Figure 1 cluster A_II, Table S2). The remaining positively correlated Erasers, lysine demethylase 2B (KDM2B) and ubiquitin specific peptidase 15 (USP15) (full-length transcripts) and the negatively correlated Erasers, lysine demethylase 4B (KDM4B) and ubiquitin specific peptidase 7 (USP7), did not show a cluster specific response to BMI but are within the sperm proteome (Table 2). This suggests that while these four Erasers are correlated with BMI, they are independent of the co-regulated clusters. The six methylation, acetylation or ubiquitin Writers included 3 full-length transcripts and 5 as present within the sperm proteome (Table 2). When considered as a function of cluster membership, cullin 4B (CUL4B, cluster B_II, Table S2) increased once obese (Figure 1, Table 2) while HECT, UBA and WWE domain containing 1, E3 ubiquitin protein ligase (HUWE1, cluster A_II, Table S2) and fibrillarin (FBL, cluster A_III, Table S2) steadily increased with BMI (Figure 1, Table 2). The remaining correlated Writers, DNA methyltransferase 1 (DNMT1), histone methyltransferase complex regulatory subunit (DPY30) and USP7 were not identified within any co-regulatory BMI responsive clusters (Table 2). Interestingly, USP7 is reported to function as a Writer of methylation [39,53], an Eraser of histone modifications [53] and a chromatin remodeler cofactor [39].

Of the 24 CREWs, 13 had binding sites for the TFs identified in this study (Table 2), with the majority sharing a MAZ, NFAT3 or NFAT4 site. TF binding to the CREW was confirmed using the Encyclopedia of DNA Elements (ENCODE) Transcription Factor Target gene set [63]. MAZ has been shown to bind, DDX5, DNMT1, NIPBL, PHF12 and SWI/SNF related, matrix associated, actin dependent regulator of chromatin subfamily C member 1 (SMARCC1) (Table 2). The Chromatin remodeler cofactors CHD2, DDX5 and SMARCC1 are full-length transcripts and within the sperm proteome. NFAT3 and NFAT4 (Table S5) are the only members of the NFAT family that are full-length transcripts. This TF family was enriched in two separate gene sets (TGGAAA_NFAT_Q4_01 and NFAT_Q4_01, Table S4), both using the same binding motif. CDH2, myosin 1C (MYO1C), HNRNPC and FBL have the potential to be bound by the NFATs (Table 2), whereas CHD2, HNRNPC and FBL have been shown to be bound by the NFAT family member NFATC1, which shares a conserved DNA-binding domain [64]. All members of the NFAT family share this DNA-binding domain, suggesting that the remaining NFAT members, including NFAT3 and NFAT4, should also bind to CHD2, HNRNPC and FBL by NFATC1. Of these 3 binding pairs with NFATC1, the Chromatin remodeler cofactor CHD2 is the only full-length transcript (Table 2) and is also a known to be bound by MAZ. Together this shows the interplay of the CREWs with TFs orchestrating the BMI response.

Discussion

Our overweight population is a growing concern with ~33% of adults in the United States classified as obese [2]. BMI is commonly used as a measure of obesity which in mice, is associated with altered small noncoding sperm RNA, DNA methylation and DNA fragmentation. Such effects may be epigenetically heritable [14,21,23,25–27,33]. Human studies have identified a series of associations of paternal BMI to birth outcome and offspring phenotype [18–20]. Similarly, common semen parameters have correlated with increasing BMI [4,5,12,13] and appear independent of total testosterone with the possible exception of Kruger’s strict morphology [65]. Few human obesity studies have assessed epigenetic alterations. Of these studies, only circulating miRNAs or sperm DNA methylation and small non-coding RNAs have been evaluated [28,29]. Here, we identified 487 SREs responsive to BMI from our long RNA libraries.

RNAs associated with some of the key hallmarks of obesity (including adipogenesis, stress and low-grade systemic inflammation) were enriched within the 487 sperm-BMI associated SREs (i.e., 421 genes) identified. Their relationships are summarized in Figure 2. Briefly, as BMI increases, adipogenic hormones such as insulin act directly in high lipid content tissues like abdominal and epidydimal adipose [66]. eWAT, similar to abdominal WAT (aWAT), responds to insulin [6] as part of a key metabolic regulatory pathway necessary for spermatogenesis [45,66]. As illustrated in Figure 3, 4. Insulin like growth factor 1 receptor (IGF1R), a key receptor involved in insulin signaling within eWAT, is present at the core of the response. It is a member of ciliary associated cluster A_II (Figure 1, Table S2). For example, IGF1R signaling is dependent on a functional cilium [67,68] and acts to regulate adipocyte differentiation [69]. This directly affects pathways related to adipogenesis, cellular stress and inflammation.

Figure 3. Network of the 24 CREWs correlated with BMI. Chromatin remodeler cofactors, RNA interactors, Erasers and Writers (CREWs) involved in chromatin organization within the 487 SREs associated with BMI acting as key regulators of the obesity related pathways. Yellow boxes highlight the location of IGF1R (center) and the Chromatin remodeler cofactor DDX5 (center) which act as a feedback loop continuing the progression of obesity altered insulin signaling. Line color indicates whether the reported interaction is validated (solid brown) or predicted (dashed light grey).

Figure 4. Network of the 12 CREWs responsive to BMI. Chromatin remodeler cofactors, RNA interactors, Erasers and Writers (CREWs) involved in chromatin organization within the 235 MFuzz BMI responsive SREs acting as key regulators of the obesity-related pathways. Yellow boxes highlight the location of IGF1R (center) and the Chromatin remodeler cofactor DDX5 (center) which act as a feedback loop continuing the progression of obesity-altered insulin signaling. Line color indicates whether the reported interaction is validated (solid brown) or predicted (dashed light grey).

The effects of obesogens and adipogenic hormones are often mediated through TFs, e.g., PPAR and RXR family members [66]. Within the 421 BMI associated genes, downstream genes responding to PPARG and RXRA [40] steadily declined with increasing BMI (Figure 1(c)). In comparison, the SREs that increased once obese are downstream from insulin signaling [45] (Figure 1(b)). As expected, cellular stress RNAs, respond to insulin and other adipogenic hormones [3,4], was identified with increasing BMI (Figure 1 clusters A_I and A_II).

As adipogenesis progresses, cell free DNA can be released which, if left unchecked, facilitates chronic low-grade systemic inflammation [2,70] mediated through the action of TLRs [71]. Interestingly, the 421 BMI associated SREs associated with the downstream response follow activation of TLR2 [GSE25146, 42], a regulator of glucose and insulin [72]. Similarly, these SREs were enriched in genes responding to the pro-inflammatory cytokine IFN-α [GSE3920, 43], are within the A cluster genes (Figure 1) that steadily increase with BMI, TGF-β1 and IL-6 [GSE39820, 44]. DDX5, a Chromatin remodeler cofactor [73–75] (Figure 3, center) belonging to cluster B_II (Figure 1), is an enriched downstream target of both IFN-α and active CD4+ Th1 cells . It increases in response to cyclooxygenase 2 (COX2) [76] and is also associated with pro-inflammatory response [77]. This is in accord with a TLR cell-free DNA inflammatory response through IFN-α, TGF-β1, IL-6 and COX2 cytokine signaling impacting sperm RNA.

The sperm’s, environment also changes utilizing CRREWs as epigenetic regulators of gene expression through chromatin accessibility. For example, the increase in the Chromatin remodeler cofactor DDX5 in response to obesity-mediated inflammation will act in a feedback loop to increase insulin signaling (Figure 3, center) thus continuing the cycle. Specifically, insulin signaling through IGF1R (cluster A_II-member Figure 1) altered during obesity promotes an increase in pro-inflammatory cytokines, including COX2. DDX5, as part of the B clusters (Figure 1), is then stimulated by COX2 to increase DDX5 production as BMI increases (Figure 3, Figure 4). Through the depression of miR-183, DDX5 action will lead to an increased expression of insulin receptor substrate 1 (IRS1), a ligand of IGF1R [76].

The interaction of CREWs and TFs in response to the effects of adipogenesis, stress and inflammation is summarized in (Figure 3). A total of 24 CREWs were contained within the 421 BMI correlated SRE associated gene names, of which 16 are full length transcripts and 21 have been reported in the human sperm proteome (Table 3, Figure 3). Chromatin remodeler cofactors were the most prevalent, followed by Writers, Erasers and finally RNA interactors. Interestingly, Readers [54] appear to be absent, within the silenced sperm genome. Twelve of these CREWs retained membership within 1 of the 8, Figure 1 clusters, primarily cluster B that directly increased when obese, in which 6 are full length transcripts and 9 CREWs are members of the sperm proteome (Table 2, Figure 4).

Using the mechanisms of obesity-induced inflammation (similar to those described in abdominal adipose tissue), this data supports the view that obesity related changes impacting eWAT will negatively affect sperm. As shown above these obesity-induced alterations connecting adipogenesis, stress and inflammatory pathways yields a BMI responsive alteration of chromatin modifiers, within the network of CREWs (Figure 2) that negatively affect the male gamete and potentially extend beyond as carried by human sperm.

Methods and materials

AMIGOS study and subset demographics

The Assessment of Multiple Gestations from Ovarian Stimulation (AMIGOS) trial [36] was designed to address use of letrozole compared to clomiphene citrate or gonadotropins on multiple gestation rate during intrauterine insemination in idiopathic infertile couples. It included a total of 900 couples from which male characteristics including BMI, age, serum hormone levels and semen parameters were obtained. A subset of 32 available semen samples with an average BMI of 28.55 kg/m² (range 19.25–40.68 kg/m²) representative of the entire AMIGOS study (average BMI 28.94 kg/m², range 15.79–63.70 kg/m²) was obtained from the Reproductive Medicine Network (RMN) biorepository (Table 1, Table S1). BMI was used as a measure of obesity in this subset of 32 samples. This included eight samples within the healthy BMI range (< 25 kg/m²), 11 within the overweight BMI range (25–30 kg/m²) and 13 within the obese I (30–34.9 kg/m²) and II (> 35 kg/m²) BMI ranges. The deidentified frozen samples were processed and analyzed in accordance with HIC protocol 095701MP2E(5R).

Sequencing and SRE discovery

Sample preparation and RNA-seq library construction was accomplished as previously described [30,78,79] followed by sequencing using the Illumina Hi-Seq 2500 platform. Sequence alignment to the human genome (hg38) and sperm RNA element (SRE) discovery was accomplished using an in-house pipeline and REDa, RNA Element Discovery algorithm [37]. An average of 97.7% of sequences passed quality control with an average of 16.5% of sequences marked as unmapped. SRE expression was calculated using Reads per Kilobase Million (RPKM). Transcript integrity was evaluated on the University of California Santa Cruz (UCSC) Genome Browser (http://genome.ucsc.edu/) [80] and SREs aligned to genomic regions (exonic, within 10kB of an exon and intronic), associated with a gene name were assigned the appropriate Ensembl ID. Genomic regions with no assigned gene were considered novel-orphan regions.

Statistical modeling

Principal component analysis prcomp() function in R [38] was used to assess semen parameters, age, BMI, hormone levels, birth outcome, smoking/alcohol use and technical sequencing statistics such as percent genomic duplicated reads and percent unmapped reads, that may impact SRE variability. A Pearson correlation coefficient > |0.45| and p-value < 0.05 was considered as a significant impact. PC1 (principal component 1) contributed approximately 98% of SRE variability but did not meet this criterion and was excluded from the final model. Three parameters, percent mitochondria, percent sperm motility and percent genomic duplicated, had a p-value < 0.05 but a Pearson correlation coefficient < |0.45|.

The generalized linear model (GLM) in the R stats package [38] was used to analyze the impact of BMI on SRE abundance using BMI as a continuous variable (Supplementary Methods). Obesity is associated with comorbidities such as diabetes and semen parameters. Additionally, within the full AMIGOS study, low total testosterone was correlated with BMI [65]. To determine if BMI was associated with any of the measured parameters in this subset, Pearson correlations were evaluated between BMI with age, hormone levels, semen parameters, birth outcome, smoking/alcohol use and technical sequencing statistics. BMI was not significantly associated (p-value < 0.05, Pearson correlation coefficient > |0.45|) with any parameter indicating no interactions between BMI and the other parameters could be tested in this subset. A total of three GLMs were used to determine the influence of BMI on SRE abundance. As no significant Pearson correlations between BMI or for PC1 were identified, the first model did not regress out any parameters. No SREs were identified as associated with BMI from this model which we hypothesized was due to a large degree of background variance resulting from a large number of lowly contributing parameters. For the second model the effects of paternal age and hormone levels (total testosterone, insulin, glucose, SHBG, FSH, LH, TSH, INHB) were removed by regression. An initial threshold of SREs > 0 RPKM in 1/3 or more samples was applied prior to modeling to remove all lowly expressed SREs. As PC1 was not regressed out of the model, SRE variance was considered by requiring a mean and median RPKM to be within 20%. This was set as the mean/median RPKM difference being between 0.8 and 1.2. Significance was met if the mean > 10 RPKM and Benjamini-Hochberg (BH) [81] adjusted p-value < 0.05. A total of 487 SREs were identified for further analysis. The third model regressed out age, hormones and the weakly correlated contributors to PC1 (percent mitochondria, percent sperm motility and percent genomic duplicated) in which the same 487 SREs were identified as associated with BMI.

Co-regulatory networks were identified with the Mfuzz algorithm [82]. Briefly, mean RPKM was calculated for each of the BMI categories; healthy (8 samples, BMI < 25 kg/m²), overweight (11 samples, BMI 25–30 kg/m²) and obese (13 samples, BMI > 30 kg/m²). Fold changes between the BMI category mean were determined and like patterns (membership score > 0.4) were assessed. Initially, a default of 10 clusters was created within the Mfuzz program. Cluster number was systematically decreased until each cluster exhibited one specific pattern. Eight clusters were used for further analysis retaining all 4 major BMI response patterns.

An inclusive list of CRREWs (Chromatin remodeler cofactors, RNA interactors, Readers, Erasers and Writers) was compiled using the WERAM 1.0 database [56], EpiFactors Database [53], GeneCards [62], SnapShot Archive: CellPress [57–61] and other literature sources [39,52,55]. This list was further filtered to include only those found within the human sperm proteome [39]. CRREW probability was determined using phyper() function [38]. Briefly, this was calculated as the probability of identifying 24 CREWs (Chromatin remodeler cofactors, RNA interactors, Erasers and Writers) within the 421 BMI associated genes or the 12 CREWs from the 224 BMI responsive genes and their occurrence within the sperm proteome.

Gene ontology (functional) analysis

Gene set enrichment analysis was performed on the 487 BMI associated SREs, 235 BMI responsive SREs and BMI responsive patterns (Table S2). Enriched GO terms relating to Biological Processes and MeSH (Medical Subject Headings) terms were analyzed using Genomatix Pathway Systems [GePS, 83]. Enriched signaling pathways identified in GePS were supplemented with the Molecular Signature Database [MSigDB, 47, 84]. Enrichment was evaluated using the MSigDB curated datasets compiled from online publication pathway databases, publications in PubMed and from knowledge domain experts; immunological signature gene sets and TF binding motif gene sets. Gene sets that corresponded to the enriched signaling pathways identified in GePS were examined further. Harmonizome, the ENCODE Transcription Factor Target gene set [63] was used to identify those TF: gene binding motifs that have been validated.

Gene network visualization of enriched Biological Processes and pathways was accomplished using Cytoscape [85]. The BiNGO module in Cytoscape [85,86] was used for identification of the broad GO annotation terms that encompass the specific enriched pathways. CREWs and upstream interactors were defined first by 1) the specific type of function performed (e.g. histone methylation) and 2) the created term Chromatin remodeler cofactor, RNA interactors, Writer, Eraser and upstream interactor of Writers/Erasers.

Data availability

The data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus (Edgar et al., 2002) and are accessible through GEO Series accession number GSE134580 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE134580).

Author contribution

SAK and GS designed the research; GS, RG processed sperm RNA samples; GS and ME performed and analyzed the data; GS and SAK wrote the manuscript. All authors contributed to the editing and review of the manuscript.

Acknowledgments

The work was supported in part by National Institutes of Health (NIH)/Eunice Kennedy Shriver National Institute of Child Health and Human Development (NICHD) Grants U10 HD039005. Support from the Charlotte B. Failing Professorship and Wayne State University Grants Boost award 19-0381 to SAK is gratefully acknowledged. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NICHD or NIH.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by the NIH/NICHD [U10 HD039005]; Wayne State University Charlotte B. Failing Professorship and Grants Boost 19-0381 award.