Bioinformatics Analysis and High-Throughput Sequencing to Identify Differentially Expressed Genes in Nebulin Gene (NEB) Mutations Mice.

BACKGROUND High-throughput sequencing of the pathological tissue of 59 patients with thyroid cancer was compared with the normal population. It was found that the mutation frequency of the Nebulin gene (NEB) at amino acid 1133 locus of thyroid cancer patients was much higher than that of the normal population, suggesting that NEB mutation may be related to thyroid cancer. Therefore, we constructed the NEB mutant mice for further investigation. MATERIAL AND METHODS The RNA extracted from the thyroid of wild-type and NEB mutant mice was analyzed by high-throughput sequencing, and the differential expression was analyzed by edgeR software. Several differentially expressed genes were selected for quantitative real-time PCR (qRT-PCR) verification, and these genes were analyzed with Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. RESULTS A total of 624 genes were significantly enriched. Analysis of GO function and pathway significant enrichment showed that differentially expressed genes were enriched in thyroid cancer, myocardial contraction, and autoimmune thyroid disease. The qRT-PCR results were consistent with the high-throughput sequencing results. CONCLUSIONS Our data indicate that the expression of some cancer-driving genes and cancer suppressor genes are significantly changed in NEB mutant mice compared to wild-type mice, which suggests that NEB function plays an important role in regulating the expression of cancer-related genes in the thyroid gland.

Background

Nebulin, encoded by NEB, is abundantly expressed in the myocardium and is the basic component of striated sarcomere filaments. It is one of the largest proteins in vertebrates (600–900 kDa). NEB contains 183 exons and covers 249 kb in chromosomal region 2q23(NM_001271208.1). Its translation starts from exon 3 and continues until the last exon 183 [1,2]. The smallest NEB exon is 42 bp (exon 4), the largest is 596 bp (exon 183), and most of the exons are between 93 and 312 bp [2,3]. Comparison of NEB sequences between mice and humans reveals highly conserved regions with the same phosphorylation motif and SH3 domain [4].

NEB mutations are a common cause of nemaline myopathy (NM). In addition, core-rod myopathy and distal myopathy are also related to NEB mutations [1,5-7]. About 50% of nemaline myopathy (NM) cases are associated with mutations in NEB, and NEB is one of the largest and most complex genes associated with neuromuscular diseases [8]. Currently, the largest mutation described in the NEB mutation is a 2.5 kb deletion in exon 55 in the Ashkenazi Jewish population [9]. NEB mutations can cause a variety of different phenotypes, manifested in the distribution and severity of muscle weakness [10]. To study the mechanism of muscle weakness caused by NEB mutations, a mouse model with NEB exon 55 deletion was constructed. This model has important phenotypic characteristics of patients and has severe muscle weakness caused by filament dysfunction [9,11]. In patients with NEB-associated myopathy, there are significant differences in the distribution and severity of their disease weaknesses, which may be related to the type of NEB mutation and its expression of nebulin [12,13].

We performed high-throughput sequencing of the pathological tissues of 59 patients with thyroid cancer. Compared with the normal population, we found that the NEB was mutated at amino acid 1133. The mutation frequency of the normal population is 0.1053, while the somatic mutation frequency of thyroid cancer patients is 0.2373. NEB mutation frequency of tumor samples is relative higher than that of the normal population, suggesting that the mutation site is involved in tumorigenesis. To detect the function of NEB in thyroid cancer, we constructed a NEB mutation mouse model. We performed high-throughput sequencing and qRT-PCR validation to investigate the effects of NEB mutation in the mouse thyroid.

Material and Methods

NEB point mutant mouse construction

