Whole blood transcriptomic analysis reveals PLSCR4 as a potential marker for vaso-occlusive crises in sickle cell disease

Sickle cell disease, a common genetic blood disorder, results from a point mutation in the β-globin gene affecting the configuration of hemoglobin, predisposing to painful vaso-occlusive crisis (VOC) and multi-organ dysfunctions. There is a huge variation in the phenotypic expressions of SCD and VOC owing to genetic and environmental factors. This study aimed to characterize the whole blood gene expression profile using Microarray technology in Bahraini patients with SCD determining the differentially expressed genes in steady-state (n = 10) and during VOC (n = 10) in comparison to healthy controls (n = 8). Additionally, the study intended to identify potential genetic marker associated with hemolysis. The analysis identified 2073 and 3363 genes that were dysregulated during steady-state and VOC, respectively, compared to healthy controls. Moreover, 1078 genes were differentially expressed during VOC compared to steady state. The PLSCR4 gene was almost 6-fold up-regulated in microarray, 4-fold in polymerase chain reaction, and a mean protein concentration of 0.856 ng/ml was observed in enzyme-linked immunosorbent assay during VOC compared to steady-state (0.238 ng/ml) (p < 0.01). Amongst these genes, PLSCR4 is involved in erythrocyte membrane deformity thus, predisposing to hemolysis, adhesion, and thrombosis. In conclusion, PLSCR4 may serve as a potential biomarker for VOC and future large-scale validation are recommended.

www.nature.com/scientificreports/ This study aimed to compare SCD patients using gene expression analysis by microarray to improve our knowledge of the disease pathophysiology and to detect genetic markers associated with the VOC. Identifying such markers may aid in the development of targeted therapy to treat SCD and prevent VOC.

Results
Characteristic of participants. Twenty sickle cell disease patients and their characteristics are listed in Table 1. Majority of the study participants were male with a mean age of 33 years. The frequency of VOC among steady-state group was 9.7 ± 6. 16 per year in which 34% needed hospital admission, while the VOC frequency was 8.3 ± 5.38 per year in the VOC group from which nearly half of the attacks (44.6%) ended up with hospital admission.
The statistical analysis showed no significant differences between both groups in baseline characteristics (p-value > 0.05). However, there was a statistically significant difference among the two groups (steady-state, VOC) in the level of hyperbilirubinemia in form of direct bilirubin (10.1 ± 3.45, 20.7 ± 10.76 in µmol/L, p = 0.03) and indirect bilirubin (20.7 ± 10.76, 35.02 ± 17.35 in µmol/L, p = 0.04). Determination of the differentially expressed genes. A total of 2073 genes were dysregulated in which 736 genes were up-regulated (p < 0.05) with a fold change of > 2, and 1337 genes were down-regulated (p < 0.05) with a fold change of < − 2 in SCD patients in steady-state compared to healthy controls (Fig. 1a). Whereas, in SCD patients in VOC compared to healthy controls, 3363 genes were differentially regulated including 1080 genes were up-regulated (p < 0.05) with a fold change of > 2 and 2283 genes were down-regulated (p < 0.05) with a fold change of < − 2 (Fig. 1b). In addition, 1078 genes were differentially expressed including 410 up-regulated genes and 668 down-regulated genes in SCD patients in VOC compared to steady-state (p < 0.05) with a fold change of > 2 and < − 2, respectively (Fig. 1c).
Furthermore, at a p-value of < 0.001 and a fold change of > 4 a total of 292 genes were up-regulated in SCD patients in steady-state compared to healthy controls, while in SCD patients in VOC compared to healthy controls 428 genes were up-regulated, Table 2. Additionally, 31 genes were up-regulated in SCD patients in VOC compared to steady-state.
Moreover, the IFIT1B (Interferon Induced Protein with Tetratricopeptide Repeats 1B) gene shows the highest fold change among the differentially regulated genes as it was up-regulated in SCD patients in steady-state compared to healthy controls with a fold change of 70.57 (p = 1.7e−9), as well as in SCD patients in VOC compared to healthy controls at a fold change of 313.47 (p = 3.65e−10) (Fig. 2). www.nature.com/scientificreports/ Potential genetic marker for vaso-occlusive crisis. To identify a potential genetic marker for VOC, the thirty-one up-regulated genes were assessed to define the "True up-regulation" against those genes that were up-regulated at a P-value of < 0.001 with 4-fold change in SCD patients in VOC compared to SCD patients in steady-state and in SCD patients in VOC compared to healthy controls, but not more than 2-folds downregulated in SCD patients in steady-state compared to healthy controls at a P-value of < 0.05. A total of 13 upregulated genes remained, Table 3. Following the data mining, an up-regulated gene PLSCR4 (Phospholipid Scramblase 4) was selected to be assessed by qRT-PCR and ELISA and further studied in association with hemolysis and inflammation. The PLSCR4 gene showed almost 6-folds up-regulation in SCD patients in VOC compared to SCD patient in steadystate (p = 1.63e−6) (Fig. 3). PCR showed statistically significant up-regulation in SCD patients in VOC compared to SCD patient in steadystate with a 4-fold changes (p = 0.00017) (Fig. 4a,c). Moreover, the PLSCR4 gene expression was confirmed by measuring the protein concentration using ELISA. The assay showed a significant increase in PLSCR4 expression in SCD patients in VOC (0.856 ng/ml) compared to SCD patient in steady-state (0.238 ng/ml) (p = 9.072e−6) (Fig. 4b).

