GWAS of genetic factors affecting white blood cell morphological parameters in Sardinians uncovers influence of chromosome 11 innate immunity gene cluster on eosinophil morphology

Abstract

Few genome-wide association studies (GWAS) analyzing genetic regulation of morphological traits of white blood cells have been reported. We carried out a GWAS of 12 morphological traits in 869 individuals from the general population of Sardinia, Italy. These traits, included measures of cell volume, conductivity and light scatter in four white-cell populations (eosinophils, lymphocytes, monocytes, neutrophils). This analysis yielded seven statistically significant signals, four of which were novel (four novel, PRG2, P2RX3, two of CDK6). Five signals were replicated in the independent INTERVAL cohort of 11 822 individuals. The most interesting signal with large effect size on eosinophil scatter (P-value = 8.33 x 10⁻³², beta = −1.651, se = 0.1351) falls within the innate immunity cluster on chromosome 11, and is located in the PRG2 gene. Computational analyses revealed that a rare, Sardinian-specific PRG2:p.Ser148Pro mutation modifies PRG2 amino acid contacts and protein dynamics in a manner that could possibly explain the changes observed in eosinophil morphology. Our discoveries shed light on genetics of morphological traits. For the first time, we describe such large effect size on eosinophils morphology that is relatively frequent in Sardinian population.

Introduction

Large-scale genome-wide association studies (GWAS) have revealed thousands of loci affecting the levels of white blood cell subpopulations (1,2). However, previous GWAS scans have tended to focus on common genetic variants, while rare coding variants were often not directly assessed. The use of dense coverage genetic maps thus accesses more of genetic landscape and can also increase the likelihood to identify causal variants or proxies for the genotyped single-nucleotide polymorphisms (SNPs) (3).

The range of studies can also be extended from standard diagnostically used parameters to relatively non-classical parameters, enabling the exploration of the mechanism of regulation of white blood cell levels and their effects on related diseases. Previous GWAS studies using data on white blood cells derived from hematology analyzers have focused mainly on cell counts. Novel technologies of hematology analyzers provide additional detailed cell morphology analysis. WBC differential technology biophysically characterizes blood cells and detects abnormal cells through flow cytometric analysis with multiple light-scatter angles and digital conductivity. Cells are unaltered by staining. For example, DxH Coulter Cellular Analysis VCS technology measures cell volume (V, estimate cell size), conductivity (C, that is used to collect information about the permeability of cells, internal structure of the cell, including chemical composition and nuclear volume) and light scatter (S, that gives information about internal complexity like cellular granularity, nuclear lobularity and cell surface structure) in neutrophils, lymphocytes, monocytes and eosinophils (4).

Standard blood cell diagnostic evaluations include levels of different cell populations. For red blood cells, many morphological features are also analyzed, and platelet volume has also been usually reported. By contrast, morphologic features of white blood cells are rarely included in previous GWAS blood cell evaluations (2). However, previous GWAS on polychromatic flow cytometry–fluorescence-activated cell sorting (FACS) analyses evaluated some of these parameters (5). To date, only one study comprehensively analyzed such parameters assessed by a hematology analyzer (6). However, they are relevant in diagnosis and research and can be determined by appropriate laser settings and light refraction. Such non-classical cell parameters, are reported based on histogram data. They may be related to clinical conditions, assessed by statistical co-localization of genetic association data with other types of data. The possible function of the associated variants can also be further analyzed to look for protein folding effects of nonsynonymous variants providing insights into the function of rare coding variants.

Here, we report a GWAS of morphologic non classical white blood cell traits such as cell volume, conductivity and light scatter (7). We replicated five signals in the Interval cohort. Of two rare Sardinian specific variants in PRG2 and HBB, the second was already described in our FACS previous study (5). Using molecular dynamics (MD) analysis we infer that the previously undetected mutation in PRG2 could interfere with the mature PRG2 non-toxic granules formation, possibly influencing the crystallization/aggregation processes that mediate PRG2 toxicity.

Results

GWAS analysis