The mouse NEB is located on mouse chromosome 2 and has a total of 157 exons. We created a C577BL/6 mouse model with point mutation at mouse NEB. We selected the 1133rd locus on exon 33 as a target, and a homologous directed point mutation was performed to mutate the A at the point to G, thus its encoded amino acid is mutated from Asparagine (Asn) to Aspartic acid (Asp). At the same time, we introduced a silent mutation (CAC to CAT) downstream to prevent the binding and re-cleavage of guide RNA to sequences following homologous directed mutagenesis (Figure 1A). NEB mutant mice were provided by Cyagen Biosciences (GenBank accession number: NM_ 010889.1; Ensembl: ENSMUSG00000026950). The mice were reared at the Animal Experiment Center of the Beijing Institute of Biotechnology. Animal experiments were approved by the Animal Care and Use Committee of the Beijing Institute of Biotechnology.

Figure 1

Construction of NEB point mutant mouse model. (A) Mutation site of NEB mutant mouse model. (B) Frequency of NEB mutations in exon sequencing of thyroid cancer patients and normal population. (C) The differential gene heatmap. Rows represent different genes and columns represent different samples; red indicates upregulation of gene expression and green indicates downregulation.

RNA extraction and sequencing

Two-month-old mice were sacrificed and total RNA of the thyroid glands was extracted using an RNA extraction kit (RNeasy Mini Kit 250, QIAGEN, Valencia, CA) and DNA was digested (FSQ-301, TOYOBO, Osaka, Japan). To ensure the accuracy of data interpretation and analysis, 3 groups of biological replicates were established in mutant mice and wild-type mice. The extraction was carried out in strict accordance with the standard operating procedure manual provided by the kit manufacturer and the total sample was extracted.

Total RNA was assayed for quality using an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA) and total RNA was quantified using a Qubit® 3.0 Fluorometer (Invitrogen, CA) and a NanoDrop One spectrophotometer (Thermo Fisher Scientific, Waltham, MA).

The reagents were prepared for sequencing as shown in the Illumina NovaSeq 6000 User Guide manual and the flow cell of the cluster was carried. Double-ended (PE) sequencing was performed using the paired-end program. The sequencing process was completely controlled by the data collection software provided by Illumina (Illumina, San Diego, CA) and the sequencing result data were analyzed in real time. The sequencing was completed by Shanghai Whale Boat Gene Technology Co.

Screening of differentially expressed genes

We use the concept of FPKM (Fragments Per Kilobase Million, or Fragments Per Kilobase of exon model per Million mapped fragments) to characterize the expression levels of different genes. In the application, StringTie [14,15] software was first used to compare the fragments in each gene segment; this was normalized using the trimmed mean of M values algorithm, and the FPKM of each gene was calculated. The obtained FPKM value was regarded as the expression level of each gene. For each sample, FPKM was used to calculate the expression level of each gene, and this value was used for comparisons between different samples.

Based on the comparison results, the differential expression of genes in NEB mutant mice and wild-type mice was analyzed using edgeR software. The differential expression multiple was calculated according to the FPKM value. The Log2 fold Change (log2FC) was calculated. The fold change (FC) was greater than or less than 2 times the differential gene.

Functional analysis of differentially expressed genes

To clarify the biological functions and involved signaling pathways of genes in vivo and in cells, we annotated each gene based on the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. GO enrichment analysis was performed on differential genes, and a hypergeometric test was used to discover the significant enrichment of GO entries in differentially expressed genes to determine the main biological functions.

Validation of differentially expressed genes

The mice were sacrificed with carbon dioxide and the RNA of the thyroid glands was extracted and purified using an RNA extraction kit (RNeasy Mini Kit 250, QIAGEN, Valencia, CA) and cDNA was inverted (FSQ-301, TOYOBO, Osaka, Japan). The differentially expressed genes were verified by qPCR using a Real-Time Quantitative Analyzer (QuantStudio 6 Flex Real-Time PCR System) and the results were analyzed by 2−ΔΔCT. Primers used are shown in Table 1. The data were analyzed by one-way ANOVA using GraphPad Prism5 software. The data are expressed as mean±SEM. Comparisons of 2 sets of data were performed using the t test, and P<0.05 was considered statistically significant.

