Identification of differential gene expression profile from peripheral blood cells of military pilots with hypertension by RNA sequencing analysis

Elevated blood pressure is an important risk factor for cardiovascular disease and is also an important factor in global mortality. Military pilots are at high risk of cardiovascular disease because they undergo persistent noise, high mental tension, high altitude hypoxia, high acceleration and high calorie diet. Hypertension is the leading cause of cardiovascular disease in military pilots. In this study, we want to identify key genes from peripheral blood cells of military pilots with hypertension. Identification of these genes may help diagnose and control hypertension and extend flight career for military pilots. We use RNA sequencing technology, bioinformatics analysis and Western blotting to identify key genes from peripheral blood cells of military pilots with hypertension. Our study detected 121 up-regulated genes and 623 down-regulated genes in the peripheral blood mononuclear cells (PBMCs) from hypertensive military pilots. We have also identified 8 important genes (NME4, PNPLA7, GGT5, PTGS2, IGF1R, NT5C2, ENTPD1 and PTEN), a number of gene ontology categories and biological pathways that may be associated with military pilot hypertension. Our study may provide effective means for the prevention, diagnosis and treatment of hypertension for military pilot and extend their flight career.


Background
Elevated blood pressure (BP) is an important risk factor of cardiovascular disease (CVD) and is also an important factor in global mortality, leading to about 9.4 million deaths per year [1]. The morbidity and mortality of CVD are associated with degrees of increased BP. For every increase of 20 mmHg in the systolic BP above 115 mmHg, the incidence of CVD risk will double [2]. In general, the threshold of hypertension is set to systolic pressure of 140 mmHg (150 mmHg for older adults) or diastolic pressure of 90 mmHg based on BP measurement in quiescent condition [3]. One-third American adults over 18 years old of age have hypertension, and 54% of old people (55-to 64-year-olds) have high BP. Among those over 75 years old of age in the United States, nearly 80% have hypertension [4].
Military pilots are at high risk of CVD because they are undergo persistent noise, high mental tension, high altitude hypoxia, high acceleration and high calorie diet [5]. Hypertension is the leading cause of CVD in military pilots. Although hypertension itself will not cause sudden disability in flight, but it is a main risk factor for disability in flight career and it is also one of the major reasons to cause the pilot grounded [6]. Wenzel et al. reported that the incidence of hypertension in Brazilian Air Force is about 22% [7]. Grossman et al. found 2.4% of the pilots had moderate or higher blood pressure in a 7.5-year follow-up study of Israeli Air Force pilots [8]. The hypertension prevalence rate was 9.7% in Chinese Air Force pilots, and the grounded rate of pilots was 21.7% among students in the flight academy because of hypertension or increased blood pressure [9]. Essential hypertension is a disease caused by complex, multifactorial and multigenic changes, and it is the result of both gene regulation and environmental impact [10]. In this study, we use RNA sequencing technology, bioinformatics analysis and Western blotting to identify key genes from peripheral blood cells of military pilots with hypertension. Identification of these genes may help diagnose and control hypertension and extend flight career for military pilots.

Study subject
For RNA-Seq, six samples of peripheral blood cell from military pilot (3 hypertensives and 3 normotensives) were collected. For quantitative RT-PCR and Western blotting analysis, another 8 samples of peripheral blood cell from fighter pilot (4 hypertensives and 4 normotensives) were collected. All samples collected with the help of doctors from Lintong Aviation Medical Evaluating and Training Center of Air Force, Xi'an, China. The average systolic blood pressure (SBP) of these hypertension pilots was above 160 mmHg and the average diastolic blood pressure (DBP) was above 100 mmHg. The average SBP of these normotensive pilots was below 135 mmHg and the average DBP was below 85 mmHg.

Peripheral blood mononuclear cells (PBMCs) isolation
Five milliliter whole blood collected was transferred to a 15 ml sterile centrifuge tube, 5 ml phosphate buffer saline (PBS) was added to dilute the whole blood. Five milliliter lymphocyte separation medium was added to another 15 ml sterile centrifuge tube. The diluted whole blood was transferred to the lymphocyte separation medium gently to avoid mixing. Then the blood was centrifuged at 2000 rpm for 20 min, room temperature. The white membrane cells of PBMSs were sucked into another sterile 15 ml centrifuge tube. Add 5 ml PBS and centrifuge at 1000 rpm for 10 min. Aspirate supernatant and add 2 ml ACK buffer to the tube, suspend the cells and keep standing for 5 min. Add 5 ml PBS and centrifuge at 1000 rpm for 10 min. Aspirate the supernatant, and the cells were quickly frozen in liquid nitrogen and used for further analysis.

