ACE2 and TMPRSS2 variants and expression as candidates to sex and country differences in COVID-19 severity in Italy

ACE2

For most X-chromosome genes, the double allelic dosage in females is balanced by the epigenetic silencing of one of the X chromosomes in early development [16]. However, the X-chromosome inactivation (XCI) is incomplete in humans and up to one third of genes are expressed from both alleles, with the degree of XCI escape varying between genes and individuals [17]. ACE2 is one of the genes escaping X inactivation, but it belongs to a subgroup of X-chromosome genes escaping XCI showing an uncharacteristically heterogeneous pattern of male-female expression, with higher expression in males in several tissues [13]. Specifically concerning the lung, a recent analysis on published expression data, reported a substantial similar level of ACE2 transcript in males and females [14], however, another study, using single-cell sequencing, found a higher expression of ACE2 in Asian males [15]. Figure 1 reports data on ACE2 mRNA expression levels in the lung as retrieved from the largest datasets available in the literature; no substantial differences were found between males and females, nor between younger and older females, thus confirming what already observed by Cai and colleagues [14].

Another possible sex-related effect might be due to the fact that males are hemizygous for the gene, therefore, in the presence of an ACE2 allelic variant increasing disease susceptibility or severity, males will have all cells expressing the risk variant. Based on this hypothesis, we looked into the genetic variation in ACE2. A recent manuscript explored this same topic in different populations using data from public databases [18]. However, a specific analysis of the Italian population is Iacking.

We have therefore exploited the available data on 3,984 exomes obtained from an Italian cohort representative of the whole country [19, 20] to extract the variants in exons and splice junctions of ACE2. Variants were filtered for quality and classified according to their predicted effect at protein level and on splicing. Concerning rare variants (i.e. those with a minor allele frequency, MAF, <1%; to be used in burden tests), we considered only null variants, abolishing or significantly impairing protein production (nonsense, out-of-frame ins/dels, and splicing variants), and missense variants predicted to be deleterious or possibly deleterious by all the 5 prediction algorithms used (see Supplementary Methods, paragraph “Definition of disrupting variants and statistical analysis”). Concerning common variants (i.e., MAF>5%), all were retained for comparing their frequency with those of the European (non Finnish) and East Asian populations, retrieved from the GnomAD repository.