Table 1

Primers used for the quantitative real-time polymerase chain reaction (qRT-PCR).

Primer	5′-3′ sequence
Actin-f	GAGACCTTCAACACCCCAGC
Actin-r	ATGTCACGCACGATTTCCC
NCOA4-f	ACCAACCCACAGGACTGGCT
NCOA4-r	TTGCCCCAGGCATCGCTGAA
ETV5-f	ACAAACATTTGCGGTCCCCCG
ETV5-r	ATGGGCTCTGACATCTGCCGGT
GPX3-f	ATGTTCGACCAGGTGGGGGCTT
GPX3-r	ATGGGTTCCCAAAAGAGGCGGC
TXNIP-f	TCAGGGACTTGCGCATCGTG
TXNIP-r	ATGCTGGATGTCCGGCTGCT
MYC-f	GCCAGCCCTGAGCCCCTAGT
MYC-r	GCGGAGGTTTGCTGTGGCCT
Gadd45b-f	CCGCTGTGGAGTGTGACTGCAT
Gadd45b-r	TCATCAGTTTGGCCGCCTCGT
CACNA1S-f	CCGCGAGTGGAAGAAGTACGAGT
CACNA1S-r	AGGGAAACACTACAAAGTACACCACG
ATP1B2-f	ACTGGCCACACCAGGCTTGA
ATP1B2-r	TATCGCCCTGGACGGCAGACAT
MLF1-f	CCGGATGCTGAGCAGCTTT
MLF1-r	TCGCCATCATCACGTTCTCGACG
UTY-f	TTGCAACCAACCCCAGGATGCC
UTY-r	GCTCTGCGGGTATTGGTAGGCT

Results

Exon sequencing

We performed whole-exome sequencing on the pathological tissue samples of 59 patients with thyroid cancer, and found that there was a mutation at the amino acid 1133, from A to G, causing the amino acid mutation to change from asparagine to aspartic acid. The frequency of NEB mutation in patients with thyroid cancer is 0.2373, which is much higher than the frequency of 0.1053 in the normal Chinese population (Figure 1B).

Sequencing data analysis

Sequencing obtained 36 GB of raw data and 35.3 GB of pre-processed, effective transcriptome data. The knockout mouse sequence data was compared to the wild-type mouse genome (Figure 1C).

Gene function annotation and classification

By GO function classification, the number of all genes in the biological process, the molecular function, and the cellular component were 3460, 716, and 2439, respectively (Figure 2A). Based on sequence homology, annotations were made using the National Center for Biotechnology Information Clusters of Orthologous Groups of proteins database for a total of 53 classifications. Among them, Cell had the largest number of genes (412), followed by Cell Part and Cellular Process, with 411 and 386, respectively. By KEGG functional classification, there were more differentially expressed genes in signal transduction, immune system, and endocrine system (Figure 2B).

Figure 2

Distribution of differentially expressed genes. (A) Differential gene-related GO function distribution map. (B) Functional distribution map of differential gene-related pathways based on KEGG database. (C) Scatter plot display of differentially expressed genes. (D) Differential gene volcano map. The abscissa is the fold change (FC) and the ordinate is −log10 (P value). Red point represents the upregulated differential gene and blue is the downregulated differentially expressed gene. Gene ontology (GO). Kyoto Encyclopedia of Genes and Genomes (KEGG).

GO analysis of differentially expressed genes

The differentially expressed genes were subjected to GO significant enrichment. The calculated P value was corrected by a multiple hypothesis test and Q value ≤0.05 was used as the threshold. The GO term satisfying this condition was defined as significantly enriched in differentially expressed genes. Correlation analysis showed that the specific gene expression levels of wild-type mice and NEB point mutated mice were well correlated (r=0.964), indicating that the RNASeq data was highly reproducible (Figure 2C). A total of 624 GO entry differential genes were significantly enriched, of which 343 genes were upregulated and 281 genes were downregulated (Figure 2D). Table 2 shows the data of 20 genes that were significantly downregulated or upregulated in NEB point mutant mice.

