Differentially Expressed Genes in Cardiomyocytes of the First Camelized Mouse Model, Nrap^c.255ins78 Mouse

1. Introduction

Camels are unique mammals that have adapted remarkably to extreme environmental conditions, thriving in harsh habitats like deserts. Numerous studies have highlighted the physiological, behavioral, and genetic mechanisms camels have developed to survive under such conditions [1,2,3]. However, our understanding of the genetic variations that drive environmental adaptation in camels remains incomplete.

Our previous research identified that exon 4 of the NRAP gene in camels plays a significant role in cold resistance. Specifically, the Nrap^c.255ins78 variant in mice exhibits unique expression in the heart, where it enhances inflammatory cytokine production under cold stress, thereby contributing to temperature resilience [4]. Cytokine expression is known to change in mammals exposed to cold [5,6,7]. For instance, cold exposure has been reported to increase interleukin (IL)-1β and IL-6 levels in cold-resistant humans [8], while in mice, cold stress elevates tumor necrosis factor (TNF)-α and IL-6 levels [9].

These findings suggest a greater upregulation of inflammatory cytokines in Nrap^c.255ins78 mice compared to wild-type mice. To investigate these cytokine mechanisms, this study analyzes changes in the heart transcriptome of Nrap^c.255ins78 mice under non-stress conditions using RNA sequencing (RNA-seq). RNA-seq is a widely used tool for analyzing transcriptomic changes in organs influenced by internal and external environments [10,11,12,13]. This approach, comparing the differentially expressed gene set under normal conditions, may provide insights into the cytokine expression pathways observed in the heart under cold stress. Ultimately, this study could further contribute to understanding temperature adaptation mechanisms in camels and environmental stress adaptation in various mammals.

2. Materials and Methods

2.1. Ethics

All animal experiments including standards of euthanasia for this study were performed according to the Korean Food and Drug Administration (KFDA) guidelines and the relevant legal guidelines and conducted in the GEM division of Macrogen Inc. (Seoul, Republic of Korea). The protocols were reviewed and approved by the IACUC (MS-2022-01).

2.2. RNA-seq

To investigate gender-based differences in gene expression, twelve-week-old Nrapc.255ins78 homozygous mice (2 males and 2 females) were used, along with twelve-week-old wild-type mice (2 males and 2 females) as controls. The excised hearts were photographed to confirm phenotypic characteristics, and total RNA was extracted using Trizol reagent (Sigma-Aldrich, St. Louis and Burlington, MA, USA) following the manufacturer’s instructions. Detailed RNA sequencing procedures were previously reported [4].

2.3. Differentially Expressed Gene (DEG) Analyses

The quality of raw read data for each sample was assessed using FastQC software v0.11.7. Adaptor trimming was conducted with Trimmomatic v0.38, based on quality results. The trimmed reads were then aligned to the reference genome (GRCm39) from the Ensembl genome browser using HISAT2 v2.1.0. Raw counts for each library were calculated based on the exons in Mus musculus GTF v110 (Ensembl) using the featureCounts function of the Subread package v1.6.3. Differentially expressed gene (DEG) analysis was conducted using edgeR v3.26.5, with raw counts normalized via the TMM (Trimmed Mean of M-values) method. DEGs were identified in the hearts of Nrap^c.255ins78 mice compared to wild-type mice, using a false discovery rate (FDR) of < 0.05 and an absolute log2 fold-change (FC) threshold of ≥ 1. A multidimensional scaling (MDS) plot analysis was also performed to illustrate the sample clustering.

2.4. Gene Ontology (GO) Functional Enrichment Analysis

The identified DEGs were annotated with Gene Ontology (GO) terms using DAVID (Database for Annotation, Visualization, and Integrated Discovery) v2024q1 [14]. GO annotations were conducted across the categories of Biological Processes (BP), Cellular Components (CC), and Molecular Functions (MF), using thresholds of P-value < 0.05 and counts ≥ 2. Enriched GO terms were grouped with similar terms and visualized in a bubble plot, displaying -log10 P-value and fold enrichment.

2.5. Interaction Network Construction

A network based on the DEGs was constructed using the Search Tool for the Retrieval of Interacting Genes (STRING) in Cytoscape v3.10.2. The interaction score was set to 0.4, representing a medium confidence level.

2.6. RT-qPCR Validation

Reverse transcription-quantitative polymerase chain reaction (RT-qPCR) was performed using the SYBR Green PCR mix (Qiagen, Hilden, Germany) to validate the RNA-seq results. Β-actin was used as a control for normalization. Relative quantification of mRNA expression was calculated using the 2-ΔΔCt method, and results were presented as the average relative fold change. The primer sequences and amplification temperatures are provided in Table S1.

2.7. Statistics

All statistical values are presented as mean ± SEM. RNA-seq data analysis and visualization were conducted using R statistical software v4.3.3. The qPCR experiments were independently repeated three times. Statistical significance was evaluated using a paired t-test, with a P-value of less than 0.05 considered statistically significant.

3. Results

3.1. Data Processing and Transcriptomes

The average overall mapping rate was 99.15%, and the average unique mapping rate was 72.60% (Table S2). Multiclustering was observed in the MDS analysis, with the transcriptomes of each sample showing differentiation by sex and genotype (Figure 1a). DEG analysis was conducted by comparing expression levels, which were visualized in a volcano plot (FDR < 0.05, absolute log2FC ≥ 1) (Figure 1b). In total, 25 DEGs were identified in heart tissue, with 7 upregulated and 18 downregulated (Table S3).

3.2. Functional Annotations

