Sex differences in viral entry protein expression and host transcript responses to SARS-CoV-2

Epidemiological studies suggest that men exhibit a higher mortality rate to COVID-19 than women, yet the underlying biology is largely unknown. Here, we seek to delineate sex differences in the gene expression of viral entry proteins ACE2 and TMPRSS2, and host transcriptional responses to SARS-CoV-2 through large-scale analysis of genomic and clinical data. We first compiled 220,000 human gene expression profiles from three databases and completed the meta-information through machine learning and manual annotation. Large scale analysis of these profiles indicated that male samples show higher expression levels of ACE2 and TMPRSS2 than female samples, especially in the older group (>60 years) and in the kidney. Subsequent analysis of 6,031 COVID-19 patients at Mount Sinai Health System revealed that men have significantly higher creatinine levels, an indicator of impaired kidney function. Further analysis of 782 COVID-19 patient gene expression profiles taken from upper airway and blood suggested men and women present distinct expression changes. Computational deconvolution analysis of these profiles revealed male COVID-19 patients have enriched kidney-specific mesangial cells in blood compared to healthy patients. Together, this study suggests biological differences in the kidney between sexes may contribute to sex disparity in COVID-19.

261,530 total deaths (as of December 2020). Such sex disparity exists in all age groups 5 . Sex-specific clinical characteristics analysis of patients in Wuhan, China showed that men with comorbidities presented a higher risk of being critically ill than men without comorbidities 6 , and a similar analysis of patients in Urban New York found that chronic kidney disease is higher in deceased men and obesity is higher in deceased women 3 . The compelling epidemiological evidence of sex disparity in COVID-19 and the urgency of this pandemic propelled the investigation of sex steroid hormone drugs in clinical trials (Estradiol: NCT04359329, Progesterone: NCT04365127, Degarelix: NCT04397718), although more preclinical evidence is needed. To underpin the biology of sex differences, hypotheses pertaining to the expression of viral entry proteins and immune systems have been actively explored 7,8 .
SARS-CoV-2 engages the receptor ACE2 (angiotensin-converting enzyme 2) for entry into the target cell through its spike protein 9 . Its internalization requires priming of the spike protein by the cellular protease TMPRSS2 (transmembrane protease, serine 2) in the host cell 10 , thus co-expression of ACE2 and TMPRSS2 on the target cell surface is required for virus entry. The high mortality in patients with COVID-19 may be partially driven by the strong affinity of the virus to ACE2 and the facilitation from TMPRSS2. For example, the higher expression of ACE2 in the nasal epithelium of children compared to adults may be linked to a lower case rate in children 11 . Plasma concentrations of ACE2 were observed to be higher in men than in women in cohorts of patients with heart failure 12 and higher expression of ACE2 and TMPRSS2 was also observed in men with asthma 13 . To perform tissue-wide expression analyses of ACE2 and/or TMPRSS2 between sexes, a few studies [14][15][16] exploited bulk or single-cell RNA-Seq samples primarily profiled from healthy individuals. However, many of the hospitalized COVID-19 patients have an underlying illness which increases mortality risks, and are in older age range (>60 years) 17 . Furthermore, the number of patients in these studies [14][15][16] was relatively small.
Sex differences in immune responses in COVID-19 were examined in a recent study, where blood samples from female patients were found to present more robust T-cell activation than male patients during SARS-CoV-2 infection 18 . Transcriptomic comparison of samples from nasopharyngeal swabs revealed that male patients had reduced B cell-specific and NK cell-specific transcripts and an increase in inhibition of nuclear factor kappa-B signaling 19 . Neutralizing autoantibodies that neutralize the ability of the corresponding type I IFNs to block SARS-CoV-2 infection were also found to express more in old men 20 . In addition to the postulation of immune responses, the genomic profiles produced in these studies and others offer a great opportunity to perform a holistic examination of the difference of cellular responses to the virus.
In this study, we repurpose existing large gene expression profiles to examine sex differences before and after SARS-CoV-2 infection. First, we leverage public gene expression profiles covering a wide range of ages, tissues, and disease conditions, to quantify expression differences of virus entry proteins, and later utilize Electronic Medical Records (EMR) to validate findings. Next, we harness the emerging COVID-19 patient gene expression profiles to characterize transcriptional response differences between sexes in the upper airway and blood.

