A Systems Biology Approach to Investigating Sex Differences in Cardiac Hypertrophy

All mice were treated and cared for in accordance with the Guide for the Care and Use of Laboratory Animals (National Institutes of Health, revised 2011), and protocols were approved by the National Heart, Lung, and Blood Institute institutional animal care and use committee. Male and female C57BL/6 mice (12–14 weeks old, obtained from the Jackson Laboratory) were given angiotensin II at 1.5 mg/kg per day or saline (the vehicle) via Alzet minipumps for 2 or 3 weeks. In some studies, GW6471, a potent inhibitor of PPARα,²⁰ was also given at 4 mg/kg per day. GW6471 has an IC₅₀ of 240 nmol/L and has been shown to function as an antagonist in mice within the range of 2 to 10 mg/kg per day.^{21, 22} Following treatment, echocardiography was performed on the mice. Mice were then anesthetized and euthanized. Their heart weights and tibia lengths were recorded, and harvested hearts were snap frozen and stored in liquid nitrogen.

Transthoracic echocardiography was performed using a high‐frequency linear array ultrasound system (Vevo 2100, VisualSonics) and the MS‐400 transducer (VisualSonics) with a center operating frequency of 30 MHz, broadband frequency of 18 to 38 MHz, axial resolution of 50 μm, and footprint of 20×5 mm. M‐mode images of the left ventricle were collected from the parasternal short‐axis view at the midpapillary muscles at a 90° clockwise rotation of the imaging probe from the parasternal long‐axis view. From the M‐mode images, the left ventricle systolic and diastolic posterior and anterior wall thicknesses and end‐systolic and ‐diastolic internal left ventricle chamber dimensions were measured using the leading‐edge method. Left ventricle functional values of fractional shortening and ejection fraction (EF) were calculated from the wall thicknesses and chamber dimension measurements using system software. Mice were lightly anesthetized with isoflurane delivered via a nose cone. The mice were imaged in the supine position while placed on a heated platform equipped with ECG leads.

0.5 mL TRI reagent was added to heart tissue along with Precellys (Bertin Technologies) homogenizing beads. Homogenization was carried out in a Precellys homogenizer (Bertin Technologies) chilled with liquid nitrogen. The samples were spun twice at 5000 rpm for 30 seconds for each cycle. RNA was isolated according to the TRI reagent protocol provided by Life Technologies. The isolated RNA samples were further cleaned with the miRNeasy Kit, according to the protocol provided by Qiagen. Samples were then treated with DNase (Ambion) and further cleaned with another Qiagen miRNeasy column. RNA concentration was determined by optical density at 260 nm.

The RNA sequencing libraries were constructed using a TruSeq Stranded Total RNA Sample Preparation Kit (Illumina), according to the manufacturer's protocols. Briefly, the ribosomal RNA was removed using Ribo‐Zero (Illumina) rRNA removal beads. The resulting RNA was then fragmented using divalent cations under elevated temperature. The RNA fragments were copied into first‐strand cDNA using reverse transcriptase and random hexamers. After second‐strand synthesis, double‐stranded cDNAs were ligated with Illumina adaptors. The final RNA sequencing library was enriched by low‐cycle polymerase chain reaction and sequenced with paired 50‐bp reads on an Illumina HiSeq 2000.

The raw data in fastq format were aligned to the mouse reference Ensembl GRCm38 genome using TopHat2²³ (tophat/2.0.13, Bowtie1/2.2.3, and samtools/0.1.19) with default settings except for the parameter –g 1. For transcript‐level analysis, the raw counts of the transcripts in the mm10_refSeq.bed, which were produced by the software RSeQC/2.6, were used as the input for the Bioconductor edgeR package. TMM (Trimmed mean of M values) algorithms were used to normalize read counts across all 24 samples. The lowly expressed transcripts were discarded by requiring a count per million >1 in at least 3 samples. The euclidean distance metric was used for principal component analysis. The normalized and log₂ transformed output count‐per‐million values were compared with a generalized linear model among 4 conditions. The differentially expressed transcripts were defined as ≥2‐fold changes with a 10% false discovery rate (FDR).

For gene‐level analysis, read counts were generated for each gene by HTSeq software²⁴ using the UCSC RefSeq annotation downloaded from iGenome in the “union” model. The statistical significance of differentially expressed genes was evaluated using edgeR.²⁵ The sex difference in hypertrophy was tested with edgeR's interaction (sex×disease) generalized linear model.