Interrogating 18 425 755 variants, either directly genotyped or successfully imputed, we conducted a GWAS on 12 white-blood cell morphological traits including measures of cell volume (V), conductivity (C) and scatter (S) in four white cell subpopulations (eosinophils, lymphocytes, monocytes, neutrophils), in 869 volunteers who are part of the SardiNIA general population cohort (Supplementary Material, Table S1, Figs S1–S3) (8). This analysis yielded seven association signals with the assessed traits (Table 1) (P < 4.17 x 10⁻⁹), three known (LYZ, FAM117B, HBB) and four novel (PRG2, P2RX3, two associations with CDK6).

Table 1

Associated loci with morphological features of white blood cells

Trait	Chr.	Position hg37	Effect allele/other allele	MAF Sardinia	MAF EUROPE	PVALUE	BETA	SE	Gene	Snp142	Consequence	Novel/known	R2	h2
EOS	11	57 156 106	G/A	0.031	6.2E-5	8.33E-32	-1.651	0.1351	PRG2	rs769591668	p.Ser148Pro NP_002719.3	novel	0.1470	0.1638
EOS*	11	57 110 622	A/G	0.114	0.142	1.37E-14	-0.5928	0.0757	P2RX3	rs80116434	Intronic	novel	0.0661	0.0710
EOV	7	92 408 370	T/C	0.183	0.096	9.68E-31	0.7048	0.0588	CDK6	rs445	Intronic	novel	0.1422	0.1485
NEV	7	92 408 370	T/C	0.183	0.096	3.94E-15	0.4902	0.0613	CDK6	rs445	Intronic	novel	0.0688	0.0710
MOS	12	69 744 014	A/C	0.082	0.040	1.54E-12	-0.6264	0.0873	LYZ	rs1800973	p.Thr88Asn NP_000230.1	known	0.0561	0.0591
MOV	2	203 591 785	CT/C	0.505	0.291	2.76E-09	0.2969	0.0494	FAM117B	n.a.	Intronic	known	0.0399	0.0441
LYV	11	5 248 004	A/G	0.058	6.2E-4	2.13E-09	0.6515	0.1077	HBB	rs11549407	p.Gln40Ter NP_000509.1	known	0.0405	0.0464

EOS: eosinophil scatter; EOV: eosinophil volume; NEV: neutrophil volume; MOS: monocyte scatter; LYV: lymphocyte volume.

MAF Europe – source gnomAD browser (gnomad.broadinstitute.org), positions according to the build GRCh37.

In cursive known signals.

R  ² represents the percentage of trait variance explained by the associated variant; h² represents the narrow heritability explained by the associated variants.

^aAfter conditional analysis.

The fraction of trait variance explained by each associated locus ranges from 0.0399 to 0.147 while the heritability from 0.0441 to 0.1638 (Table 1).

Novel loci

In particular, PRG2:p.Ser148Pro, rs769591668 mutation, in exon 4 of the proteoglycan 2 gene (PRG2) gene was associated with reduced eosinophil scatter (EOS, cell granularity), P-value = 8.33 x 10⁻³², effect = −1.651 (Table 1, Fig. 1, Supplementary Material, Fig. S2A). PRG2 encodes a cytotoxic polypeptide contained in the eosinophil granules (9). PRG2 148Pro variant is rare or absent outside but relatively common in Sardinia (Minor Allele Frequency [MAF] Europe = 6.2 x 10⁻⁵, MAF Sardinia = 0.031, respectively in gnomAD browser) (Table 1). PRG2 represents the predominant constituent of the crystalline core of the eosinophil granule protein, and enhances natural killer cell activity, providing microbicidal activity against fungi and gram-negative and gram-positive bacteria (10). It may also be involved in antiparasitic defense mechanisms as cytotoxin and helminthoxin and in immune hypersensitivity reactions as well as in non-cytolytic histamine release from human basophils. After conditional analysis to identify independent association signals in the same locus, we discovered another association in the same region. Allele A of rs80116434 was associated with reduction of EOS. This variant is located within the P2RX3 gene that encodes for the purinergic receptor P2X P-value = 1.37 x 10⁻¹⁴, effect = −0.5928 (Table 1, Fig. 1, Supplementary Material, Fig. S2B). This variant is described to reduce eosinophil count (1). Interestingly, rs80116434 is located in the eosinophil regulation innate immunity genes cluster (P2RX3, PRG3, PRG2) on chromosome 11.