RNA sequencing
RNA sequencing was entrusted to Novel Bioinformatics Co., Ltd., Shanghai, China. Using high-throughput Life technologies Ion Proton Sequencer, the transcript with poly(A)-containing RNA of Human were analyzed.

Quality control
Fast -QC software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was employed to evaluate the overall quality of sequencing data.

Gene expression calculation
The expression of genes is mainly calculated by RPKM (Reads Per Kb per Million reads) method. The formula is: RPKM ¼ 10 6 C NL=10 3 RPKM is gene expression, C is the unique number of reads on the gene, and N is the unique number of total reads in the reference gene, and L is the base number of the gene coding region. RPKM method can eliminate the effects of gene length and sequencing on the expression of genes, which can be directly used to compare the gene expression differences between different samples.

GO analysis
The difference of gene analysis based on the database from BP, MF, CC GO annotation in the three dimensions and all GOs were obtained. Each GO significance level was obtained by using the Fisher test (P Value) and so gene enrichment significant difference GOs were screened out [11,12].

Pathway analysis
The differential expression genes filter out were annotated in KEGG database (http://www.genome.jp/kegg/) for Pathway annotations, and all the Pathway Terms of different genes involved were got. The Pathway of significance level is obtained by using the Fisher test (P Value), and the significant Pathway Term of differential expression gene enrichment was screened out.

Gene-act-network
The construction of gene interaction is to sort out the regulation of all the genes, and through the construction of signal transduction network, we could easily find the vein of gene signal transduction. Based on KEGG database, we could get gene interactions. So, we constructed gene and adjacency matrix using gene interactions [13][14][15][16][17].

Overview of RNA sequencing data
We used Fast-QC online software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) to assess the quality of sequencing results (Results were not shown). The results showed that the sequencing data of all samples were qualified. Total raw reads of these six samples were about 12-15 million. GC content of these six samples was about 53%. The average reads mapped to human genome sequence were about 1.28 ± 0.78 × 10 7 reads (96% of the total reads) in the six samples (Additional file 1: Table S1) and about 92% of the total reads (1.24 ± 0.74 × 10 7 reads) were mapped to human genome sequence uniquely. MapSplice was employed to map the reads. According to the mapping results, about 1 × 10 7 reads were mapped to the transcript exon, 7 × 10 6 reads were mapped to CDS, 2 × 10 6 reads were mapped to intron, 2 × 10 6 reads were mapped to the UTR regions and the rest reads (about 5% or less) were mapped to intergenic regions, TSS (transcription start site) and TES (transcription end site) (Fig. 1a). We also detected the distribution on chromosomes of these sample mapped reads (Fig. 1b). The results showed that the most reads were aligned to chromosome 1 (about 10% or more) and the least reads were aligned to chromosome Y (less than 0.2%).

Differential gene expression profiles between PBMCs from hypertensives and normotensives military pilots and GO analysis
The six samples were screened for difference gene expression with 3:3. The total number was 26, lower than the minimum standard of data analysis. Two discrete samples (Con4 and Hyp3) were deleted according to results of principal component analysis, and 744 differential genes were acquired. This number was enough for further analysis, so we adopted difference gene expression screened with 2:2.
To characterize the gene expression changes between hypertensives and normotensives military pilots' PBMCs, differentially expressed genes were screened by the following standard: Log 2 FC > 0.585 or Log 2 FC < − 0.585, FDR < 0.05 and P value < 0.05. We found 744 differentially expressed genes between hypertensives and normotensives military pilots' PBMCs. Of these genes, 121 genes were up-regulated in hypertensives military pilots' PBMCs and 623 genes were down-regulated. PCA (principal components analysis) cluster analysis was used to compare differential expression of these two groups. The differential gene expression patterns of these samples were showed by gene thermal map (Fig. 2). Gene ontology analysis was used to seek the functions of these differentially expressed genes. The results showed that there were 337 genes belonged to protein binding and 64 genes belonged to ATP binding for molecular function (MF). There were 91 genes belonged to signal transduction and 70 genes belonged to small molecule metabolic process for biological process (BP). As for cellular component (CC), there were 311 genes belonged to membrane and 248 genes belonged to cytoplasm. Results from GO term analyzing showed that inflammatory response, protein binding and phagolysosome were the most significant for BP, MF and CC respectively (Fig. 3).

