Growth hormone receptor disruption causes different gene changes in fat tissue of adult males and females

Mice with a C57BL/6J genetic background carrying loxP sites flanking exon 4 of the Ghr gene were previously produced by the Knockout Mouse Project (KOMP) [17]. C57BL/6J mice expressing a ubiquitous Cre recombinase gene driven by the ROSA26 gene promoter/enhancer (B6.129-Gt(ROSA)26Sortm1(cre/ERT2)Tyj/J mice) were purchased from The Jackson Laboratory [15, 18]. Homozygous mice for both the floxed Ghr gene and the Cre recombinase gene were used to generate the 6mGHRKO mice as previously described [18]. Briefly, mice of 6 months of age were used to generat the 6mGHRKO and control mice by injecting Tamoxifen or peanut oil, respectively, for 5 consecutive days [18]. Mice were housed in a temperature and humidity-controlled room at 22 °C under a 14-hour light, 10-hour dark cycle, with 3–4 mice per cage, and with ad libitum access to water and standard laboratory chow (ProLab RMH 3000). All mouse protocols were approved by Ohio University's Institutional Animal Care and Use Committee.

A Bruker Minispec NMR analyzer (Bruker Corp., The Woodlands, TX) was used to measure body composition at the time of dissection, (12 months of age, suppl. Table 1), as previously described [18]. Longitudinal body composition and body weight determined every month starting 1 day before Tam treatment until the time of dissection was described in a previous publication [18]. Body-length was measured at the time of dissection from the tip of the nose to the anus.

Twelve-month old mice were euthanized by cervical dislocation, and adipose samples were dissected and weighed. For RNA isolation, tissue was flash frozen in liquid nitrogen and stored at −80 °C until further processing. Frozen tissues were homogenized using a Precellys 24-Dual homogenizer, and RNA isolated using the QIAGEN RNeasy Lipid Tissue Mini Kit following manufacturer's instructions. The quantity and quality of total RNA were measured with the NanoDrop 2000c (Thermo Scientific). The RNA used for subsequent experiments had a 260/280 and a 260/230 ratio ≥ 1.8 to ensure purity and quality of the RNA. The quantity and quality of the isolated RNA were determined by the Agilent 2100 Bioanalyzer from Agilent Technologies at the Genomics Facility (Ohio University). The RNA used for RNA-seq had an RNA Integrity Number ≥ 8.

Male and female adipose mRNA samples isolated from three 6mGHRKO mice of each sex and three control littermates of each sex were used to create cDNA libraries. First, mRNA was isolated from the total RNA samples with Takara Clontech SMART technology, according to the manufacturer's protocol. Subsequently. the samples were indexed and sequenced on an Illumina NovaSeq S4. Basecalls and demultiplexing were performed with Illumina's bcl2fastq software and a custom python demultiplexing program with a maximum of one mismatch in the indexing read for pooled samples. Alignment of the RNA-seq reads was performed using Ensembl top-level assembly with STAR. Subread: featureCount was used to generate gene counts derived from the number of uniquely aligned unambiguous reads. Sequencing performance was assessed for the total number of aligned reads, total number of uniquely aligned reads, and features detected. RSeQC was used to quantify the ribosomal fraction, known junction saturation, and read distribution over known gene models. It is important to note that RNA from each individual mouse (n = 3 per sex/genotype) was sequenced independently. No physical pooling of RNA occurred at the experimental stage. The term "pooling" in this study refers exclusively to the statistical integration of data during downstream bioinformatics analysis (edgeR/limma with voom normalization), where replicate samples were modeled together to increase statistical power and reduce variability.

All gene counts were then imported into the R/Bioconductor package EdgeR and TMM normalization size factors were calculated to adjust the samples for differences in library size. Ribosomal genes and genes not expressed in the smallest group size minus one sample greater than one count-per-million were excluded from further analysis. The TMM size factors and the matrix of counts were then imported into the R/Bioconductor package Limma. Performance of the samples was assessed with Spearman correlations, a MultiDimensional Scaling plot, and hierarchical clustering. Weighted likelihoods based on the observed mean-variance relationship of each gene and sample were then calculated with the voomWithQualityWeights. The performance of all genes was assessed with plots of the residual standard deviation of each gene to their average log-count with a robustly fitted trend line of the residuals. Differential expression analysis was then performed to analyze differences among conditions, and the results were filtered for only those genes with Benjamini-Hochberg false-discovery rate (FDR) adjusted p-values less than or equal to 0.05.