Functional enrichment analysis of differentially expressed transcripts was carried out using GOstats (https://www.bioconductor.org/packages/release/bioc/html/GOstats.html), and the results were summarized and visualized using REVIGO²⁶ (http://revigo.irb.hr/).

Transcription factor families whose binding sites were overrepresented in the promoter regions of the genes of interest were identified using Genomatix (https://www.genomatix.de/).

Coexpression networks shared in all 24 samples were constructed using the WGCNA package (weighted correlation network analysis, https://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/).²⁷ WGCNA identifies robust clusters of highly correlated genes that serve as modules in coexpression networks and allows users to further examine correlation of network modules with experimental conditions and external traits. Global coexpression network construction is computationally intensive; therefore, to improve efficiency, we included only genes that showed differential expression at an FDR <10% and |log₂ fold change|>log₂(1.2), in at least 1 of the 4 comparisons (baseline male and female, male control versus hypertrophy, female control versus hypertrophy, and male and female hypertrophy), totaling 4422 genes.

To examine the relevance to experimental conditions (sex and hypertrophy) of the resultant coexpression network modules, we constructed four 24‐dimension condition‐specific vectors (for male control, male hypertrophy, female control, and female hypertrophy, respectively). In the vectors, each element represents 1 sample, with a value of 1 if the sample belongs to the condition and 0 otherwise. The Pearson correlation between eigengenes of network modules and the condition vectors was then calculated for association estimation.

Protein–protein interaction annotation was downloaded from STRING (version 9),²⁸ a database of known and predicted protein–protein interactions (http://string-db.org/). Only interactions with a confidence score >700 were kept. Protein–protein interaction subnetworks with gene expression variations significantly associated with a factor (sex, hypertrophy, and the interaction between sex and hypertrophy) were identified using jActiveModulesTopo, a software package for trait‐relevant subnetwork identification that takes network topology into consideration.²⁹

We developed a 3‐step pipeline that determines regulatory relationships between microRNA (miRNA) and mRNA by integrating sequence‐based prediction and experimental condition‐dependent correlation in expression of miRNA and mRNA. First, TargetScan (http://www.targetscan.org/mmu_61/) was used to predict miRNA–mRNA regulatory relationships. TargetScan identified likely target mRNAs for each given miRNA by searching mRNAs for the presence of conserved sequence sites that matched the miRNA's seed regions. Second, correlation in expression variations of mRNA and miRNA were integrated with TargetScan predictions through Lasso (least absolute shrinkage and selection operator) regression to further narrow the list of candidate regulatory pairs. An R package developed by Lu et al³⁰ was used to carry out this step. Briefly, let y_i denote the expression level of the ith mRNA and x_k,i the expression level of its k‐th out of a total number (n) of candidate targeting miRNA in the simple linear regression model for y_i is given by:In equation (1), β represents the regression coefficients, and ε_i is a random error term. The lagrangian form of Lasso regression³¹ of this model is as follows:In equation (2), λ is the lagrangian multiplier, and the regression is to find coefficients {β_k,i} that minimize the value of equation (2). Lasso performs the regression with the sum of the absolute value of the regression coefficients constrained to balance between improving prediction accuracy and avoiding overfitting.

Third, the predicted miRNA–mRNA pairs were further filtered, and only those showing differential expression in 1 of the 3 statistical models—sex with hypertrophy controlled, hypertrophy with sex controlled, and the interaction between sex and hypertrophy—were retained. FDRs of 20% and 40% were used as the thresholds for mRNA and miRNA, respectively.

Although previous studies have reported sex differences in hypertrophy,^{6, 7, 8, 32, 33, 34, 35} the mechanistic basis for this difference is poorly understood. In this study, we measured global changes in mRNA and used a systems biology approach¹⁸ to determine pathways and gene networks responsible for these sex differences. We first confirmed previous studies showing sex differences in hypertrophy. We initially treated male and female mice with angiotensin II for 3 weeks. As shown in Figure 1, after 3 weeks of angiotensin II treatment, male mice exhibited significantly more hypertrophy than female mice (Figure 1A). Furthermore, after 3 weeks of treatment, EF in female mice was not different compared with the control value (54% versus 56% at baseline). In male mice, however, EF dropped from the baseline level of 55% to 37% (P<0.05), which was significantly lower than EF in female mice after 3 weeks of treatment (Figure 1B). We also measured the percentage increase in aortic velocity (mean and peak) following angiotensin II treatment, and we observed similar increases in both male and female mice (36% increase in females and 29% increase in males).

We next used RNA sequencing to examine sex differences in mRNA at baseline and after 2 weeks of angiotensin II treatment. We measured mRNA after 2 weeks of treatment because at that time, sex differences in EF had not occurred (65±9% in females versus 61±9% in males, P=0.53, n=5–7). Thus, the changes in mRNA that occur at 2 weeks of treatment are more likely to be a cause rather than an effect of sex differences. Data were analyzed using pairwise comparisons of the conditions. As shown in Figure 2, principal component analysis revealed not only clustering differences between baseline and hypertrophy but also discrete clustering for mRNA of male and female mice at baseline that persisted following hypertrophy. The full data set of mRNA, with significant changes, is provided in Table S1. The volcano plots in Figure 3 show differences in expression of mRNA between the groups. As illustrated in Figure 3A, comparing male and female mice at baseline, 174 transcripts showed a significant sex difference (using an FDR <10% and log₂ fold change >1). Male mice tended to have more upregulated than downregulated genes compared with female mice (127 versus 47, respectively). Hypertrophy led to differential expression of mRNA in both male and female mice (Figure 3B through 3D). The mRNA changes that occur with hypertrophy largely overwhelm the baseline sex difference such that with hypertrophy, fewer mRNAs show a sex difference (compare Figure 3A and 3B). These differences are shown in more detail in the Venn diagram in Figure 4A. At an FDR <10% and |log₂(fold change)| >1, there were 174 mRNAs (141 unique to baseline +33 overlapping with the list from hypertrophy) that exhibited a sex difference at baseline; however, following hypertrophy, only 71 mRNAs showed a sex difference, with 33 of these mRNAs showing a sex difference at both baseline and following hypertrophy. For the 141 transcripts that showed a sex difference only at baseline, 61 did not display a change with hypertrophy in either sex compared with their corresponding baselines, 10 displayed a hypertrophy‐induced change in both sexes, 63 changed only in female hypertrophy, and 7 changed only in male hypertrophy.

If we examine differences between baseline and hypertrophy as a function of sex, we find that 328 mRNAs changed with hypertrophy in female mice, whereas only 174 mRNAs changed with hypertrophy in male mice. There are 117 mRNAs that changed in both male and female mice with hypertrophy; thus, 211 mRNAs are altered only in females with hypertrophy compared with 57 mRNAs that change only in males (see Figure 4B). Taken together, these data suggest that the process of hypertrophy overwhelms the sex differences observed at baseline.

To confirm the validity of our model, we confirmed an elevation in natriuretic peptide B (Figure 4C) and skeletal α‐actin 1 (Figure 4D) with hypertrophy. Interestingly, although there were sex differences in the development of hypertrophy and in cardiac function following angiotensin II treatment, consistent with the RNA sequencing data, both male and female hearts showed significant increases in brain natriuretic peptide (Figure 4C). We found a similar discrepancy in a transaortic constriction model of hypertrophy in which we also observed a sex difference in hypertrophy but found no sex difference in natriuretic peptide B levels in heart.⁶ We also confirmed the sex differences that were found in the transcriptomic analysis for metallothionein 1 and tissue metalloproteinase inhibitor 4 (see Figure S1).

Coexpression network analysis identifies gene sets showing significant and robust coregulations under the conditions of interest (hypertrophy and sex), manifested as modules in coexpression networks. Compared with the conventional pathway enrichment analysis that examines gene sets belonging to predefined pathways, it has the advantage of uncovering condition‐specific gene signatures. In this study, we used WGCNA²⁷ to construct the coexpression network of the 4422 genes showing significant variations in response to sex, to hypertrophy, or to the interaction of sex and hypertrophy. Using all 24 samples, 20 coexpressed modules were identified. Figure 5 shows the module–condition correlation. The complete list of 4422 genes and their module memberships are given in Table S2. Coexpression network construction using gene lists filtered at other statistical stringencies or by variance in expression yielded similar module composition (data not shown).

The light yellow module with 25 mRNAs showed the most significant condition‐specific differences along with the highest correlation and P value. At baseline, genes in this module were highly expressed in male but not in female mice. With hypertrophy, the sex difference diminished, and the genes were moderately expressed in both sexes. Consequently, this module displayed a sex difference in changes with hypertrophy. The module contained several mRNAs involved in metabolism including PPARα, Ccrn4l (carbon catabolite repression 4‐like), Fitm1 (fat storage‐inducing transmembrane domain 1), Acot1 (acyl–coenzyme A thioesterase 1), and Sik1 (salt inducible kinase 1). Several transcription factors (in addition to PPARα) and splicing factors were also contained in this module (Per1 [period circadian clock 1], SRSF3 [serine rich splicing factor 3], and Siah2 [siah E3 ubiquitin protein ligase 2]).

To dissect the contributions from sex and hypertrophy and their interactions, we identified genes that showed significant changes in hypertrophy when sex as a cofactor was controlled, significant sex‐dependent variation when hypertrophy as a cofactor was controlled, and significant sex‐dependent differential hypertrophy‐associated changes. This analysis was done at the gene level to avoid complications due to multiple splicing variants. At an FDR <40%, this approach gives 379 genes with significant interactions. The complete lists of genes with the statistics and pathway enrichment analysis results are available in Table S3. The pathway results were further simplified by removing redundancy in the tree‐structured Gene Ontology terms and visualized using REVIGO (Figure 6. Four main categories were identified: cell‐cycle and growth‐related biosynthesis, cellular development process, extracellular matrix organization, and response to organic substances.

To gain further insight into the mRNA networks that played a role in the sex‐dependent differences in hypertrophy‐induced gene expression changes, we selected genes that were significant at an FDR <40% for the interaction and mapped them to the protein–protein interaction network constructed using STRING data. This provided a network with 178 genes and 313 interactions. We used an in‐house software, jActiveModulesTopo^{29, 36} to identify the subnetworks that were most relevant to the sex–hypertrophy interaction (ie, connected sets of genes with high levels of sex–hypertrophy difference), using a simulated annealing method and setting the search depth at 2, and the result is given in Figure 7, with 17 genes and 22 interactions. This subnetwork is centered on PPARα, a clear hub with 8 interactions with other members, whereas the interactions for other genes ranged from 1 to 4. Note that PPARα is also implicated in the coexpression network analysis, being a member of a module that exhibited the most significant sex difference at baseline but not at hypertrophy (Figure 5, light yellow module).

PPARα has been shown previously to be involved in regulating hypertrophy.^{37, 38, 39, 40, 41, 42} Esrrg (estrogen‐related receptor γ), a regulator of mitochondrial function that has been shown to play a role in regulating hypertrophy, is also found in this subnetwork. Intriguingly, a number of genes involved in circadian rhythms were also shown to have a sex and disease interaction; these include Per1, Arntl (aryl hydrocarbon receptor nuclear translocator‐like protein 1; also known as Bmal [brain and muscle ARNT]), and Ccrn4ls. Among them, Per1 and Ccrn4l were also members of the light yellow module shown in Figure 5. The expression levels of the key genes in this subnetwork that showed a sex bias are given in Figure S2. We also confirmed a sex difference in PPARα level with quantitative polymerase chain reaction (Figure S1).

It is known that mRNAs are regulated in groups or networks by common transcriptional regulators, such as miRNAs; therefore, we performed an analysis of miRNAs on sham‐ and angiotensin II–treated male and female hearts after 2 weeks of treatment, using the same hearts that were used for mRNA measurement. The volcano plots illustrating miRNA differences between male and female hearts at baseline, male control and hypertrophy hearts, female control and hypertrophy hearts, and male and female hearts with hypertrophy are shown in Figure S3. The full data set is available in Table S4. Consistent with previous data in the literature and supporting the validity of the model, we found increases in miRNAs 15b, 21, 34, 199, 208b, and 214, which have all been reported previously to change with hypertrophic stimuli.^{43, 44, 45, 46, 47}

We further performed an integrative analysis of mRNA–miRNA interaction. As described in the Methods, a Lasso regression model was used for identification of miRNA–mRNA targeting relationships that combine sequence‐based prediction and experimental condition–dependent correlation in miRNA/mRNA expression variations. This was designed to overcome the problem of high false‐positive rates in sequence‐based predictions. This filtering with dynamic information of biological context‐specific interactions improves reliability.

In total, TargetScan predicted 72 311 miRNA–mRNA targeting pairs, and less than one‐third (23 323 pairs) remained after filtering with the Lasso regression. We then used cytoscape/Partek to generate a network of miRNA–mRNA targeting pairs most relevant to sex and hypertrophy interaction, as shown in Figure 7B. This analysis demonstrates sex and hypertrophy regulation of miRNAs 208b and 124. TargetScan shows miRNA 124 as a regulator of PPARα.

PPARα was found to be at the center of the network of mRNAs that were significantly different based on sex and hypertrophy (Figure 7); therefore, we examined whether PPARα contributes to sex differences in cardiac hypertrophy. To test whether PPARα regulates sex differences in hypertrophy, we examined whether inhibition of PPARα would block the sex differences observed with hypertrophic stimuli.

Cardiac hypertrophy was induced by treating male and female mice for 3 weeks with angiotensin II with and without GW6471, an inhibitor of PPARα. Consistent with the data in Figure 1, a sex difference in hypertrophy was observed with angiotensin II treatment (Figure 8). Interestingly, this sex difference in cardiac hypertrophy was blocked by treatment with GW6471 (Figure 8). These results are consistent with the hypothesis that PPARα contributes to sex differences in cardiac hypertrophy.