Discussion
In SCD, numerous biomarkers were identified to be involved in the disease pathophysiology 14 but their clinical utility was questioned owing to influence of several factors are contributing to inflammation, oxidative stress, adhesion and coagulation 15,16 . Few studies were conducted in detecting the transcriptomic biomarkers of SCD and their results are encouraging; however, these studies were carried out in West African population 17,18 . To the best of our knowledge, this is the first study to explore the use of microarray technology in determining the whole blood gene expression profiles in Arabs with SCD.
Recently, a gene expression meta-analysis combined with a genome-wide association study data analysis identified that major pathways involved in SCD include innate immunity, hemostasis, response to stress, hemopoiesis, heme biosynthesis and apoptosis 19 . Several up-regulated genes identified in the previous transcriptomic studies were also found to be up-regulated in our study, notably the genes related to erythrocyte development and innate immune system. Also, the meta-analysis discovered some of protein coding genes that were not previously studied in relation to SCD including RUNDC3A, TMCC2, OSBP2 and IFI27 20 similar to the present study.
However, characteristically, only in our study population, the IFIT1B gene was the highest up-regulated gene in SCD patients in both states with a higher fold change in patients with VOC. This gene was not previously studied in relation to SCD. According to the GeneCards database, the IFIT1B is a protein coding gene, which has a part in the Toll-like Receptor signaling pathway 21 and belongs to the IFIT family which exhibits an Table 3. Differentially regulated genes in SCD patients in VOC compared to SCD patients in steady-state. www.nature.com/scientificreports/ antiviral activity 22 . Moreover, the GeneHancer revealed a relation between IFIT1B gene and phenotype SNPs in leukocyte count, C-reactive protein measurement and myocardial infarction 23 . Nevertheless, IFIT1B gene is uncharacterized yet and its function and expression need to be identified 22 .