Table 2

Top 20 down/upregulated genes between mutant mice and wild mice.

Down regualted			Up regualted
Gene	Log2FC	P-value	Gene	Log2FC	P-value
UTY	−8.85	4.15E-51	GPR83	8.17	1.72E-21
LRRC30	−6.91	6.99E-07	IFI44L	6.91	3.80E-20
ODF3L2	−5.36	1.27E-04	IGLC2	5.49	3E-04
MYOD1	−5.18	2.25E-04	IGHG2B	5.35	3.21E-04
MSS51	−4.98	1.19E-07	NR1D2	1.92	3.33E-09
KCNJ4	−4.90	4.44E-14	PER3	1.61	2.13E-07
GM12240	−4.57	4.35E-04	MYC	1.57	2.34E-04
CMYA5	−4.10	3.21E-04	NCOA4	1.50	1.76E-07
MLF1	−3.04	3.15E-04	GPX3	1.05	2.61E-04
NEXN	−2.85	2.42E-0.5	ETV5	1.04	0.0387

Differentially expressed genes were mainly enriched in myofibril assembly and response to interferon-alpha in biological process components (Figure 3A). Immunoglobulin receptor binding, antigen binding, T cell receptor binding, and titin binding were significant enrichment terms in the molecular functional component (Figure 3B). Contractile fiber and myofibril were a type of significant enrichment in cellular component components (Figure 3C).

Figure 3

Enrichment results of differentially expressed genes. (A) GO domain: Biological process. (B) GO domain: Molecular functions. (C) GO domain: Cellular components. (D) KEGG enrichment results of differentially expressed genes. The ordinate is the specific path name. The color of the point indicates the significance of the path (Q value) and the size of the point characterizes the number of genes mapped to the pathway. Gene ontology (GO). Kyoto Encyclopedia of Genes and Genomes (KEGG).

KEGG analysis of differentially expressed genes

Pathway significant enrichment analysis was performed based on the differentially expressed KEGG functional annotation results. Five significant enrichment pathways were obtained, involving thyroid cancer, oxytocin signaling pathway, cardiomyocyte adrenergic signaling pathway, mitogen-activated protein kinase (MAPK) signaling pathway, and autoimmune thyroid disease (Figure 3D). Among them, MAPK signaling pathway plays an important role in cellular physiology and pathological processes, which was significantly convergent with the GO function of differentially expressed genes in the immune system, development, and other biological processes.

The thyroid gland is an important endocrine organ in the body. When it is cancerous, it not only changes the cancer-driving and tumor-suppressing genes, but also hormone synthesis. Therefore, we selected relevant differentially expressed genes for real-time PCR (ABI 7500 Real-time PCR instrument, Life Technologies, US), including MYC, TXNIP, and CACNA1S for cancer-driven genes (Figure 4A); Gadd45b and NCOA4 for thyroid cancer(Figure 4B); MLF1, UTY, and ETV5 for transcriptional disorders in cancer (Figure 4C); and GPX3 and ATP1B2 for thyroid hormone synthesis (Figure 4D). The results showed that the expression of Gadd45b, GPX3, MYC, TXNIP, CACNA1S, NCOA4, and ETV5 in mutant mice was decreased and the expression of MLF1, UTY, and ATP1B2 was increased compared to wild mice, which was consistent with the results of transcriptome analysis.

Figure 4

The quantitative real-time polymerase chain reaction (qRT-PCR) validation of the differentially expressed genes using RNA extracted from wild mice and NEB mutant mice. (A) Cancer-driven genes. (B) Genes involved in thyroid cancer. (C) Transcriptional disorders in cancer. (D) Genes involved in thyroid hormone synthesis. P<0.05 was considered statistically significant in all statistical analyses. The t test was used to compare 2 sets of data. *, **, *** in the graph indicate P<0.05, P<0.01, and P<0.001, respectively.