Pathways analysis of differential expression genes
The differential expression genes were involved in multiple GO, so we constructed functional relation network with significant GO-Term (p-value< 0.05) to reveal relationship between genes clearly based on hierarchical structure of GO (Additional file 2: Figure S1). Of these pathways, protein phosphorylation, toll-like receptor signaling pathway and cell surface receptor signaling pathway were in the core position (Additional file 2: Figure S1). Then, pathway-analysis was carried out to detect significant and important pathways of these differential expression genes. Influenza A and osteoclast differentiation were the most significant (Fig. 4a). Also, the top 20 of pathway enrichment was displayed in Fig. 4b. P-Value and gene number was indicated as circle size and color. Next, the pathways interaction network was built to analysis deeply. The analysis results showed that the most important pathways were apoptosis, Jak-STAT signaling pathway, toll-like receptor signaling and cytokine-cytokine receptor interaction (Additional file 3: Figure S2). Because these four pathways are located at the center of the all significant pathways and have the most arrowheads around, these four pathways are likely to be most important in the elevated blood pressure of military pilots. This result suggested that differential expression genes related to apoptosis, Jak-STAT signaling pathway, toll-like receptor signaling and cytokine-cytokine receptor interaction may have important role in the occurrence and development of elevated blood pressure of military pilots.

Gene act network of differentially expressed genes
Although we got four important pathways related to elevated blood pressure of military pilots, we know that one gene may be involved in multiple signal transduction pathways at the same time. So next, we built gene act network based on the relationships between the differentially expressed genes including expression, binding, inhibition, activation and compound. This analysis method can form the corresponding regulation relationship between gene and gene and is easier to find important related genes under the intervention measures. By analysis the gene act network, we found that NME4, PNPLA7, GGT5, PTGS2, IGF1R, PLCG2, NT5C2, ENTPD1, PIK3CD and PTEN these ten genes were located at the center of the all significant genes and have the most arrowheads around (Fig. 5). And also, these ten genes were Western blotting showed that NME4 and PNPLA7 these two genes were up-regulated significantly and GGT5, PTGS2, IGF1R, NT5C2, ENTPD1 and PTEN these six genes were down-regulated significantly in PBMCs of military pilots with hypertension (Fig. 6). Although the expression change of PLCG2 and PIK3CD was not significant between hypertensives and normotensives, there was downward trend in PBMCs of hypertensive military pilots (Fig. 6). The results indicated that the expression change of NME4, PNPLA7, GGT5, PTGS2, IGF1R, NT5C2, ENTPD1 and PTEN could be as sign of elevation of blood pressure for military pilots.

