Expression array data and high throughput sequencing data regarding “COVID-19”, “SARS-CoV-2”, and “male infertility” were screened and downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The GSE164805 dataset contained transcriptome data of peripheral blood mononuclear cells taken from COVID-19 patients as well as healthy controls. The GSE160749 dataset contained transcriptome data of sperm samples taken from patients with male infertility and healthy controls. To further improve our data’s validity and authenticity, we used the GSE147507 dataset (COVID-19) and GSE26881 dataset (male infertility, including oligospermia, asthenozoospermia, and teratozoospermia) as independent validation sets.

To acquire the normalized expression matrixes, we used the R package “limma” (14) to convert the probe names, normalize the expression of each sample, and apply log2 transformation (if the original data was not transformed). For the RNA sequencing data (Count), we used the R packages “edgeR” (15) and “limma”. For all datasets, we set |log2 fold-change| (|log2FC|) > 0.5 and P value < 0.05 as the screening criteria to obtain CORGs. The display of DEGs was used by the R packages “ggplot2” (16) and “heatmap” (https://stat.ethz.ch/R-manual/R-devel/library/stats/html/heatmap.html). When we need to take the intersection of different gene sets, the R package “Venn diagram” (17) was utilized.

DAVID online tool (https://david.ncifcrf.gov/) was explored to provide functional classification, biochemical pathway maps, and conserved protein domain architectures implicated in the various physiological and pathological processes. Herein, we analyzed the significant signaling after uploading the CORGs. Three types of Gene Ontology (GO) terms, including biological process (BP), cellular component (CC), and molecular function (MF) were obtained (18). The requirements of screening the significant GO terms included gene number>2 and Ease < 1. For pathway analysis, we utilized the KEGG, REACTOME, and WikiPathway databases to acquire the most enriched pathways.

To comprehensively understand the interrelationship of CORGs, the protein-protein interaction (PPI) network construction was carried out using the Search Tool for the Retrieval of Interacting Genes (STRING) database (http://string-db.org). We removed CORGs with low connectivity according to the criteria of medium confidence > 0.4. The visualization of the PPI network was achieved with the help of Cytoscape software (version 3.7.1) (19). Module analysis was performed to obtain the gene clusters with more important values in the gene set of CORGs. We first calculated the module score by using the Molecular Complex Detection (MCODE) plugin in the Cytoscape software. Then, we selected the top three modules with the highest scores and downloaded the module composition of each module. GO and pathway analysis were used to obtain enriched biological processes and signaling alterations in each module using the Metascape online tool (http://metascape.org). Our cutoff in the Metascape was a P value < 0.01, a minimum count of 3, and an enrichment factor > 1.5.

To explore pharmacologic agent candidates that predominantly target and influence CORGs, we used the DrugBank database (Version 5.0) to predict the pharmacologic agent candidates. The interaction network of CORGs and candidates were constructed in the NetworkAnalyst online tool (https://www.networkanalyst.ca/) according to the connectivity of each pharmacologic agent candidate. The visualization of the network was implemented via the Cytoscape software.

With the use of the cytoHubba plugin in the Cytoscape software, different connectivity ranking calculation was applied to CORGs. The algorithms included MCC, EPC, BottleNeck, EcCentricity, Closeness, and Radiality. Only the significant CORGs with high scores were used for further investigation.

Receiver operating characteristic(ROC)analysis is a well-established statistical method that illustrates the diagnostic ability of variables and assesses the accuracy of model predictions. Herein, we used ROC analysis to determine the diagnostic ability of significant CORGs in the COVID-19 dataset, MI dataset, and validation sets. The R package “timeROC” (https://cran.r-project.org/web/packages/timeROC/) (20) was used to perform the ROC analysis. We first selected the CORGs with the area under the curve (AUC) > 0.9 in the COVID-19 and MI datasets. Subsequently, we performed ROC analysis once again in the two independent validation sets (GSE147507 and GSE26881), only the CORGs with AUC > 0.75 in the validation sets were regarded as hub CORGs.

To mine the immunological influence in the progression of COVID-19 and MI, we explored the CIBERSORT tool (21) to analyze the composition level of different immune cells. With the use of the R packages “e1071”, “parallel”, and “preprocessCore”, we obtained the ratio of immune cells in each sample in both the COVID-19 and MI datasets. The differences between the two groups (COVID-19 vs. control, and MI vs. control) were evaluated by Wilcoxon’s test. The relevance of immune cells and hub CORGs was analyzed using Spearman’s test. The threshold for statistical significance was set at P < 0.05.

We performed the network analysis of hub CORGs and transcription factor (TF) (the transcriptional level), and hub CORGs and miRNA (the post-transcriptional level). The TF-CORGs interactions and miRNA-CORGs interactions were obtained from the RegNetwork repository (22). The visualization of the above interactions was achieved by the Cytoscape software.

The interacting genes of hub CORGs were discovered using the web program GeneMANIA (http://genemania.org/) to create the hub CORGs-gene interaction network. These interacting genes should meet one or more of the following relationships: physical interaction, co-expression, predicted, co-localization, genetic interactions, pathway, and shared protein domains.

To explore which diseases hub CORGs are predominantly associated with, we analyzed the most enriched diseases for each gene in the Comparative Toxicogenomics Database (CTD, http://ctdbase.org/) (23). We focused on male reproductive system diseases, and the threshold for screening was set at an inference score > 40.

In our previous findings (24), we recruited a cohort of male COVID-19 patients and healthy controls attending Tongji hospital. After reviewing their medical records, we found seven healthy controls and seven patients with male infertility which manifested no pregnancy after unprotected sexual intercourse for more than one year, and different degrees of abnormal semen parameters (oligospermia, asthenozoospermia, and teratozoospermia) based on the WHO laboratory manual (25). For the diagnosis of COVID-19, real-time reverse transcriptase-polymerase chain reaction (RT-PCR) analysis of pharyngeal swab specimens was used to demonstrate that the patients were SARS-CoV-2 positive. Then, all patients were regarded as in the recovery state according to relieving symptoms and two consecutively negative testing for SARS-CoV-2. The New Coronavirus Pneumonia Prevention and Control Program (7th edition) issued by the National Health Commission of China, served as the basis for the diagnosis of COVID-19. The study received permission from Huazhong University of Science and Technology’s institutional research ethics committee. Each participant was provided with written informed consent, and all procedures followed the Helsinki Declaration.

The clinical and demographic information was recorded at the participant’s hospital visit. For patients, their semen examinations were performed at the first visit during the recovered period, and the recovered time for the seven patients was all ≤ 100 days. The serum Follicle−stimulating hormone (FSH), luteinizing hormone (LH), total T, and estradiol (E₂) concentrations were measured using a chemiluminescence immunoassay (Beckman Coulter, Fullerton, USA). Following the three to seven days of advised abstinence, freshly obtained semen was processed for 30 to 60 minutes at room temperature to liquefy it. Sperm examination was performed in strict accordance with the WHO laboratory manual (25). Purified sperm samples were obtained by centrifuging them for 20 minutes at 1000 g while using a Percoll gradient (Solarbio, Beijing, China). After washing them with PBS twice, the purified sperm samples were used for subsequent experiments.

Following the manufacturer’s instructions, we extracted RNA from purified sperm samples using RNA-easyTM Isolation Reagent (Vazyme, Nanjing, China) and synthesize cDNA using the Hifair^® II 1st Strand cDNA Synthesis Kit (Yeasen, Shanghai, China). Herein, we diluted the cDNA to the proper concentrations, and performed quantitative real-time polymerase chain reaction (qRT-PCR) as previously described (26). Expression levels of hub CORGs were calculated after normalization of the reference gene GAPDH. The detailed information on primers used in the study is shown in Supplementary Table 1 .

We retrieved the immunohistochemistry data from the Human Protein Atlas database (HPA) (https://www.proteinatlas.org/), and accessed the positive region of hub CORGs to identify the expression signature in testes.

The clinical and demographic data of participants and mRNA levels of hub CORGs were presented as mean ± standard deviation (SD). The comparisons were performed using Student’s t test in built-in algorithms of Graphpad software (LLC, San Diego, California USA). The threshold for statistical significance was set at P < 0.05.

The study strategy was depicted in Figure 1 . With the threshold of log2FC > 0.5, we obtained the DEGs taken from the COVID-19 dataset and MI dataset. The results showed that there were 944 upregulated DEGs and 312 downregulated DEGs between the COVID-19 and healthy control samples ( Figures 2A, C ), and there were 9372 upregulated DEGs and 8667 downregulated DEGs between the MI and control samples ( Figures 2B, D ). To obtain the common genes implicated in both the disease of COVID-19 and MI, we took the intersection of COVID-19 DEGs and MI DEGs, and a total of 460 overlapped CORGs were discovered in both DEGs sets ( Figure 2E ).

To figure out which signaling the CORGs predominantly participated in, we used four databases to perform GO analysis (GO database) and pathway analysis (KEGG, REACTOME, and WikiPathways databases). As shown in Figure 3A , for biological process, GO terms were mainly enriched in protein ubiquitination, cell division, cellular response to DNA damage stimulus, apoptotic process, cell cycle, and DNA repair. GO terms for cellular component were tightly related to cytosol, nucleus, cytoplasm, and nucleoplasm. GO terms for molecular function were focused on protein binding, ATP binding, RNA binding, GTPase activator activity, and ubiquitin-protein transferase activity.

In another three pathway databases, the results showed the overlapped CORGs were notably enriched in DNA damage and repair-associated (DNA repair pathways/full network, DNA IR-damage and cellular response via ATR, DNA Damage Bypass, Recognition of DNA damage by PCNA-containing replication complex, and DNA Repair), cell cycle-associated (Cell Cycle, Cell Cycle Mitotic, Cell Cycle Checkpoints), ubiquitination-associated (Ubiquitin mediated proteolysis, Antigen processing: Ubiquitination & Proteasome degradation), and coronavirus-associated (SARS-CoV Infections, Viral Messenger RNA Synthesis, and Host-pathogen interaction of human coronaviruses-autophagy) pathways ( Figures 3B–D ). Combining the above results, DNA damage and repair-associated, cell cycle-associated, and ubiquitination-associated signaling were most enriched, and these signaling may play an important role in the initiation and development of COVID-19 and MI.

In the next step, we aimed to investigate the interaction features of CORGs by using the STRING database. A total of 738 interactions were identified as shown in Figure 4A . Among all the interactions, we found that three genetic modules had the top highest density score, indicating the stronger level of interaction of CORGs in the three modules ( Figure 4B ). The enriched GO functions for Module one were involved in intraflagellar transport, cilium assembly, and intraciliary transport. GO functions for Module two had a close relationship with meiotic cell cycle and DNA biosynthetic processes. GO functions for Module three captured purine metabolism and pyrimidine metabolism ( Figures 4C–E ). The functional annotation confirmed the role of cell cycle and DNA metabolism once again.

To explore the pharmacologic agent candidates which may target the overlapped CORGs, we utilized prediction algorithm in the DrugBank database, and we harvested 135 pharmacologic agent candidates and 137 interactions with CORGs ( Figure 5A ). Among these candidates, we mined their therapeutic effects in Pubmed. Only four candidates (ascorbic acid, biotin, caffeine, and L-cysteine) had favorable effects on male infertility with in vitro or in vivo evidence, and their two-dimensional structures were shown in Figure 5B . It was implied that the four pharmacologic drugs have a good medical value in COVID-19-associated MI.

Considering the large number of CORGs, we used connectivity ranking analysis and ROC analysis to narrow down the gene number. After implementing six connectivity ranking algorithms, we found 153 CORGs with high scores in Figure 5C . Moreover, ROC analysis showed that there were 100 CORGs with AUC > 0.9 in the COVID-19 datasets, and there were 18 CORGs with AUC > 0.9 in the MI datasets ( Supplementary Table 2 ). A total of twelve significant CORGs (RNF7, UBE2G1, ATG7, AKT3, STIL, PDE3B, ENTPD6, EIF3B, ZEB1, RBP4, C2orf69, CIB1) showed AUC > 0.9 in both datasets, suggesting that these genes have higher research potentials ( Figure 5D ).

To add authority and persuasion of significant CORGs, we examined the diagnostic value of twelve CORGs in other independent transcriptome data. For further validation, we first obtained 3129 DEGs between the COVID-19 and control samples in the COVID-19 dataset (GSE147507). Similarly, 441 DEGs were identified between the MI and control samples in the MI dataset (GSE26881) ( Supplementary Figure 1 ). Next, we performed ROC analysis to validate the diagnostic effectiveness of CORGs in the above validation sets. As shown in Figure 6 , it was revealed that three hub CORGs (ENTPD6, CIB1, and EIF3B) had good diagnostic performance to discriminate COVID-19/MI and control samples (AUC > 0.75). These hub CORGs were well verified in external transcriptome data and may have the potential as biomarkers for COVID-19-associated MI.

Activation of immune response is inevitable in patients with COVID-19, thus, we mined the immune infiltration landscape in COVID-19 and MI based on genomic data. We first focused on the difference in immune composition between the patient and control groups. It is found that there were significant differences in the levels of CD8+ T cells, monocytes, macrophages M0, dendritic cells activated, mast cells activated, eosinophils, and neutrophils ( Figure 7A ). Furthermore, the level of regulatory T cells (Tregs) in the MI group was remarkedly different from the control group ( Figure 7B ). In the subsequent step, we analyzed the correlation degree between the immune cells and hub CORGs. Correlation results indicated that hub CORGs (ENTPD6 and CIB1) were positively related to the level of significant immune cells, while hub CORG (EIF3B) was negatively related. The detailed interrelationship was as follow: the level of CD8+ T cells (ENTPD6, r = 0.7, p = 0.0049; CIB1, r = 0.74, p = 0.0022; EIF3B, r = 0.86, p < 2.2 e^-16), the level of monocytes (ENTPD6, r = 0.54, p = 0.041; CIB1, r = 0.65, p = 0.01; EIF3B, r = 0.74, p = 0.0025), and the level of macrophages M0 (ENTPD6, r = -0.61, p = 0.018; CIB1, r = -0.7, p = 0.0049; EIF3B, r = -0.66, p = 0.0075) ( Figure 7C ). These data suggested that hub CORGs may be implicated in affecting the sperm parameters in COVID-19 patients via disturbance of CD8+ T cells, monocytes, and macrophages M0.

In this part, we investigated the upregulation and downregulation of hub CORGs. Therefore, we characterized the profile of interacting TFs and miRNAs. In Figure 8A , 24 TF-CORGs interactions and 13 miRNA-CORGs interactions were shown. Among all the interacting TFs, ARNT, EGR1, MAX, MYC, and USF1 were predicted to target at least two hub CORGs. For interacting miRNAs, 13 miRNAs were shown to have the potential to regulate the post-transcriptional level of hub CORGs, including hsa-miR-24, hsa-miR-181b, hsa-miR-181c, hsa-miR-130b, hsa-miR-485-5p, hsa-miR-497, hsa-miR-181d, hsa-miR-588, hsa-miR-608, hsa-miR-645, hsa-miR-423-5p, hsa-miR-874, and hsa-miR-923, which deserved further exploration in COVID-19 and MI.

With the use of the GeneMANIA program, we constructed the hub CORGs-gene interaction network ( Figure 8B ). Moreover, we found that this network was tightly relevant to the viral translation, viral gene expression, and site of polarized growth (false discovery rate < 0.05). It suggested that hub CORGs may mediate cellular viral replication and damage. For disease association, three hub CORGs strongly impacted diseases of the male reproductive system. Oligospermia, testicular diseases, and male infertility were highlighted in the diseases of the male reproductive system with an inference score > 40 ( Figure 8C ), which suggested the close association of hub CORGs in MI-associated illnesses.

As shown in Supplementary Table 3 , seven patients with male infertility had varying degrees of oligospermia, asthenozoospermia, and teratozoospermia. During their recovery period, they had normal serum hormone levels and no medical history of varicocele, genitourinary tract inflammation and infection, scrotal trauma, and so on. In comparison to the control group, the sperm concentration, motility, and progressive motility were significantly lower in the patient group (P < 0.05, Figures 9A–C ).

To enhance authority and stringency, qRT-PCR was utilized to measure the mRNA levels of hub CORGs between the COVID-19 patients with MI and healthy controls. In Figures 9D–F , the results of qRT-PCR suggested that notable changes were observed in the mRNA levels of CIB1 and EIF3B in the patient group than those in the control group (P < 0.01 and P < 0.001, respectively). However, no pronounced alteration was observed in the mRNA of ENTPD6 between the two groups. Furthermore, immunohistochemical data showed that it is hard to detect the positive area of CIB1 in seminiferous tubules and interstitial region of testes. As for EIF3B, moderate positive blots were seen in spermatogenic and Sertoli cell in seminiferous tubules and Leydig cell in interstitial region ( Supplementary Figure 2 ). It is noted that two hub genes (CIB1 and EIF3B) had the significant change pattern in the sequencing findings and qRT-PCR results, emphasizing the critical genetic links of CIB1 and EIF3B between COVID-19 and MI.