ID Gene symbol Chromo-some p-values
In our study, the PLSCR4 was significantly up-regulated in SCD patients in VOC compared to SCD patient in steady-state and this was further confirmed by qRT-PCR and ELISA. These results suggest that PLSCR4 gene may serve as a potential genetic marker for the development of VOC in SDC patients.
The PLSCR4 was not previously identified in studies investigated gene expression in SCD patients and was not previously assessed for its role at transcriptional level in SCD patients.
PLSCR4 is a protein coding gene located on chromosome 3q24 and belongs to the scramblase family 21 . The gene function facilitates fast bidirectional transbilayer transportation of phospholipids including PS which triggered by Ca2 + binding and independent of ATP resulting on loss of the normal asymmetry of plasma membrane 24,25 . Additionally, it has a crucial role in fibrin clot formation initiation, mast cells activation, detection and removal of damaged cells by the reticuloendothelial system [26][27][28][29][30] .
A study on murine models of SCD found that majority of sickle mice RBCs exhibit phospholipid scramblase that was linked to phosphatidylserine (PS) externalization, in addition, the results reveal a rapid peripheral destruction of RBCs which was strongly correlated with the degree of scramblase activation suggesting that scramblase activity and PS exposure serve as markers of hemolysis 31 . Moreover, PS-exposing cells act as a target for interaction with proteins and enzymes in plasma such as secretory phospholipase A2 which is a potent lipid mediator in inflammation as it hydrolyzes the lipids in PS-exposing RBCs producing lysophospholipids and free www.nature.com/scientificreports/ fatty acids that affects vascular integrity contributing to vascular damage, as well as secretory phospholipase A2 was identified to be involved in the development of acute chest syndrome 32,33 . Likewise, scramblase function has a major role in SCD pathogenesis especially in initiating VOC as it facilitates cell adhesion, aggregation, increasing hemolysis, vascular dysfunction, promoting coagulopathy and inflammation which is the hallmark of microvascular occlusion, yet its mechanism in SCD is still not fully explored and need to be further studied 32 . Furthermore, in a study done by Francis et al. to characterize the function of human phospholipid scramblase 4 (hPLSCR4) which is an isoform of the scramblase family through cloning a recombinant hPLSCR4, it confirms that hPLSCR4 is a Ca2 + -binding protein and identified a point mutation (Asp290 → Ala) resulting in almost 50% reduction in scramblase activity in the presence of Ca2 + 34 .
In addition, Kenneth Dobie invented a patent antisense oligonucleotide (US 10/673,523) targeting nucleic acids encoding phospholipid scramblase 4 and modulating its expression for treatment of diseases such as Scott syndrome 35 . Such a discovery will open the path for treating diseases associated with increased expression of PLSCR4. And based on our study finding this may help as a targeted therapy for SCD patient to reduce the risk of developing VOC, chronic anemia and other disease complications.
In conclusion, our study characterized the transcriptomic signature in Bahraini SCD patients with highlights on the role of hemolysis and inflammation in disease state. Analyzing the transcriptional changes in SCD during two status (steady-state and VOC) resulted in the discovery of new genes to be associated with SCD for the first time and were significantly differentially expressed. Amongst these genes, PLSCR4 is involved in causing RBC membrane deformity thus, predisposing to hemolysis, adhesion, and thrombosis. Further validation in a larger sample size is recommended and its pathway needed to be further studied in relation to the disease pathogenesis as it may serve as a potential genetic biomarkers and aids in the discovery of novel therapeutic target. Finally, the study yielded a transcriptomic database of SCD patients from Arab's ethnicity that may help future studies in further understanding the disease heterogeneity and facilitating the development of personalized medicine and targeted treatment in managing SCD patients.