Expression of ACE2 and TMPRSS2 in a diverse and comprehensive set of human samples
We first compiled three large independent expression datasets consisting of 220,835 samples from diverse tissue types and patient populations (healthy and disease conditions) and completed their meta-information, including sex, age group (younger: 0-19, middle: 20-59, and older: >60), and tissue of origin (14 main tissues), through machine learning and manual annotation ( Figure 1). To minimize batch effects, all the samples in each dataset were profiled under the same platform and processed using the same pipeline. The first dataset was compiled from the Treehouse project (T), where 17,654 RNA-Seq samples primarily from consortium projects including TCGA, GTEx, and TARGET were processed through the Toil pipeline 21  Leveraging their expression features, we built machine learning models (deep multi-task neural network and XGBoost) that completed metadata for the majority of samples with high confidence (Table S1). The average accuracy of the predictions in A and G are 93.6% (Sex), 80.8% (Age group), and 87.3% (Tissue) ( Table S2).
All of these predictions were further manually inspected based on unstructured sample meta-information available in the source files when possible. In total, we ended with 220,835 samples for the following analysis ( Figure 1B, detailed distribution for individual datasets in Figure S1).   Barrett's esophagus, trachoma, and ichthyosis have overall higher ACE2 expression in disease samples compared to control (Student's t-test, P < 0.05, Figure S3). The small sample size for each disease obstructed the subsequent sex difference analysis. In summary, although expression difference of entry proteins between sexes was not observed in the Healthy dataset, higher ACE2 expression was found in men, especially in older men, in the Merged dataset.   (Table S3). Notably, the kidney is the only tissue showing a remarkable difference in ACE2 expression between sexes in both A and G ( Figure 2B). After adjusting for age and data source, the OR is 1.45 [1.26-1.67] (P < 0.001) ( Table S3). The kidney is also the only tissue showing a significant difference in TMPRSS2 expression between sexes in both A and G ( Figure S2B).
We were able to further map 28% of those samples with high ACE2 expression to their diseases based on sample meta-information. The top mapped diseases are clear cell/renal cell carcinoma (60.8%), renal interstitial fibrosis (9.1%), acute kidney injury (7.5%), glomerulosclerosis (6.7%), nephritis   23 . Wilcoxon rank-sum test was used to test the difference and GPL570 data were employed for both figures.
Besides, steroid hormone receptors regulate the renin-angiotensin-aldosterone-system, where ACE2 is an essential component 24  0.26, P < 0.001) in the kidney ( Figure S4A). The genes regulated by AR also highly overlap with the genes positively co-expressed with ACE2 in the healthy adrenal gland (P = 3.08E-5, Figure S4B), suggesting that ACE2 expression might be associated with androgen receptor activity. Androgen signaling was identified to be a key modulator of ACE2 levels and treatment with antiandrogenic drugs reduced ACE2 expression and protected hESC-derived lung organoids against SARS-CoV-2 infection 25  To find clinical evidence of these findings, we analyzed 6,031 COVID-19 patients (4,621 inpatients and 1,410 outpatients with available labs) for serum creatinine levels measured in five member hospitals at Mount Sinai Health System up to May 10, 2020. We observed that men have significantly higher serum creatinine levels than women after normalizing to sex-specific reference ranges and adjusting for age and race (Inpatients OR: and clinical data analysis suggest that sex difference in the kidney is not specific to the older group (Table   S3 and Extended Data). Recent studies reported that acute kidney injury is common in patients hospitalized with COVID-19 and is associated with increased mortality 26,27 . Our accompanying study indicates chronic kidney disease is more common in men than in women who died from COVID-19 3 . Together, the expression difference of entry proteins in the kidney between sexes might be a factor contributing to sex differences in COVID-19 susceptibility.

Sex stratified analysis of host responses to SARS-CoV-2
After contracting the virus, individuals may have varying responses to the infection. We then sought to dissect the differences through leveraging the emerging COVID-19 patient profiles. We searched GEO and SRA to obtain COVID-19 patient RNA-Seq samples and reprocessed raw sequence data when possible. We compiled four datasets with gender information (one from upper airway nasopharyngeal, one from upper airway naso/oro-pharyngeal, one from blood PBMC, and one from blood leukocytes), totaling 782 samples (Table S5). In each dataset, the ratio of the number of samples between sexes is close to 1. For each large upper airway dataset, we stratified samples into an older age group (>60 years) and a middle age group (20-60 years). We enumerated all the possible comparisons for each dataset (i.e., female control vs. female patient, male control vs. male patient, female control vs. male control, and female patient vs. male patient), with each comparison using the same thresholds to select differentially expressed genes.
In comparing female and male samples either in the control group or in the COVID-19 group, only a few sexspecific genes were dysregulated between women and men. Neither ACE2 nor TMPRSS2 was differentially expressed between sexes in either group. As SARS-CoV-2 likely infects nasal epithelial cells to enter into the body 28 Figure 5B). The enrichment of mesangial cells was detected in a substantial number of ICU COVID-19 blood, although the absolute value is low ( Figure 5C, 5D).
Together with the higher expression of ACE2 and higher creatine levels in men, this analysis implies that impaired kidney function could be one source of sex differences. Patients were stratified into two groups (<= 60 and > 60) whenever age information was available and sample size in each group is greater than 10. The percentage of common DE genes was computed as the ratio between the number of common DE genes and the total number of distinct genes (Male_CT_vs_SARS2 plus Female_CT_vs_SARS2).

