Comprehensive analyses of human transcriptome across tissues were initiated by downloading on 02/22/2019 public human RNA-seq data (GTEx Analysis v7; dbGaP Accession: phs000424.v7.p2, release date: June 30, 2017) deposited in the GTEx portal²) (Figure 1). Gene expression values in the data were derived from non-diseased normal tissues of postmortem human donors (aged 21–70) and normalized using transcripts per million (TPM). Before processing the data, tissues from less than 30 donors and non-tissue samples such as cultured transformed fibroblasts and EBV-transformed lymphocytes were excluded from the current study. Then, tissue samples with high quality RNA (RNA integrity number, RIN, of 6.0 or higher) were selected. Finally, the filtered data included 47 human tissues, and the sample size (i.e., donor numbers) for each tissue ranged from 35 to 560 (Supplementary Table S1). In addition, 56,202 human protein-coding and non-coding genes were filtered by general thresholds of at least 20% of samples having TPM > 0.1³) and median TPM > 0.5 in the heart (Ahn et al., 2019), leaving 16,480 and 14,911 protein-coding and non-coding genes for data analysis for the atrial appendage (AA) and left ventricle (LV), respectively, (Figure 2A).

Enrichment analysis on randomly selected genes was conducted as follows. From the above 16,480 (AA) and 14,911 (LV) genes, the top 500 genes were selected by the rank of relative median value. From either 16,480 (or 14,911) genes or the top 500 genes, 5 to 500 genes were randomly selected by adding one gene per new selection, using the sample_n function of the dplyr package (v0.8.5⁴). The same random selection procedure was also applied for AA-specific/enhanced genes (5 to 81 genes) and LV-specific/enhanced genes (5 to 79 genes). The kyoto Encyclopedia of Genes and Genomes (KEGG) pathway was selected as a reference database, and we ran an Enrichr algorithm⁵ (Chen et al., 2013; Kuleshov et al., 2016) to query to the KEGG 2019 human database (308 terms, gene coverage of 7,802, and 92 genes per term) as presented in a previous publication (Zachariou et al., 2018). In the current study, a Python script (enrichr-cli.py, W. De Coster⁶) was used to run the algorithm and obtain Enrichr combined scores for input gene lists, which are derived from a correction to the Fisher’s exact test p-value by the z-score of deviation from the expected rank (Chen et al., 2013).

Heart-specific genes were defined as genes with fold changes of the median expression levels (FCMs) higher than 5.0 in the heart versus all other tissues. Heart-enhanced genes were defined as genes having FCMs higher than 5.0 in the heart versus all other tissues, with an exception of 1 or 2 other tissues, and a higher median expression in the heart in the case of exceptions. That is, we selected a trend in which the heart showed a greater expression than any other analyzed tissues even in the case of heart-enhanced genes as well as heart-specific genes.

Using lmFit function in the R/Bioconductor limma package (Ritchie et al., 2015), a linear model was fitted to our data and a contrast matrix was created for all 45 pair-wise comparisons between each tissue. In limma, the empirical Bayes (eBayes) method (Smyth, 2004) followed by the decideTests function with “global” option was used to test t-statistic for each contrast and F-statistics for all 45 contrasts simultaneously, along with multiple testing correction by the Benjamini–Hochberg’s procedure (Benjamini and Hochberg, 1995). As a result, genes with at least a 5-fold higher median expression in the heart and whose expression is significantly different between tissues, formulated as: [FDR of i versus j pairwise comparison] < 0.01, where i = A1 (heart-atrial appendage) or A2 (heart-left ventricle) and j = B, C, D, ⋯, or AT (non-heart), were selected.

For grouping protein-coding and non-coding genes (81 genes in AA and 79 genes in LV) related to long intergenic non-coding RNA (linc RNA), pseudogene, antisense RNA, and sense intronic transcripts, the Ensembl Gene IDs retrieved from the GTEx portal (hg19/GRCh37) were reannotated using the human reference genome assembly (GRCh38.p12⁷). Heat maps were generated to visualize relative heart-specific and heart-enhanced expressions using the R package “heatmap3”⁸ (Zhao et al., 2014).

