Occurrence of Orientia Tsutsugamushi Genotypes in Areas of Union Territory of Puducherry and Tamil Nadu State, India

Scrub typhus, also called ‘tsutsugamushi disease’, is an acute febrile illness caused by Orientia tsutsugamushi and transmitted by trombiculid mites; predominantly Leptotrombidium deliense, which accounts for about 20% of pyrexias in endemic regions [1,2]. Although scrub typhus was restricted to the so-called called ‘Tsutsugamushi triangle’, which included most of the South East Asian countries, it is recently reported from areas outside this region; most notable among these being Dubai, Africa and South America [3-5]. In India the disease is reported from wide geographic areas such as Jammu and Kashmir, Rajasthan, Himachal Pradesh North Eastern states and Southern parts of the country [6-8].

Introduction
Scrub typhus, also called 'tsutsugamushi disease' , is an acute febrile illness caused by Orientia tsutsugamushi and transmitted by trombiculid mites; predominantly Leptotrombidium deliense, which accounts for about 20% of pyrexias in endemic regions [1,2]. Although scrub typhus was restricted to the so-called called 'Tsutsugamushi triangle' , which included most of the South East Asian countries, it is recently reported from areas outside this region; most notable among these being Dubai, Africa and South America [3][4][5]. In India the disease is reported from wide geographic areas such as Jammu and Kashmir, Rajasthan, Himachal Pradesh North Eastern states and Southern parts of the country [6][7][8].
The most common clinical presentations of scrub typhus are high fever, headache and myalgia. However, it can also present with varied nonspecific symptoms often leading to a clinical dilemma during the course of infection. Scrub typhus can also result in severe complications that may prove fatal. Although O. tsutsugamushi has been the only known species of the genus, recently a new species, Orientia chuto has been reported from Dubai [3,5]. Based on the variation in the variable domain of the immunodominant surface 56 kDa antigen, prototypes and serotypes of O. tsutsugamushi have been identified. Among these, Karp, Gilliam and Kato are considered as protoypes of Orientia, while there are genotypes such as Kuroki, Shimokoshi, Irie, Hirano, Kawasaki, Boryong etc. [9][10][11]. Information on the clinical presentations of these prototypes/ serotypes/genotypes is scanty. A few reports indicate that the clinical presentation and severity of scrub typhus varies with the infecting strain of O. tsutsugamushi [12]. A study from Korea highlighted the differences in clinical manifestations between Boryong and Karp genotypes of O.
tsutsugamushi. Patients infected with the Boryong genotype were found to present with significantly more generalized weakness, conjunctival injection and skin rash than in the group infected with Karp genotype [10]. Studies conducted in mice model also demonstrated differences in the virulence among strains of O tsutsugamushi, and the strains were classified into a high virulence group (Karp, Kato and KN-3 genotypes), an intermediate virulence group (Gilliam genotype) and a low virulence group (Kuroki, Kawasaki and KN-2) [13,14]. In view of the differences in clinical manifestations with the infecting strain, it is important to investigate the distribution of genotypes across endemic areas and their association with clinical manifestations [10,12,15]. Although scrub typhus is reported from several geographic areas of India, there is a paucity of information on the occurrence of common genotypes of O. tsutsugamushi and the associated clinical presentations. Increasing incidence of scrub typhus with varying clinical manifestations has been reported recently in Puducherry, a coastal city in southern peninsular India, and its adjoining areas of Tamil Nadu state [8,16]. We investigated the occurrence of different genotypes of O. tsutsugamushi among

Materials and Methods
Puducherry is located between 11°56' North and 79°50' East, in south-eastern coastal peninsular India. All suspected cases of scrub typhus reporting to a tertiary care centre (Pondicherry Institute of Medical Sciences) in Puducherry, between March 2011 to March 2013 were recruited to the study. The inclusion criteria were fever for a period of 5 days or more for which the cause was not known or fever with a clinical suspicion of Scrub typhus according to WHO guidelines, in which the cases had an eschar, headache, lymphadenopathy or conjunctival injection. The exclusion criteria were cases with fever due to other common infections in this region like malaria, typhoid and dengue. A sample size of 145 was arrived at, assuming a prevalence rate of 43% with 8% margin of error and 95% confidence interval [17]. The approval of Institutional Ethics Committee was obtained prior to the start of the study (No. IEC/RC/12/54). Blood samples were collected from 145 patients, after obtaining written informed consent. Serum was separated and aliquots of sera and blood were stored at -70°C for further use. Other possible causes of febrile illnesses like Malaria, Typhoid, Dengue, Leptospirosis and Filariasis, were ruled out based on appropriate laboratory investigations as follows; Malaria was ruled out based on microscopic examination of stained peripheral blood smear and Parascreen Rapid test for Malaria (Zephyr Biomedicals, India), Leptospirosis and Dengue were ruled out based on Leptospira IgM and IgG card test (SD Bioline, India) and Dengue day 1 test (J. Mitra, India), respectively. Typhoid was ruled by both conventional blood culture and Widal test, which was done using stained Salmonella antigens (Span, India), and lymphatic filariasis by stained peripheral blood smear and BinaxNow Filariasis card test (Alere, USA).

Detection of IgM antibodies
A commercially available kit (Inbios International Inc. Seattle, USA) was used for detecting anti O. tsutsugamushi IgM antibodies. The kit utilizes recombinant 56kDa antigen of O. tsutsugamushi which was coated on wells of microtiter plate and the ELISA was performed as per the instructions of the manufacturer. An O.D value of the sample more than the cut-off (as per kit insert) was considered as positive for IgM antibodies.

Genotyping of samples
DNA was extracted from blood clots or whole blood samples using the Blood Genomic Extraction kit (Sigma-Aldrich, USA) and stored at -20 o C until further use. The DNA was used in nested PCR targeting 56 kDa antigen of O. tsutsugamushi, employing the method described by Furuya et al. [18] and modified by Saisongkorh et al. [19].
Amplicons of 483 bps obtained were purified using PCR clean-up gel extraction kit (Macherey-Nagel, Germany) and the product was subjected to cycle sequencing reaction with either the forward or reverse primer. This product was further sequenced; using ABI automated Genetic analyzer 3130 XL (ABI-Applied Biosystem, CA, USA). The amplicons were sequenced in forward direction, and if the sequence was not readable, sequencing was performed in the reverse direction also.

Analysis of sequences
The fragment of 56 kDa gene located between variable domains II and III was sequenced from all the PCR positive samples. The nucleotide sequences obtained were edited using BioEdit (Version 7.0.0) and analyzed using Basic local alignment search tool (BLAST) in the NCBI website. Phylogenetic analysis was carried out along with reference sequences obtained from GenBank, in MEGA 6 software. After alignment using BioEdit Sequence alignment editor, 110 nucleotides (nt20 to nt130), which aligned in a similar manner were used for construction of phylogenetic tree. The tree was constructed by Maximum Likelihood method using Tamura-Nei model, with bootstrapping (1000 iterations) and gap positions were included in the analysis. The nucleotide sequences obtained have been submitted to Genbank (Accession numbers KT 970944-KT 970979 and KU 000002-KU 000004).

Clinical analysis
The demographic parameters and clinical presentations of 145 patients were documented by physicians. The clinical presentations were grouped into undifferentiated fever, hepatitis, Acute Respiratory Distress Syndrome (ARDS), meningitis, Multiorgan Disorder Syndrome (MODS) and Disseminated Intravascular Coagulation (DIC) [20,21].

