ACE2 gene variants may underlie interindividual variability and susceptibility to COVID-19 in the Italian population

INTRODUCTION

In December 2019, a new infectious respiratory disease emerged in Wuhan, Hubei province, China (1–3). An initial cluster of infections likely due to animal-to-human transmission was rapidly followed by a human-to-human transmission (4). The disease was recognized to be caused by a novel coronavirus (SARS-CoV-2) and termed coronavirus disease 19 (COVID-19). The infection spread within China and all over the world, and it has been declared as pandemic by the World Health Organization (WHO) on 2^nd March 2020. The symptoms of COVID-19 range from fever, dry cough, fatigue, congestion, sore throat and diarrhea to severe interstitial bilateral pneumonia with a ground–glass image at the CT scan. While recent studies provide evidence of a high number of asymptomatic or paucisymptomatic patients who represent the main reservoir for the infection progression, the severe cases can rapidly evolve towards a respiratory distress syndrome which can be lethal (5). Although age and comorbidity have been described as the main determinants of disease progression towards severe respiratory distress, the high variation in clinical severity among middle-age adults and children would likely suggest a strong role of the host genetic asset.

A high sequence homology has been shown between SARS-associated coronavirus (SARS-CoV) and SARS-CoV-2 (6). Recent studies modelled the spike protein to identify the receptor for SARS-CoV-2 and indicated that angiotensin converting enzyme 2 (ACE2) is the receptor for this novel coronavirus(7,8). Zhou et. al. conducted virus infectivity studies and showed that ACE2 is essential for SARS-CoV-2 to enter HeLa cells (9). Although the binding strength between SARS-CoV-2 and ACE2 is weaker than that between SARS-CoV and ACE2, it is considered as much high as threshold necessary for virus infection. The spike glycoprotein (S-protein), a trimeric glycoprotein in the virion surface (giving the name of crown -corona in latin-), mediates receptor recognition throughout its receptor binding domain (RBD) and membrane fusion (10,11). Based on recent reports, SARS-CoV-2 protein binds to ACE2 through Leu455, Phe486, Gln493, Ala501 and Tyr505. It has been postulated that residues 31, 41, 82, 353, 355, and 357 of the ACE2 receptor map to the surface of the protein interacting with SARS-CoV-2 spike protein (12), as previously documented for SARS-CoV. Following interaction, cleavage of the C-terminal segment of ACE2 by proteases, such as transmembrane protease serine 2 (TMPRSS2), enhances the spike protein-driven viral entry (13,14). Thus, it is possible, in principle, that genetic variability of the ACE2 receptor is one of the elements modulating virion intake and thus disease severity. ACE2 is located on chromosome X. Although it is one of the genes escaping X inactivation several lines of evidence suggest that a different degree of X-chromosome inactivation is present in distinct tissues (15).