To access to PubMed⁹ and download abstracts of publications on the heart and heart-related medical conditions, the PubMed literature database was queried on May 26, 2020 using the R package “easyPubMed”¹⁰. Briefly, the PubMed database was queried without restriction on publication date and using search terms independently for each of the 57 heart- specific/enhanced protein-coding genes as follows: ‘heart[TIAB] OR atrial[TIAB] OR ventricle[TIAB] OR cardio[TIAB] OR cardiac[TIAB] OR coronary[TIAB] OR cardium[TIAB] AND (genesymbol [TIAB])’, where TIAB stands for title/abstract and genesymbol was replaced with each different gene symbol. Abstracts were downloaded in XML format. Extraction of biomedical entities such as genes, diseases, species, and PubMed identifier (PMID) from the downloaded abstracts was performed using the R package “pubmed.mineR”¹¹ (Rani et al., 2015). When no abstract was retrieved for a certain gene by easyPubMed, it was referred to as a previously “unreported gene” in PubMed (hereafter unreported gene). When an abstract was retrieved and biomedical entities were extracted, it was referred to as a “reported gene” in PubMed (hereafter, a reported gene). For some cases in which abstracts were retrieved but related diseases were not extracted, a manual review on abstracts and full texts under each extracted PMID was performed as a part of data mining strategy (Pospisil et al., 2006) (e.g., research on non-pathological conditions such as heart rate was regarded as a reported case).

The Ingenuity Pathway Analysis (IPA) was used to perform pathway enrichment analysis as follows. IPA Comparison Analysis was conducted on two groups of 20 tissues (i.e., either AA and 19 other tissues or LV and 19 other tissues) to examine similarities and differences in pathway enrichment. The selection of 20 tissues was due to the maximum analyzing capacity of IPA Comparison Analysis of 20 observations (or experimental groups), and tissues with greater size and larger sample numbers were chosen. Among the top 500 genes, all heart-specific/enhanced genes were analyzed, and their log ratios (i.e., log₂RMVs) were used to rank them. For those heart-specific/enhanced genes, adjusted p-values of differential expression between tissue i versus rest, where i is one tissue from 20 tissues and rest is 19 other tissues, were obtained from the limma package using a contrast matrix of one versus 19 followed by the eBayes method and t-statistic. Using log₂RMVs and adjusted p-values, IPA Core Analysis was performed on each of the 20 tissues and then IPA Comparison Analysis was conducted using 20 Core Analysis results. Enrichment was considered significant when p < 0.01. A cutoff set for activation or inhibition of the pathways predicted by IPA was |activation z-score| ≥ 2.

Read count matrix from RNA-seq of atrial appendage and left ventricle was downloaded from the GTEx portal on 02/22/2019. Read counts ≥ 6 in at least 20% of samples¹² was satisfied with 23,288 genes in AA and 22,137 genes in LV, and the union (23,610 genes) of the two gene sets were retained. For differential expression analysis, the R/Bioconductor “DESeq2” package (Love et al., 2014) was used to normalize samples for sequencing depth and calculate DEGs between AA and LV. To obtain significant DEGs, combined criteria of FDR < 0.01 and absolute value of log₂(fold change) > 3 were used, where a fold change is defined as the expression in samples of AA divided by that of LV. IPA Core Analysis was conducted to analyze pathway enrichment of those significant DEGs, using log₂(fold change) and adjusted p-values as input values. The p < 0.01 was considered a significant enrichment. The |activation z-score| ≥ 2 was used to predict activation or inhibition of the pathways.

With the purpose of identifying heart-specific genes and their associations with diseases and traits, our integrative workflow was designed (Figure 1). To initiate the workflow, the human RNA-seq-based transcriptome data (GTEx Analysis v7) were downloaded from the GTEx portal¹³. In each tissue, the medians of gene expression values (TPMs) were obtained in virtue of their robustness to outliers. Then, the median TPM value in the heart was divided by an average of the rest of the median TPMs from other tissues, producing ‘relative median values (RMVs)’, and distribution of these values were plotted against the number of genes (Figure 2A). The plots were relatively symmetric and displayed about one hundred and fifty values with more than 10-fold, indicating the presence of exclusive or abundant expression in the two heart tissues: atrial appendage and left ventricle.