No significant differences in the burden of rare deleterious variants were observed comparing the Italian population with Europeans and East Asians (Table 1A). Concerning common exonic variants, the only striking difference, as also noticed by Cao and colleagues [18], was observed for the single nucleotide polymorphism (SNP) rs2285666 (also called G8790A), with the frequency of the rare A allele being 0.2 in Italians and Europeans, and 0.55 in East Asians (uncorrected P=2.2*10-16 for difference in Italians vs East Asians; corrected P=7.9*10-15; Table 1B). This variant was extensively studied as a potential risk factor for hypertension, type 2 diabetes, and coronary artery disease [21, 22], hence possibly constituting a predisposing factor also for the comorbidities observed in COVID-19 patients. A single paper reports the association of the three rs2285666 genotypes with ACE2 protein level measured in serum by ELISA, with the A/A genotype having an expression level almost 50% higher than the G/G genotype, while heterozygous G/A individuals had intermediate levels [23]. Given the position of the variant, at nucleotide +4 in the donor splice site of intron 3 (c.439+4G>A), we calculated the predicted effect on splicing and indeed the substitution of G with an A is predicted to increase the strength of the splice site of about 9.2% (calculation made through the Human Splicing Finder v.3.1 webtool, http://www.umd.be/HSF/), consistently with the higher level of ACE2 protein in serum. It would be crucial to compare the frequency of this variant with ACE2 expression in the lung and with susceptibility to viral infection and severity of COVID-19 manifestations. Of note, no eQTL for ACE2 in the lung has been described so far in the GTEx database, and investigations on this topic are recommended.

TMPRSS2

TMPRSS2 is a gene well known to oncologists as genetic rearrangements producing a fusion between TMPRSS2 and ERG (or, more rarely, other members of the ETS family) are the most frequent genetic lesions in prostate cancer patients [24]. As TMPRSS2 is an androgen responsive gene, the fusion results in androgen dependent transcription of ERG in prostate tumor cells. Therefore, we can hypothesize that males might have higher TMPRSS2 expression also in the lung, which might improve the ability of SARS-CoV-2 to enter cells by promoting membrane fusion. Looking into GTEx and GEO data, the overall expression of TMPRSS2 in the lung is only slightly increased in males (P=0.029; Figure 2A). However, TMPRSS2 expression is also promoted by estrogens [12], and therefore the situation might be different when considering individuals above 60 years, who are at higher risk of fatal events due to COVID-19, as in this group females will all be postmenopausal. According to this hypothesis, we checked the expression of the gene in lungs of males and females at different ages, but no substantial differences emerged between males and females (neither below, nor above 60 years of age; data not shown).

Finally, we explored genetic variation in TMPRSS2 in search of variants, possibly already annotated as eQTL in the lung, which might have an impact on the serine protease expression as well as on its catalytic activity. Again, we used the available Italian exome data, as well as data deposited in GnomAD [25].

Firstly, we looked at the overall burden of deleterious rare variants, using the variant classification described above. Italians had a nominally significant decrease in the burden of deleterious variants compared to Europeans (uncorrected P=0.039, not significant after correction for multiple testing; Table 2A). This decrease was even more evident for the East Asian population (corrected P=8.6*10-4); however, in this case, we must consider that the number of individuals over 65 years of age in Italy is more than double the one in the Hubei province (22.7 vs. 10%, respectively) and this is a major determinant of disease lethality.

Focusing specifically on common exonic variants, 4 SNPs showed significantly (P<2.2*10-16) different frequencies when comparing the Italian population with East Asians (and with Europeans) (Table 2B); 3 of them are synonymous variants, whereas one is the missense substitution p.Val160Met, which impacts on a residue far from the serine protease catalytic triad. This variant was previously found significantly associated with genomic rearrangements involving TMPRSS2, with the risk of prostate cancer [26] and with shorter time to prostate cancer diagnosis for high-risk patients [27].

Concerning eQTLs, a number of variants significantly impacting on TMPRSS2 expression in the lung (GTex data) are reported in the 3’ region of the gene (Figure 2B). In Table 2C, a list of the most significant (P<1*10-8), together with their GnomAD frequencies in the East Asian and European populations, are reported. As for the Italian frequencies, we took advantage of the genome-wide association study (GWAS) performed on the above-described cohort (for a total of 3,284 individuals) [28]; in this case, we had to infer genotype frequencies by an imputation approach (for details, see Supplementary Methods). Interestingly, all these eQTLs appear to have extremely different frequencies among populations. In particular, 2 different haplotypes can be inferred from frequency data:

1. A frequent “European” haplotype (composed at least of SNPs rs463727, rs34624090, rs55964536, rs734056, rs4290734, rs34783969, rs11702475, rs35899679, and rs35041537), which is totally absent in the Asian population. Interestingly, this haplotype has been functionally linked to another eQTL (rs8134378), located at a known androgen-responsive enhancer for TMPRSS2, 13 kb upstream of the TMPRSS2 transcription start site [29] (Figure 2B). Hence, this haplotype is expected to up regulate TMPRSS2 gene expression in an androgen-specific way.

2. A second haplotype, predicted to be associated with higher TMPRSS2 expression, is characterized by 3 SNPs (rs2070788, rs9974589, rs7364083), whose MAF is significantly increased in Europeans (9% increase in Italians respect to East Asians, corrected P<6.8*10^-9). Importantly, a small-scale GWAS, comparing the distribution of genetic variants in severe and mild cases of patients with A(H1N1)pdm09 influenza, identified rs2070788 as being associated with increased risk to both human A(H7N9) and severe A(H1N1)pdm09 influenza [11]. Of note, also in A(H7N9) influenza, the proportion of male patients was more than double that of female patients [30].

Gene expression data

Expression data for ACE and TMPRSS2 genes were obtained through the: 1) genotype-tissue expression (GTEx) database (https://gtexportal.org/home/), which was also used to extract quantitative trait loci (eQTLs) for the two genes (all data based on RNAseq experiments); and 2) Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/geo/). In particular, two GEO datasets were extracted and analyzed: 1) GSE66499, reporting microarray data on 152 normal lung samples from Caucasian individuals; 2) GSE19804, reporting microarray data on 60 normal lung samples from Taiwanese females (see also Supplementary Methods, paragraph “Datasets and statistical power estimations”).

Genetic data

Genetic data for general European and East Asian populations were retrieved through the GnomAD repository, which contains data on a total of 125,748 exomes and 71,702 genomes (https://gnomad.broadinstitute.org/).

As for Italians, details on whole-exome sequencing (on 3,984 individuals) and genome-wide microarray genotyping (on 3,284 individuals) of the analyzed cohort are specified elsewhere [19, 20, 28], as well as in Supplementary Methods (paragraphs “Sequencing” and “Datasets and statistical power estimations”). Imputation procedures are detailed in Supplementary materials (paragraph “Dataset imputation>”).

Statistical analysis

Expression levels were compared by using either the Kruskal-Wallis test (RNAseq data) or the student t test (microarray data). Allele frequencies were compared using the chi square test. All calculations were performed using the R software (https://www.r-project.org/). P values are presented as non-corrected for multiple testing, but the Bonferroni-corrected threshold of significance is indicated below each set of comparisons presented in Tables. Power calculations have been described in Supplementary Methods (paragraph “Datasets and statistical power estimations”).