Discussion
Military pilots are in a state of persistent noise, high mental tension, high altitude hypoxia, high acceleration and their high calorie diet, so hypertension is a very common disease in this group. To identify the key genes that related to hypertension of military pilots is very necessary for prevention, diagnosis and treatment of this disorder. In this study, we used RNA sequencing of PBMCs to identify differential gene expression profile between the hypertensive and normotensive military pilots. Six samples were sequenced and the difference gene expression was compared with 3:3. Because two samples (Con4 and Hyp3) were discrete, so we took them out and adopted different gene expression screened with 2:2. 121 up-regulated genes and 623 down-regulated genes were identified as different expressed genes between hypertensive and normotensive military pilots. We selected 10 important and significant genes according to results of gene act network analysis. Western blotting was employed to validate the expression change of the 10 genes above. The results showed that expression change of NME4, PNPLA7, GGT5, PTGS2, IGF1R, NT5C2, ENTPD1 and PTEN were consistent with the results of RNA sequencing.
NME4 (also known as NDPK-D, Nm23-H4) belongs to the Nm23 family, which includes 10 isoforms (NME1 to NME10) [18]. The classical function of NME4 as a group I isoform is its NDP kinase activity [19]. PNPLA7 (also known as NTE-R1 or NRE) is a member of the PNPLAs family, and its encoded protein is very conservative in mice, rats and human. PNPLA7 plays an important role in the hydrolysis of triglycerides, energy metabolism, lipid formation and adipocyte differentiation [20]. Up-regulation of PNPLA7 indicates elevated blood lipids, which has some correlation with hypertension.
The main expression of Gamma-glutamyl transferase 5 (GGT5) is on the surface of cell membrane, and the role of GGT5 is to hydrolyze the gamma-glutamyl bond glutathione [21]. There is no phenotypic abnormality in GGT5 knockout mice under normal conditions. But GGT5 gene knockout mice were unable to metabolize LTC4, resulting in diminished potential of neutrophils infiltrating into the peritoneum. The expression and function of GGT5 in human are rarely reported [22]. PTGS2(also known as COX-2) could promote carcinogenesis and  [23,24]. Park et al. reported that high fat diet could reduce PTGS2 expression [25], indicating that expression change of PTGS2 may be caused by noise exposure, high altitude hypoxia and high fat diet. IGFR is a member of the receptor tyrosine family and can form homodimers with insulin receptor (InsR) to identify and bind to the ligand of insulin receptor IGF1 and IGF2 [26]. Heterozygous deficiency of Igf1r reduced postnatal growth and develop age-dependent insulin resistance. Old-aged Igf1r +/− mice had increased adiposity and exhibited increased adipogenesis [27], indicated that reduced expression of IGF1R may have a correlation with hypertension. NT5C2 plays an important role in purine metabolism. Some papers reported that there were some relationships between somatic mutations of NT5C2 and T-acute lymphoblastic leukemias (T-ALL) [28,29]. ENTPD1 (also known as CD39) is a plasma membrane protein and its role is to hydrolyze extracellular ATP and ADP to AMP. Helenius et al. reported that the expression of ENTPD1 was significantly reduced in small arterial endothelial cells in patients with pulmonary arterial hypertension (PAH). The attenuation function of ENTPD1 is closely related to vascular dysfunction, suggested that ENTPD1 may be a novel target for PAH therapy [30]. Their results strongly suggested that down-regulation of ENTPD1 may be a sign of hypertension, consistent with our results. PTEN signaling pathway is one of the most important signaling pathways that regulate  [32]. These results suggested that abnormal expression of PTEN could lead nervous system abnormalities and hematologic disease, maybe have some relationship with hypertension. Some of the above eight genes are associated with a high-fat diet, some are associated with noise and high-altitude hypoxia, some may be associated with high mental tension. We hypothesized that the expression change of these genes may affect the level of insulin-like growth factor and small vessel remodeling to induce hypertension under the long-term impact of flight environment. Therefore, we should pay attention to the expression changes of these eight genes in the prevention and control of military pilots' hypertension. Early prevention and early treatment according to the expression change of these eight genes, may extend flight career of military pilots.
At last, we acknowledge that the lack of non-pilot controls is a limitation of the study. The eight differential genes we identified may not be pilots specific. Therefore, our results may lead to some errors in the prevention, diagnosis and treatment of pilot hypertension. But our study still has good guidance for the prevention, diagnosis, and treatment of hypertension in pilots, because the samples used in RNA sequencing and Western blotting experiments in this study are all from military pilots. We will employ non-pilot controls including hypertensives and normotensives to identify differential genes more accurate in the future.

Availability of data and materials
The datasets used and analyzed during the current study are available from the corresponding author on reasonable request.
Authors' contributions XCZ and ZBY conceived and designed the study. XCZ, SHY, YQY, LZ, BJ and SJ performed the experiments. SHY and XZ collected the blood sample. YQY, LZ and SJ extracted and quantified protein sample. XCZ, SHY, XZ and BJ did Western blotting and data analysis. XCZ, YQY, LZ and SJ analyzed the sequencing data. XCZ interpreted the data and wrote the paper. ZBY revised the manuscript. All authors read and approved the final manuscript.
Ethics approval and consent to participate Ethics approval of this study was obtained from committee of Xijing Hospital and all procedures performed involving human participants were in accordance with the ethical standards of this committee. Informed consent was obtained from all individual participants as verbal form included in the study. All individual participants agreed to choose verbal form instead of written form and this was approved by the committee of Xijing Hospital.

Consent for publication
Not applicable.