Taking advantage of the Network of Italian Genomes (NIG), a consortium established to generate a public database (NIG-db) containing aggregate variant frequencies data for the Italian population (http://www.nig.cineca.it/), here we describe the genetic variation of ACE2 in the Italian population, one of the newly affected countries by the SARS-CoV-2 outbreak causing COVID-19. Three common c.2158A>G p.(Asn720Asp), c.77A>G p.(Lys26Arg) and c.631G>A p.(Gly211Arg) variants and 27 rare missense variants were identified, 9 of which had not previously been reported in public databases. We show that p.(Asn720Asp), which lies in a residue located close to the cleavage sequence of TMPRSS2, likely affects the cleavage-dependent virion intake. Along with the other two common variants, this substitution is represented in the Italian and European populations but is extremely rare in the Asian population. We also show that two rare variants, namely, c.1051C>G p.(Leu351Val) and c.1166C>A p.(Pro389His) are predicted to cause conformational changes impacting RBD interaction. As the uncertainty regarding the transmissibility and severity of disease rise, we believe that a deeper characterisation of the host genetics and functional characterization of variants may help not only in understanding the pathophysiology of the disease but also in envisaging risk assessment.

RESULTS

ACE2 variants identification

The extent of variability along the entire ACE2 coding sequence and flanking intronic stretches was assessed using 6930 Italian WES, out of which 4171 males and 2759 females which sum up to 9689 alleles. Identified variants and predicted effects on protein stability are summarized in Table 1, Table S1 and Table 2, and represented in Fig.1. Three more common variants, c.2158A>G p.(Asn720Asp), c.77A>G p.(Lys26Arg), c.631G>A p.(Gly211Arg), were identified. The c.2158A>G p.(Asn720Asp) substitution was estimated to have a frequency of 0.011 (103/9689 alleles), which is in line with the frequency of the variant reported in the gnomAD database (0.016), and is lower than the frequency reported in gnomAD for the European non-Finnish population (0.025, 2195/87966 analyzed alleles). Given the ACE2 localization on X chromosome we focused our attention on the females alleles. All analysed females (2759 out of 6930) belonging to the Italian population, were heterozygotes for the variant. Notably, this variant has not been reported in the Eastern Asia population (13,784 exomes). The c.77A>G p.(Lys26Arg), c.631G>A p.(Gly211Arg) variants were found with a frequency of 0.0011 (lower than the frequency in the European-non Finnish population, 0.0058) and 0.0012 (European-non Finnish population frequency, 0.0019), respectively. Out of ~92708 analyzed alleles in the European non-Finnish population, one homozygous female has been reported for the c.77A>G p.(Lys26Arg) while no homozygous females were reported for the c.631G>A p.(Gly211Arg). According to gnomAD, the allele frequency of the c.77A>G p.(Lys26Arg) variant in the Eastern Asia population was 6×10^-5, while the c.631G>A p.(Gly211Arg) has not been reported in 14.822 exomes. In addition to these variants, 28 rare missense variants were identified, out of which 10 had not previously been reported in GnomAD database and 9 truncating variants that had not been reported in gnomAD database (Table 1, Supplementary Table 1). Out of these variants, two fall in the neck domain, which is essential for dimerization and one in the intracellular domain. Many of them truncate the protein in different positions of the Protease domain embedded in the extracellular domain, which contains the receptor binding site for SARS-CoV-2. Only three truncating variants have been previously described for ACE2 likely due to a low-tolerance for loss-of function variants. In line with this evidence, all these variants were very rare and no homozygous females were detected for the identified variants. Three missense changes c.1517T>C p.(Val506Ala), c.626T>G p.(Val209Gly), c.1129G>T p.(Gly377Glu) were predicted to have destabilizing structural consequences (Table 2); among these, c.1517T>C p.(Val506Ala) is indeed the only amino acid change reported in the European non-Finnish population (rs775181355; allele frequency 1.40×10^-5, CADD 27,2) and is predicted as probably damaging for the protein structure by Polyphen and deleterious by SIFT. Similarly, c.1051C>G p.(Leu351Val) and c.1166C>A p.(Pro389His), which affect a highly hydrophobic core, were predicted to induce conformational changes influencing the interaction with spike protein. The amino acid substitution c.1166C>A p.(Pro389His) (rs762890235, European-non Finnish population allele frequency: 2.45 10^-5, CADD 24,8) was predicted to be probably damaging by Polyphen and deleterious by SIFT. Moreover, this rare variant has never been reported in Asian populations.

View this table:

• View inline
• View popup
• Download powerpoint

View this table:

• View inline
• View popup

Table 1. Missense, frameshift and truncating variants identified in ACE2 gene

The table reports the genomic position, the nucleotidic and protein change of exonic ACE2 identified variants. The genomic reference sequence is NM_021804.2 (hg19). CADD_phred scores are reported for the missense variants. When available, dbSNP rs number and the genome aggregation database gnomAD allele frequency are reported. For all variants are reported the number of individuals hemizygous, heterozygous or homozygous.

View this table:

• View inline
• View popup

Table 2. Predicted changes in ACE2 protein stability as consequence of residues changes

DUET program results that display predicted change in folding free energy upon ACE2 missense variant (∆∆G in kcal/mol). In the first three columns are reported single missense variants with specific position on ACE2 protein. The residues in the first column highlighted in grey are involved in N-glycosylation pattern NxT/S, therefore those missense variants determine the loss glycosylation of Asparagine 53 and 90 respectively. In the fourth column is reported G analysis predict effects of missense variants on protein stability using an integrated computational approach. The column “Interaction Network around (5 Å)” shows for each single missense variant the residues around 5Å. In this column, we highlight in green residues involved in spike SARS-CoV protein interaction, in yellow residues involved in Zinc coordination and finally in magenta residues of Asn involved in N-glycosylation. The last column defines the outcome of protein stability for each single missense variant. An increasing negative value for the ∆∆G is correlated with a higher destabilizing effect, while a positive value is associated with a variant predicted as stabilizing.

• Download figure
• Open in new tab

Figure 1. ACE2 Crystal structure with PDB ID 1R42.

Surface and cartoon representations of protein in grey. In blue stick are represented each single mutated positions, cartoon region in yellow represent segment between amino acid 30-41, cartoon in green represent segment between amino acid 353-357 and cartoon in red represent segment between amino acid 82-84 that are protein regions responsible of interaction with 2019-nCOv spike glycoprotein

ACE2 variants likely affect protein stability and SARS-CoV-2 binding

Molecular Dynamic analysis provides bona-fide simulations of protein structural changes caused by missense variants effects. Yet, its computationally expensive procedure led us to perform MD simulation for only a selection of representative candidate variants. Indeed, we selected the following five variants and corresponding effects: c.1517T>C p.(Val506Ala) which has the higher destabilizing effect, c.77A>G p.(Lys26Arg) and c.631G>A p.(Gly211Arg) with higher allele frequency along with c.1051C>G p.(Leu351Val) and the c.1166C>A p.(Pro389His) with a predicted effect on spike protein interaction. To analyse differences in protein structure between wild type and mutants, we performed 100ns MD simulations. The comparison was performed by Root Mean Square Deviation (RMSD) analysis. The global effects of the residue substitutions on flexibility and global correlated motion of ACE2 protein are represented in Figure 2 and the simulation is provided in Supplementary Video S1, S2, S3, S4, S5. While a similar trend for wild-type, c.77A>G p.(Lys26Arg) and c.1517T>C p.(Val506Ala) was observed with a steady course in the RMSD value, which stabilizes at an average of 0.2 nm, 0.25 nm and 0.3 nm respectively (Figure 3, panel A), the c.1166C>A p.(Pro389His) and c.1051C>G p.(Leu351Val) variants show a difference in comparison with the native protein with a gradual increase in RMSD value, which stabilizes at an average of 0.5 nm (Figure 3, panel A). Finally, the c.631G>A p.(Gly211Arg) shows a bigger difference with a higher increase in RMSD value, which stabilizes at an average of 0.6 nm (Figure 3, panel A). Structural analysis between WT and mutant c.1517T>C p.(Val506Ala) MD simulations showed that the c.1517T>C p.(Val506Ala) forms a hydrophobic centre together with Leu456, Leu503 and Phe516 with minimum differences in protein rearrangements when the residue is mutated in Ala as reported in Figure 2, Supplementary Video S1. The c.77A>G p.(Lys26Arg) is located at the N-terminus and the sidechain engages a hydrogen bond with Asn90 thus determining a minimal destabilizing predicted effect as shown in Table 2 and confirmed by RMSD analysis. The c.1166C>A p.(Pro389His) and the c.1051C>G p.(Leu351Val) variants, located in the region for the spike protein interaction, are characterized by a partially overlapping destabilizing effect. The c.1166C>A p.(Pro389His) variant sidechain being more bulky causes the shift of ACE2 (30-40) helix involved in spike protein interaction which being freer to move engages an interaction with Gln96 (Figure 2, Supplementary video S5). The c.1051C>G p.(Leu351Val) shorter sidechain is enable to interact with the hydrophobic core composed by Trp349 and Leu45 with a consequent rearrangement of the protein conformation. Finally, while c.631G>A p.(Gly211Arg) has theoretically a smaller destabilizing effect because of an external sidechain which is not involved in particular interaction network, as shown by MD simulation, it confers a wide flexibility to this region because the polar sidechain is able to engage different interactions with vicinal amino acid residues (Figure 2, Supplementary Video S2). During MD simulations, we have also investigated the surrounding region of ACE2 WT and previously selected variants by calculating change in Solvent Accessibility Surface Area (SASA). Differences in average SASA value would suggest for the native protein a wider surface exposed to solvent and subsequently a different ability to interact with spike SARS-CoV-2 in comparison with the studied variants (Figure 3, panel B).

• Download figure
• Open in new tab

Figure 2. ACE2 wild type and variants superimposed structures after 100ns MD simulation.

Cartoon representation of ACE2 wild type (orange) and variants (green) in blue sticks the wild type residues while in red the corresponding variants. In cyan and pink sticks residues interacting with each specific position.

• Download figure
• Open in new tab

Figure 3: Structure superimposition snapshot between wild type protein and variant proteins

Panel A. Root Mean Square Deviation (RMSD) trends for the backbone of ACE2 WT (Black line) and some selected variants (colored lines, see legend) during 100 ns of simulation. The molecular dynamics simulation shows a good stability for all systems with exception of G211R mutants. (RMSD) is a parameter used to define the stability of an element. Wild type shows a steady course in the RMSD value, stabilizing at an average of 0.2 nm, while, the G211R variant shows a gradual increase in RMSD value, stabilizing at an average of 0.6 nm. Panel B. SASA graphical representation of ACE2 WT (Black line) and ACE2 variants (colored lines, see legend).

Differences in ACE2 allelic variability in COVID-19 patients compared to controls

In order to shed light on the role of ACE2 variants on interindividual variability and susceptibility to COVID-19 in Italian population we performed WES analysis on a cohort of 131 patients and 258 controls who agreed in participating to the study (see Materials and Methods). Data analysis of ACE2 variants identified a different distribution of variants in controls compared to patients (Figure 4) with the c.2158A>G p.(Asn720Asp) variant being present in two hemizygous male patients (allele frequency 0,012) compared to 7 heterozygous female and 4 hemizygous male controls (allele frequency 0,028). A silent variant the c.2247G>A p.V749V, was also detected in 26 control individuals (allele frequency 0,069) compared to 5 COVID-19 patients (allele frequency 0,030). Although any single variant was not statistically significantly enriched in one cohort compared to the other, a cumulative analysis of the identified variants detected a statistically significant higher ACE2 allelic variability (P value<0,029) in the control group compared to the patient cohort.

• Download figure
• Open in new tab

Figure 4. Differences in ACE2 variants in COVID-19 patients compared to controls

The figure shows the variants located in the ACE2 protein domains. The variants present in controls are shown in black while the variants in cases are shown in red. The number of patients carrying the variant is shown in brackets.

DISCUSSION

According to recent reports, ACE2 is essential for SARS-CoV-2 to enter cells. Recent single-cell RNA studies have also shown that ACE2 is expressed in human lung cells (30). The majority of ACE2-expressing cells are alveolar type 2 (AT2) cells. Other ACE2 expressing cells include alveolar type 1 (AT1) cells, airway epithelial cells, fibroblasts, endothelial cells, and macrophages although their ACE2-expressing cell ratio is low and variable among individuals. The expression and distribution of the ACE2 receptor can thus justify the route of infection and the main localization at the alveolar level. Although the different density of ACE2 receptors in the upper respiratory tract among individuals can partially give reason for the clinical variability, which ranges from asymptomatic/paucisymptomatic patients to severely affected ones, it could not be the only reason for such variability. In addition, recent works did not observe significantly different viral loads in nasal swabs between symptomatic and asymptomatic patients (31). Italy has been the first European country that experienced the COVID-19 outbreak with a rapid increase in the positive cases in a very short-time period and a morbidity and lethality (~10%) definitely higher in comparison with Asian countries, such as China (4%) and South Korea (1.2%) (How deadly is COVID-19? A rigorous analysis of excess mortality and age-dependent fatality rates in Italy (32). These considerations raise the possibility of a predisposing genetic background accounting for or contributing to the wide inter-individual clinical variability, as well for the differential morbidity and lethality observed among different countries, population awareness and constrictive measures apart.

We integrated genomic WES data produced by five Italian centers (Siena, Naple, Turin, Bologna, Rome) interconnected by the Network for Italian Genome (NIG) in the attempt to identify variation encompassing the ACE2 gene, which could account for a difference in SARS-CoV-2 spike binding affinity, processing or internalization. Previous studies showed that the residues near lysine 31, and tyrosine 41, 82-84, and 353-357 in human ACE2 are important for the binding of S-protein to coronavirus (12). In line with previous reports (33), we did not find polymorphism or rare variants in these residues in the Italian population. However, we identified three variants namely c.2158A>G p.(Asn720Asp), c.1166C>A p.(Pro389His) and c.1051C>G p.(Leu351Val), one of which polymorphic c.2158A>G (p.(Asn720Asp), moderately expressed in the Italian and European-non Finnish populations and with a very low allele frequency or not occurring in the Eastern Asia population. These variants which surround residual essentials for the SARS-CoV-2 spike protein binding were predicted to likely affect the cleavage dependent virion intake, such as the polymorphic c.2158A>G p.(Asn720Asp) (allele frequency 0.011) which lies 4 amino acids from the cleavage sequence of TMPRSS2 or to have a substantial impact on protein structure and spike protein interaction by MD simulation (Figure 3, Panel A). The relatively frequent c.631G>A p.(Gly211Arg) (allele frequency 0.0012, 12 /6930 individuals) was predicted to confer a wide flexibility to the region because of the ability to engage different interactions with the nearby amino acid residues. Along with these more common variants we also identified very rare variants such as the c.1166C>A p.(Pro389His) and the c.1051C>G p.(Leu351Val), some of which only described in the non-Finnish European population, that could give reason for a different affinity for the SARS-CoV-2 spike protein (Figure 2, Figure 3 panel A, Supplementary Video S4). Interestingly all the studied variants affect residues highly conserved among species (Supplementary Figure S1). Given their rarity in other populations, we cannot exclude that these variants can partially account for the clinical outcome observed in the Italian population. WES data generated from a wide cohort of COVID-19 Italian patients revealed a statistically significant (P<0,029) higher allelic heterogeneity for ACE2 in controls compared to patients with a higher chance to find at least one ACE2 variant in the cohort of controls compared to the cohort of patients. Therefore, it is plausible to think that the effect of allelic variability on ACE2 conformation would at least partially account for the interindividual clinical differences and likely modulate clinical severity. This finding reinforces the hypothesis that at least some of the identified variants or the cumulative effect of few of them confer a different susceptibility to virus cell entry and consequently to disease onset and progression. We cannot exclude that also silent variants such as the c.2247G>A (p.Val749Val) with no effect on the protein could play a role because of an unpredictable impact at a posttranscriptional level.

Notably, morbidity and lethality have been reported definitely higher in men compared to women (~70% vs 30%, 20th March 2020 ISS report). Although several parameters have been brought to case to explain this difference, i.e. smoking, differences in ACE2 localization and/or density in alveolar cells, hormonal asset, it is noteworthy that ACE2 is located on chromosome X and that given the low allele frequency of the identified variants the rate of homozygous women is extremely low (see Results section). The X-chromosome inactivation (XCI) is incomplete in humans and some genes show a degree of XCI escape which vary between individuals and tissues (34). ACE2 is one of the genes escaping X inactivation, but it belongs to a subgroup of X-chromosome genes showing a higher expression in men in several tissues thus mostly suggesting that ACE2 gene XCI is present although different in distinct tissues (15). Therefore, the impact of X-inactivation on the alternate expression of the two alleles would guarantee, in the affected tissues, a heterogeneous population of ACE2 molecules, some of which protective towards the infection until the point of a complete or almost complete protection in the case of a X-inactivation skewed towards the less SARS-CoV-2 -binding prone allele. This hypothesis would justify the high rate of asymptomatic or paucisymptomatic patients. However, the presented data does not allow to confirm a clear cause-effect relationship and, since most of the identified variants have very low frequencies, further functional studies are needed to validate these results. ACE2 is definitely one of the main molecules whose genetic heterogeneity can modulate infection and disease progression; however, a deeper characterisation of the host genetics and functional variants in other pathway-related genes may help in understanding the pathophysiology of the disease opening up the way to a stratified risk assessment and to tailored preventive measures and treatments.

Ethical Approval

The NIG study was approved by the Ethical Committee of each center, Prot Name NIG Prot n 10547_2016, approved CEAVSE on 18.07.2016 n 88/16. The COVID-19 patients study was approved by the University Hospital (Azienda Ospedaliera Universitaria Senese) Ethical Committee, Siena, Italy (Prot n. 16929, dated March 16, 2020).

Conflict of interest

The authors declare that they have no competing interests.

Author’s contribution

EB, RT, OS, AMP and AR have made substantial contributions to conception and design, acquisition of data, analysis and interpretation of data and have been involved in drafting the manuscript. RA, GB, AB, AB, AT, GD, AG, FM, TP, AT, AT and FV has made substantial contributions to acquisition and analysis of the data. MB, MB, AC, SF, FM, GM, VN, MS and MT have made substantial contributions to interpretation of data and clinical evaluation. All authors have been involved in drafting the manuscript; have given final approval of the version to be published and agree to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.