Results
Scrub typhus cases reported in the study came from 22 different rural and urban locations in the Union territory of Puducherry and surrounding areas of Tamil Nadu state, India. About 25% of the cases were from coastal hamlets of Puducherry such as Kalapet, Pillaichavady and Keezhputhupet, while the remaining cases were from neighboring inland areas of Tamil Nadu state (Figure 1). O. tsutsugamushi specific IgM antibodies were detected in 135 (93%) cases out of 145 suspected scrub typhus cases. The remaining 10 IgM negative cases were diagnosed as scrub typhus based on the clinical presentations suggestive of the disease and the presence of the characteristic eschar. Nested PCR targeting 56 kDa gene was done on all the 145 samples of which 40 samples were positive. The amplicons were sequenced and sequences obtained were confirmed to be that of O. tsutsugamushi upon BLAST search.
Seasonal variation in the incidence of scrub typhus cases was observed, with a rise of cases in relatively cooler months, peaking in September numbering 25 (17.2%) and another peak in December and January with 23 and 21 cases respectively (15.8 and 14.4%). The period between June and July months recorded lesser number of cases ( Figure  2). The analysis of age profile of 145 suspected cases of scrub typhus showed an adult preponderance with 65.5% of the cases being adults. Among adults with scrub typhus, 49 were male and 46 were female. The youngest cases were three 1 year infants, while the oldest case was an 86 years old male, who succumbed to the infection due to complications.

Clinical manifestations and laboratory parameters of the cases
The clinical manifestations and laboratory parameters of cases are presented in Figure 3 and Table 1. About 84% of the cases had highgrade intermittent fever of less than 2 weeks duration. Majority of the patients presented with nausea and headache which was seen in 50% and 36.5% of the cases, respectively. The other common features seen were hepatomegaly and abdominal pain. Hepatomegaly, eschar and lymphadenopathy were seen in 29%, 21.3% and 24.8% of the cases, respectively. Manifestations such as chills, rigor and abdominal pain were also observed in some cases. The laboratory parameters found to be elevated included SGPT (48.9%), and bilirubin (51.7%). Anemia was observed in 40.6% of the cases. The clinical classification of cases is given in fig. 4. Although 67% of the cases had an acute febrile illness with no other complications, hepatitis, Acute Respiratory distress The clinical manifestations found in majority of the cases such as nausea/vomiting, headache, chills, rigor, and lymphadenopathy are non-specific symptoms that could be found in other case of febrile illnesses. Characteristic eschar, which occurs at the site of mite bite and reported with varying frequencies, was seen only in a minor proportion of the cases. Thirty one cases (21.3%) had an eschar, predominantly in the thigh, infra axillary or inguinal region. Out of these cases with an eschar, 21 were IgM positive.  Gastrointestinal symptoms like abdominal pain, nausea, vomiting and loose stools, which are reported to be common among scrub typhus cases, were also found to occur in cases in the present study.
Lymphadenopathy is another common feature seen in this study. Elevated ESR (44%), anemia (40.6%) and thrombocytopenia (30.3%) were the other significant laboratory findings. Burning micturition was found to occur in 10.3% of the cases and all these cases had negative urine culture, which is not reported so far in scrub typhus cases. The symptoms more common in children when compared to adults were nausea, presence of an eschar, cough and lymphadenopathy.
Elevated levels of bilirubin, SGPT and thrombocytopenia were seen more among adults while anemia was more common in children.
Liver damage as indicated by simultaneous high levels of serum alkaline phosphatase, SGPT, bilirubin and gamma-glutamyltransferase and ARDS was the most common complication seen in patients included in the study. A single case had co-infection with Chikungunya, 2 with Malaria and 11 cases with Dengue. Most of the patients responded to treatment with doxycycline (100 mg bid for 7 days), however three cases, who reported to the hospital in later stages of the disease and were not on any antibiotics expired due to myocarditis, severe liver dysfunction and shock.

Phylogenetic analysis and clinical features of genotypes:
Thirty five nucleotide sequences of 56 kDa gene of O. tsutsugamushi generated in the present study were aligned along with those of O. tsutsugamushi and Rickettsia species retrieved from the NCBI and used for constructing phylogenetic tree ( Figure 5). The tree broadly formed two major branches. One of these major branches formed two major sub-branches, one comprising the Karp prototype and 15 sequences of the present study and another containing 6 sequences of the present study and a sequence reported from Vellore, India (Vellore_JF421231). This major branch also had Kuroki, Boryong, Chuto and Kato genotypes as independent minor branches. The second major branch divided into two sub-branches; one sub-branch with Gilliam prototype (DQ485289), and the second sub-branch with 12 sequences of the present study as one minor branch. In addition to this, one sequence (SS6) formed a distinct minor branch aligning with an isolate from Japan, and another isolate (SS271) fell under another minor branch along with sequences from Taiwan and Thailand. Thus the O. tsutsugamushi encountered Leukopaenia (<3,500/cu.mm) 6.2 Table 1: Laboratory parameters associated with Scrub typhus cases (n=145), the major marker being elevation of liver enzymes. in the study belonged to several different prototypes/genotypes, with 15 isolates belonging to Karp prototype and interestingly 12 isolates forming an independent group from other genotypes.
An analysis of clinical features among various genotypes of the present study showed that each genotype is associated with a different set of clinical manifestations. The isolates that formed an independent group along with Taiwan isolates had wider spectrum of manifestations involving multiple organs and more severe manifestations like ARDS and hypotension, whereas 6 isolates forming a sub-branch with Vellore isolate were characterized by comparatively less severe and narrow spectrum of clinical manifestations ( Figure 5). Strains that closely aligned with Karp prototype were associated with manifestations like hepatitis, meningitis and MODS.