To examine whether the degree of enrichment changes by selecting genes with high RMVs, the number of enriched pathways were queried by increasing the input gene set cumulatively (Figure 2B). At first, five genes were randomly selected from two groups: one is the gene set with all valid genes (i.e., pre-selected 16,480 genes for AA or 14,911 genes for LV) and another is with the top-ranked 500 genes (more than 2.51- and 2.32-fold in AA and LV, respectively). For the cumulative increase, one gene was added into the gene set each time after random selection from either group. The KEGG 2019 human database was selected as a reference pathway database for this test query owing to its comprehensiveness and a wide range of usage (Kanehisa et al., 2017). Upon running through the Enrichr algorithm, cumulative numbers of total enriched pathways with Enrichr combined scores, which outperform other measurements on enrichment (Kuleshov et al., 2016; Zachariou et al., 2018), were plotted for each gene list with respect to numbers of randomly selected genes (Figure 2B).

As shown in Figure 2B, in either case, with less than about 50 to 70 genes, the number of total enriched pathways was higher in gene sets selected from the top 500 genes compared to the all genes. However, for larger gene sets, the number of enriched pathways for the top 500 gene set became substantially smaller than the all genes, and this pattern was continued in both cases of AA and LV. It indicated that the extent of enrichment was affected by RMVs and choosing genes with high RMVs increased the detection sensitivity. Taken together, highly expressed genes in the heart could be detected by sorting RMVs, and selection of the top 500 genes improved enrichment results.

In order to identify more specific genes in the heart, we further narrowed down those top 500 genes to “heart-specific genes” and “heart-enhanced genes” using rigorous criteria as described in the section “Materials and Methods.” As a result, we selected a total of 81 protein coding and non-coding genes in AA and 79 in LV (all FCMs > 5 with up to 2 exceptions, FDR < 0.01; Figure 2C and Supplementary Table S1). Among them, the number of protein coding genes were 45 and 44 in AA and LV, respectively, and there were 32 overlapping genes. In addition, with the same random selection strategy used in the case of the top 500 genes (Figure 2B), those AA- and LV-specific/enhanced genes were subjected to the cumulative pathway enrichment test and showed a more refined enrichment than all genes similar to the top 500 genes (Figure 2D). It also indicated that the detection sensitivity increased with AA- and LV-specific/enhanced genes because the refined enrichment occurred with smaller gene sets, compared to the top 500 genes, and these specific/enhanced genes were further studied below.

Among those overlapping genes, three genes (RD3L, FBXO40, and SMCO1) have not been documented in PubMed (i.e., unreported genes) regarding their related diseases and functions in the heart in humans based on our literature search. Genes whose function has been investigated in non-pathological conditions such as CCDC141 on heart rate (van de Vegte et al., 2019), as well as disease conditions, were classified as reported genes. Among non-overlapping genes, 13 protein-coding genes were AA-specific/enhanced and among them the following four genes were unreported genes: SBK2, PRR32, and SBK3 showed exclusive expression (AA-specific), and expression of CHRNE was more than 5-fold compared to other tissues, except pituitary (AA-enhanced). In addition, 11 non-overlapping LV-specific/enhanced protein coding genes were found and all of them were reported. Heat maps displayed z-scores of RMVs of the heart-specific/enhanced protein-coding genes in various tissues (i.e., each row) with extreme specificity or abundancy to AA and LV (Figure 3). Those heat maps were extended for 16,480 genes in AA and 14,911 genes in LV and divided by each chromosome, displaying specificity or abundancy of the heart-specific/enhanced protein-coding and non-coding genes in each chromosome compared to the rest of the genes (Supplementary Figures S1, S2).

Given the extreme specificity and abundancy of the heart-specific/enhanced genes, pathological or functional implications of those genes were explored. Using IPA software, enriched “Diseases and Biological Functions” were identified in the following three groups of gene sets: 58 common-specific/enhanced protein-coding and non-coding genes (α), 81 AA-specific/enhanced protein-coding and non-coding genes (α + β), and 79 LV- specific/enhanced protein-coding and non-coding genes (α + γ) (Figure 4A). Based on the aforementioned cumulative enrichment test, more than 50 specific/enhanced genes were selected to produce enrichment results with IPA.