Figure 1

Step-wise conditional analyses of the region of PRG2 on chromosome 11 associated with EOS. (A) EOS association with the top variant 11:57156106; (B-C) first and second conditional analysis association. After the first conditional analysis the association of the second variant region 11:57110622 with EOS has been showed. In each plot, the y-axis shows the association strength (−log₁₀(P-value)) versus genomic position (hg19/GRCh37) of the most significant SNP on the x-axis. Other SNPs are colored to reflect their LD with the top SNP. Genes and positions of exons are shown below.

The associated variant in CDK6 included rs445, with allele T increasing the volume of myeloid lineage eosinophils and neutrophils: (EOV) P-value = 9.68 x 10⁻³¹, effect = 0.7048; (NEV), P-value = 3.94 x 10⁻¹⁵, effect = 0.4902 (Table 1, Supplementary Material, Figs S2C, D and S3A and B). Notably, rs445 had already been associated with increased volumes of granulocytes and monocytes in a previous study that used flow cytometry FACS data (5). In addition, it has been previously associated with reduced counts of myeloid cell line subpopulations (1,2). Here, we discovered that it is associated with eosinophil and neutrophil volume, pointing to a role for this variant in the myeloid cell lineage.

Known loci

As we previously reported in flow cytometry FACS data, the nonsynonymous LYZ: p.Thr88Asn mutation, rs1800973 located in the lysogene (LYZ) gene (allele 88Asn) is associated with reduced monocyte scatter (MOS) P-value = 1.54 x 10⁻¹², effect = −0.6264 (Table 1, Supplementary Material, Figs S2E and S3C). This variant was classified as benign by ClinVar. Moreover, this variant increases monocyte counts, reduces granulocyte and lymphocyte count and side scatter in monocytes in our recent study (5). It creates an eQTL and pQTL in blood, augmenting both transcript and protein levels, as reported in Open Targets data (11). Thus, it has been associated with increased LYZ protein in blood plasma and augmented gene expression (12).

In addition, we detected an intronic variant rs149760691 within the FAM117B gene. Allele CT at this locus augmented monocyte volume (MOV) P-value = 2.76 x 10⁻⁹, effect = 0.2969 (Table 1, Supplementary Material, Supplementary Figs S2F and S3D). This variant has been previously associated with side scatter in flow cytometry data (5).

We also replicated another of our recent discovery in flow cytometry FACS data of a founder Sardinian mutation p.Gln40Ter (rs11549407) within HBB gene. Allele 40Ter was associated with reduced lymphocyte volume (LYV; Table 1, Supplementary Material, Figs S2G and S3E) (5).

Replicated loci in INTERVAL study

We replicated five signals in the INTERVAL cohort (Table 2) (6,13,14). For that study, the cohort of 11 822 whole-genome sequenced samples was used. In the INTERVAL study, 63 non classical phenotypes were measured, of which 11 summarize measurements of intracellular complexity/granularity (SSC), 16 summarize measurements of cell nucleic acid content/membrane lipid content (SFL) and 15 summarize measurements of cell morphology/volume (FSC). About 20 of the traits were platelets phenotypes, 10 were red cell phenotypes and 33 were white cell phenotypes (6). We considered a replication threshold of P-value = 0.05 adjusted for multiple variant testing using the Bonferroni method (0.05/7 variants to be replicated, thus P-value = 0.007). At this threshold, we replicated five signals, of which three novel. The remaining two signals are rare, Sardinian specific variants in PRG2 and HBB (MAF = 6.2E-5 and 6.2E-4 in Europeans in GnomAD browser, gnomad.broadinstitute.org, respectively), and were not replicated in the INTERVAL cohort study (Tables 1 and 2).