Discussion
COVID-19-related death is mainly associated with being male, older age, and comorbidities 30 . The virus first enters the nose and throat, and then travels down to attack lung, which likely causes substantial respiratory pathology including acute respiratory distress syndrome. Its reach can extend to many other organs like blood vessels, liver, kidney, heart, and brain 31 . ACE2 and its partner proteins are the key facilitators of virus entry into different organs. Therefore, their expression levels could be associated with disease severity and further mortality, and their expression differences between sexes could partially explain the higher mortality in men.
Previous efforts did not detect the difference 14,15 , likely because only healthy tissues were examined. In our Healthy dataset, we did not observe the difference either. Therefore, a comprehensive dataset covering a wide range of tissue samples including diseased samples is necessary to solve the puzzle.   An additional difference resulted from host responses after the infection of SARS-CoV-2. In addition to reported immune response differences, we found vast differences in non-immune cells such as mitochondria functions, phagocytosis, and cholesterol biosynthesis, suggesting men and women have unique response trajectories after infection and latency. Both large-scale genome-wide CRISPR screenings confirmed that cholesterol synthesis is a key host response factor and targeting cholesterol synthesis using small molecules like 27-hydroxycholesterol could reduce SARS-CoV-2 infection in vitro [32][33][34] . Of note, independent of sex, younger people (<=60) have much more differentially expressed genes than older people (>60); independent of age, women tend to have more differentially expressed genes than men. In two independent upper airway datasets, the subtle difference of host response between older male COVID-19 patients and older male healthy individuals might suggest the attenuated transcriptional responses to viral infection with increasing age, although the causal relation between host responses and disease susceptibility needs further exploration.
Our work has the following limitations. First, in the quantitative analysis of ACE2 differences, the batch effect between studies, variations among the three datasets, and other unknown confounding effects do exist; however, the main conclusions are robust, as we compared the samples that were produced using the same platform and processed under the same pipeline, used three relatively large independent datasets, and analyzed both categorical and continuous data. Such remedies may not adequately address the following confounding effects. A non-functional isoform of ACE2 (dACE2) has recently been discovered to be upregulated in the lung, gastrointestinal, and urinary tract in response to interferon activity and viral infections 35 . dACE2 might be absent or weakly expressed in normal tissues, but could be induced by the inflammatory tissue microenvironment. Although the survey of TCGA tumors confirmed that ACE2, not dACE2, was most expressed in kidney tumors 35 , it remains possible the expression of dACE2 in various tissues may confound the analysis. The analysis could be also confounded by the disease state. The annotation of disease states for all the samples is cumbersome, if not impossible, given the urgency of this project. Third, when studying cell composition of bulk COVID-19 patient samples, although the computational tool Xcell we employed has been widely used, the exact cell fractions could be quantified more precisely using single-cell technologies.
Furthermore, because of heterogeneous cell compositions in the samples, the DE genes might primarily reflect the difference of various cell types rather than that within infected cells. Lastly, as more and more COVID-19 patient clinical data are accumulated, it would be of great interest to further associate creatinine levels with more comorbidities as well as assess efficacy of steroid sex hormones based on real-world evidence. Nevertheless, our analysis of extensive genomic, and clinical data quantified sex differences before and after virus infection.