With common-specific/enhanced genes (α), out of a total of approximately 900 terms listed in “Diseases and Biological Functions,” 18 terms had z-scores in at least one tissue among 20 tissues. Among them, 5 functional terms (Function of Muscle, Cardiac Contractility, Contractility of Cardiac muscle, Cardiogenesis, and Contraction of Heart) were predicted to be activated (z-score ≥ 2) and 4 disease terms (Congestive Heart Failure, Fibrosis, Failure of Heart, and Organismal Death) were predicted to be inhibited (z-score ≤ −2), both when common-specific/enhanced genes were compared between AA and other tissues (Figure 4B) and between LV and other tissues (Supplementary Figures S3A,B). Most other tissues showed opposite prediction of activation or inhibition, and prediction was not significant for skeletal muscle (for all of the above 9 terms), testis (8 terms), spleen (5 terms) and coronary artery (2 terms) showing these diseases and biological functions were not oppositely inhibited or activated in those tissues (Figure 4B and Supplementary Figure S3B). Nonetheless, none of the other tissues showed the same pattern of prediction as the heart, indicating the exclusive role of common-specific/enhanced genes in promoting those biological functions and inhibiting those diseases.

With AA-specific/enhanced genes (α + β), 6 terms were predicted to be activated (z-score ≥ 2) and four terms were predicted to be inhibited (z-score ≤ −2) in AA (Supplementary Figure S3C), and with LV-specific/enhanced genes (α + γ), 7 terms were predicted to be activated (z-score ≥ 2) and 5 terms were predicted to be inhibited (z-score ≤ −2) in LV (Supplementary Figure S3D). Those terms were all related to activation in cardiac muscle contraction and cardiogenesis (functional terms) and inhibition of heart failure, fibrosis and organismal death (disease terms), as shown in the case of common-specific/enhanced genes. It suggested essential roles of common-specific/enhanced genes in promoting contraction of the heart and cardiac development and preventing dysfunction and failure of the heart.

Specific canonical pathways were further explored to compare common-specific/enhanced genes with AA- or LV-specific/enhanced genes. From IPA, prediction of activation for “Factors Promoting Cardiogenesis in Vertebrates” were found to be higher in AA (z-score 2.45) and LV (z-score 2.65), compared to common-specific genes (z-score 2.24). The AA- or LV-specific/enhanced genes, that are not common-specific/enhanced and might increase those z-scores were BMP10 in AA and MYH7 (or β-MHC) and MYL2 (or MLC2v) in LV (Figures 4C,D). In this study, the BMP10 gene showed the highest specificity (Figure 3) but the expression ratio was slightly changed when comparing 20 tissues for IPA ranking BMP10 the second in Figure 4C. The MYH7 and MYL2 genes were LV-enhanced genes (Figure 3).

In addition, in relation to cardiac muscle contraction, prediction of activation for another canonical pathway “Actin Cytoskeleton Signaling” was higher in AA (z-score 2.24) and LV (z-score 2.65), compared to common-specific genes (z-score 2.00). The AA- or LV-specific/enhanced genes, that are not common-specific/enhanced and might increase those z-scores were FGF23 and MYL4 in AA and MYL3, MYH7 and MYL2 in LV (Figures 4E,F).

For disease terms, prediction of inhibition was increased the most for “Organismal Death” with AA-specific/enhanced genes (z-score −3.05) and LV-specific/enhanced genes (z-score −3.07) compared to common-specific/enhanced genes (z-score −2.38). Also, prediction of inhibition for “Failure of Heart” (z-score −2.219) with common-specific/enhanced genes was increased with AA-specific/enhanced genes (z-score −2.43), and “Fibrosis” (z-score −2.09) with common-specific/enhanced genes was predicted to be inhibited more with LV-specific/enhanced genes (z-score −2.31). It suggests that those above canonical pathways related to cardiogenesis and cardiac muscle contraction and AA- and LV-specific/enhanced genes might be related to suppressing these pathological conditions and promoting organismal survival and function.