Table 2

The replication of associations in interval study

Trait	Chr.	Position Hg37	Effect allele/other allele	MAF Sardinia	MAF Interval	P-value	BETA	SE
EOSa	11	57 156 106	G/A	0.031	NA	NA	NA	NA
EOS	11	57 110 622	A/G	0.114	0.135	9.29E-72	−0.3704	0.0205
EOV	7	92 408 370	T/C	0.183	0.099	1.61E-18	0.2063	0.0235
NEV	7	92 408 370	T/C	0.183	0.099	1.03E-27	0.2511	0.0230
MOS	12	69 744 014	A/C	0.082	0.060	2.78E-20	−0.2626	0.0284
MOV	2	203 591 785	CT/C	0.505	0.299	0.0038	0.0430	0.0149
LYVa	11	5 248 004	A/G	0.058	NA	NA	NA	NA

^aFounder Sardinian mutations.

Molecular dynamic analysis of PRG2 mutation

To examine possible consequences of the PRG2:p.Ser148Pro mutation on molecular structure and conformational changes of PRG2 (also known as eosinophil granule major basic protein or MBP-1) we performed molecular modeling analysis based on the X-ray crystal structure (PDB ID 1H8U) (15). Analysis of the trajectories revealed that the mutation of Ser148-to-Pro148 mainly affects the atomic fluctuations (RMSF) of the loop 164–171 and the surroundings of position 148 (Fig. 2A-C). The presence of a Pro residue at position 148 restricts the atomic fluctuations of the neighboring residues. The stiffness induced by proline is transmitted through the contiguous beta-strand to the residues located on loop 164–171, whose positions and fluctuations are modified (Fig. 2D-E).

Figure 2

Consequences of the PRG2:p.Ser148Pro mutation on molecular structure and conformational changes of PRG2. (A) Root mean square deviation (RMSD), (B) radius of gyration (RG) and (C) root mean square fluctuations (RMSF) of the simulated systems with WT (black) and S148P (red) PRG2. Region corresponding to C-end of helix 138–152 and loop 164–171 is shadowed in gray, and the serine residue that is mutated is labeled by dashed line. (D) Average ribbon representation of the above systems after MD simulation were calculated from the trajectories of the last 50 ns of simulation. The mutated position (S148P) is colored in blue. (E) Detailed view of the vibrational transmission from extended, rigid C-end Pro148-containing α helix (red) to loop 164–171 (purple) through contiguous β strands (blue) of S148P PRG2.

Moreover, an in-depth analysis of the trajectories and contacts between residues showed that substitution of Ser148 by Pro148 also affects the contact network of Tyr109 (Supplementary Material, Fig. S4). Tyr109 establishes close contact with Ser148 along the trajectory of WT PRG2. Pro148 impairs this contact and induces the rotation of Tyr109 sidechain that forms a different contact network compared to the WT PRG2. This change in the sidechain configuration could contribute to the transmission of the perturbation induced by proline to the loop 164–171. In addition, Ser148-by-Pro148 mutation also affects the electrostatic surface potential of residue located nearby loop 164–171, as Asp176 is more exposed to solvent (Supplementary Material, Fig. S5).

Functional annotation of the novel signals

We performed comprehensive queries of the Open Target PHEWAS catalog for novel variants (genetics.opentargets.org). Multiple associations with white blood subpopulations were obtained for rs445 including reduced counts of neutrophils and eosinophils, and for rs80116434 with reduced eosinophil counts, corroborating their involvement in these cells function (Supplementary Material, Table S2).

According to the Human Protein Atlas (proteinatlas.org) RNA single cell analysis confirms that PRG2 is expressed in granulocytes and in flow sorting data that include RNA-seq data set that is expressed in eosinophils. In the same analysis, CDK6 is ubiquitously expressed gene in different tissues that includes granulocytes in RNA single cell data set and eosinophil and neutrophil white blood cells in flow cytometry RNA data set.

