- 1Department of Cardiovascular Medicine, The Second Affiliated Hospital of Nanchang University, Nanchang, China
- 2Department of Science and Technology, The Second Affiliated Hospital of Nanchang University, Nanchang, China
- 3Jiangxi Key Laboratory of Molecular Medicine, The Second Affiliated Hospital of Nanchang University, Nanchang, China
Hypertrophic cardiomyopathy (HCM) is an autosomal dominant disease and mitochondria plays a key role in the progression in HCM. Here, we analyzed the expression pattern of nuclear-encoded mitochondrial genes (NMGenes) in HCM and found that the expression of NMGenes was significantly changed. A total of 316 differentially expressed NMGenes (DE-NMGenes) were identified. Pathway enrichment analyses showed that energy metabolism-related pathways such as “pyruvate metabolism” and “fatty acid degradation” were dysregulated, which highlighted the importance of energy metabolism in HCM. Next, we constructed a protein-protein interaction network based on 316 DE-NMGenes and identified thirteen hubs. Then, a total of 17 TFs (transcription factors) were predicted to potentially regulate the expression of 316 DE-NMGenes according to iRegulon, among which 8 TFs were already found involved in pathological hypertrophy. The remaining TFs (like GATA1, GATA5, and NFYA) were good candidates for further experimental verification. Finally, a mouse model of transverse aortic constriction (TAC) was established to validate the genes and results showed that DDIT4, TKT, CLIC1, DDOST, and SNCA were all upregulated in TAC mice. The present study represents the first effort to evaluate the global expression pattern of NMGenes in HCM and provides innovative insight into the molecular mechanism of HCM.
Introduction
Hypertrophic cardiomyopathy (HCM) is an autosomal dominant genetic disease that is mainly characterized by ventricular hypertrophy with asymptomatic or serious complications such as sudden cardiac death (SCD), heart failure, and thrombosis (Marian and Braunwald, 2017). The prevalence of HCM in the general population was estimated to be 1/500 (Gersh et al., 2011), which was underestimated due to the limited HCM diagnostic technology. HCM is considered a leading cause of SCD in younger people and the leading cause of heart failure in cardiac diseases originating primarily from the myocardium (Weissler-Snir et al., 2019).
Normal myocardial energy metabolism from mitochondria is also an important material basis for keeping the normal heart tissue structure and the internal environment stable. Cardiac function will inevitably be impaired by mitochondrial dysfunction. Clinical and experimental studies have shown that the myocardial energy source switching from fatty acid oxidation to glycolysis is a common event in HCM (Tian, 2003). Mutations in a wide spectrum of nuclear-encoded mitochondrial genes (NMGenes) have been reported to be able to cause HCM characterized by impaired mitochondrial function (Marin-Garcia and Goldenthal, 2002b). For example, mutations in ELAC2 (ElaC ribonuclease Z 2) encoding a short form of RNase Z were found to be associated with HCM (Saoura et al., 2019). Mitochondrial function depends on proteins encoded by both mitochondrial DNA (mtDNA) and nuclear DNA (nDNA). The mitochondrial proteome has been estimated to contain approximately 1000–1500 proteins, more than 99% of which are encoded by nuclear DNA (nDNA), while mtDNA refers to only 13 protein-coding genes (Pfanner et al., 2019). Considering the importance of mitochondria in HCM and the fact that functional proteins in mitochondria are encoded mainly by nDNA genes, exploring the function of NMGenes in HCM would help us better understand the novel role of mitochondria in the development of HCM.
With the development of genetic studies, high-throughput omics technologies (such as DNA microarrays and next-generation sequencing) that investigate gene function and expression at the genome-wide level have been widely used in basic research, clinical diagnosis, drug research and other fields. As a powerful technique, gene expression microarray-based bioinformatics analyses have also been widely used to identify HCM-related genes or noncoding RNAs, possible molecular mechanisms, and biological pathways (Lim et al., 2001; Yang et al., 2015; Hu et al., 2019; Li et al., 2019; Liu et al., 2019). For example, microarray analysis was performed to explore the expression pattern of lncRNAs (long noncoding RNAs) and mRNAs (messenger RNAs) in HCM, which identified hundreds of differentially expressed lncRNAs and genes (Yang et al., 2015). A recent study systemically analyzed RNA-seq data from 28 HCM patients and 9 healthy controls and identified 43 potential pathogenic variants in 19 genes and four subnetworks with significant roles in the progression of HCM (Gao et al., 2020). Although previous studies have highlighted the importance of integrative gene expression analysis in exploring the molecular mechanism of HCM, a systemic analysis of the expression pattern of NMGenes in HCM patients has never been reported.
To investigate the potential role of NMGenes in the pathogenesis of HCM, in this study, we performed a computational systems biology analysis based on large-scale HCM-related transcriptional data. A total of 316 differentially expressed NMGenes (DE-NMGenes) were identified. Based on these DE-NMGenes, gene ontology (GO) and pathway enrichment analyses were performed, and 17 KEGG-dysregulated pathways were identified. We also constructed a PPI (protein-protein interaction) network that consisted of 215 DE-NMGenes and 440 interactions. Finally, a total of 17 TFs (transcription factors) were predicted to potentially regulate the expression of the 316 DE-NMGenes. We provided a systematic view of the roles of mitochondrial genes in HCM and revealed some available candidates for future experimental verification.
Results
Nuclear-Encoded Mitochondrial Genes Are Significantly Changed in HCM
The normalized gene expression dataset GSE36961 was downloaded from the GEO (Gene Expression Omnibus) database1, which included 107 HCM samples and 40 control samples (Clough and Barrett, 2016). Differentially expressed genes (DEGs) between HCM and the corresponding control samples were detected using the “Limma” package from R software (Ritchie et al., 2015). By keeping genes with a BH (Benjamini-Hochberg)-corrected p-value less than 0.01 and fold change (FC) larger than 1.2, we obtained 2927 DEGs, 1499 of which were upregulated and the remaining 1428 were downregulated (Figure 1A and Supplementary Table 1). To explore the expression pattern of NMGenes in HCM, we collected 1943 mitochondrial genes from the MitoCarta (Calvo et al., 2016), MitoMiner (Smith and Robinson, 2019), IMPI and UniProt databases (UniProt, 2019) (Figure 1B, see section “Materials and Methods” for details). After removing 13 mtDNA-encoded genes, 1930 NMGenes were retained for further analysis. Among these genes, 1562 genes were detected on microarray, and 316 genes were differentially expressed. Compared with the overall genes detected on the microarray, the proportion of DEGs in NMGenes was significantly higher (Figure 1C, Fisher’s exact test, p-value < 2.20 × 10–16). The extensive expression changes of NMGenes in HCM indicate that mitochondria play critical roles in the progression of HCM. Table 1 lists the top ten upregulated and downregulated NMGenes in HCM. Among these genes, four upregulated genes [namely, PDK4 (thpyruvate dehydrogenase kinase isozyme 4), STAT3 (Signal Transducer and Activator of Transcription 3), HCLS1 (Hematopoietic Cell-Specific Lyn Substrate 1) and FKBP11 (FKBP Prolyl Isomerase 11)], and four downregulated genes [namely, GATM (Glycine Amidinotransferase), ATPIF1 (ATP Synthase Inhibitory Factor Subunit 1), CPT1B (Carnitine Palmitoyltransferase 1B), and GJA1 (Gap Junction Protein Alpha 1)] have already been proven to play important roles in pathological hypertrophy (summarized in Table 1).
Figure 1. The expression pattern of NMGenes in HCM. (A) A volcano plot shows DEGs in HCM detected from GSE36961. The vertical lines correspond to 1.2-fold up and down, respectively, and the horizontal line represents a BH-corrected p value of 0.05. Therefore, red and blue dots represent upregulated and downregulated DEGs, respectively. Triangular nodes and circle nodes represent NMGenes and the other genes, respectively. (B) A four-way Venn diagram shows the number of NMGenes collected from the MitoCarta, MitoMiner, IMPI, and UniProt databases. (C) A bar plot shows the proportion of DEGs in all genes on microarray and NMGenes. The percentage of DEGs in NMGenes was significantly larger than the percentage of DEGs in overall genes detected on microarray (Fisher’s exact test, p-value < 2.20 × 10–16). *** represents P-value < 0.001.
Downregulated DE-NMGenes Are More Functionally Diverse Than Upregulated DE-NMGenes
GO biological process (BP) and KEGG pathway enrichment analyses for 316 DE-NMGenes were performed using DAVID (Database for Annotation, Visualization, and Integrated Discovery) (da Huang et al., 2009). Although the numbers of upregulated and downregulated DE-NMGenes were similar, downregulated DE-NMGenes were more functionally diverse than upregulated DE-NMGenes. By keeping terms with BH-corrected p-values less than 0.05, we obtained 4 GO BP terms and 4 KEGG pathways for 141 upregulated DE-NMGenes and 16 GO BP terms and 17 KEGG pathways for 175 downregulated DE-NMGenes (Figures 2A,B and Supplementary Table 2). The top 3 enriched GO terms in downregulated DE-NMGenes were “oxidation-reduction process,” “branched-chain amino acid catabolic process” and “fatty acid beta-oxidation.” The GO terms “oxidation-reduction process” and “translation” were both enriched in 141 upregulated and 175 downregulated DE-NMGenes. KEGG pathway enrichment analysis showed that downregulated DE-NMGenes were significantly enriched in energy metabolism-related pathways such as “Carbon metabolism” (15 genes, BH-corrected p-value = 1.94∗10–9), “Pyruvate metabolism” (7 genes, BH-corrected p-value = 1.25∗10–4), “Fatty acid metabolism” (7 genes, BH-corrected p-value = 3.24∗10–4), and “Citrate cycle” (5 genes, BH-corrected p-value = 4.35∗10–3). For upregulated DE-NMGenes, the KEGG pathways “Biosynthesis of antibiotics” (15 genes, BH-corrected p-value = 1.75∗10–5) and “Biosynthesis of amino acids” (7 genes, BH-corrected p-value = 8.88∗10–3) were significantly enriched.
Figure 2. GO and KEGG enrichment results of 141 upregulated DE-NMGenes and 175 downregulated DE-NMGenes. Bar graph of the significantly enriched BP (A) and KEGG (B). Blue and red bars represent the enriched BP terms or KEGG pathways for upregulated and downregulated DE-NMGenes, respectively.
A Group of 215 DE-NMGenes Are Biologically Connected to Form a Network
The 316 DE-NMGenes were analyzed together to construct a PPI network. Consequently, a PPI network including 440 interactions and 215 nodes was obtained by using STRING (Search Tool for the Retrieval of Interacting Genes/proteins database) (Szklarczyk et al., 2019), with parameters including a minimum required interaction score larger than 0.7 (high confidence) and only query proteins being displayed. Thus, 215 out of the 316 DE-NMGenes were included in the final PPI network (Figure 3A). The 316 DE-NMGenes had significantly more interactions than would be expected (p-value < 2.2∗10–16) from a randomly chosen set of proteins of the same size drawn from the genome. In a PPI network, highly connected nodes are called hubs, which are expected to play an important role in understanding the biological mechanism of disease (Barabasi and Oltvai, 2004). Then, we calculated the degree for each node and selected genes with the degree ranked in the top 5% as hubs. Of the 215 nodes in the PPI network, 13 nodes were ranked in the top 5% and selected as hubs (Table 2). DLD (dihydrolipoamide dehydrogenase) was the hub gene with the largest degree and interacted with 19 proteins in the PPI network.
Figure 3. Network analysis of 316 DE-NMGenes. (A) The PPI network consists of 440 edges and 215 nodes. Red and blue nodes represent upregulated and downregulated DE-NMGenes, respectively. (B) The network representation of modules detected from the PPI network. Red and blue nodes represent upregulated and downregulated NMGenes, respectively.
The MCODE (Molecular Complex Detection) plugin in Cytoscape was used to detect network modules from the PPI network (Bader and Hogue, 2003). A module is a group of closely related proteins that act in concert to perform specific biological functions through a PPI network that occurs in time and space (Lin et al., 2015). A total of 12 modules were extracted from the PPI network, of which five modules (modules 1–5) had nodes ≥ 5 (Figure 3B and Supplementary Table 3). The associated BPs for module 1, module 2, module 4 and module 5 were “mitochondrial respiratory chain complex I assembly,” “mitochondrial translation,” “protein N-linked glycosylation via asparagine” and “folic acid-containing compound biosynthetic process,” respectively, but no-GO term was significantly enriched in module 3.
TFs Potentially Regulating DE-NMGenes Play Key Roles in HCM
iRegulon (Janky et al., 2014), available as a Cytoscape plugin, was used to predict TFs potentially regulating the 316 DE-NMGenes. A total of 17 TFs were obtained with the minimum normalized enrichment score >3 and the FDR on motif similarity <0.001 (Table 3). The TF with the largest number of targets is PBX3 (pre-B-cell leukemia transcription factor 3), which regulates 145 DE-NMGenes. Eight of the 17 TFs, namely, BACH1 (BTB Domain and CNC Homolog 1), ATF3 (Activating Transcription Factor 3), XBP1 (X-Box Binding Protein 1), KLF4 (Kruppel Like Factor 4), MEF2C (Myocyte Enhancer Factor 2C), JUND (JunD Proto-Oncogene), MYC (MYC Proto-Oncogene) and YY1 (YY1 Transcription Factor), have already been proven to play important roles in pathological hypertrophy. The remaining nine genes with unknown roles in HCM were good candidates for further experimental verification.
Validation of the Differentially Expressed NMGenes in vivo
To validate the identified genes in vivo, the samples were extracted from control and transverse aortic constriction (TAC) mice to identify whether the mRNA levels of the top five genes that have not been proven to play important roles in cardiac hypertrophy were consistent with the bioinformatic analysis. In the TAC group, MYH7, ANP, and BNP expression levels were increased (Figures 4A–C), indicating that pressure overload successfully induced cardiac hypertrophy in the mouse TAC model. Interestingly, the expression of DDIT4, TKT, CLIC1, DDOST, and SNCA in the mouse TAC model were all increased compared with the sham operation group (Figures 4D–H).
Figure 4. Validation of the differentially expressed NMGenes in TAC mice. (A–C) mRNA of the genes indicating cardiac hypertrophy increased. (D–H) mRNA of the differentially expressed NMGenes also increased. Data are presented as the mean ± SD. *P < 0.05, **P < 0.01, ***P < 0.001 (Student’s t-test); n = 5 samples per group. TAC, transverse aortic constriction.
Discussion
Although many studies have been conducted to explore the pathogenesis of HCM, the role of mitochondria in HCM development and progression remains largely unknown. More than 99% of mitochondrial proteins are encoded by nDNA, so NMGenes are more responsible for mitochondrial function (Ferramosca and Zara, 2013). Exploring the expression pattern of NMGenes in HCM will help us better understand the molecular mechanism of mitochondria in HCM. Therefore, we performed a comprehensive comparative analysis of NMGenes in HCM by comparing transcriptional data in HCM patients and normal healthy controls.
Differential expression analysis showed that the proportion of differentially expressed genes in NMGenes was significantly higher than the proportion of overall genes detected on the microarray, highlighting the importance of NMGenes in HCM. For the top 10 NMGenes with the highest fold change, four upregulated genes (i.e., PDK4, STAT3, HCLS1, and FKBP11) and four downregulated genes (i.e., GATM, ATPIF1, CPT1B, and GJA1) have already been shown to play important roles in pathological hypertrophy (Table 1). Importantly, the other genes with undetermined roles in HCM are good candidates for further experimental verification. Consistent with the bioinformatic analysis, DDIT4, TKT, CLIC1 and DDOST mRNA expression increased in TAC mice compared with the sham operation group, suggesting that these genes may play an important role in promoting pathological hypertrophy. However, in contrast with the bioinformatic analysis, its expression at the mRNA level increased at 4 weeks after TAC. We speculate that it would decrease in the earlier or later time of TAC, as the duration of the pressure overload can affect the expression of associated genes (Souders et al., 2012). These DE-NMGenes provide a new perspective on the mechanisms in HCM. For example, CLIC1, as a metamorphic protein, is abundantly expressed in the heart (Ponnalagu et al., 2016); however, its function in the heart is far from fully understood. Direct evidence has shown that CLIC1 plays a significant role in ischemia-reperfusion (IR) injury by regulating reactive oxygen species (ROS) generation (Gururaja Rao et al., 2020). Previous studies have demonstrated that abnormal production of ROS in cardiomyocytes is closely related to the occurrence and development of HCM (Hafstad et al., 2013; Brown and Griendling, 2015). We speculate that CLIC1 is involved in the progression of HCM by regulating the generation of ROS and might be a potential therapeutic target for cardiac hypertrophy. DDIT4 is an inhibitor of mTOR signaling, which plays a key regulatory role in cardiovascular pathology (Sciarretta et al., 2014). It is possible that DDIT4 is involved in the progression of HCM by regulating mTOR signaling.
KEGG pathway analysis showed that abnormal expression of metabolically related pathways such as pyruvate metabolism and fatty acid metabolism (Figure 2B) may contribute to the pathogenesis of HCM. In the normal heart, mitochondrial fatty acid oxidation is the main (70–80%) source of energy, and the remaining 20–30% of ATP production derives largely from glucose oxidation (Sacchetto et al., 2019). Fatty acid metabolism disturbances are common in HCM patients, and mutations in the fatty acid oxidation pathway can result in HCM (Marin-Garcia and Goldenthal, 2002a). Fatty acid oxidation involves two key steps: fatty acid transfer and ββ-oxidation. Our results showed that genes involved in fatty acid transport (CPT1B and CPT2) or β-oxidation (ACADSB, ACADM, ACADL, and HADH) were all significantly downregulated (Supplementary Table 1). CPT1B was one of the top 10 downregulated DE-NMGenes and its deficiency could cause heart lipotoxicity, leading to exacerbated cardiac pathology (He et al., 2012). ACADM and HADH were two hubs in the PPI network. Rats with hypertrophic myocardium showed impaired fatty acid oxidation and decreased expression of ACADM and ACADL (Doenst et al., 2010). In rats with cardiac hypertrophy caused by left ventricular volume overload, HADH activity was significantly reduced (Lachance et al., 2014). Pyruvate metabolism is a key step in glucose oxidation. Compared with normal hearts, glucose oxidation was actually lower in hypertrophied hearts (Allard et al., 1997). Glucose oxidation and fatty acid oxidation are under fine regulation during disease progression, although there is still controversy, allowing us to consider the treatment of HCM from the perspective of energy metabolism. New treatments include inhibiting enzymes related to fatty acid oxidation and directly increasing the oxidation of glucose and pyruvate, which may bring light to patients with cardiomyopathy in the future.
Genes rarely act alone and usually perform their functions in connection with other genes. Moreover, genes with relatively small but significant changes in expression can also contribute to the phenotypes of interest. However, differential expression analysis focuses only on individual gene expression without considering its close connection with other genes. PPI network-based analysis might largely overcome these limitations by combining gene expression and connections. In the present study, we performed PPI network analysis and obtained 440 interactions among 215 DE-NMGenes. We found that compared with random gene sets from the genome, DE-NMGenes formed significantly more interactions, which indicated that DE-NMGenes were biologically connected to form a group. Our results identified five closely connected modules that might contribute to the development of HCM. In addition, we highlighted 13 hub genes with a high level of network connectivity but relatively modest changes in expression. Hubs DLD, CAT, ACADM, and HADH have already been proven to be involved in the progression of HCM or pathological hypertrophy, and the role of the remaining hubs in HCM deserves further investigation. The top 2 hub gene, ACLY, is an essential cytosolic enzyme for generating acetyl-CoA, a key metabolite for glucose, fatty acid, and amino acid catabolism. A mendelian randomization study found ACLY to be a promising target for cardiovascular protection (Ference et al., 2019). In addition, we also noticed that five of the 13 hubs were genes encoding mitochondrial ribosomal proteins (MRPs), which assist protein synthesis within mitochondria. These MRPs were grouped into module 2 in the network module analysis (Figure 3B). Mutations in MRPL3, MRPS14, MRPS22, and MRPL44 could cause HCM accompanied by multiorgan diseases (Smits et al., 2011; El-Hattab and Scaglia, 2016). Therefore, these MRPs are functionally connected, and the consistent downregulation of MRPL46, MRPL53, MRPL1, and MRPL40 may cause mitochondrial translation deficiency, which would result in a severe phenotype in HCM.
Generally, gene expression is under the fine turn regulation of TFs. Among the 17 TFs predicted in this work, more than half have been shown to be associated with pathological hypertrophy (Table 3). The remaining 9 TFs (i.e., ATF4, NFYA, PBX3, NFYC, GATA1, MYBL2, GATA5, and RARA were good candidates for further experimental verification. NFY (nuclear transcription factor Y) is a heterotrimeric TF complex consisting of three subunits, NFYA, NFYB and NFYC. In this work, NFYA and NFYC were predicted to regulate 75 and 117 DE-NMGenes, respectively. By analyzing the targets of NFYA and NFYC in DE-NMGenes, we found that they were both enriched in the GO term “negative regulation of apoptotic process” with p-values of 4.1∗10–7 and 2.5∗10–4, respectively. Although the role of NFY in cardiovascular disease has not been reported, NFY is involved in cancer by regulating apoptosis (Gurtner et al., 2010). Moreover, NFYA and NFYC were both significantly differentially expressed in HCM. We speculate that NFYA and NFYC may be involved in the pathogenesis of HCM by regulating apoptosis, which provides us with a new perspective to understand the relationship between NFY and HCM.
The GATA TF family comprises six members (named GATA1-6) that are involved in the regulation of growth, differentiation, survival and maintenance of body function. Previous studies have underscored the pivotal roles of the GATA family in cardiac hypertrophy (Pikkarainen et al., 2004). Mutations in GATA2, GATA4, and GATA6 were identified in patients with HCM (Alonso-Montes et al., 2017). Overexpression of either GATA4 or GATA6 could induce cardiac hypertrophy both in vitro and in vivo (Liang et al., 2001). GATA5 and GATA1 are closely related to cardiomyopathy diseases such as dilated cardiomyopathy, although their role in HCM has not yet been reported (Zhang et al., 2015). In this work, GATA1 and GATA5 were predicted to regulate 19 and 13 DE-NMGenes, respectively (Table 3). Given that the functional characteristics of GATA5 and GATA1 overlap at least partly with those of other GATA TFs and that GATA1 and GATA5 regulate DE-NMGenes, it is reasonable to speculate that GATA1 and GATA5 may contribute to HCM.
Conclusion
The present study was the first effort to evaluate the global expression pattern of NMGenes in HCM. Based on differential expression analysis, we found that NMGenes were significantly changed and identified 316 DE-NMGenes. Further GO enrichment analysis showed that downregulated DE-NMGenes were more functionally diverse. These DE-NMGenes participated in 10 significant pathways, and nine of these pathways were metabolically related. PPI network analysis showed that 13 DE-NMGenes with high node connectivity were selected as hubs. Finally, a total of 17 TFs were predicted to potentially regulate the expression of the 316 DE-NMGenes, and TFs (such as ATF4, NFYA, NFYC, GATA1, and GATA5) might play roles in HCM. This analysis will provide valuable information for future research on the molecular mechanisms of HCM and offer clues for the discovery of novel therapeutic strategies.
Materials and Methods
Data Collection
Normalized gene expression data (GSE36961) were collected from the GEO database (Clough and Barrett, 2016). NMGenes were collected from the MitoCarta (Version 2.0) (Calvo et al., 2016), MitoMiner (Version 4.0) (Smith and Robinson, 2019), IMPI2 and UniProt databases (UniProt, 2019).
Differential Expression Analysis
To identify DEGs between HCM and normal healthy hearts, limma (Version 3.40.6), an R package in Bioconductor, was utilized (Ritchie et al., 2015). Genes with BH-corrected p-values less than 0.01 and fold changes (FCs) larger than 1.2 were selected as significantly differentially expressed. We have deposited the analysis code to a public repository3.
Functional Enrichment Analysis
GO BP and KEGG pathway enrichment analyses of DE-NMGenes were performed using DAVID, an online functional annotation tool, to understand the biological significance of a list of genes (da Huang et al., 2009). In this work, GO BP and KEGG pathways with BH-corrected p-values less than 0.05 were selected as significant.
PPI Network Construction
The PPI network of DE-NMGenes was constructed using the STRING database, and an online database provides information regarding the predicted and experimental protein interactions (Szklarczyk et al., 2019). In this work, PPIs between DE-NMGenes with interaction scores larger than 0.7 were retained.
Network Module Analysis
A network module is defined as a group of genes participating in the same biological function. In this work, we detected network modules from the constructed PPI network using MCODE (Bader and Hogue, 2003), a plugin in Cytoscape4. Given the following parameters: a degree cutoff = 2, node score cutoff = 0.2, k-score = 2 and max. depth = 100, modules with scores > 3 and a number of nodes > 5 were selected. GO BP enrichment analysis of modules was performed using DAVID, and BH-corrected p values < 0.05 were selected as significant.
Prediction of TFs Regulating DE-NMGenes
To predict TFs regulating DE-NMGenes, iRegulon (Version: 1.3), a plugin in Cytoscape, was employed (Janky et al., 2014). The iRegulon plugin uses motif and track discovery in a set of coregulated genes to identify regulons. Given the following parameters: motif collection (10 kb, 9,713 PWMs), track collection (1120 ChIP-seq tracks of ENCODE raw signals), putative regulatory region (20 kb centered around TSS), motif rankings database (20 kb region centered around TSS, 7 species), track of rankings database (20 kb centered around TSS, ChIP-seq-derived), minimum identity between orthologous genes = 0 and maximum false discovery rate on motif similarity = 0.001, TFs with the NES (normalized enrichment score) larger than 3 were selected. The higher the NES was, the more reliable the TFs were.
Animals and Surgical Procedures
All experiments involving animals were approved by the Animal Ethics and Experimentation Committee of Nanchang University and carried out according to the “Guide for the Care and Use of Laboratory Animals.” Male C57BL/6 mice, aged 8 weeks and weighing 20–25g, were purchased from the SlacJingda Experimental Animals Company [Changsha, Hunan Province, China]. A total of 20 mice were divided into two groups (ten mice per group): the sham operation group and the TAC group. TAC was performed as previously described (Zhang et al., 2020). Briefly, mice were induced with 5% isoflurane and intubated orally and then maintained at 2% isoflurane during surgery with mechanical ventilation. After a midline sternotomy, the aortic arch was exposed. Constriction was performed by tying a 5-0 silk suture around a 27-gauge needle overlying the arch between the origin of the brachiocephalic trunk and left common carotid artery. For the sham operation group, 10 mice underwent the same procedure, but the suture was withdrawn without tying. Then, the thorax and skin were closed by using 6-0 polypropylene sutures. Four weeks after surgery, the mice were euthanized, and their hearts were quickly excised for further evaluation.
Quantitative Real-Time PCR Analysis
Total RNA was extracted from mouse cardiac tissues using TRIzol reagent (Invitrogen, New York, United States), and then the quality and concentration of RNA were determined using an Agilent Bioanalyzer 2100 according to the manufacturer’s instructions. The cDNAs were generated by MMLV transcriptase (BioRAD, United States), and quantitative real-time PCR assays were performed as previously described (Yu et al., 2018). Triplicate PCR amplifications were performed for each sample, and the mRNA levels were normalized to GAPDH. The comparative threshold cycle method (2-ΔΔCT) was applied to estimate the relative gene expression of cardiac tissues between the TAC and sham operation groups. The primer sequences for CLIC1, DDIT4, TKT, DDOST, SNCA, MYH7, ANP, and BNP are listed in Supplementary Table 4. The differences in mRNA levels between the two groups were evaluated by using Student’s t-tests. A P-value < 0.05 was considered significant.
Data Availability Statement
Publicly available datasets were analyzed in this study. This data can be found here: GEO (Gene Expression Omnibus) database, Accession number GSE36961: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36961.
Ethics Statement
The animal study was reviewed and approved by the Animal Ethics and Experimentation Committee of Nanchang University.
Author Contributions
KH and ZJ were responsible for the entire project and revised the draft of the manuscript. ZT, LW, and ZJ collected the data, performed the analyses, and drafted the first version of the manuscript. All authors took part in the interpretation of the results and preparation of the final version of the manuscript.
Funding
This work was supported by the National Natural Science Foundation of China (NSFC, 32060154).
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.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.670787/full#supplementary-material
Supplementary Table 1 | 2927 DEGs between HCM and control samples
Supplementary Table 2 | Annotation results for 141 upregulated DE-NMGenes, 175 DE-NMGenes, and 4 network modules
Supplementary Table 3 | Modules detected from PPI network
Supplementary Table 4 | Real-time PCR Primer sequences.
Footnotes
- ^ https://www.ncbi.nlm.nih.gov/geo/
- ^ http://impi.mrc-mbu.cam.ac.uk/
- ^ https://github.com/ZhenhongJiang/Nuclear-encoded-mitochondrial-genes-in-HCM.git
- ^ https://cytoscape.org/
References
Allard, M. F., Henning, S. L., Wambolt, R. B., Granleese, S. R., English, D. R., and Lopaschuk, G. D. (1997). Glycogen metabolism in the aerobic hypertrophied rat heart. Circulation 96, 676–682. doi: 10.1161/01.cir.96.2.676
Alonso-Montes, C., Rodriguez-Reguero, J., Martin, M., Gomez, J., Coto, E., Naves-Diaz, M., et al. (2017). Rare genetic variants in GATA transcription factors in patients with hypertrophic cardiomyopathy. J. Investig. Med. 65, 926–934. doi: 10.1136/jim-2016-000364
Bader, G. D., and Hogue, C. W. (2003). An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics 4:2. doi: 10.1186/1471-2105-4-2
Barabasi, A. L., and Oltvai, Z. N. (2004). Network biology: understanding the cell’s functional organization. Nat. Rev. Genet. 5, 101–113. doi: 10.1038/nrg1272
Brown, D. I., and Griendling, K. K. (2015). Regulation of signal transduction by reactive oxygen species in the cardiovascular system. Circ. Res. 116, 531–549. doi: 10.1161/CIRCRESAHA.116.303584
Calvo, S. E., Clauser, K. R., and Mootha, V. K. (2016). MitoCarta2.0: an updated inventory of mammalian mitochondrial proteins. Nucleic Acids Res. 44, D1251–D1257. doi: 10.1093/nar/gkv1003
Chen, L., Zhao, L., Samanta, A., Mahmoudi, S. M., Buehler, T., Cantilena, A., et al. (2017). STAT3 balances myocyte hypertrophy vis-a-vis autophagy in response to Angiotensin II by modulating the AMPKalpha/mTOR axis. PLoS One 12:e0179835. doi: 10.1371/journal.pone.0179835
Clough, E., and Barrett, T. (2016). The gene expression omnibus database. Methods Mol. Biol. 1418, 93–110. doi: 10.1007/978-1-4939-3578-9_5
da Huang, W., Sherman, B. T., and Lempicki, R. A. (2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57. doi: 10.1038/nprot.2008.211
Doenst, T., Pytel, G., Schrepper, A., Amorim, P., Farber, G., Shingu, Y., et al. (2010). Decreased rates of substrate oxidation ex vivo predict the onset of heart failure and contractile dysfunction in rats with pressure overload. Cardiovasc. Res. 86, 461–470. doi: 10.1093/cvr/cvp414
Duan, Q., Ni, L., Wang, P., Chen, C., Yang, L., Ma, B., et al. (2016). Deregulation of XBP1 expression contributes to myocardial vascular endothelial growth factor-A expression and angiogenesis during cardiac hypertrophy in vivo. Aging Cell 15, 625–633. doi: 10.1111/acel.12460
El-Hattab, A. W., and Scaglia, F. (2016). Mitochondrial cardiomyopathies. Front. Cardiovasc. Med. 3:25. doi: 10.3389/fcvm.2016.00025
Ference, B. A., Ray, K. K., Catapano, A. L., Ference, T. B., Burgess, S., Neff, D. R., et al. (2019). Mendelian randomization study of ACLY and cardiovascular disease. N. Engl. J. Med. 380, 1033–1042. doi: 10.1056/NEJMoa1806747
Ferramosca, A., and Zara, V. (2013). Biogenesis of mitochondrial carrier proteins: molecular mechanisms of import into mitochondria. Biochim. Biophys. Acta 1833, 494–502. doi: 10.1016/j.bbamcr.2012.11.014
Fontes, M. S., van Veen, T. A., de Bakker, J. M., and van Rijen, H. V. (2012). Functional consequences of abnormal Cx43 expression in the heart. Biochim. Biophys. Acta 1818, 2020–2029. doi: 10.1016/j.bbamem.2011.07.039
Gao, J., Collyer, J., Wang, M., Sun, F., and Xu, F. (2020). Genetic dissection of hypertrophic cardiomyopathy with myocardial RNA-seq. Int. J. Mol. Sci. 21:3040. doi: 10.3390/ijms21093040
Gersh, B. J., Maron, B. J., Bonow, R. O., Dearani, J. A., Fifer, M. A., Link, M. S., et al. (2011). 2011 ACCF/AHA guideline for the diagnosis and treatment of hypertrophic cardiomyopathy: executive summary: a report of the American College of Cardiology Foundation/American Heart Association task force on practice guidelines. Circulation 124, 2761–2796. doi: 10.1161/CIR.0b013e318223e230
Gurtner, A., Fuschi, P., Martelli, F., Manni, I., Artuso, S., Simonte, G., et al. (2010). Transcription factor NF-Y induces apoptosis in cells expressing wild-type p53 through E2F1 upregulation and p53 activation. Cancer Res. 70, 9711–9720. doi: 10.1158/0008-5472.CAN-10-0721
Gururaja Rao, S., Patel, N. J., and Singh, H. (2020). Intracellular chloride channels: novel biomarkers in diseases. Front. Physiol. 11:96. doi: 10.3389/fphys.2020.00096
Hafstad, A. D., Nabeebaccus, A. A., and Shah, A. M. (2013). Novel aspects of ROS signalling in heart failure. Basic Res. Cardiol. 108:359. doi: 10.1007/s00395-013-0359-8
He, L., Kim, T., Long, Q., Liu, J., Wang, P., Zhou, Y., et al. (2012). Carnitine palmitoyltransferase-1b deficiency aggravates pressure overload-induced cardiac hypertrophy caused by lipotoxicity. Circulation 126, 1705–1716. doi: 10.1161/CIRCULATIONAHA.111.075978
Hilfiker-Kleiner, D., Hilfiker, A., Castellazzi, M., Wollert, K. C., Trautwein, C., Schunkert, H., et al. (2006). JunD attenuates phenylephrine-mediated cardiomyocyte hypertrophy by negatively regulating AP-1 transcriptional activity. Cardiovasc. Res. 71, 108–117. doi: 10.1016/j.cardiores.2006.02.032
Hu, X., Shen, G., Lu, X., Ding, G., and Shen, L. (2019). Identification of key proteins and lncRNAs in hypertrophic cardiomyopathy by integrated network analysis. Arch. Med. Sci. 15, 484–497. doi: 10.5114/aoms.2018.75593
Janky, R., Verfaillie, A., Imrichova, H., Van de Sande, B., Standaert, L., Christiaens, V., et al. (2014). iRegulon: from a gene list to a gene regulatory network using large motif and track collections. PLoS Comput. Biol. 10:e1003731. doi: 10.1371/journal.pcbi.1003731
Jensen, M., Muller, C., Choe, C. U., Schwedhelm, E., and Zeller, T. (2020). Analysis of L-arginine:glycine amidinotransferase-, creatine- and homoarginine-dependent gene regulation in the murine heart. Sci. Rep. 10:4821. doi: 10.1038/s41598-020-61638-3
Koren, L., Elhanani, O., Kehat, I., Hai, T., and Aronheim, A. (2013). Adult cardiac expression of the activating transcription factor 3, ATF3, promotes ventricular hypertrophy. PLoS One 8:e68396. doi: 10.1371/journal.pone.0068396
Lachance, D., Dhahri, W., Drolet, M. C., Roussel, E., Gascon, S., Sarrhini, O., et al. (2014). Endurance training or beta-blockade can partially block the energy metabolism remodeling taking place in experimental chronic left ventricle volume overload. BMC Cardiovasc. Disord. 14:190. doi: 10.1186/1471-2261-14-190
Li, J., Wu, Z., Zheng, D., Sun, Y., Wang, S., and Yan, Y. (2019). Bioinformatics analysis of the regulatory lncRNAmiRNAmRNA network and drug prediction in patients with hypertrophic cardiomyopathy. Mol. Med. Rep. 20, 549–558. doi: 10.3892/mmr.2019.10289
Liang, Q., De Windt, L. J., Witt, S. A., Kimball, T. R., Markham, B. E., and Molkentin, J. D. (2001). The transcription factors GATA4 and GATA6 regulate cardiomyocyte hypertrophy in vitro and in vivo. J. Biol. Chem. 276, 30245–30253. doi: 10.1074/jbc.M102174200
Liao, X., Haldar, S. M., Lu, Y., Jeyaraj, D., Paruchuri, K., Nahori, M., et al. (2010). Kruppel-like factor 4 regulates pressure-induced cardiac hypertrophy. J. Mol. Cell Cardiol. 49, 334–338. doi: 10.1016/j.yjmcc.2010.04.008
Lim, D. S., Roberts, R., and Marian, A. J. (2001). Expression profiling of cardiac genes in human hypertrophic cardiomyopathy: insight into the pathogenesis of phenotypes. J. Am. Coll. Cardiol. 38, 1175–1180. doi: 10.1016/s0735-1097(01)01509-1
Lin, C. Y., Lee, T. L., Chiu, Y. Y., Lin, Y. W., Lo, Y. S., Lin, C. T., et al. (2015). Module organization and variance in protein-protein interaction networks. Sci. Rep. 5:9386. doi: 10.1038/srep09386
Liu, X., Ma, Y., Yin, K., Li, W., Chen, W., Zhang, Y., et al. (2019). Long non-coding and coding RNA profiling using strand-specific RNA-seq in human hypertrophic cardiomyopathy. Sci. Data 6:90. doi: 10.1038/s41597-019-0094-6
Marian, A. J., and Braunwald, E. (2017). Hypertrophic cardiomyopathy: genetics, pathogenesis, clinical manifestations, diagnosis, and therapy. Circ. Res. 121, 749–770. doi: 10.1161/CIRCRESAHA.117.311059
Marin-Garcia, J., and Goldenthal, M. J. (2002a). Fatty acid metabolism in cardiac failure: biochemical, genetic and cellular analysis. Cardiovasc. Res. 54, 516–527. doi: 10.1016/s0008-6363(01)00552-1
Marin-Garcia, J., and Goldenthal, M. J. (2002b). Understanding the impact of mitochondrial defects in cardiovascular disease: a review. J. Card. Fail. 8, 347–361. doi: 10.1054/jcaf.2002.127774
Marino, G., Pietrocola, F., Kong, Y., Eisenberg, T., Hill, J. A., Madeo, F., et al. (2014). Dimethyl alpha-ketoglutarate inhibits maladaptive autophagy in pressure overload-induced cardiomyopathy. Autophagy 10, 930–932. doi: 10.4161/auto.28235
Mito, S., Ozono, R., Oshima, T., Yano, Y., Watari, Y., Yamamoto, Y., et al. (2008). Myocardial protection against pressure overload in mice lacking Bach1, a transcriptional repressor of heme oxygenase-1. Hypertension 51, 1570–1577. doi: 10.1161/HYPERTENSIONAHA.107.102566
Mori, J., Basu, R., McLean, B. A., Das, S. K., Zhang, L., Patel, V. B., et al. (2012). Agonist-induced hypertrophy and diastolic dysfunction are associated with selective reduction in glucose oxidation: a metabolic contribution to heart failure with normal ejection fraction. Circ. Heart Fail. 5, 493–503. doi: 10.1161/CIRCHEARTFAILURE.112.966705
Odievre, M. H., Chretien, D., Munnich, A., Robinson, B. H., Dumoulin, R., Masmoudi, S., et al. (2005). A novel mutation in the dihydrolipoamide dehydrogenase E3 subunit gene (DLD) resulting in an atypical form of alpha-ketoglutarate dehydrogenase deficiency. Hum. Mutat. 25, 323–324. doi: 10.1002/humu.9319
Olson, A. K., Ledee, D., Iwamoto, K., Kajimoto, M., O’Kelly Priddy, C., Isern, N., et al. (2013). C-Myc induced compensated cardiac hypertrophy increases free fatty acid utilization for the citric acid cycle. J. Mol. Cell Cardiol. 55, 156–164. doi: 10.1016/j.yjmcc.2012.07.005
Pereira, A. H., Clemente, C. F., Cardoso, A. C., Theizen, T. H., Rocco, S. A., Judice, C. C., et al. (2009). MEF2C silencing attenuates load-induced left ventricular hypertrophy by modulating mTOR/S6K pathway in mice. PLoS One 4:e8472. doi: 10.1371/journal.pone.0008472
Pfanner, N., Warscheid, B., and Wiedemann, N. (2019). Mitochondrial proteins: from biogenesis to functional networks. Nat. Rev. Mol. Cell Biol. 20, 267–284. doi: 10.1038/s41580-018-0092-0
Pikkarainen, S., Tokola, H., Kerkela, R., and Ruskoaho, H. (2004). GATA transcription factors in the developing and adult heart. Cardiovasc. Res. 63, 196–207. doi: 10.1016/j.cardiores.2004.03.025
Ponnalagu, D., Gururaja Rao, S., Farber, J., Xin, W., Hussain, A. T., Shah, K., et al. (2016). Molecular identity of cardiac mitochondrial chloride intracellular channel proteins. Mitochondrion 27, 6–14. doi: 10.1016/j.mito.2016.01.001
Reibel, D. K., O’Rourke, B., Foster, K. A., Hutchinson, H., Uboh, C. E., and Kent, R. L. (1986). Altered phospholipid metabolism in pressure-overload hypertrophied hearts. Am. J. Physiol. 250, H1–H6. doi: 10.1152/ajpheart.1986.250
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43:e47. doi: 10.1093/nar/gkv007
Sacchetto, C., Sequeira, V., Bertero, E., Dudek, J., Maack, C., and Calore, M. (2019). Metabolic alterations in inherited cardiomyopathies. J. Clin. Med. 8:2195. doi: 10.3390/jcm8122195
Saoura, M., Powell, C. A., Kopajtich, R., Alahmad, A., Al-Balool, H. H., Albash, B., et al. (2019). Mutations in ELAC2 associated with hypertrophic cardiomyopathy impair mitochondrial tRNA 3′-end processing. Hum. Mutat. 40, 1731–1748. doi: 10.1002/humu.23777
Sciarretta, S., Volpe, M., and Sadoshima, J. (2014). Mammalian target of rapamycin signaling in cardiac physiology and disease. Circ. Res. 114, 549–564. doi: 10.1161/CIRCRESAHA.114.302022
Shany, E., Saada, A., Landau, D., Shaag, A., Hershkovitz, E., and Elpeleg, O. N. (1999). Lipoamide dehydrogenase deficiency due to a novel mutation in the interface domain. Biochem. Biophys. Res. Commun. 262, 163–166. doi: 10.1006/bbrc.1999.1133
Smith, A. C., and Robinson, A. J. (2019). MitoMiner v4.0: an updated database of mitochondrial localization evidence, phenotypes and diseases. Nucleic Acids Res. 47, D1225–D1228. doi: 10.1093/nar/gky1072
Smits, P., Saada, A., Wortmann, S. B., Heister, A. J., Brink, M., Pfundt, R., et al. (2011). Mutation in mitochondrial ribosomal protein MRPS22 leads to Cornelia de Lange-like phenotype, brain abnormalities and hypertrophic cardiomyopathy. Eur. J. Hum. Genet. 19, 394–399. doi: 10.1038/ejhg.2010.214
Souders, C. A., Borg, T. K., Banerjee, I., and Baudino, T. A. (2012). Pressure overload induces early morphological changes in the heart. Am. J. Pathol. 181, 1226–1235. doi: 10.1016/j.ajpath.2012.06.015
Sucharov, C. C., Dockstader, K., and McKinsey, T. A. (2008). YY1 protects cardiac myocytes from pathologic hypertrophy by interacting with HDAC5. Mol. Biol. Cell 19, 4141–4153. doi: 10.1091/mbc.E07-12-1217
Szklarczyk, D., Gable, A. L., Lyon, D., Junge, A., Wyder, S., Huerta-Cepas, J., et al. (2019). STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 47, D607–D613. doi: 10.1093/nar/gky1131
Tan, C. Y., Wong, J. X., Chan, P. S., Tan, H., Liao, D., Chen, W., et al. (2019). Yin Yang 1 suppresses dilated cardiomyopathy and cardiac fibrosis through regulation of Bmp7 and Ctgf. Circ. Res. 125, 834–846. doi: 10.1161/CIRCRESAHA.119.314794
Tian, R. (2003). Transcriptional regulation of energy substrate metabolism in normal and hypertrophied heart. Curr. Hypertens. Rep. 5, 454–458. doi: 10.1007/s11906-003-0052-7
UniProt, C. (2019). UniProt: a worldwide hub of protein knowledge. Nucleic Acids Res. 47, D506–D515. doi: 10.1093/nar/gky1049
Wang, X., Deng, Y., Zhang, G., Li, C., Ding, G., May, H. I., et al. (2019). Spliced X-box binding protein 1 stimulates adaptive growth through activation of mTOR. Circulation 140, 566–579. doi: 10.1161/CIRCULATIONAHA.118.038924
Weissler-Snir, A., Allan, K., Cunningham, K., Connelly, K. A., Lee, D. S., Spears, D. A., et al. (2019). Hypertrophic cardiomyopathy-related sudden cardiac death in young people in Ontario. Circulation 140, 1706–1716. doi: 10.1161/CIRCULATIONAHA.119.040271
Yang, K., Long, Q., Saja, K., Huang, F., Pogwizd, S. M., Zhou, L., et al. (2017). Knockout of the ATPase inhibitory factor 1 protects the heart from pressure overload-induced cardiac hypertrophy. Sci. Rep. 7:10501. doi: 10.1038/s41598-017-11251-8
Yang, W., Li, Y., He, F., and Wu, H. (2015). Microarray profiling of long non-coding RNA (lncRNA) associated with hypertrophic cardiomyopathy. BMC Cardiovasc. Disord. 15:62. doi: 10.1186/s12872-015-0056-7
Yu, P., Hu, L., Xie, J., Chen, S., Huang, L., Xu, Z., et al. (2018). O-GlcNAcylation of cardiac Nav1.5 contributes to the development of arrhythmias in diabetic hearts. Int. J. Cardiol. 260, 74–81. doi: 10.1016/j.ijcard.2018.02.099
Zang, R., Tan, Q., Zeng, F., Wang, D., Yu, S., and Wang, Q. (2020). JMJD1A represses the development of cardiomyocyte hypertrophy by regulating the expression of catalase. Biomed. Res. Int. 2020:5081323. doi: 10.1155/2020/5081323
Zhang, N., Zhang, Y., Qian, H., Wu, S., Cao, L., and Sun, Y. (2020). Selective targeting of ubiquitination and degradation of PARP1 by E3 ubiquitin ligase WWP2 regulates isoproterenol-induced cardiac remodeling. Cell Death Differ. 27, 2605–2619. doi: 10.1038/s41418-020-0523-2
Keywords: hypertrophic cardiomyopathy, microarrays, bioinformatics analysis, nuclear-encoded mitochondrial genes, transcription factors
Citation: Tan Z, Wu L, Fang Y, Chen P, Wan R, Shen Y, Hu J, Jiang Z and Hong K (2021) Systemic Bioinformatic Analyses of Nuclear-Encoded Mitochondrial Genes in Hypertrophic Cardiomyopathy. Front. Genet. 12:670787. doi: 10.3389/fgene.2021.670787
Received: 22 February 2021; Accepted: 08 April 2021;
Published: 12 May 2021.
Edited by:
Chunjie Jiang, University of Pennsylvania, United StatesReviewed by:
Renxiang Yan, Fuzhou University, ChinaZhou Zhou, Fuwai Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College, China
Copyright © 2021 Tan, Wu, Fang, Chen, Wan, Shen, Hu, Jiang and Hong. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zhenhong Jiang, anh6aGppYW5nQDE2My5jb20=; Kui Hong, aG9uZ2t1aTg4QDE2My5jb20=
†These authors have contributed equally to this work