In order to gain functional and mechanistic insights, the above canonical pathways were further investigated. In relation to “Factors Promoting Cardiogenesis in Vertebrates,” bone morphogenetic protein 10 (BMP10) is a peptide growth factor and both SMAD-dependent and independent pathways were predicted to be activated via phosphorylation (Figure 5). A subsequent binding of SMAD1/5/8 and co-SMAD (SMAD4) and activation of NKX2.5, the determinant of cardiac development and the common-enhanced gene in our case, through both SMAD complex and GATA-4 (Brown et al., 2004) might be a key determination step of cardiogenesis in atrial appendage. The role of BMP10 in the adult heart remains elusive, but as reported recently BMP10 may play multifunctional roles not only in promoting cardiac development in prenatal stages but also in inhibiting cardiomyocyte apoptosis and cardiac fibrosis in adults (Qu et al., 2019). In the downstream of this canonical pathway, cardiac transcription factors were predicted to activate target genes including AA-specific and LV-enhanced genes: natriuretic peptide B (NPPB) and natriuretic peptide A (NPPA). These paralogous genes are transcriptionally regulated as a cluster and encode peptide hormones that are secreted by cardiomyocytes (Man et al., 2018), and was predicted to induce cardiomyocyte differentiation (Figure 5). The expression of the other two target genes, sodium voltage-gated channel alpha subunit 5 (SCN5A) that encodes an integral membrane protein, a sodium channel, for transmitting electric signals to maintain normal heart contraction (Zaklyazminskaya and Dzemeshkevich, 2016) and myosin heavy chain 6 (MYH6) which encodes the sarcomeric filament protein alpha-myosin heavy chain (Posch et al., 2011), were relatively lower than NPPB and NPPA. In LV, additional genes include myosin heavy chain 7 (MYH7 or β-MHC) that encodes beta-myosin heavy chain and myosin light chain 2 (MYL2 or MLC2v) that encodes regulatory myosin light chain that are up-regulated as LV-enhanced genes.

In addition, regarding “Actin Cytoskeleton Signaling” which was predicted to be activated with heart-specific/enhanced genes, one of the genes that were additionally involved in this enriched pathway in AA was fibroblast growth factor 23 (FGF23), which is an AA-enhanced gene (Figure 6). Fibroblast growth factor 23 (FGF23) is a circulating hormone that is associated with cardiac hypertrophy and also acutely elevates intracellular Ca²⁺ and increases cardiac contractility (Touchberry et al., 2013). Binding of FGF23 to receptor tyrosine kinases leads to tyrosine phosphorylation and stimulation of Ras/MARK and PI3K/Akt signaling (Grabner et al., 2015). Another additional gene in AA was myosin light chain 4 (MYL4), a previously known atrial-specific myosin light chain, whose mutation causes familial atrial fibrillation (Orr et al., 2016). In LV, MYH7, MYL2, and myosin light chain 3 (MYL3) were additionally involved in this canonical pathway. Although both AA-specific/enhanced genes and LV-specific/enhanced genes were predicted to promote actin polymerization and cytoskeleton reorganization, those different cardiac sarcomere protein genes might give rise to distinct contractile and electrophysiological properties between AA and LV.

A total of 341 DEGs (Log₂(fold change) ≥ 3, FDR < 0.01) were identified, among which 236 and 105 DEGs showed up-regulated expression in AA and LV, respectively, (Figure 7A and Supplementary Table S2). Protein-coding DEGs were 170 and 55 in atrial appendage and left ventricle, respectively, (Figure 7B). By doing IPA Core Analysis on 341 DEGs, twenty networks regarding diseases and functions were identified. Among them, the top ten networks were scored 19 to 42 [p-score = −log₁₀(p-value)] with 13 to 23 focus molecules (Supplementary Table S3). The first-ranked network was “Cardiovasuclar System Development and Function, Embryonic Development, Organ Development” with 23 focus molecules (p-score = 42). In this network, a subnetwork with the most significant enrichment was “Morphogenesis of Cardiac Muscle” (p = 8.96 × 10^–7) with 4 molecules: BMP10, HAND1, MYL3, and XIRP2. Among those four genes, BMP10 was upregulated in AA (log₂FC = 12.38), and the other three genes were downregulated in AA (HAND1 log₂FC = −3.03, MYL3 log₂FC = −4.50, and XIRP2 log₂FC = −3.05) and upregulated in LV. Another significant subnetwork was “Hypertrophy of Left Ventricle” (p = 4.08 × 10^–5), and related molecules were IRX4 (log₂FC = −7.33), MYL3 (log₂FC = −4.50), NPPA (log₂FC = 6.05), and REN (log₂FC = 7.20).