Data collection
Treehouse OCTAD (T): We downloaded the processed TPM data and phenotype data from the Treehouse project 21 . We used the same pipeline TOIL to process additional samples and integrated them into a single dataset OCTAD 36 . This dataset has been used for drug repositioning. In this work, we took the subset of OCTAD where samples have tissue, sex, and age information. The subset includes samples from healthy normal, primary cancer, and adjacent normal tissue samples. This dataset has complete information of sex and age, as well as a fairly complete annotation of tissues. In Figure 1, we only analyzed healthy normal Each batch was then normalized with the Affy package using RMA. Selected batches were normalized with justRMA to maintain large batch size. Median was used to merge expression of multiple probes. We included all the samples profiled under this platform. The dataset covers human tissue sample, cultured human cell samples, and those treated with various perturbagens. In order to perform unbiased analysis, we did not exclude cell line samples frequently used at the bench. Sample metadata such as title, source name, and characteristics were pulled out from GEOmetadb.

Meta-data curation and imputation
For Treehouse OCTAD (T), since sex and age information were relatively complete, we used only labeled data to perform the analysis, which also served as a reference for other datasets. For ARCHS4 GPL11154 (A) and GEO GPL570 (G), only less than 1/3 samples were labeled (Table S1). For samples taken from sexspecific tissues such as prostate and testis, we manually imputed their sex as male. However, the missing rate was still high after this process. Therefore, we further utilized state-of-the-art machine learning models to impute missing labels. Given a dataset, for each target of interest, i.e., sex, age, or tissue, we extracted the gene expression of target-specific (enriched) genes 37,38 for all the samples and used labeled data to build a prediction model. For sex and tissue, we built a multi-task deep neural network 39 ( Figure S6) to predict them together since sex and some tissues are dependent on each other (e.g., only males have prostate/testis samples while breast samples are more likely from females; building models independently for sex and tissue may result in futile predictions such as female-testis). For age prediction, when adding age prediction into a multi-task learning framework, sex and tissue tasks were dominant over age prediction which yielded unsatisfactory results. Therefore, we treated age prediction as a single task and used the extreme gradient boosting (XGBoost) 40 to build the model. All the hyper-parameters were tuned using 5-fold cross-validation and the detailed experiment settings can be found in Supplementary Text. After the classifiers were built, we applied it to unlabeled data and predicted their labels. Only predictions with high confidence were kept for later use. For sex, we considered prediction with probability of <0.1 (female) and >0.9 (male) as confident.
After that, human inspections were also involved to correct the potentially mislabeled samples based on the characteristics of each sample from the raw source file. Age prediction using gene expression data remains challenging, especially with the context of disease states, so we only predicted three age groups (0-19, 20-59, and >60). We observed that ages 40s and 50s are more likely to be misclassified. In the analysis, we only chose highly confident predictions and focused on the comparison of the younger and older group. In summary, we explored multiple options to ensure that each step produces optimal results. The complete labeled percentage for sex, age, and tissue before and after imputation can be found in Table S1. The model accuracy for all the tasks across all datasets can be found in Table S2, and the visualizations for learned representations are shown in Figure S7.