Methods
Study population. Twenty Bahraini patients with SCD (10 at steady state and 10 with VOC) and eight healthy participants from Salmaniya Medical Complex were recruited in this cross-sectional study during September 2019. Approval from the Research and Ethics Committee of the Arabian Gulf University, and the Research Technical Support Team of the Ministry of Health, Kingdome of Bahrain were obtained. Written informed written consents were obtained from each study participant. All experiments and methods were performed in accordance with relevant guidelines and regulations.
The healthy volunteers were confirmed to have Hemoglobin AA genotype by High Performance Liquid Chromatography (HPLC), while all patients with SCD were confirmed to have HbSS genotype. SCD patients were divided into two groups of ten participants each: SCD patients in steady state defined as participants without any history of VOC that required neither evaluation in an emergency department nor hospital admission 12 weeks prior to the study enrollment; and SCD patients during VOC defined as participants with a history of acute, severe pain at the time of enrollment (self-rated score of ≥ 7 out of 10 on a Numerical Rating Scale (NRS)) 36,37 . All SCD patients were not under treatment with hydroxyurea.
Sample collection. For all subjects, blood samples were collected in two separate tubes. For VOC group, the samples were collected within the first 48 h of the crisis. First, 5 ml of venous blood were collected in serumseparating tube and kept for 30 min in room temperature for clot formation and then centrifuged at 3500 rpm for 15 min. The separated serum was stored at − 80 °C until the analysis. In the second container, 2.5 ml of venous blood were collected in PAXgene® Blood RNA Tube (PreAnalytiX GmbH, Hombrechtikon, Switzerland) for immediate stabilization of intracellular RNA and was kept for minimum 2 h at room temperature to allow for complete lysis of blood cells, and then stored at 4 °C until the analysis that was carried out within 3 days.
RNA extraction and gene expression analysis. RNA was extracted from whole blood samples using PAXgene® Blood RNA kit (PreAnalytiX GmbH, Hombrechtikon, Switzerland) following the manufacturer's instructions. The quantity and purity of RNA samples were determined using the NanoDrop 1000 Spectrophometer (Thermo Fisher Scientific, Inc., Waltham, MA, USA) and the acceptable RNA purity of A260/A280 was 1.8-2.2. The RNA integrity was assessed using 1.2% agarose gel electrophoresis. All RNA samples were stored at − 80 °C until further analysis.
The assessment of gene expression was carried out using Affymetrix ClariomTM S Assays for human and GeneChip™ WT PLUS Hybridization, Wash and Stain Kit (Applied Biosystems™, California, USA) according to the manufacture's protocol. In brief, reverse transcription 100 ng of total RNA of each sample was converted to double-stranded cDNA using the T7 promoter sequence primer. Followed by synthesis and amplification of cRNA by an in vitro transcription of the second-stranded cDNA using T7 RNA polymerase. Then through reverse transcription of cRNA, the second cycle of single-stranded cDNA was synthesized, which contains dUTP. After hydrolyzing the RNA, 5.5 μg of purified single-stranded cDNA was fragmented using uracil-DNA glycosylase and apurinic/apyrimidinic endonuclease 1. Next, by terminal deoxynucleotidyl transferase the fragmented cDNA was labeled with DNA Labeling Reagent which binds to biotin. The fragmented and biotin-labeled singlestranded cDNA samples were hybridized to GeneChip™ WT PLUS for sixteen hours in Affymetrix GeneChip® Hybridization Oven 645. Followed by washing and staining using the Affymetrix GeneChip® Fluidics Station 450 and Affymetrix® GeneChip® Command Console™ (AGCC) software. Finally, the arrays were scanned using Affymetrix GeneChip® Scanner 3000 7G.  The sequences of the primers for PLSCR4 and GAPDH were as follows: PLSCR4 forward primer, 5′-CAT  GGG TCT CTG GCG TTT CT-3′, and PLSCR4 reverse primer 5′-AGT TTG TAC GGT GCC CT-3′; GAPDH forward  primer 5′-TCC CTG AGC TGA ACG GGA AG-3′, and GAPDH reverse primer 5′-GGA GTG GGT GTC GCTGT -3′. The reaction was carried out in 20 µL capillaries and incubated in Light Cycler® 2.0 (Roche). The used Light-Cycler run protocol was as the follows: Denaturation program at 95 °C for 10 min, amplification and quantification program repeated 45 times at 95 °C for 10 s, 60 °C for 30 s and 72 °C for 30 s, and finally a cooling step at 40 °C for 30 s. The accumulation of PCR products during each cycle was determined by observing the rise in fluorescence of DNA-binding SYBR Green. Afterwards, the crossing point of each sample was detected and normalized to the expression of housekeeping gene. Then the fold change of expression was calculated using the 2^-ΔΔCt method.
Enzyme-linked immunosorbent assay. The protein produced by PLSCR4 gene was measured by Enzyme-Linked Immunosorbent Assay (ELISA) (MyBioSource, California, USA) according to the manufacturer's instructions. The assay was performed by assigning duplicated wells for all standards and samples on plates, then 100 µL per well of standard or serum sample were pipetted into the assigned well. After incubation and washing, a 100 µL of Biotin-Conjugate was added, incubated, and washed. After that, a 100 µL of Streptavidin-HRP was added, incubated, and washed. Then, a 100 µL of Substrate Solution was pipetted to each well and incubated at 37 °C for 15-20 min. After getting the desired blue color intensity, the reaction was terminated by adding a 50 μL of Stop Solution to each well. Immediately, the optical density (OD) at 450 nm was measured for each well using Synergy™ HTX Multi-Mode Microplate Reader (BioTek®) (BioTek Instruments, Inc., Winooski, VT, USA) and analyzed by Gen5 2.07.17.

Statistical analysis.
The demographic data were analysed for their differences in SCD patients in steady state and during VOC. The values of continuous data were analysed by Student's two-sided unpaired t-test and presented in mean ± standard deviation (SD). The categorical variables were presented in numbers (percentage) and were analysed using Fisher's exact probability test. p-value of < 0.05 was considered significant.
The Transcriptome Analysis Console (TAC) software version 4.0.0.25 by Thermo Fisher Scientific were used to define the differential expression profile within the different groups, performs statistical analysis and provides a list of differentially expressed genes. Genes with a fold change of > 2 or < − 2 and with a t-test or ANOVA P-value of < 0.05 were considered significantly altered between the conditions of each group.
The Light Cycler Software version 4.1.1.21 was used for the analysis of qRT-PCR results by identifying the crossing points for the target and the reference gene in each sample. Average of crossing points for each target gene was calculated in relative to the housekeeping gene GAPDH in all the groups using relative Mono-Color Relative Quantification assay. Then the fold change was calculated using delta Ct (2^-ΔΔCt) method.
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.