Discussion

It has long been recognized that NEB mutations are linked to a variety of muscle disorders [12,13]. Our study indicated that the frequency of NEB mutations in patients with thyroid cancer was much higher than that in the normal population, suggesting that NEB mutation may be related to thyroid cancer. We constructed a NEB mutant mouse model to investigate whether there was an association between NEB and thyroid.

High-throughput sequencing analysis was used to compare the differentially expressed genes of wild-type mice and NEB mutant mice. The expressions of 624 genes were significantly changed in mutant thyroid. GO and KEGG analysis showed that the upregulated genes are mainly enriched in dilated cardiomyopathy (DCM), hypertrophic cardiomyopathy (HCM), and myocardial contraction. These pathways are mainly related to nebulin in muscle-related functions. The downregulated genes are mainly enriched in thyroid cancer, autoimmune thyroid disease, iron death, and PPAR signaling pathways, which are mainly related to the immune system and endocrine system. For example, Gadd45b can limit tumor growth and participate in DNA repair, cell survival, and aging or death [16]. Research in Gadd45b knockout mice proved that it was important for tumor immunity [17]. The genetically engineered mouse breast cancer model showed that most genes inhibiting cell proliferation, including Gadd45b, were decreased [18]. MYC is a cancer-driven gene and plays a vital role in tumor transformation [19]. It has been reported that MYC regulates the expression of 2 immune checkpoint proteins on the surface of tumor cells, and its downregulation can enhance the anti-tumor immune response [20-22]. Our data showed that the expressions of Gadd45b and MYC in NEB mutant mice were decreased, which may indicate the occurrence and development of tumors.

Ferroptosis genes were also enriched in our data. Ferroptosis is a new form of programmed cell necrosis discovered in recent years. When iron death occurs, intracellular lipid reactive oxygen species (ROS) are generated [23,24]. Glutathione peroxidase (GPXs) can remove ROS. GPX3 is significantly downregulated in differentially expressed genes. GPX3 is an antioxidant enzyme distributed in glandular epithelial cells, such as kidneys and thyroid glands. It can use glutathione as a reducing matrix to remove H2O2, soluble lipid hydroperoxides, phospholipids, and hydroperoxide to protect kidneys, thyroid glands, and other glands from oxidative damage. It has been reported that GPX3 expression is reduced in prostate cancer model mice and that the risk of cancer is increased in GPX3 knockout mice [25]. This suggests that a significant decrease of GPX3 in NEB mutant mice leads to an imbalance of lipid reactive oxygen species in the glands, which then leads to tumorigenesis.

In addition, thioredoxin interacting protein (TXNIP) can increase the production of ROS and oxidative stress induced apoptosis, and has been identified as a tumor suppressor gene [26]. It has been reported that TXNIP has strong growth inhibition, metastasis inhibition, and pro-apoptotic effects [27,28]. Defects in TXNIP were found to cause cancer in Txnip-deficient mice [29]. The significant reduction of TXNIP in NEB mutant mice suggests that our mice have an increased risk of cancer. Recent studies found that TXNIP also plays a key role in type 1 and type 2 diabetes [28]. This suggests that NEB mutations affect not only the immune system, but also the endocrine system.

Our study may only show the tip of the iceberg in NEB related to thyroid function, and functional research is also needed. More experiments and analysis are necessary to further confirm our conclusions.

Conclusions

In this study, NEB mutant model mice were constructed and compared with wild-type mice by high-throughput sequencing analysis. Differentially expressed genes are mainly concentrated in thyroid cancer, autoimmune thyroid disease, iron death, and myocardial contraction. Our qRT-PCR validation is consistent with our sequencing results. Our data suggest that the mouse NEB mutation has an effect on the mouse thyroid gland.