Moreover, according to the UCSC genome browser (genome.ucsc.edu) rs445 that relay in the second intron of CDK6 gene is located in the regulatory region marked by the H3K27Ac site in GM12878 cell line, in the DNase I hypersensitivity cluster and in the enriched of transcription factor binding sites region (Supplementary Material, Fig. S6). Also, GM12878 cells chromatin state segmentation by HMM from ENCODE/Broad in UCSC genome browser shows strong enhancer feature in the rs445 region (Supplementary Material, Fig. S6). We used assay for transposase-accessible chromatin using sequencing (ATAC-seq) data, which interrogate chromatin accessibility that depicted enhancer presence in the rs445 region (Supplementary Material, Fig. S7). We showed that rs445 shows significant chromatin accessibility in all hematopoietic linages. This suggests that rs445 is located within active enhancer in all hematopoietic lineages. Moreover, chromatin regulation 3D datasets show a contact between the region that includes rs445 and the promoters of genes CDK6 and AC002454.1 (Supplementary Material, Fig. S8A). The gene expression data in different blood cells showed that CDK6 is significantly co-expressed with AC002454.1 (antisense located in the 5’UTR of CDK6) in hematopoietic precursors that include myeloid cell line precursors, showing potential regulatory character of rs445 in those cells (Supplementary Material, Fig. S8B).

Discussion

A population sequencing GWAS for non-classical morphologic parameters of white blood cells enabled the detection of a coding variant in the PRG2 gene affecting eosinophil scatter (EOS)-reduced cell granularity. After conditional analysis a second variant in the P2RX3 gene that reduced cell granularity in eosinophils (EOS) was discovered. Furthermore, a variant in the CDK6 gene, previously associated with myeloid cell counts and intracellular complexity, was found here associated with neutrophil and eosinophil volume (EOV, NEV) (Supplementary Material, Fig. S3).

PRG2 novel locus

