We evaluated differential expression by re-mining an aggregated dataset of 7,142 RNA-Seq samples 42 modified from the normalized and batch-corrected data from the GTEx and TCGA projects 43 . The Genotype-Tissue Expression (GTEx) project provides data representing "non-diseased" conditions from diverse tissues. The well-curated TCGA project is the largest project available with easily accessible metadata on the races of the individuals who contributed samples representing diseased tissues (tumors) of multiple origins. These large data provide a unique opportunity to evaluate differences in gene expression across populations in multiple tissues in diseased and normal states.
Re-mining existing RNA-Seq data and metadata has several caveats. Because race assignments are self-reported, many of the individuals sampled will be from admixed populations 41, 44 . We are terming those self-reporting as "Black or African American" as "African Americans" and "White" as "European Americans". The ancestry of the preponderance of African Americans is Western Africa 44 , thus our results for African Americans would mostly reflect more specifically Western Africans. Those self-reporting as White are presumed to be predominantly European Americans, but this group also would include individuals of other populations, including Indians and admixed Hispanic individuals, depending on how these individuals chose to self-report.
Also, we were limited to comparison of differences between gene expression in African American and European American populations, because even in these large studies, the sample numbers for the other three major population groups (Asian, Native American, and Pacific Islanders) 45 were too low for robust statistical assessment. Even between these two populations, not all conditions (cancers or "non-diseased" tissues) had sufficient samplings of African Americans for robust statistical assessment (Table 1; Supplementary Table S1 ).
Finally, although the GTEx project analyzes non-diseased tissues, many of the individuals who donated tissues were severely ill or postmortem from varied causes, which would likely effect gene expression in these tissues.
We analyzed the data and metadata with MetaOmGraph (MOG), software that supports interactive exploratory analysis of large, aggregated data 42 . Exploratory data analysis uses statistical and graphical methods to gain insight into data by revealing complex associations, patterns or anomalies within that data at different resolutions 46 .
Differentially expressed (DE) genes among samples from European Americans and African Americans were identified in a tumor-specific or tissue-specific manner using MOG (Mann-Whitney U test). Of the tumor types in the TCGA data, BRCA, COAD, KIRC, KIRP, LUAD, LUSC, THCA, and UCEC had sufficient numbers of samples for DE analysis (Table 1) . GTEx normal tissues analyzed were: breast, colon, esophagus, lung, liver, prostate, stomach, thyroid, and uterus (Table 1) . We define a gene as DE between two groups if it meets each of the following criteria:
1. Estimated fold-change in expression of 2-fold or more (log fold change, |logFC| ě 1 where logFC is calculated as in limma 47 .)
2. Mann--Whitney U test is significant between the two groups (BH corrected p-value ă 0.05)
Supplementary Tables S2-S25 contains the full list of DE genes between African American and European American populations for each condition. The numbers of DE genes vary depending on the condition the samples were obtained from. Many genes follow a similar DE trend in diseased as in non-diseased tissues, however, in each case, the fold-change difference of expression among the DE genes was larger in the cancers than in the corresponding non-diseased tissues. We investigate the distribution of the gene expression values with two additional statistical analyses. We performed the twosample Kolmogorov-Smirnov test (KS test) to assess whether there is a significant difference in distribution of gene expression between African Americans and European Americans. Hartigans' dip test was used to test whether a given distribution shows bi or multi-modality. Bi-or multi-modal distributions indicate there may be hidden or unknown covariates affecting the gene expression. Within a race, this could imply presence of sub-population structure. These analyses indicate the not only are genes differentially expressed between the two populations, but the overall distribution in expression values often differed between the populations, for an obvious example, GSTM1 expression in KIRC (Figure 4 . Also, in many cases, one or both populations have a bimodal distribution, for example, GSTM1 expression in BRCA has a bimdodal distribution in European Americans but not in African Americans ( Figure 4 ).
GO terms related to the biological processes of infection, inflammation and immunity are overrepresented among the genes differentially expressed between African Americans and European Americans. Table 2 lists these GO terms. Supplementary  Table S29 provides the overrepresented Gene Ontology (GO) terms among DE genes for each diseased or non-diseased tissue.
We drew from in silico studies 48 and experimental analyses especially human responses to infection by SARS-CoV-2 and other coronaviruses [49] [50] [51] [52] to identify ten genes implicated in cellular responses to SARS-CoV-2 infection that are among those differentially expressed in African American and European American populations. Molecular functions of these genes include receptor kinases, cytokines, other signal transduction molecules, and antioxidants. These genes are integral to central cellular processes that affect pathogenesis by SARS-CoV-2, including endosomal development, autophagy, immunity and inflammation 6, 51, 53 . Molecular functions of these genes include receptor kinases, cytokines, other signal transduction molecules, and antioxidants.
A number of genes of the immune response are differentially expressed between African American and European American populations. Expression of IL6ST, a component of the cytokine receptor complex that acts as signal transducer for cytokine interleukins IL6 and IL7, is 2-fold lower in African Americans than European Americans in BRCA ( Figure 1 ). Circulating chemokines CXCL9 and CXCL10 are also differentially expressed in African Americans as compared to European Americans in several cancers. CXCL9 expression is 2-fold greater in KIRP and over 2-fold lower in COAD and KIRC; CXCL10 expression is 1.5-fold higher in BRCA, and over 2-fold lower in COAD, KIRC and TCHA ( Figure 2 ).
The small inducible chemokine CCL3L3 is upregulated in African Americans by 2-to 3-fold in BRCA, COAD, and KIRP ( Figure 3 ) and in several non-diseased tissues (Supplementary Table S16-S25) .
Carcinoembryonic Antigen-related Cell Adhesion Molecules CEACAM5 and CEACAM6 are both downregulated 2 to 3-fold in BRCA in African Americans. (Supplementary Table 2 ).
Expression of GSTM1, a key enzyme involved in oxidative stress, differs between African Americans and European Americans. Expression is up to to 6-fold higher in African Americans in the cancers we evaluated ( Figure 4 ). Distribution of expression in European Americans, but not African Americans is bimodal. GSTM1 expression is also DE in non-diseased esophagus and thyroid gland (Supplementary Table S18, S24) .
Endocytosis and autophagy are intimately interrelated with Covid-19 54 . One little-studied player implicated in early endosome motility 55 and hence the endocytotic pathway and autophagy, is the seven tetratricopeptide-like repeat F8A/HAP40 (HAP40) protein 56 . In human genomes, three genes, F8A1, F8A2, and F8A3, encode the HAP40 protein 56 . The F8A genes are located on the X chromosome. F8A1 is within intron 22 of the coagulation factor VIII gene, which has a high frequency of mutations 57 ; F8A2, and F8A3 are located further upstream.
F8A1, F8A2 and F8A3 are each differentially expressed in African Americans versus European Americans. F8A1 is more highly expressed by about 2-fold in European Americans in every cancer analyzed ( Figure 5 ) and by 2-fold in non-tumor colon (Supplementary Table S17 ). Conversely, F8A2 and F8A3 are more highly expressed in African Americans in all cancer types. Expression of F8A2 in African Americans is up to 24-fold greater; expression of F8A3 is up to 6.6-fold greater. In LUSC, F8A2 and F8A3 are the only DE genes (Supplementary Table S7 ). F8A2 and F8A3 follow a similar trend in non-diseased tissues, being more highly expressed in African Americans by up to 4-fold in colon, esophagus, and thyroid (Supplementary Table S17-S18). Distribution of F8A2 and F8A3 expression is bimodal in European Americans for most cancers. Thus, part of the difference in levels of F8A2 and F8A3 expression between the two populations is due to their distribution, with low levels of expression in a subset of the European American population.
Because of the paucity of literature on HAP40 58 and because, to our knowledge, the relationships among F8A1, F8A2, and F8A3 genes have not been described, we investigated further the sequences, sequence variants, and the expression patterns of these genes.
The sequences of the HAP40 proteins of F8A1, F8A2, and F8A3 are identical to each other in human reference genome GRCh38.p13 (https://www.ncbi.nlm.nih.gov/assembly/GCF_000001405.39). Allele variants of HAP40 proteins encoded by F8A1, F8A2, and F8A3 were mined from The Genome Aggregation Database (gnomAD) 59 , a open database, which contains sequences of over 140,000 exomes and genomes from individuals of diverse populations as categorized by clustering of genetic features (rather than being self-reported). Individuals are assigned to one of the five major populations and to sub-populations within these. The search of gnomAD identified no variants in the HAP40s encoded by F8A2 or F8A3, and a single very rare variant (<1/1000) of F8A1 found only in European (non-Finnish) populations. The variant encodes a missense mutation (https://gnomad.broadinstitute.org/gene/ENSG00000197932?dataset=gnomad_ r2_1). No structural variants were identified for HAP40 of F8A1 or F8A3; F8A2 has a rare duplication of 54 aa.
To our knowledge, F8A2 and F8A3 gene expression has not been described. We analyzed coexpression of the three F8A genes in the context of the other 18,212 genes represented in the full TCGA-GTEx dataset, using the "Pearson correlation" function in MOG. Although the F8A genes are proximately located on the X chromosome, the genes are not highly coexpressed with each other (|Pearson Corr.| ă 0.46). Furthermore, no F8A gene is coexpressed with any other of the 18,212 genes represented in the dataset (|Pearson Corr.| ă 0.46) (Supplementary Table S30 ). Of the 18,212 genes, the expression of F8A1 is most negatively (anti-) correlated with those of F8A2 and F8A3, with Pearson Correlations of´0.45 and´0.24, respectively (Supplementary Table S30 ). Mutual information (MI) can detect linear as well as complex non-linear associations, whereas Pearson's correlation measure quantifies linear dependencies. Using the MI function in MOG, we found that F8A2 and F8A3 genes are most associated with F8A1 (Supplementary Table S30 ), indicating that there may be a negative interaction among these genes.
Genetics of human populations contribute to the propensity and severity of diseases 37, 39, 41, 41, 45, [60] [61] [62] [63] [64] [65] . Sometimes the contribution is straightforward; a single allele variation found in Ashkenazi Jews, causes the vast majority of Tay-Sachs disease 66 . Sometimes it is more complex; for example, hypertension, which more prevalent in African American than European American populations 60 in part due to detrimental APOL1 mutations that are more frequent in West African populations 62 . Despite the paucity of studies focused on Western African populations, the propensity and severity of several other diseases among this population have been attributed to genetics 41, 62, 67 .
The individual's immune system is key to fighting viral infections. However, conversely, many COVID-19 deaths have been attributed to a cyclic over-excitement of the innate immune system. This latter process, often termed a cytokine storm, results in a massive production of cytokines and the body attacking itself, rather than specifically destroying the pathogen-containing cells 14 . Thus, people with comorbidities, the elderly, and immunosuppressed individuals, may be at a greater risk for COVID-19 morbidity and mortality because they may not respond to infection with sufficient immune response 50 and/or because they may be more likely to develop a cytokine storm 14 . We focused on ten genes that are differentially expressed between African Americans and European Americans are implicated in these biological processes.
The most dramatic differences in gene expression were in expression of the F8A genes, F8A1 being upregulated in European Americans, and F8A2 and F8A3 being upregulated in African Americans. Each F8A gene encodes an identical HAP40 protein. HAP40 function has been researched only in the context of its critical role in early endosome maturation in Huntington's disease 56 . In Huntington's, HAP40 forms a bridge between the huntingtin protein and the regulatory small guanosine triphosphatase, RAB5; formation of this complex reduces endosomal motility by shifting endosomal trafficking from the microtubule to the actin cytoskeleton 55 . High F8A1 expression has been reported in several conditions: Huntington's 68 ; a SNP variant for type 1 diabetes risk 69 ; cytotrophoblast-enriched placental tissues from women with severe preeclampsia 70 ; and mesenchymal bone marrow cells as women age 71 . Although its non-disease biology has been little explored, because of its role in early endosome motility in Huntington's, HAP40 is considered a potential molecular target in therapy of autophagy-related disorders 72 .
Endosome motility and development play an important but complex role in the innate immune response, which can either promote or hinder the battle between SARS-CoV-2 and its human host 51, 54, 73, 74 . Coronaviruses including SARS-CoV-2 mainly enter host cells via binding to the ACE2 receptor followed by endocytosis 20, 29, 51 . Nascent early endosomes are moved along the microtubule cytoskeleton, fusing with other vesicles; varied molecules can be incorporated into the membrane or the interior 51, 54, 73, 74 . This regulated development enables diverse fates. For example, in the context of SARS-CoV-2, endosomes might release viral RNA or particles; they might merge with lysosomes and digest their viral cargo; or might fuse with autophagosomes (autophagy) and subsequently with lysosomes that digest the cargo 51, 54, 73, 74 . SARS-CoV-2 might reprogram cellular metabolism, suppressing autophagy and promoting viral replication, 75 . The cell might modify autophagy machinery to decorate viral invaders with ubiquitin for eventual destruction, activate the immune system by displaying parts of the virus, or catabolize excess pro-cytokines. Autophagy might induce cytokine signaling, which could promote protective immune response or engender a destructive storm of cytokines, inflammation and tissue damage 54 .
Cytokines and other immunomodulatory molecules, including CCL3L3, CXCL9, CXCL10, CEACAM5 and CEACAM6, were differentially expressed between African Americans and European Americans. CCL3L3 is a member of the functionallydiverse C-C motif chemokine family. It encodes CCL3, which acts as ligand for CCR1, CCR3 and CCR5 recruits and activates granulocytes; it also inhibits HIV-1-infection 76 . CCL3L3 is upregulated in younger and impoverished white males 77 . Circulating chemokines CXCL9 and CXCL10 initiate human defenses, and potentially instigate autoimmune and inflammatory diseases, by activating G protein-coupled receptor CXCR3 [78] [79] [80] . CEACAM5 and CEACAM6 are members of the C-Type Lectin Domain Family. This gene family encodes a diverse group of calcium (Ca2+)-dependent carbohydrate binding proteins, several of which, including CEACAM5 and CEACAM6 have been implicated as having specific cell adhesion, pathogen-binding and immunomodulatory functions 81 . CEACAM5, a driver of breast cancer 82 and modulator of inflammation in Crohns Disease 83 , and CEACAM6, an inhibitor of breast cancer when coexpressed with CEACAM8 84 , are both downregulated 2 to 3-fold in BRCA in African Americans.
Reactive Oxygen Species (ROS) generated in the mitochondria promote the expression of proinflammatory cytokines and chemokines, thus playing a key role in modulating innate immune responses against RNA viruses [85] [86] [87] [88] [89] . Mitochondrially-targeted glutathione S-transferase, GSTM1, which was more highly expressed in African Americans than European Americans, is a key enzyme in the metabolism of ROS, as well as xenobiotics including pharmaceuticals 88 . GSTM1 is induced by nuclear factor erythroid 2-related factor 2 (Nrf2), a transcription factor that integrates cellular stress signals [90] [91] [92] [93] . Low expression of GSTM1 can lead to increased mitochondrial ROS, which may ultimately result in a cytokine storm that triggers inflammation and/or autoimmune disease. Conversely, if GSTM1 is too highly expressed, pharmaceuticals may be metabolized and thus rendered inactive, and ROS may be metabolized too rapidly to maintain a sufficient signaling role in the immune system. Allele frequencies of GSTM1 vary among Asian, African and European populations 94 ; the biological significance of these alleles is being investigated 95, 96 .
CEACAM5 and CEACAM6, both downregulated 2 to 3-fold in BRCA in African Americans, are members of the C-Type Lectin Domain Family. This gene family encodes a diverse group of calcium (Ca2+)-dependent carbohydrate binding proteins; CEACAM5 and CEACAM6 have been implicated as having specific cell adhesion, pathogen-binding and immunomodulatory functions 81 . CEACAM5 is a driver of breast cancer 82 and modulator of inflammation in Crohns Disease 83 , and CEACAM6 is an inhibitor of breast cancer when coexpressed with CEACAM8 84 .
By revealing differential expression of genes implicated in COVID-19 morbidity and mortality between African Americans and European Americans, we emphasize the importance of integrating gene expression data into the mix of factors considered in studying this pandemic. Our study indicates that, under both diseased and non-diseased conditions, many genes involved in infection, inflammation, or immunity are differentially expressed between African Americans and European Americans. One contributing explanation of the finding that disease-related genes are overrepresented among DE genes is that the selection pressure due to disease is very strong on both (ancestral) regions, but these regions have very different complements of pathogens. Humans living in Europe and those living in Western Africa would have had to evolve the ability to resist the prevalent local pathogens.
Archived expression data has tremendous value. However, studies such as this one are hampered by several factors. For example, obtaining adequate sample sizes for statistical analysis of populations is very important but difficult or expensive to address. Ethnic bias and practical factors (such as subject availability) often result in insufficient numbers of subjects from many populations to be represented in medical studies; this lack of representation prevents the development of precise prognosis or therapy based on genetics 39, 97 . Similarly, diverse socioeconomic contexts may not be well represented among the individuals sampled. Yet these are clearly a factor in disease 98 . The statistical predictor provided for the U.S. by the Robert Wood Johnson Foundation (https://www.rwjf.org/en/library/interactives/whereyouliveaffectshowlongyoulive. html) reflects the concept that "your zip code can be greater than your genetic code" (although unfortunately there is no accompanying genetic information).
Another factor that would greatly increase the utility of archived expression data is improving and extending the metadata for all (future) studies; this would be relatively simple to implement. Well-constructed metadata is key to the usefulness of data. Among the vast body of human RNA-Seq data being deposited, fields for age and gender are typically represented and available 9/14 in the metadata. However, except in specialized studies, metadata on the race and ethnic heritage of the sampled individuals are often not included, or are very difficult to access. The same is even more true of fields that would provide socio-economic information, such as postal code or risk factors such as occupation. Because of the absence of socioeconomic metadata, even in GTEx and TCGA, the arguably most comprehensive RNA-Seq datasets to date, we were unable to distinguish genetic effects from environmental causes of the differences in gene expression. Without routine inclusion of diverse metadata for human 'omics samples, data re-mining is hampered, and important information is lost.
Multiple genes implicated in COVID-19 are differentially expressed in African American and European American populations. The differential expression is evident despite the fact that race is self-reported in and metadata, and that many Americans are racially admixed 41 .
Gene expression represents the interaction of genetic and environmental factors. Routine inclusion of information on ethnicity, race, postal code (as a proxy for socioeconomic condition), and profession in the metadata for each individual sampled would empower large-scale data-driven approaches to dissect the relationships between race, socio-economic factors, and disease.
By highlighting the wide-ranging differences in expression of several disease-related genes across populations, we emphasize the importance of harvesting this information for medicine. Such research will establish prognostic signatures with vast implications for precision treatment of diseases such as COVID-19.
The MOG tool was used to interactively explore, visualize and perform differential expression and correlation analysis of genes. We downloaded the precompiled MOG project http://metnetweb.gdcb.iastate.edu/MetNet_MetaOmGraph. htm 42 , created using the data processed by Wang et. al., in which expression values have been normalized and batch corrected to enable comparison across samples 43 . This MOG_HumanCancerRNASeqProject contains expression values for 18,212 genes, 30 fields of metadata detailing each gene, across 7,142 samples representing 14 different cancer types and associated non-tumor tissues (TCGA and GTEX samples) integrated with 23 fields of metadata describing each study and sample.
Since the data was normalized and batch corrected, we used Mann-Whitney U test, a non-parametric test, to identify differentially expressed genes between two groups. An R script was written to perform KS and dip tests, and create the violin plots and executed via MOG interactively.
Pearson correlation values were computed, after data was log 2 transformed within MOG, in MOG's statistical analysis module.
Information on how to reproduce the results are available at https://github.com/urmi-21/COVID-DEA.
We subscribe to FAIR data and software practices 99 
Supplementary data are available at bioRxiv.
This work is funded in part by the Center for Metabolic Biology, Iowa State University, and by the National Science Foundation award IOS 1546858, Orphan Genes: An Untapped Genetic Reservoir of Novel Traits.