Clinical data analysis
Using the Mount Sinai Data Warehouse, we compiled de-identified electronic medical records (EMR) data including age, sex, race, and creatinine levels for inpatients and outpatients confirmed with COVID-19 at Mount Sinai Health System up to 05/10/2020. We normalized the creatinine levels by using x/1.2 for males and x/1.1 for females (x was the creatinine level of the corresponding patient) since the reference range for males is 0.6~1.2 mg/dL and for females is 0.5~1.1 mg/dL. For logistic regression, we set the normalized creatinine > 1 as 1 for both in-and outpatients. In logistic regression, we set sex-female and race-white as the reference. Patient summary statistics can be found in Extended Data and more characteristics of the patients were reported elsewhere 41,42 Statistics and data analysis For the comparison of sex differences in the proportion of samples with high expression of virus entry proteins, we employed logistic regression with the adjustment of age, tissue, and data source whenever possible.
Since smoking status is emerging as an inconclusive factor accounting for sex differences in COVID-19, we adjusted smoking status in a small set (non-smoker: 859, smoker: 1,542), and found men still have higher ACE2 expression than women (data not shown). Because of the limited size, we did not include smoking status in the main analysis. In tissue analysis, we chose 14 main tissues based on their sample counts and their importance; unknown refers to samples either belonging to other tissues or with low-confidence tissue prediction.
The expression of virus entry proteins is not normally distributed across samples: a majority of them with undetectable expression (TPM < 0.5), and a considerable number of samples with high expression of these entry proteins. We chose the top 10% most highly expressed samples as the high group (label 1) and the remaining as the normal group (label 0). We also explored the thresholds 75%, 80%, 85%, and 95% (the TPM is less than 0.5 at the threshold of 75%). There are variations of ORs among these thresholds, but the conclusions did not change ( Figure S8). Such grouping enabled the incorporation of three datasets, the inclusion of samples with undetectable expression, and the mitigation of batch effect between studies. In addition, we applied a similar analysis to individual studies (i.e., GEO GSEs) (Supplementary Text).
The sex difference in continuous gene expression of target genes was evaluated by two metrics: the ratio of average gene expression between female and male groups and the difference of median gene expression between female and male groups. The 95% CI of ratio estimate between female and male groups is obtained by bootstrapping: sampling with replacement and calculating the ratio of average gene expression between female and male for 1000 times, with 2.5% and 97.5% quantile recorded. The difference of median gene expression between women and men is computed using a two-sided Wilcoxon rank-sum test with 95% CI and p-value reported. Since low abundance genes in RNA-Seq samples (TPM < 0.5) are often noisy, these samples were removed accordingly.

COVID-19 patient sample processing
We searched GEO and SRA using keywords "COVID" and "SARS-CoV-2" and chose those datasets with >3 samples in each sex group (e.g., male healthy, male COVID). We ended with four datasets in the following analysis (Table S4). Two datasets were from upper airways, one dataset was from PMBC and one was from leucocytes. The dataset from leucocytes includes patient severity (ICU vs. non-ICU). Reprocessed raw sequence data available at SRA were downloaded and mapped to the human Hg38 transcriptome using the ENSEMBL GRCh38.p3 annotation using STAR aligner (https://github.com/alexdobin/STAR). The read count mapped on transcriptome was used for DE analysis. All possible comparisons between male and female samples (female-CT vs. male-CT, female-SARS-Cov2 vs. male-SARS-cov2, female-CT vs. female-SARS-Cov2, and male-CT vs. male-SARS-Cov20) were performed. The absolute log2foldchange >=1 with a false discovery rate < 0.01 computed from the EdgeR 43 package was chosen to identify DE genes. The DE genes of these comparisons were further compared via vennDiagram R package to find out common genes and specific genes for each comparison. Furthermore, DE genes specific to Control vs. SARS-cov2 in male and Control vs. SARS-cov2 in female were applied to identify enriched pathways through the enrichR API 44 . The sequence alignment, DE computation, and pathway enrichment were implemented in the OCTAD package 36 . XCell was employed to infer cell composition 29 .  ACE2 expression across (A) age groups and (B) tissues in three datasets Treehouse (T), ARCHS4 (A), and GPL570 (G). In T, only normal and adjacent normal tissue samples were selected. Due to the limited sample size (n), the 0-19 age group, and the blood and bone tissues were not included in (A) and (B), respectively. M/F represents the ratio of mean expression between males and females. The 95% con dence interval of ratio was obtained using bootstrapping, i.e., sampling with replacement and calculating the ratio for 1000 times with 2.5% and 97.5% quantile reported. P-value was computed from a two-sided Wilcoxon rank-sum test on the difference of median expression between sexes. * indicates P < 0.001. The bar and the dashed line show the mean of expressions in each group and the mean of all expressions, respectively. Wilcoxon rank-sum test was used to test the difference and GPL570 data were employed for both gures.