Functional enrichment analysis was conducted based on Gene Ontology (GO) terms in heart tissue and visualized with a bubble plot (Figure 1c). The most significantly enriched Biological Processes (BP) included cell division (GO:0051301; P-value = 8.78E-10; FDR = 1.21E-07) and cell cycle (GO:0007049; P-value = 1.30E-9; FDR = 1.21E-07). A comprehensive list and detailed values for each GO term are provided in Table S4.

3.3. Expression Pattern and Validation

We visualized the expression patterns of the 25 DEGs across all samples in a heatmap (Figure 2a) and performed a co-expression network analysis to examine interactions among genes identified in enriched GO terms (Figure 2b). The network analysis revealed interactions based on the GO terms “cell division” and “cell cycle.” Within this network, the Kif11 gene displayed a text-mining association with Cib3, which was identified as an upregulated gene

To confirm the reliability of the RNA-seq data, we selected 3 upregulated genes (Cib3, Aldob, Rtn4r), 4 downregulated genes (Ccnb1, Kif11, Aspm, Ncapg), and 2 non-significant genes (Nrap, Spkh2) for validation. The expression levels of all genes closely matched the transcriptome analysis results (Figure 2c).

4. Discussion

In mammals, the heart is a critical organ for regulating body temperature and circulating blood to accommodate sudden physiological changes [15,16,17]. In this study, we grouped the heart transcriptomes of Nrap^c.255ins78 and wild-type mice by genotype and sex. Transcriptomic differences were observed in each group through MDS visualization (Figure 1a), suggesting that Nrap^c.255ins78 exhibits distinct phenotypes depending on gender.

Functional annotation of the 25 DEGs showed significant associations with the “cell cycle” and “cell division” GO terms (Figure 1c). Notably, all 12 genes (Kif11, Anln, Ccnb1, Kntc1, Kif23, Ncapg, Aspm, Mki67, Pimreg, Knl1, Top2a, Nr4a3) associated with cell cycle and division terms were consistently downregulated in Nrap^c.255ins78 mice (Figure 2a, b). Reduced expression of Nr4a3 has been linked to an increased population of cells susceptible to infection [18]. In T. cruzi-infected human fibroblasts, cytokines are secreted as part of the immune response; however, cell cycle downregulation due to infection occurs independently of cytokine signaling [19]. This suggests that inflammatory cytokine regulation and cell cycle function may operate separately, potentially explaining why cell cycle downregulation was not observed under cold exposure conditions.

The DEG network illustrates interactions among 12 genes involved in the cell cycle and division, with Cib3 (calcium and integrin-binding protein 3) showing an association with Kif11 (Figure 2b). Cib3 was the most significantly upregulated gene (Table S3) and belongs to the CIB family, which includes Cib1, Cib2, and Cib4 [20]. CIB proteins are small EF-hand calcium-binding proteins that interact with the cytoplasmic domain of integrin αIIbβ3, playing a role in hemostasis. Among these, Cib3 is expressed at relatively low levels in the heart and muscle tissues [21,22]. Notably, prior research has demonstrated that overexpression of Cib1, which shares strong homology with Cib3, inhibits cell proliferation [23]. These findings align with the results of our study.

Maintaining appropriate calcium levels in the heart is essential for preserving calcium homeostasis in the cytoplasm through ATPase pumps. Although our analysis alone cannot determine the causal relationship between cell proliferation and Cib3 upregulation, the observed inhibition of cell cycle and division genes in Nrap^c.255ins78 mouse hearts suggest a regulatory mechanism to control abnormal calcium flow [24,25]. This altered calcium flow can potentially affect intracellular signaling pathways, and the significant upregulation of Cib3, which is expressed in a calcium-dependent manner, can be interpreted as a response to the calcium flow in cardiomyocytes [26].

This transcriptome analysis did not identify differential expression of the Nrap (nebulin-related anchoring protein) gene, but qPCR results indicated a slight increase in its expression (Figure 2c). The Nrap is a Z-disc protein with nebulin-like super repetitive sequences, playing a crucial role in myofibril formation [27,28]. Nebulin can bind to calmodulin (CaM) and regulate calcium release in skeletal muscle [29]. The predictions of the domains and protein structure of Nrap^c.255ins78 suggest that there be structural and functional alteration in the protein [4]. Although further in-depth analysis of the protein structure is necessary, this change in the domain structure of the Nrap^c.255ins78 protein may lead to altered interactions with proteins like CaM, potentially affecting calcium levels in cardiomyocytes [30].

Despite the downregulation of genes involved in cell cycle and division, the hearts of Nrap^c.255ins78 mice did not significantly differ in size from those of wild-type mice (Figure S1). In our findings, downregulated genes such as Angptl4, Ccnb1, Mki67, and Top2a are known to be influenced by YAP/TAZ signaling [31,32,33,34]. Consequently, YAP/TAZ inactivation could lead to alterations in extracellular matrix (ECM) stiffness [35]. Although the precise mechanisms by which ECM stiffness regulates YAP/TAZ remain unclear, compensatory effects in actin cytoskeleton stability might have helped maintain organ size under these conditions [36].

Future investigations should incorporate histological staining, heart rate measurements, and blood composition analyses, including calcium concentration assessments, to clarify physiological impacts. Additionally, integrating transcriptome analyses of other organs, such as the lungs and kidneys, could offer a more comprehensive understanding of genetic influence on temperature stress adaptation.

5. Conclusions

The heart transcriptomes of Nrapc.255ins78 mice under non-stress conditions revealed consistent downregulation of cell proliferation-related genes alongside significant upregulation of the Cib3 gene. These findings suggest a homeostatic response within cardiomyocytes to altered calcium dynamics, ultimately indicating heightened cytokine sensitivity compared to wild-type. This study provides foundational insights into temperature adaptation mechanisms, broadening the understanding of the NRAP gene’s exon 4 function beyond camels.