To visualize spatial differentiation among genotypes and sexes, R with RStudio was used with edgeR and limma libraries to create a principal component analysis (PCA). Non-zero values were filtered from the dataset, and then two standardized and centered PCAs were performed, clustered by genotype and sex separately, and graphed using ggplot2. The gene expression profiles in Figs. 1 and 2 were adapted from code published by Duncan et al. [19]. Heatmaps were generated using ComplexHeatmap [20].

The gene ontology (GO; http://geneontology.org/↗) enrichment of DEGs was implemented using the hypergeometric test, in which p-value is calculated and adjusted as a q-value. Hierarchical GO terms with q < 0.05 were considered significantly enriched. Ingenuity Pathway Analysis (IPA) from QIAGEN was used to perform the functional enrichment pathway and upstream transcriptional regulator analysis https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/↗). All differentially expressed genes (log₂ fold change > 1.5, FDR < 0.05) were used for the analysis. The z-score was used to infer the activation state which was calculated using the log2 fold change values of each gene.

The genotype contrast revealed more DEGs in the males than in the females. Specifically, 1187 DEGs (224 upregulated, 963 downregulated) were found in male 6mGHRKO mice versus controls, whereas female 6mGHRKO mice showed only 138 DEGs (2 upregulated, 136 downregulated), an 8.6-fold higher DEG count in males. In the sex comparisons, 1368 genes (876 upregulated, 492 downregulated) were differentially expressed between control females and males, compared to 940 genes (763 upregulated, 177 downregulated) in the 6mGHRKO female vs. male comparison, representing a 1.4-fold reduction in sex-associated DEGs following GHR deletion. Principal component analysis (PCA) illustrated clear segregation of groups based on both genotype and sex (Figs. 1b and 2b), highlighting distinct expression profiles.

To further explore genotype-related gene expression differences in both sexes, DEGs were categorized into upregulated (Figs. 1c–f), downregulated (Figs. 1g–j), and divergent (Figs. 1k–m) subsets in 6mGHRKO relative to controls. The top 20 genes with the highest or lowest log₂ fold change (LFC) values in 6mGHRKO males and females versus controls are displayed in Fig. 1c, e and g, and 1i, with complete gene lists, corresponding LFCs, and false discovery rates (FDR) detailed in Supplementary Tables 2–5.

In the genotype comparisons, male mice displayed a greater number and proportion of uniquely expressed genes in AT following GHR deletion than females. Specifically, 223 genes (representing 99% of all upregulated DEGs) were uniquely upregulated in males, whereas only one gene -Corneodesmosin, Cdsn- (0.44%) was uniquely upregulated in females (Fig. 1d; Suppl. Tables 2 and 3). Similarly, 942 genes (87% of downregulated DEGs) were exclusively downregulated in males, while females showed 115 unique downregulated genes (11%) (Fig. 1h; Suppl. Tables 4 and 5). Notably, only one upregulated gene (0.44% of the total) (Fig. 1b and f) and 21 downregulated genes (1.9%) (Fig. 1h and j) were common to both sexes. Interestingly, a single gene, Immunoglobulin Lambda Variable 2 (Iglv2), displayed opposite regulation—upregulated in males and downregulated in females—following GHR ablation (Fig. 1l and m).

To evaluate how GHR ablation affects sex-specific gene expression, DEGs in females versus males were analyzed within each genotype group (control and 6mGHRKO), focusing on upregulated, downregulated, shared, and divergent genes. Notably, the proportion of unique DEGs between female and male AT was lower in the 6mGHRKO group than in controls. Specifically, female controls showed 388 uniquely upregulated genes (34% of total upregulated DEGs), while 6mGHRKO females had 275 (24%) when compared to their respective male counterparts (Fig. 2a and d; Suppl. Tables 6 and 7). Similarly, female controls exhibited 437 unique downregulated genes (71% of total), compared to only 122 (20%) in 6mGHRKO females (Fig. 2a and h; Suppl. Tables 8 and 9). The top 20 DEGs exclusive to either control or 6mGHRKO females in comparison to males are highlighted in Fig. 2c, e and g, and 2i.

Among the DEGs shared between control and 6mGHRKO females versus males, 488 genes (42% of upregulated) and 55 genes (9% of downregulated) were commonly regulated (Fig. 2d, f and h, and 2j). To identify divergent expression patterns, DEG lists from both genotypes were compared (Fig. 2k and n). Only four genes showed inverse regulation—being downregulated in control females but upregulated in 6mGHRKO females when compared to males (Fig. 2k and l). Additionally, one gene, Iglv2, was uniquely upregulated in control females and downregulated in 6mGHRKO females (Fig. 2m and n).

To assess the biological and functional implications of the differentially expressed genes (DEGs), we performed Gene Ontology (GO) enrichment analysis on the DEG datasets derived from both genotype and sex comparisons (Suppl. Tables 2–7). Genes were categorized under three GO domains: biological process, molecular function, and cellular component. Figure 3 highlights the top five enriched terms for each category across both comparisons.

In the comparison between male 6mGHRKO and control mice, the most significantly enriched terms were 'G protein-coupled receptor signaling pathway' under biological processes, and 'Transmembrane signaling receptor activity' under molecular functions (Fig. 3a, yellow and red bars). Conversely, the comparison of 6mGHRKO vs. control female mice yielded fewer DEGs (see Table 1), which translated to lower GO enrichment significance (x-axis, Fig. 3b). Notably, the enriched terms in this group were predominantly associated with ion signaling, including processes such as 'Sodium ion transmembrane transport', 'Inorganic cation import across plasma membrane', and 'Iodide transport', as well as related functions like 'Sodium ion transmembrane transporter activity' and 'Iodide transmembrane transporter activity' (Fig. 3b, yellow and red bars). Additionally, GO terms such as 'Negative regulation of testosterone' and 'Estrogen response element binding' were notably enriched in female samples.

In sex-based comparisons within genotypes, GO analysis revealed enrichment in processes related to extracellular matrix (ECM) remodeling, immune function, and ion transport in both control and 6mGHRKO mice. For control animals, prominent biological processes included 'Anatomical structure development' and 'Complement activation' (Fig. 3c, yellow bars), while enriched molecular and cellular components featured 'Calcium ion binding', Extracellular matrix structural constituent', and 'Complement C1 complex' (Fig. 3c, red and blue bars). In the 6mGHRKO sex comparison, 'Developmental process' and 'Protein binding' emerged as the top enriched biological process and molecular function, respectively (Fig. 3d, red and blue bars). This group also displayed enrichment in ECM and ion-related terms such as 'Sodium ion transport', Calcium ion binding', 'Metal ion transmembrane transporter activity', 'Collagen trimer', and 'Cell junction' (Fig. 3d). Notably, the term 'Cell periphery' was consistently the most enriched cellular component across all GO analyses (Fig. 3a–d, blue bars).

Pathway analysis was used to investigate the pathways that were significantly activated, repressed, or modulated by the DEGs of both genotype (Fig. 4a and b) and sex (Fig. 4c and d) pairwise comparisons. Interestingly, the pathways of the genotype comparison showed either a repressed or neutral status. Specifically, 14 out of the 20 most significant altered pathways showed a repressed status in the male 6mGHRKO mice when compared to controls (Fig. 4a). Furthermore, the suppressed pathways are associated with either extracellular matrix remodeling (ECM) and fibrosis or vesicle formation. For example, while the most significantly altered pathway was 'Hepatic Fibrosis/Hepatic Stellate Cell Activation', some of the pathways that showed repressed status and are associated with ECM remodeling are 'Tumor Microenvironment Pathway' and 'Wound Healing Signaling Pathway' (Fig. 4a). Similarly, the pathways that showed repressed status and are associated with vesicle formation are 'Phagosome Formation', 'SNARE Signaling Pathway', and 'Actin Cytoskeleton Signaling' (Fig. 4a). As expected, because of the small number of DEGs between female 6mGHRKO and control mice, the pathway analysis did not show a strong significance (-log (p-value)) or clear activation status. Nevertheless, 10 of the 20 most significantly altered pathways in 6mGHRKO female mice were associated with amino acid or nucleotide metabolism; some of these pathways are: 'Arginine Degradation I (Arginase Pathway)', 'Urea Cycle Arginine Degradation VI (Arginase 2 Pathway)', 'Tryptophan Degradation to 2-amino-3-carboxymuconate Semialdehyde', 'Adenine and Adenosine Salvage III', 'Purine Ribonucleosides Degradation to Ribose-1-phosphate', 'NAD biosynthesis II (from tryptophan)', 'γ-glutamyl Cycle', 'Adenosine Nucleotides Degradation II' and 'Purine Nucleotides Degradation II (Aerobic)' pathway (Fig. 4b).

Pathway analyses for the sex comparisons showed that female controls and female 6mGHRKO mice when compared to male controls and male 6mGHRKO, respectively, showed activation of pathways associated with ECM remodeling and fibrosis (Fig. 4c and d). That is, pathways such as 'Pulmonary Fibrosis Idiopathic Signaling Pathway', 'Wound Healing Signaling Pathway', 'Hepatic Fibrosis/Hepatic Stellate Cell Activation', 'Cardiac Hypertrophy Signaling (Enhanced)' and 'Molecular mechanisms of cancer' are made up mostly by genes involved in ECM remodeling (Fig. 4c and d). Other pathways that are also active in females compared to males are 'IL15 production' and 'GP6 signaling pathway', while 'RhoGDI signaling' pathway' was repressed and' WNT/B-catenin signaling' and 'Phagosome formation' pathways were inversely active in female vs. male comparison in both genotypes (Fig. 4c and d).

To further assess differential gene expression among the WT male, 6mGHRKO male, WT female, and 6mGHRKO female groups, we generated a gene expression profile highlighting genes significantly altered in either the control male vs. 6mGHRKO male comparison or the control female vs. 6mGHRKO female comparison (Fig. 5a). These genes were categorized into four distinct expression pattern groups (Fig. 5b), followed by overrepresentation analysis (ORA) to identify enriched biological functions. Representative pathways associated with each gene group are shown in Fig. 5c. We also performed differential gene expression for the sex comparison (females vs. males of controls and 6mGHRKOmice). We generated a gene expression profile highlighting genes significantly altered in either the control female s vs. control male comparison or the 6mGHRKO females vs. 6mGHRKO males comparison. These genes were categorized into four distinct expression pattern groups, followed by overrepresentation analysis to identify enriched biological functions. Representative pathways associated with each gene group are shown in Fig. 5.

In the alternate gene expression analysis on the genotype comparison, the largest cluster of genes (544 genes) had high expression in both control groups and low expression in both 6mGHRKO groups, and these genes were associated with the "olfactory transduction" and "biogenic amine synthesis" pathways. In contrast, a cluster of 244 genes that has high expression specifically in control males was associated with the "Biosynthesis of DHA-derived SPMs" pathway, while the other two clusters did not have significantly enriched pathways.

The alternate gene expression analysis on the significant genes from the sex comparison yielded a similar result, with a cluster of genes that are specifically high in control males being associated with the "Biosynthesis of DHA-derived SPMs" pathway, but also with the Complement and coagulation cascades" pathway. A cluster of 334 genes that had high expression in both male groups was also associated with the "Complement and coagulation cascades" pathway, but also with the "Innate immune system" and "Microglia pathogen phagocytosis" pathways. The two remaining clusters were high in both female groups, and were associated with "Axon guidance" and "Ephrin signaling" pathways or the "ECM receptor interaction", "Signaling by receptor tyrosine kinases", and "Alpha 6 beta 4 integrin signaling" pathways.

Ingenuity Pathway Analysis (IPA) was utilized to pinpoint upstream transcriptional regulators potentially responsible for the differential gene expression observed in both genotype and sex-based analyses. These upstream regulators are defined as molecules (such as transcription factors) that influence the expression of other genes. Table 2 highlights the top five most significantly impacted upstream regulators for each comparison. In the genotype analysis, progesterone emerged as the most significantly affected regulator in male 6mGHRKO mice compared to controls. In contrast, beta-estradiol exhibited a predicted inhibition pattern in female 6mGHRKO mice relative to controls. Notably, beta-estradiol also stood out in the sex-based comparison, showing a consistent activation in females regardless of genotype. Additionally, several of the regulators predicted to differ between sexes—such as TGFB1, TNF, and ADAMTS18—are linked to extracellular matrix (ECM) remodeling processes (Table 2).