Discussion
Scrub Typhus is being increasingly reported from South India, with a definitive seasonal occurrence [22]. The present study also documented a seasonal trend, beginning from the end of monsoon (September), to the cooler months of December and January. Cases declined during the hotter season. Land cover with secondary or transitional vegetation is an important risk factor in determining the occurrence of the disease [23]. Post-monsoon, the rich growth of vegetation favors the breeding of mites, unlike the dry season, leading to increased man-vector contact and hence scrub typhus cases. Pondicherry Institute of Medical Sciences is situated in the outskirts of Puducherry town and patients who present to the hospital are mostly from Union Territory of Puducherry and bordering areas of Tamil Nadu state. Hence, the cases were found to be predominantly concentrated in the coastal areas of Puducherry such as Kalapet, Pillaichavady and Keezhputhupet and scrub typhus has been reported from coastal regions earlier [24]. More number of cases was also found to occur in adjoining areas of Tamil Nadu such as Panruti, Villupuram, Tiruvannamalai and Thirukoilur which are characterized by the presence of scrub forests and cashew groves [25].
Adult preponderance seen among cases in this study, could be related to outdoor exposure and work habits, although occupation related exposure was seen only in 4 cases, who were farmers. Also, there is encroachment of farm areas for residential purposes. The eschar, which is said to be a pathognomonic sign of scrub typhus and reported with varying frequencies in earlier studies [7,8,14] was seen only in a minor proportion of cases. The areas in the body where these eschars were seen suggest that the mite clings to tight clothing in these areas as observed in other studies [26]. The unique observation of the present study is burning micturition seen in a small proportion of cases, which needs further clinical evaluation. Some clinical symptoms were more common in children than in adults indicating that the clinical course could be different in different age groups. Co-infections that were seen in the study like Dengue, Malaria and Chikungunya can be expected, as the vectors of these diseases also breed in high densities during the post-monsoon season, increasing the possibilities of coinfections. Mortality, which was observed in three cases, could be due to the organism gaining virulence due to extremes of age, or due to underlying disorders like cirrhosis of liver, which was seen in one case who died.
Scrub typhus is emerging as a major public health problem in India with reports of its occurrence from several wide geographic areas. But, the information on the occurrence of genotypes/prototypes is sketchy. In the present study, genotyping based on 56kDa gene, showed that majority of the study isolates clustered into two groups viz., the Karp prototype and an independent group; probably a new variant. A few belonged to minor independent groups; one sequence (SS6) distinctly branched with isolates from Japan, and another sequence (SS271) with a Thailand isolate. Varghese et al. reported the occurrence of two prototypes viz., Kato and Karp from South India, similar to the former group in the present study [27]. Another study, reported by the same group, genotyped samples from 3 scrub typhus endemic geographic regions viz., in southern, northern and northeastern India and found Kato, Karp, Gilliam, Ikeda and also Neimeng-65 genotype like strains in the South and Northeast regions [28]. These findings indicate that all the prototypes and several genotypes are prevalent in different regions of India. However, in the present study Karp prototype and a new variant genotype were found in addition to other genotypes that have been reported from other endemic areas in South Eastern Asia.
It may be noted that clinical manifestations vary with the genotypes of O. tsutsugamushi. Earlier studies have shown that patients infected with Boryong strain had significantly more generalized weakness, eschar, skin rash, conjunctival congestion, high albumin levels, elevated ESR and fibrinogen levels compared to patients infected with the Karp strain. Response to treatment was also significantly slower in the Karp group as compared to Boryong group [10]. In an experiment done in BALB/c mice using O. tsutsugamushi isolates, remarkably different fatality rates were observed between Boryong strains and it was concluded, that fatality rates due to these strains, correlated with differences in both serotypes and virulence genes. In another study in BALB/c mice, Karp strain was found to multiply intraperitoneally within macrophages and kill the animal, whereas Kuroki strain invaded macrophages but were eliminated from the cells, allowing mouse survival [29,30]. The clinical presentations of the genotypic groups identified in the present study varied considerably and the group which formed an independent cluster along with Taiwan and Thailand isolates presented a broad range of clinical symptoms which ranged from mild features like headache, cough, abdominal pain and vomiting to more severe manifestations like hypotension and hepatitis. Two new genotypes have been reported from Himalayan region of India in an earlier study [31,32]. Yang et al. reported 4 prototype/genotype strains (Karp, Kato, Kawasaki and Gilliam) and eleven different Taiwanese genotypes (Taiwan-A, B, C, D, E, G, H, J, N, O and P) based on RFLP analysis of the 56-kDa type-specific antigen (TSA) gene [33]. Taiwan-H was the major genotype in eastern Taiwan and it was designated as the Taiwan Gilliam-variant (TG-v) cluster. Further, there are now evidences to show that recombinant strains of this bacterium are emerging [10,28,29].
One of the limitations of the study was that, only patients presenting to only one tertiary care hospital were included. Similar cases presenting to neighbouring hospitals could not be enrolled in the study due to logistic reasons. Nevertheless the findings of the present study suggest that there is a great genetic diversity and hence strain diversity among O.tsutsugamushi within and between different scrub typhus endemic regions. A wider and complete networking of hospitals to identify cases and the genotypes infecting amongst them, would give a more comprehensive view of the predominant genotypes. Genotypic variations will greatly impact the diagnosis of scrub typhus as the 56 kDa gene of O. tsutsugamushi and the protein coded by it are targets for PCR and immunological assays. These variations could be the reason for variation in the sensitivity of the 56 kDa protein based ELISA reported by different workers [34][35][36][37]. Similarly, development of vaccine for scrub typhus is beset with problems of polymorphism among O. tsutsugamushi strains. Hence, it is worthwhile to study the genetic diversity of O. tsutsugamushi in different regions and countries where the disease is endemic.


s2fieldsofstudy: Biology Chemistry

text:
QSAR Modeling useful in Anti-Cancer Drug Discovery: Prediction of V600EBRAF-Dependent P-ERK using Monte Carlo Method