The network between 23 focus molecules under the first-ranked “Cardiovasuclar System Development and Function, Embryonic Development, Organ Development,” including the top two subnetwork molecules (the above seven molecules), was displayed (Figure 7C). In addition, the relationship of this network to canonical pathways was examined and the “Calcium Signaling” pathway showed a significant enrichment (p = 3.31 × 10^–5) and the highest prediction of activation (z-score = 2.24) with this experimental group (the combination of up-regulated DEGs in AA and down-regulated DEGs in LV) (Figure 7C). Eight molecules (MYL7, MYL1, MYL3, ERK1/2, Mlc, MYL4, MYL2, and Rlc) were involved in this calcium signaling pathway. On the other hand, other canonical pathways that were related to AA- and LV-specific/enhanced genes were not significantly enriched or not predicted to be activated or inhibited: “Factors Promoting Cardiogenesis in Vertebrates” (p = 0.06, z-score = 0) and “Actin cytoskeleton signaling” (p = 5.37 × 10^–3, z-score = −0.447). It suggests that there may be functional differences between heart-specific/enhanced genes and DEGs.

In order to find heart-specific/enhanced genes that are differentially expressed between AA and LV and are previously unreported in PubMed, we first compared heart-specific/enhanced genes with DEGs. We found that 9 AA-specific/enhanced protein-coding genes were up-regulated DEGs in AA. Thus, these genes (BMP10, MYBPHL, CHRNE, NPPA, SBK2, MYL4, OR10P1, PRR32, and MYL7) were “AA-specific/enhanced DEGs” meaning that they were specific/enhanced to AA compared to 45 non-heart tissues as well as up-regulated in AA compared to LV (Supplementary Table S2). Also, there were 3 LV-specific/enhanced protein-coding genes that were up-regulated DEGs in LV (MYL2, MYL3, and MYH7), and they were named ‘LV-specific/enhanced DEGs’ which were specific/enhanced to LV compared to 45 non-heart tissues as well as up-regulated in LV compared to AA.

Among both AA- and LV-specific/enhanced DEGs, three AA-specific/enhanced DEGs [cholinergic receptor nicotinic epsilon subunit (CHRNE), SH3 domain binding kinase family member 2 (SBK2), and proline rich 32 (PRR32)] were “unreported” genes in PubMed regarding the heart. On the other hand, all the LV-specific/enhanced DEGs are known myosin genes. From the IPA Network Analysis for DEGs, networks for two out of the three AA-specific/enhanced DEGs were found (Supplementary Table S3). The CHRNE gene was one of the focus molecules under “Developmental Disorder, Hereditary Disorder, Immunological Disease” (p-score = 37 with 21 focus molecules) (Figure 8). Also, the PRR32 gene was one of the focus molecules under “Cell Morphology, Embryonic Development, Hair and Skin Development and Function” (p-score = 19 with 13 focus molecules) (Figure 9).

CHRNE is a transmembrane acetylcholine receptor (AChR) and its mutations lead to congenital myasthenic syndrome which is responsible for perturbation on calcium signaling and post-synaptic Ca²⁺ accumulation (Zayas et al., 2007; Shen et al., 2018). Along with calcium voltage-gated channel alpha subunits (CACNA1D, CACNA1G, and CACNA2D2 in Figure 8), the CHRNE expression may induce normal calcium flux. Subsequently, activation of AMPK via phosphorylation by calcium sensitive kinase CAMKK2 (Herzig and Shaw, 2018) may lead to prevention of cardiomyocyte injury in the atrial appendage, as shown in AMPK-mediated prevention of cellular injury (Zhao et al., 2013; Bi et al., 2015).

PRR32 is one of the proline-rich proteins that have been recognized as tanning-binding salivary proteins which counteracts a negative effect of a phenolic compound, tannin, in mammalian herbivores (Shimada, 2006). Although PRR32 expression was AA-specific in this current study, the role of PRR32 is yet to be investigated in the heart, and also the relationship between PRR32 and ESR1 has not been established according to the network analysis (Figure 9). Estrogen receptor 1 (ESR1) is a nuclear receptor of 17beta-estradiol and ESR1 activation in animal models exhibited anti-inflammatory effects (Puzianowska-Kuznicka, 2012). ESR1 shared one toxicity-related list, “Cardiac Hypertrophy,” with proinflammatory cytokine tumor necrosis factor (TNF) (Figure 9). Because a relationship between ESR1 and TNF has also not been established (Figure 9) but expression of TNF in the heart leads to cardiac pathogenesis, the role of ESR1 and PRR32 will need to be further investigated. Additionally, network analyses were conducted for three unreported common-specific/enhanced genes, and results were presented in Supplementary Figures S4–S6. Mechanistic insights on these unreported genes were provided in the “Discussion” section.