The most interesting result is the association of EOS with the Sardinian-specific nonsynonymous p.Ser148Pro variant within proteoglycan 2 (PRG2), which is relatively frequent in Sardinia. p.Ser148Pro was also associated with a null/very low phenotype of eosinophil counts (almost complete absence or low of eosinophils in four waves data; unpublished own data). Thus, individuals homozygous for rs769591668-C were characterized by absence of eosinophils or, in some instances, very reduced eosinophil levels. Furthermore, eosinophils, when detected, were also characterized by reduced granularity. Importantly, Prg2^tm1Nal/Prg2^tm1Nal mice have abnormal inflammatory responses along with respiratory system inflammation and abnormal eosinophil morphology (MGI; http://www.informatics.jax.org), in agreement with our observations. Moreover, in another study, apparently healthy MBP1^−/−/EPX^−/− double knockout mice displayed significantly fewer circulating peripheral blood eosinophils along with disruption of eosinophilopoiesis (16). However, the presence of reduced scatter in our data on this cell population suggests that there are present eosinophils that are missed in the classical measurement due to the absence of granules, similar to the observations in the mouse knockout model.

The relatively high frequency of PRG2 148Pro variant in Sardinia was unexpected, in light of recent findings that innate immunity genes globally evolved to remove loss-of-function variants through purifying selection, which could explain some auto-inflammatory and hypersensitivity reactions now observed in humans (17,18). This relatively high frequency of morphological changes in eosinophils further agrees with the lower levels of eosinophils reported recently in children exposed to dust in traditional Amish farming (19).

Replication in independent cohort

We replicated five of the seven associated signals in the independent INTERVAL cohort. Notably, the two not replicated loci PRG2 and HBB are Sardinian-specific, founder mutations that are almost absent outside of Sardinia. Sardinians are one of the oldest European population where the distinct evolutionary history and genetic structure, genetic drift and natural selection caused the higher allele frequencies that are rare in the world. However, the second signal in HBB locus was previously described in our FACS data. With these results we demonstrated that our data are robust and warranted followed-up with functional analyses.

Molecular dynamic analysis of PRG2 mutation

PRG2 toxicity is believed to function by disrupting membranes of parasites and bacteria. PRG2 (MBP-1) is packed in the specific granules of human eosinophils as a distinctive non crystalline structure, permitting the storage of otherwise toxic proteins. PRG2 toxicity is restrained via crystallization in eosinophil secretory granules with tinctorial properties typical of amyloid-like structures (20). Soragni and coworkers have shown that residues 114–119, 131–143, 146–159 and 194–202 of PRG2 have the highest aggregation propensity. Among them, the stretch between residues 131 and 159 form a nearly continuous solvent-exposed aggregation-prone area. Our computational analysis showed that substitution of Ser148 by Pro148 not only alters the dynamic properties of the loop between residues 164–171 but also affects the contact network of Tyr 109. The changes observed in the eosinophil morphology could be explained based on these findings. Furthermore, because Ser148 is located in the middle of the solvent-exposed aggregation-prone area (residues 131–159), we propose that mutation by Pro can interfere with the formation of PRG2 non-toxic granules, tuning the crystallization/aggregation processes that mediate PRG2 toxicity.

Functional annotation of the novel signals

We observed that all three novel variants regulate eosinophil morphology and in parallel two of them (rs445, rs80116434) were associated with reduced eosinophil counts in public data pointing on the genetic regulation of those traits. Additionally, consistently with CDK6 function, the rs445 regulated neutrophil volume and was associated with reduced counts of different myeloid cell subpopulations. This is in agreement with its ubiquitous expression in different tissues that included neutrophils and eosinophils. Using public data we showed the regulatory character of rs445 that is located on active enhancer in all hematopoietic lineages. In particular, rs445 region is included in the contact of promoters of genes CDK6 and nearby located in 5’UTR antisense AC002454.1 . In line with rs445 function in myeloid cell line and underlying its regulatory character, those two genes were co-expressed in myeloid cell line precursors.

Materials and Methods

Cohorts and sample description

The SardiNIA project is a longitudinal cohort study of about 8000 general population volunteers (57% females), ranging from 18 to 102 years recruited in the Lanusei Valley, Sardinia. All volunteers are deeply genetically characterized and phenotyped for thousands of quantitative trait. In this study, we used data from 869 individuals. All participants gave informed consent to the study protocols, which were approved by the Sardinian local research committee Ethical Committee of ASSLL of Sassari (2171/CE).

The INTERVAL cohort includes blood donors aged 18 years or older from donor centers of NHS Blood and Transplant across England, UK (13,14). For this study, 11 822 whole-genome sequenced samples were used.

Phenotype measurement

We studied 12 non-diagnostic blood cell phenotypes: cell volume (V), conductivity (C) and light scatter (S) in four white blood cell subpopulations (eosinophils, lymphocytes, monocytes, neutrophils). These traits were measured with the UniCel DxH Coulter Cellular Analysis system of Beckman Coulter according to the manufacturer’s instructions. Thus, we studied associations in eosinophils (volume EOV, conductivity EOC, scatter EOS), lymphocytes (volume LYV conductivity LYC, scatter LYS), monocytes (volume MOV, conductivity MOC, scatter MOS) and neutrophils (volume NEV, conductivity NEC, scatter NES) (Supplementary Material, Table S1).

In the INTERVAL study, the phenotypes were acquired with the Sysmex XN-1000 instrument. The fluorescence intensity of light emitted by the dye measures the nucleic acid content and membrane permeability of the cell, while the intensities of diffracted light scattered sideways (SSC) and forward (FSC) by the cell measure intracellular organelle complexity and cell size, respectively (6).

Genotyping and imputation

Genetic analyses were performed on 869 samples genotyped with OmniExpress Illumina array as previously described for the entire SardiNIA cohort (21). Imputation was performed on a genome-wide scale using the Sardinian sequence-based reference panel of 3514 individuals and the Minimac software on pre-phased genotypes. After imputation, only markers with imputation quality (RSQR) > 0.3 for MAF ≥ 1% or RSQR > 0.6 for MAF < 1% were retained for association analyses (21), yielding over 18 425 755 variants (17 209 121 SNPs and 1 216 634 indels) useful for analyses.

Association analysis

Before performing GWAS, the 12 immune phenotypes were both inverse-normalized and adjusted by sex, age, age² as covariate. Additive effects were searched using EPACTS version 3.2.6 (22). The software implements a linear mixed model adjusted with a genomic-based kinship matrix calculated using all quality-checked genotyped SNPs with MAF > 1%. The advantage of this model is that the kinship matrix encodes a wide range of sample structures, including both cryptic relatedness and population stratification. As a proof of appropriate adjustment of all confounders, the genomic inflation factor (lambda GC) for 12 GWAS was on average 1.005 (range 0.990–1.016). To identify independent signals, a conditional GWAS analysis for each trait was performed by adding the leading SNPs found in the primary GWAS as covariates to the basic model. A SNP reaching the standard genome-wide significance threshold (P-value < 5 x 10⁻⁸) adjusted using Bonferroni multiple testing method (P-value < 4.17 x 10⁻⁹) was considered significant.

Genotyping and quality controls for the INTERVAL study was described in Astle, Elding, Jiang et al. (2).

The replication P-value threshold was calculated by adjusting the nominal P-value (0.05) for multiple variant testing using the Bonferroni method (0.05/7 variants to be replicated), obtaining a P-value = 0.007.

The proportion of trait variance explained by each variant was calculated as the R² of the single variant association model in EPACS. The heritability explained by each variant was obtained using the following formula:

|${\textrm{Effect}}^2\times 2\times \textrm{MAF}\times \Big(1\hbox{--} \textrm{MAF}\Big)$|

Molecular modeling analysis

Initial structure for Homo sapiens PRG2 was based on the crystallographic structure of PRG2 (PDB ID 1H8U) (15). Sulfate and glycerol molecules present in the crystallographic structures were removed prior to the simulations. PRG2S148P variant was built using WT protein as a template with UCFS Chimera (23). MD simulations of WT and S148P proteins were performed with OpenMM software (24) in a CUDA platform with AMBER suite ff19SB force field (25), particle mesh Ewald for long-range electrostatic interactions with and Ewald summation cutoff of 10A, and constraint of the length of all hydrogen bonds. Each system was neutralized with chlorine and sodium counter-ions according to the total charge of the protein and solvated with optimal-3-charge, 4-point rigid water model (OPC) water molecules (26). Firstly, each system was subjected to 2500 steps of energy minimization at 298 K. Langevin thermostat was used to control the temperature with a friction coefficient of 1 ps⁻¹ and a step size of 0.002 ps. Each system was subjected to 250 ns MD simulation and snapshots of the trajectories were saved every 1 ps to evaluate the motions of residues and atoms. Finally, the trajectory analysis was performed with the CPPTRAJ module of AMBER suite (27). The coulombic electrostatic potentials of both proteins at physiological ionic strength were calculated with UCFS Chimera program from the average structures after the MD simulation. Graphic representations and further analysis were performed with OriginPro 2019b and molecular graphics with UCFS Chimera.

Acknowledgements

We thank all the volunteers who generously participated in this study.

Conflict of Interest statement. C.S. and P.A. are currently the employee of Regeneron Pharmaceuticals and the beneficiary of stock options and grants in Regeneron. The other authors declare no competing interests.

Funding

Intramural Research Program of the National Institute on Aging (N01-AG-1-2109 and HHSN271201100005C); National Institutes of Health (NIH); by research grants from the Ministry of Science and Innovation (PGC2018-096049-B-I00); European Regional Development Fund (FEDER); Andalusian Government (BIO-198, US-1254317, US-1257019, P18-FR-3487 and P18HO-4091, US/JUNTA/FEDER, UE), University of Seville (VI PPIT) and the Ramón Areces Foundation. G.P.-M. was awarded a PhD fellowship from the Spanish Ministry of Education, Culture and Sport (FPU17/04604).

Author notes

Michele Marongiu and Gonzalo Pérez-Mejías authors contributed equally to this work.