Quantitative structure−activity relationship (QSAR) modeling is one of the major computer aided modeling employed in medicinal chemistry. It is used for developing relationships between the effects (e.g. activities and properties of interest) of a series of molecules with their structural properties. The aim of this work was to develop QSAR models to predict the inhibition of V600EBRAF-dependent extracellular regulated kinase (ERK) phosphorylation in WM266.4 melanoma cell lines (IC50p ERK) using the CORAL software (http://www.insilico.eu/coral). These models make use of descriptors based on simplified molecular input-line entry system (SMILES), optimized with the Monte Carlo method. The statistical quality of the newly built models was satisfactory on three random splits of data.

Introduction
Mitogen-activated protein kinase (MAPK) modulates extracellular signals to control several cell functions (from proliferation and survival to differentiation and senescence) [1,2] . One of the most studied MAPK pathways is the extracellular signal-regulated kinase (ERK) pathway. ERK is a subgroup of MAPKs that is activated by external factors such as growth factors and mitogens [3] . The cascade of ERK signalling is triggered by the activation of receptor tyrosine kinases and flows through RAS GTPase [4] . RAF kinases are key components of this pathway: once RAS is turned on, it recruits and activates proteins necessary for the propagation of growth factor and other receptor signals, such as RAF. There are three RAF proteins in mammals, ARAF, BRAF, and CRAF, and they can all activate MAP kinase (MEK) just upstream of ERK [5] . The over-expression or mutations of the components of the RAS-RAF-MEK-ERK-MAP kinase pathway has been found to occur in up to 30 % of human cancers [6] . In different cancer types, the mutation of RAS or RAF causes the hyper activation of ERK followed by un-2 most studied inhibitor for this class of kinases. It was approved by FDA in 2005 for the treatment of renal cell carcinoma and in 2007 for the treatment of hepatocellular carcinoma, and is still undergoing multiple clinical trials in other types of cancer [11,12] .
The experimental measurement of the inhibition activity of chemicals is difficult, expensive and time-consuming. QSAR based analysis can be used as a tool to screen or filter anti-cancer drug candidates, before they are subjected to more intensive calculations, such as docking, or to experimental in vitro measurement of activity and finally under in vivo conditions [13] . In this study, a model was built to predict the half maximal inhibitory concentration (IC 50 ) of V600 EBRAF-dependent ERK phosphorylation in WM266.4 cells. This model made use of descriptors based on simplified molecular input-line entry system (SMILES) [14][15][16][17] , calculated with the CORAL software. Then, chemical information obtained from the calculation of these descriptors was used for drawing up hypothesis of molecular design of Sorafenib derivatives.
The aim was to obtain molecules able to inhibit 50% of ERK phosphorylation in WM266.4 melanoma cells at lower concentrations compared to the template molecule.

Method Data
The dataset of 142 chemical structures used in this study was taken from the literature [18,19] . All these compounds shared a common scaffold ( Figure 2). The dataset contained chemical structures represented as SMILES strings, and experimental values relative to the inhibition of the phosphorylation of extracellular signal-regulated kinases (p-ERK). The selected endpoint was IC 50 of V600 EBRAF-dependent [20] ERK phosphorylation in WM266.4 melanoma cells (IC 50 , pERK) [1,14] . IC 50 values were expressed as negative decimal logarithm (pIC 50 ). In this study ACD ⁄ ChemSketch [21] , was used to build up SMILES of each structure.

CORAL and Optimal descriptors
CORAL is freely available standalone application software for building up regression or classification QSAR models based on the Monte Carlo technique [22] . The calculations are repeated several times for various splits into training, calibration and external validation sets. Optimal descriptors for constructing QSAR models are based on graph or SMILES. The molecular graph includes hydrogen suppressed graph (HSG), hydrogen filled graph (HFG), and graph of atomic orbital's (GAO).
The SMILES-based optimal descriptors include a group of local and one of global attributes. Sk, SSk, and SSSk are local descriptors represented by a sequence of atoms and bonds present in the SMILES string [23] .
A more detailed explanation of CORAL descriptors has been provided by [24] . All the SMILES-based local descriptors and BOND, NOSP and HALO global attributes were selected to build the regression models. Data were randomly split three times into training (60%), calibration (20%) and validation sets (21%). Table 1 The training and calibration sets were used for building the QSAR model and the validation set was used as external set for estimating the predictive potential of the developed model [25] . In order to develop a model with a good predictive potential, the preferable parameters of the Monte Carlo optimization, the threshold (T*) and the number of epochs (N*) that give the best statistics for the calibration set should be defined.
The threshold is a tool for classifying codes as both rare (and thus likely less reliable features, probably introducing noise into the model) and common features, which are used by the model and labelled as active. The optimal descriptors are calculated with the correlation weights (CW) only of active features, excluding those related to rare ones.
The Nepoch is the number of cycles (sequence of modifications of correlation weights for all codes involved in model development) for the optimization [26] . For modelling activity of potential p-ERK inhibitors the descriptor of correlation weights (DCW) was calculated as follows: DCW (Threshold, Nepoch) =∑ CW (Sk) + ∑ CW (SSk) + ∑ CW (SSSk) + ∑ CW (BOND) + ∑ CW (HALO) + ∑ CW (NOSP) (1) where CW(SAk) is correlation weight for a molecular or structural feature (SAk) extracted from k-th SMILES; The endpoint pIC 50 is a function of DCW according to the following equation: Endpoint= C0 + C1 x DCW (T, Nepoch) (2) Where C0 and C1 are the intercept and the slope for the training and calibration set, respectively.
The role of molecular features or structural attributes  (4) can be used to classify the active attributes, and the Defect SMILES defines the domain of applicability for SMILES [26] .

Results and discussion
The models for three random splits are the following: Split 01 pIC 50 = 1.4942 (± 0.0237) + 0.0379 (± 0.0001) * DCW (3, 3) (07) Split 02 pIC 50 = -1.0822 (± 0.0222) + 0.0342 (± 0.0001) *DCW (5, 4) (08) Split 03 pIC 50 = -1.8219 (± 0.0236) + 0.05103 (± 0.0002) *DCW (2, 6) (09) Figure 3(a), Figure 3     The values of the parameters of statistical analysis of p-ERK inhibition models built on the three random distributions of data into training, calibration, and validation sets are collected in Table 2. n is the number of compounds in each set; r² is the coefficient of determination; q² is cross validated r 2 ; MAE is mean absolute error; s is the root-mean-square error; F is the Fischer F ratio; T and N are the preferable values for Threshold and Number of epochs, respectively. CR_p^2 indicates if the models developed are obtained by chance or not, based on the Y-randomization test. For an acceptable model, CR_p^2 should be greater than 0.5 [27,28] . r 2 range from 0.78 to 0.92, and q 2 range from 0.75 to 0.92 (for training, calibration and test sets). According to the criteria defined by Tropsha [26][27][28][29] , a model has high predictive power if the following conditions are fulfilled: q² > 0.5 r² > 0.6 (r² -r²0)/r² < 0.1 or (r² -rʹ0²)/r² < 0.1 0.85 ≤ k ≤ 1.15 or 0.85 ≤ kʹ ≤ 1.15 |r0 2 −rʹ0 2 | < 0.3 Where r0² and rʹ0² are squared correlation coefficients for regression through the origin, calculated between predicted versus experimental values and between experimental versus predicted values, k and kʹ are the slopes in the former and later cases respectively; q² is calculated for the training sets, while all other criteria are calculated for the validation sets [26] . Table 3 contains the values for the criteria iii, iv, and v. Y-randomization test for all models showed that these are not chance correlations, since the CR_p^2 is larger than 0.5 ( Table 2). All the models successfully fulfill the criteria proposed by Tropsha for predictive ability. 86 out of 142 dataset's compounds assigned to the training set are used to develop the QSAR models, 28 compounds are included in the calibration, and 28 compounds, not involved in the model development, are kept out for the external validation . Based on r² and q² values, the three models of equations (07), (08) and (09) have similar statistical behaviours. Indicating that the random selection does not affect the results. The predictability of the models was also checked calculating the statistical parameters r m 2 , average r m 2 , and ∆r m 2 on the calibration sets (Table  3), as proposed in the literature [29,30] . According to the literature, a model has predictive potential if r m 2 > 0.5, average r m 2 > 0.5, and ∆r m 2 < 0.2 [31] . The presented results reveal that the predict ability of all models is very good, especially for the model of equation (09), with q² and r m 2 values equal to e 0.92 and 0.88, respectively (Table 3). The numerical values of (r 2 -rʹ0 2 )/r 2 ; (r 2 -r0 2 )/r 2 ; k; kʹ and |r0 2 -rʹ0 2 |.

Split
(r²-r0 2 )/r 2 (r 2 -rʹ0 2 )/r 2 |r0 2 -rʹ0 2 | k kʹ  The list of all SAk, with the correlation weights for the three splits from the Monte Carlo optimization process of the built QSAR model is given in Table S2,  The mechanistic interpretation obtained from the analysis of SAk can be used in the search and computer aided design of novel p-ERK inhibitors candidates with desired pIC 50 values. Figure 4 shows the structures of proposed kinase inhibitors derivatives, the structure Sorafenib was used as a template for molecular design by added the SMILES notation c...S......., BOND10000000 and c...O....... defined as promoters of increase.
The molecule E has a higher value of pIC 50 in comparison to other molecules.   50 values calculated by using model of split 2 for kinase inhibitors derivatives designed by using the results of QSAR model obtained in this study. Supplementary materials section contains technical details of the models for splits examined in this work.

Conclusion
In this work three QSAR models were developed to predict V600 EBRAF-dependent ERK phosphorylation of 142 molecules as kinase inhibitors, using SMILES based optimal descriptors. The models showed acceptable predictive capability on three random split of data. The analysis of the structural features (obtained from SMILES) with their positive and negative correlation weights allowed designing possible pERK inhibitors with an increased activity compared to Surafenib.
The activity modulating effect of the structural features may serve to draw up preliminary hypothesis of novel anticancer drug candidates, and also to provide a better understanding the drugs mechanisms of action. Table S1. Experimental and calculated of pIC50 p-ERK inhibitors for three splits; training (Tr); calibration (C); and validation (V) sets.


extfieldsofstudy: Biology Medicine

text:
Integrated Analysis of DNA Methylation and mRNA Expression Profiles to Identify Key Genes in Severe Oligozoospermia

Severe oligozoospermia (SO) is a complex disorder, whose etiology is the combined effect of genetic factors and epigenetic conditions. In this study, we examined DNA methylation and mRNA expression status in a set of testicular tissues of SO patients (n = 3), and compared methylated data with those derived from obstructive azoospermia (OA) patients (n = 3) with normal spermatogenesis phenotype. We identified 1,960 differentially methylated CpG sites showing significant alterations in SO vs. OA using the Illumina Infinium HumanMethylation450 bead array. By integrating above DNA methylation data and mRNA expression results, we totally identified 72 methylated CpG sites located in 65 genes with anti-correlation between DNA methylation and mRNA expression. Integrated pathways analysis indicates that these genes are involved in response to hormone stimulus, activation of protein kinase activity, and apoptotic process, among others. We also observed some genes with inversely correlated difference is novel in male infertility field, including PTPRN2, EPHX1, SERPINB9, SLIT3, etc. Our results lay a groundwork for further biological study of SO. Moreover, we generated a workflow for integrated analysis of DNA methylation and mRNA expression, which is expandable to other study types.

INTRODUCTION
Fertility is a major worldwide health problem and affects about 15% of couples, 40-50% of which are attributed to male factors (Anderson et al., 2009). Genetic defects contribute to 15-30% cases of male infertility. Recently, great attentions has been paid to the epigenetic studies due to the great achievement on the novel knowledge in the regulation of gene expression. In contrast to the definition of genetic defects that modify the DNA sequence itself, epigenetic defects are heritable alterations in the function of specific genes but not related to modification in the DNA sequence (Stuppia et al., 2015). DNA methylation is one of the most extensively investigated epigenetic mechanisms. DNA methylation is the stable, covalent addition of a methylgroup to 5 ′ -position of cytosine residues, mainly in CpG dinucleotides. In human, nearly 60-80% CpG are methylated in the promoter regions. Methylation of CpG within the promoter regions commonly leads to the silence of transcription process. In mammals, formation of DNA methylation pattern is one of the important epigenetic events that affect development of various tissues, including testis. The abnormal sperm DNA methylation patterns of both imprinted and non-imprinted genes are proved to be associated with infertility or subfertility in oligospermic men (Urdinguio et al., 2015). However, no previous studies of the global DNA methylation aberrations have been reported in severe oligozoospermia (SO).
SO (sperm concentration <5 million/ml) is found in around 5% of couples presenting with infertility. The pathogenesis of this disease remain largely unknown, although deletions of the Deleted in Azoospermia -like (DAZL) gene were demonstrated to be responsible for severe oligozoospermia in some extent (Fernandes et al., 2006). Moreover, due to the reported heredity of genetic defects, such as Y chromosome microdeletions, the importance of careful evaluation of associated genetic factors in male infertility prior to assisted reproduction is evident. Early studies performed in imprinted genes suggested that there was an enhanced risk of congenital imprinting disorders in children conceived through assisted reproductive technologies (ART). Aberrant DNA methylation patterns have also been proved in some non-imprinted genes involved in spermatogenesis impairment, such as MTHFR, CREM and DAZL (Urdinguio et al., 2015). It is well established that alterations of DNA methylation and differences of gene expression are causally related, with hypomethylation generally resulting in gene expression, and hypermethylation leading to gene silencing. The purpose of this study was to identify key gene determined by both regulatory mechanisms from the testicular tissues of infertile patients via integrated analysis, and to elucidate the underlying biological significance.
By using an Illumina Infinium Human Methylation 450K BeadChip DNA methylation array and Agilent SurePrint G3 Human Gene Expression 8 × 60 K platform, we systematically defined DNA methylation and mRNA expression profiles in testicular tissues of SO (n = 3) by comparing obstructive azoospermia (OA) patients (n = 3) with normal spermatogenesis (included as controls). We considered that the list of differentiated genes contributes to tissue specific functions and integrated data improve us better understanding of the mechanisms of SO. In this article, moreover, we propose an approach to common human disease that incorporates DNA methylation and mRNA expression profiles (Figure 1).

Patients
Testicular biopsy specimens were obtained from 3 patients (aged from 23 to 37 years) with SO and from 3 patients (aged from 25 to 28 years) with OA for the microarray analysis. These patients underwent testicular sperm extraction (TESE) for assisted reproduction and/or diagnostic biopsy for histological examination. For the histological evaluation, the specimens were FIGURE 1 | Schematic pipeline depicting the strategy of identification of key genes from DNA methylation and mRNA expression array integrative data. The method proceeds as follows steps: (1) Identify significantly differentially methylated genes and expressed genes in the case and control groups; (2) Retain those genes whose methylation and expression levels are highly anti-correlated; (3) Enrich the Gene Ontology, KEGG pathway, and PPI networks of anti-correlated genes; (4) Identify key genes from a combination of biological information, including interactions and relationships among genes, gene functional annotations and pathway maps. stained with hematoxylin and eosin (H&E) and analyzed by microscopy. SO was defined as a reduced number of sperm in the male ejaculate or less than 5 million sperm per milliliter. OA was defined as in the reference (Okada et al., 2008): (1) motile spermatozoa were obtained from a microsurgical epididymal sperm aspiration (MESA) or (2) a great number of mature spermatozoa were provided by TESE. The ideal controls in the study would be the normal male population with known fertility, however, the difficulties in sampling testicular tissues leads to this strategy impractical. Instead, control samples of urology males who had no history of meiotic impairment or infertility and exhibited normal spermatogenesis upon histological examination. Karyotype analysis and Y-chromosome microdeletion analysis were performed on all patients to confirm a normal karyotype. Additionally, none of the controls experienced adjuvant hormonal treatment prior to orchiectomy. The infertile male patients who visited Xiamen University Affiliated First Hospital had a routine semen examination on the basis of 2010 WHO criteria. The ethics committees of Xiamen University Affiliated First Hospital (Institutional Review Board Number: KYX-2015-001) approved the study design. All subjects provided written informed consent.

RNA Extraction
Immediately after retrieval, the samples were stored at −80 • C until further RNA process. Total RNAs (including miRNAs) from the frozen testicular tissues was extracted with the miRNeasy Micro Kit (Catalog no. 217084, Qiagen, Germany) following the manufacturer's protocol. On-column DNase digestion was performed to eliminate DNA contamination. The RNA concentration, purity, RNA integrity number (RIN) were determined using a NanoDrop ND-1000 spectrophotometer (Peqlab, Erlangen, Germany), an Agilent 2100 Bioanalyzer and an RNA 6000 NanoLabChip Kit (Agilent Technologies), respectively.

Gene Expression Microarray
For microarray hybridization, 100 ng of total RNA were prepared using the Agilent's One-Color Microarray-Based Gene Expression Analysis Low Input Quick Amp Labeling kit (Agilent 5190-2305) according to the manufacturer's protocol, and were then hybridized to the Agilent SurePrint G3 Human gene expression 8 × 60 K microarray following the manufacturer's instruction. In short, an input of 100 ng of entire RNA was used to generate cDNAs, followed by in vitro transcription and incorporation of cyanine-3 into the nascent cRNAs. The cyanine-3-labeled cRNAs were fragmented and hybridized to the array for 17 h at 65 • C in an Agilent hybridization oven. Arrays were scanned with the Agilent G2565CA microarray scanner at 5 µm resolution. The microarray data were processed with Agilent Feature Extraction 10.7.3.1 Software (Agilent Technologies).

Gene Expression Microarray Data Analysis
The raw signal intensity data were normalized with GeneSpring GX software version 12.0 (Agilent Technologies). Specifically, the signal values were log 2 transformed using the quantile method (Bolstad et al., 2003). Unequal t-test (p-values) was applied to identify the differentially expressed genes in SO from OA. The p-values were corrected by false discovery rate (FDR) of Benjamini and Hochberg (q-values). Fold change (FC) values were calculated for each gene as the difference of the mean intensity of the SO samples from OA samples. Genes with the FC value > 2 or < 1/2 and the q-value < 0.05 were considered as differentially expressed. Quality control analysis of microarray gene expression data was performed as previously described (Zahurak et al., 2007).

DNA Isolation and Preparation
Testicular tissues were dissolved in 200 µl lysis buffer of DNA Micro Kit (Catalog no. 56304, Qiagen, Germany), and incubated with proteinase K overnight at 56 • C for two nights. DNA was extracted according to the manufacturer's protocol (QIAamp DNA Micro Kit, Qiagen), and DNA concentration was determined at 260 nm using the NanoDrop ND-1000 spectrophotometer (Nanodrop Technologies Inc., Wilmington, NC, USA). Bisulfite modification of 500 ng DNA of each samples was carried out with the EZ DNA methylation kit (Zymo Research, Orange, CA).

DNA Methylation Microarray
Infinium HumanMethylation450 BeadChips (Illumina, Inc.) were used to analyze DNA methylation. With this analysis, it is possible to cover >485,000 methylation sites per sample at single nucleotide resolution. This BeadChip covers 99% of RefSeq genes with an average of 17 CpG sites per gene distributed across the promoter, 5 ′ UTR, first exon, gene body, and 3 ′ UTR regions (Bibikova et al., 2011). Additionally, it covers 96% of CpG islands, with an average of five CpG sites each, as well as the corresponding shores and shelves. After bisulfite conversion, we used 4 µl of bisulfite-converted DNA to hybridize on the array. A pair of bead-bound probes is applied to quantify the amount of thymine or cytosine in the process of hybridization. Fluorescent signal of two bead types are obtained from BeadArray Reader. The methylation status of a CpG site was calculated with the average beta-value, which was between 0 (unmethylated) and 1 (completely methylated). Beta-value is the ratio of the signal intensity of the methylated probe to the sum of the methylated and unmethylated probes.

Detection of Differentially Methylated Probes
The raw data were subject to quality control and normalization with the recommended guidelines of the bioconductor R package minfi. Wilcoxon rank test was performed to compare SO samples to OA samples. Probes with p-value was <0.05 and the delta beta ( β) >0.2 or <0.2 were considered statistically significant and differentially methylated. Delta beta value is defined as the average beta value of SO samples minus the average beta value of OA samples.

Integrating Gene Expression and DNA Methylation
To perform the integrated analysis of DNA methylation and gene expression, we used a two-step analysis process. (1) Identify differentially methylated genes in two groups. Genes differentially methylated between SO and OA groups were selected if the absolute value of delta beta was greater than 0.2 and p-value < 0.05. (2) For each differentially methylated gene, test whether there is a strong association between DNA methylated expression and transcriptional expression. Spearman rank correlation for the median of the methylation level of CpGs in on amplicon against the expression probes on the Agilent expression arrays was used to access their correlations. In this analysis, we applied correlation coefficient rho and p-value to investigate the correlation and significance between mRNA expression and DNA methylation with the function cor.test in R. Significantly negative correlation was considered if rho < 0 and p-value < 0.05, and significantly positive correlation was considered if rho > 0 and p-value < 0.05.
We next conducted unsupervised clustering analysis and displayed a different separation of SO and OA with evidence of gene clusters that are differentially expressed or methylated. Hierarchical clustering was performed using Ward linkage with Euclidean distance for samples. A clustering method available in Charm package was applied to the differentially methylated CpGs (Aryee et al., 2011).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Analysis
To find biological mechanisms underlying the genes analyzed above, we conducted GO term and KEGG pathway analysis using Cytoscape V2.7 (http://cytoscape.org/) with the ClueGo V1.3 plug-in (Bindea et al., 2009), which is able to extract biological features and annotations of anti-correlated genes. To only obtain significant enrichment categories of GO terms and KEGG pathways, the p-value was performed with right-sided hypergeometric tests and corrected by Benjamini-Hochberg adjustment (q-value). The categories with q-value < 0.05 were selected for further analysis, which were considered to be statistically significant deviation from the expected distribution. Results were visualized with the ClueGO to generate clusters of functions associated with genes.

Protein-Protein Interaction (PPI) Network Construction and Analysis
The latest experimentally confirmed human PPI data is available from human protein-protein interactions database (HPRD) (http://www.hprd.org/), which has been widely applied in human PPI network research for various disease. NetworkAnalyzer (Assenov et al., 2008) (http://www.mpi-inf.mpg.de/) is used to analyze the topological properties of biological networks, describing networks as collections of nodes and edges. Cytoscape software was used for network visualization.

Identification of Differentially Methylated Regions in SO
We performed high-throughput DNA methylation screening to compare testicular tissue-specific DNA patterns in three SO and three OA patients. To this end, we extracted genomic DNA from each sample, and made a treatment with sodium bisulfite, and calculate the DNA methylation level via the Illumina Infinium HumanMethylation450 BeadChip, which covers more than 450,000 CpG sites across the 99% of RefSeq genes. After filtering the raw data and statistical analysis, we determined that 1960 CpG sites, related to 1440 different genes, showing significant differences of DNA methylation in SO, compared with controls (absolute value of delta beta > 0.2, p < 0.05). Specifically, there were 197 hypomethylated CpG sites located in 149 genes and 1763 hypermethylated CpG sites associated with 1,317 genes. The majority of the 1,960 differentially methylated sites were located in gene bodies (42%) and within 1,500 bp from the transcription start site (TSS1500) (33%) (Figure 2A). Only 5% differentially methylated sites located in the 3 ′ UTR region, or within 200 bp from the transcription start site (TSS200). There are 4 and 11% differentially methylated sites in first Exon and 5 ′ UTR region, respectively. Moreover, based on the functional distribution standpoint, we also analyzed the amount of differentially methylated probes in CpG islands (CGI) and promoters of corresponding genes of SO vs. OA ( Figure 2B). Fort the CGI distribution, 9% of CpGs in CGI was hypermethylated, vs. 8% of CpG in CGI was hypomethylated. Additionally, 8% of CpG in promoter regions (defined as the region 1 kb upstream to 1 kb downstream of the transcription start site) hypermethylated, vs. 7% of CpG in promoter was hypomethylated. In the differentially methylated genes from the comparison between SO and OA include a number of genes with known functions in male infertility, and also some potentially novel genes for male infertility. SOX30 is one of the good examples of highly interesting genes. The SOX gene encodes a group of transcriptional factors containing a DNA binding HMG-box domain, which play important roles for the sex determination protein called SRY. In human and mouse studies, Sox30 is found to be associated with spermatogonial differentiation and spermatogenesis (Ballow et al., 2006). However, up to now, the precise function and regulatory mechanism of Sox30 have remained unknown. Recently, Han et al. found that Sox30 expression is under the control of DNA methylation status, and this expression pattern is associated with testis development in mice . Our results indicated that SOX30 is hypermethylated in SO with respect to OA, and that hypermethylation is probably related to SOX30 downexpression in patients with SO. Another interesting example in the hypermethylated gene list is PTPRN2, a receptor-like protein tyrosine phosphatase that modulate plasma membrane PI(4,5)P2 levels to promote actin remodeling. We identified that PTPRN2 gene presented highly differentially methylated difference in SO vs. AO. PTPRN2 lately was reported to be near the hypermethylated regions of sperm DNA in the patients with failure-to-conceive (Jenkins et al., 2016). Conversely, some hypermethylated genes (such as, DENND2D and LOC401127) and hypomethylated genes (such as, HLA-DPA1 and COL18A1) have not previously been associated with male infertility, which may provide novel insight into the etiology of idiopathic SO.

Integration of DNA Methylation Data with Expression Data
DNA methylation usually is considered to repress the transcriptional expression, notably when methylated sites located at the CGI of promoters. To further study the implication of DNA methylated patterns in the regulation of gene expression, we compared the detected DNA methylated profiles with our gene expression data obtained from the same sample, which could contribute to reduce the influence of potential confounding factors. After Spearman rank, our integrated analysis showed that 72 annotated CpGs and 65 genes displayed an inverse correlation with coefficient rho < 0, p-value < 0.05 (Supplemental File 1). Table 1 shown the top 20 genes with the negative correlation of DNA methylation and mRNA expression of SO vs. OA. The two hierarchical clustering of the expression levels of just these 72 methylated probes and 65 genes could completely separate out SO and OA ( Figure 3A). Of these, we identified 53 probes (located in 49 genes) that were significantly hypermethylated in SO, and 19 probes (located in 16 genes) that were significantly hypomethylated using our criteria of p-value < 0.05 and absolute delta beta > 0.2 ( Figure 3B). The genomic locations of these 65 genes are visualized using Circos. Except for Chr8, Chr13, and ChrY, these genes were distributed in other 21 chromosomes. Clusters of hypermethylated and down-expressed genes were mainly located in Chr3, Chr7, and Chr9 (Figure 4). This global view when amplified is helpful for researchers to sinvestigate the detail relationship between chromosomes, gene expression, and DNA methylation.
We then set out to determine these differentially methylated genes are involved in which biological significance relevant to SO pathogenesis. We therefore performed Gene Ontology analysis to test whether some molecular functions or biological processes were significantly associated with the genes shown the inverse difference in DNA methylation status between SO and OA. We detected some enriched GO terms showing significantly changes in SO (Figure 3C), including the following categories: activation of protein kinase activity (GO:0032147), response to alkaloid (GO:0043279), adult behavior (GO:0030534), ovulation cycle process (GO:0022602), regulation of secretion (GO:0051046), regulation of cellular localization (GO:0060341), response to organic cyclic compound (GO:0014070), exocytosis (GO:0006887), response to corticosteroid stimulus (GO:0031960), organ morphogenesis (GO:0009887). One of the hypermethylated but down-expressed genes, BRAF, plays important roles in the most significant biological process of the enrichment: activation of protein kinase activity (GO:0032147). Moreover, we also analyzed the significant pathways enriched by these 65 genes with inversely correlation between DNA methylation and mRNA expression. Signaling pathway analysis revealed functions in the regulation of actin cytoskeleton, ErbB signaling pathway, Insulin signaling pathway, neuroactive ligand-receptor interaction, chemokine signaling pathway, MAPK signaling pathway ( Figure 3D). It has been previously reported that the intracellular pathways associated  with spermatogenesis were subject to be modulated or damaged in male infertility, including MAPK signal transduction. The MAPK signaling pathway is a vital cellular signal cascade for the regulation of cell proliferation and differentiation for mammalian, which utilizes a chain of protein kinases to transduce signals from the receptors on the cell surface to the DNA in the nucleus. One gene of the series of molecules, BRAF, known to be one of components of MAPK signaling pathway, was down-expressed in mRNA level and hypermethylated in our platform. Interestingly, it has been demonstrated that BRAF-AKAP9 gene fusion leads to activate the MAPK pathway in human thyroid cancer (Ciampi et al., 2005), and AKAP9 is crucial for sertoli cell maturation and spermatogenesis in mice (Schimenti et al., 2013). This generated the hypothesis that MAPK dysregulation contributes to male subfertility might result from abnormal level of BRAF-AKAP9 gene fusion.

Correlating Network Topology with SO Mechanisms
We next examined the PPI subnetwork of 65 inversely correlated genes between mRNA expression and DNA methylation levels.
The PPI network includes 148 nodes and 285 edges ( Figure 5A).
In the PPI network the nodes with high degree are defined as hub protein, and degrees are defined to measure how many neighbors a node directly connects to. Hub proteins in the network were considered to be the key nodes that regulate molecular mechanisms of the disease. In our PPI network, the top 10 hub proteins with significant degrees contained COL13A1 (degree, 23), FYN (degree, 18), ABAT (degree, 16), etc. (Table 2). For example, FYN, a kinase which is abundantly expressed in the testis and is involved in testis and sperm function, is closely associated with proper shaping of the head and acrosome (Luo et al., 2012). Furthermore, to identify the most relevant cellular activities controlled by these 22 both inversely correlated and PPI existed genes, we analyzed GO biological process terms and KEGG pathways. The most significant GO terms were related to synaptic transmission, purine ribonucleoside triphosphate catabolic process, response to steroid hormone stimulus, cellular response to hormone stimulus, positive regulation of response to stimulus, apoptotic process, cell surface receptor signaling pathway, regulation of immune response, behavior ( Figure 5B). Three genes EPHX1, SERPINB9, and SLIT3 play common roles in the cellular response to the hormone stimulus (q-value = 0.001836) and the response to steroid hormone stimulus (q-value = 0.000691). Compared our data with those reported studies, we found that EPHX1 (Epoxide Hydrolase 1) appeared to be luteinizing hormone (LH) dependent in adult-type Leydig cells of mice (Griffin et al., 2010). In addition, genetic mutations in epoxide hydrolases affect the predisposition for oligozoospermia and asthenospermia in FIGURE 5 | Using network topology to infer elements involved in SO. (A) Interaction network for genes differentially expressed and DNA methylated in SO. The red circles represent the hub genes inversely correlated DNA methylation and gene expression. The purple circles indicate the genes closely interacted with the hub genes that are highly linked nodes in network. (B) Gene ontology analysis of the relationship between connectivity of the genes in this network. Node shown in red, green, or blue is associated with biological process, cellular component, or molecular function, respectively. The p-values were calculated using hypergeometric tests and corrected using the Benjamini-Hochberg adjustment (q-value). The q-values are expressed and presented as negative logarithms (base 10).
Han-Chinese population (Qin et al., 2012). SERPINB9, which hypopemethylation is associated with mRNA overexpression in our platform, probably protect against inadvertent release of granzyme B, and premature or unwanted activation of cell death pathways (Bird et al., 1998). Recently, Ottesen and his colleagues shown that SLIT3 is one of most strongly upregulated genes in the a spinal muscular atrophy model, which presented disturbed testis development, degenerated seminiferous tubules, decreased sperm count and low fertility (Ottesen et al., 2016). Several lines of evidence have shown that these key genes in our list are potentially involved in spermatogenesis, and may have clinical applications as novel biomarkers of male infertility.

DISCUSSION
In this study, we have identified both novel and reported dysregulated genes involved in SO at the DNA methylation level and mRNA expression level. By using the integrated analysis of the DNA methylation and mRNA expression data, we have established a pipeline for investigating the complexity of gene dysregulation in the context of SO when using primary samples. The stepwise procedure is as follows: (1) Identify significantly differentially methylated genes and expressed genes in the case and control groups; (2) Retain those genes whose methylation and expression levels are highly anti-correlated; (3) Enrich the Gene Ontology, KEGG pathway and PPI networks of anticorrelated genes; (4) Identify key genes from a combination of biological information, including interactions and relationships among genes, gene functional annotations and pathway maps. In particular, the topological properties of this approach could be applied into explaining roles of unannotated proteins in diseases. The function of unannotated proteins could be predicted by further analysis of its direct neighbors and their neighborhood in the networks.
Recently, integrated analysis between transcriptional expression and DNA methylation has been carried out in a number of different disorders, including pilocytic astrocytomas (Zhou and Man, 2016), lung adenocarcinoma (Selamat et al., 2012), melanoma (Li et al., 2013), Crohn's disease-associated fibrosis (Sadler et al., 2016). For severe oligozoospermia, few studies have previously examined the existence of genes with aberrant DNA methylation or altered mRNA expression. There was no previous systematic effort to decipher the genetic mechanisms of SO in both regulatory layers. To our knowledge this study presents the first work of integrating microarray data of SO and OA in the progress of deciphering the cause of human male infertility. Our integrated analysis not only has allowed us to confirm changes (such as, PTPRN2 and FYN) described by others, but also to identify several novel genes with aberrant DNA methylation profiles, including DENND2D, HLA-DPA1, COL18A1. We also confirmed some genes that are relevant to the SO phenotype in mice model, such as SOX30 and EPHX1 . Additionally, a small percent of methylation modifications in our SO result occurs at CGI of promoters and TSS downstream of genes. We observed that transcriptional activity is also affected by methylation alterations at gene bodies, which is corresponding to recent reports (Yang et al., 2014). We listed all above genes whose DNA methylation changes were inversely correlated with mRNA expression changes in the Supplemental Material, we also identified several genes that were coordinately hypermethylated and up-regulated, or hypomethylated and down regulated. While these two groups of genes do not fit into the classical paradigm of DNA methylation regulation, emerging evidence from recent high through-put studies show that DNA methylation regulation may be more complex. From other aspects, miRNA regulation is one of mechanisms also affects gene expressions (Zhuang et al., 2015;Li et al., 2016). Integration of the analysis of three profiles of mRNA expression, DNA methylation and miRNA expression suggests that the latter two regulations could act in the same or in opposite directions on gene expression (De La Rica et al., 2013).
In this study, a mixture of several kinds of cell was performed. Given that variations in the expression and methylation profiles might be affected by the difference in the ratio of each kind of testicular cell in the primary sample. Germ cell fractionation or laser microdissection also constitutes an approach of validating our present result. However, above two methods have their own limitations. Microdissection is time consuming for a large number of cells required in the microarray studies (Ferfouri et al., 2013). Germ cells fractionation seems impossible because a few amount of germ cells in the specimens of the SO patients. We can expect that much fewer cells will be required with the development of technology in the future.
The emergence of new technologies to analyze gene expression and DNA methylation has allowed the observation of changes at the global genome level. The primary sample is usually a good option for the research on male infertility using these high-throughput technologies. Distinct with many studies of integration analysis which used different resources of samples from the datasets, our study systematically investigated both DNA methylation and mRNA expression in order to overcome the limitations of sample size. In fact, integrative DNA methylation and gene expression analysis has been successfully used to identify a novel dysregulated biomarker discoidin domain receptor 1 (DDR1) in non-obstructive azoospermia (NOA) patients (Ramasamy et al., 2014). Given the small sample size and high-throughput data gathering techniques, the results should nevertheless be interpreted with caution. In our study, patients of two groups are selected for their similar pathological traits, including histology characteristics, hormone levels, and biochemical parameters (Supplemental File 2). Thus, the screened genes with significantly different expressions in subjects of two groups should be due to distinct pathology.


The frequency of Y chromosome microdeletions in infertile men of Vellore cohort

Many male infertility cases are because of genetic and environmental factors and most of them are idiopathic. About 10-20% of azoospermic patients are showing the microdeletion in Y-chromosome. The azoospermia factor or AZF region at the Yq11 position which consists of genes those are necessary for spermatogenesis. In Y chromosome microdeletion, deletion in AZFa region is very rare. The aim of our study is to find out the frequency of microdeletions in Y chromosome particularly in AZFa region in azoospermic men of Vellore cohort. For this study, we collected 10 azoospermic patient and 10 control men samples from the Sandhya hospital, Vellore. In this study, we mainly focused on AZFa region to analyze the frequency of microdeletions in Y chromosome using SY82 (264bp) and SY83 (275bp) STS markers. There was an absence of microdeletion in patient samples for SY82 and SY83 markers of AZFa region. We concluded that sample size should be increased to confirm our results.

Introduction
Male infertility alludes to a male's inability to bring pregnancy in a prolific female [1]. Y chromosome microdeletion (YCM) is the second most successive hereditary reason for male infertility. It is a group of the hereditary issue brought on by missing genes on the Y chromosome [2]. Many men with Y chromosome microdeletion showed no side effects and had ordinary existences. There are mainly four regions present on Y chromosome i.e. AZFa, AZFb, AZFc, and AZFd responsible for Y chromosome microdeletion. AZFa microdeletion is found very rare [3]. Its size is 1MB and it has two candidate genes USP9Y and DBY. AZFb region is located in the middle region of Yq11 where a proximal piece of AZFb comprises of vast repeated sequences sorted out in palindromes, yet a large portion of it is the sequence of single copy [4]. At the point when both markers sY254 and sY255 are deleted, a finding of complete deletion of the AZFc region can be made. A few reviews have demonstrated that the AZFc deletion example is somewhat consistent,in spite of the fact that not generally indistinguishable. The most successive deletion sort is the AZFc region deletion trailed by AZFa, AZFb and AZFbc deletion [5]. Deletions which are distinguished as AZFabc are doubtlessly identified with irregular karyotype. Yq microdeletion testing is helpful for diagnosing infertility in men.
AZFa region deletion was observed very less frequently in the men who are infertile and were also associated with highly severe testicular phenotype with completely absent germ cells. It also contains three candidate genes which are USP9Y, DBY and UTY. They represent high homology with interval or deletions which may cause spermatogenic impairment. The main role of these genes in the human spermatogenesis is not known. It is also not clear about which of them may be responsible for the phenotype of AZFa. Thecomplete AZFa deletion removes ~792 kb of the DNA sequences, including the two candidategenes [6]. AZFa genes deletion analysis in a vast number of men who are infertile is characterized by the highly defined spermatogenic alteration. Both the genes USP9Y and the DBY may also cause testiculopathies, but the gene DBY only appears to be the major candidate of AZFa.
YCM diagnosis grants the cause for the patient's azoospermia to be built up and to detail a prognosis [7]. Genetic counselling is required to provide information about the danger of conceiving a child with weakened spermatogenesis. Men with YCM will pass on their defect in their male offspring after going through ICSI treatment [8]. The testing of Yq microdeletion is done using PCR and involves the use of various STS markers for the detection of deletion occurred on Y chromosome. Therefore, in the present study we performed molecular diagnosis for AZFa region using STS-PCR.

Study population
This study was done with the 10 azoospermic patients and 10 healthy controls between the age group of 18 to 45 years by performing semen analysis according to WHO, 2010 guidelines [9]. All of those individuals who were taken in the group of control were the healthy male without any symptoms of infertility or urological disease. They were recruited from the Sandhya hospital in Vellore, Tamil Nadu. The study was further approved by the ethical committee of Vellore Institute of Technology (VIT), Vellore. The sheet of the data collection from volunteers included the information of the population such as family history, age, BMI, smoking habits and physiological features. Informed consent form was given to volunteers. Genetic counselling gave some clue that mental state of the patient which indicated that up to some extent even depression, tension or the state of unhappiness may lead to the condition of infertility.

Genotype Determination
The extraction of the DNA was carried out from the 2ml intravenous blood sample which was collected in EDTA vacutainer. And it was further stored at -20°C for molecular analysis. This extraction was done by using standard lab protocol [10]. We have quantified DNA sample using Nano drop, 2000. We also qualified those DNA samples using Agarose gel electrophoresis technique. The STS primers SY82 and SY83 were used. The product size of SY82 primer was 264bp and for SY83 product size was 275bp. The sequence for SY82 forward primer was 5´-ATC CTG CCC TTC TGA ATC TC-3´and reverse was 5´-CAG TGT CCA CTG ATG GAT GA-3´. The sequence for the SY83 primer was 5´-CTT GAA TCA AAG AAG GCC CCT-3´ and reverse was 5´-CAA TTT GGT TTG GCT GAC AT-3´. The conditions used in PCR were such as-initial denaturation phase was at 95ºC for 5 minutes, denaturation phase was at 94ºC for 1 minute, the annealing temperature was 54ºC for both primers, the extension was done at 72ºC for 1 minute and final extension was at 72ºC for 10 minutes followed by thirty-five cycles. The PCR products were then analyzed for Y chromosome microdeletion by electrophoresis on 2% agarose gels. The absence of band in any STS markers indicated microdeletion by comparing it with the control. Sample showing deletions were tested again to confirm the results.

Results and Discussion
The patient's mean age was 33±6.7 (mean ± SD) years (range: 18-45 years) at the time of inclusion and all had primaryin fertility. There was no microdeletion found in SY82 and SY83 of AZFa region in azoospermic patients. Example of SY82 PCR products in AGE was shown in figure 1. In previous studies, it is reported that the deletion of the single the frequency of AZFa deletion is almost double (11 vs 5 %) in Indians than western population [11,12]. Therefore, we assumed that more sample size should be needed to confirm our results.

Conclusion
In a summary, this study showed the absence of Yq microdeletion in AZFa region in azoospermic men of Vellore population. Our findings are disagreeing with previous reports performed in different societies and ethnic groups in India. So, it is needed to screen more sample size with Yq microdeletions in azoospermic patients during diagnosis for male infertility and it should be carried out before ART. This will help to find out the etiology of male infertility.