Endometrial gene expression differences in women with coronavirus disease 2019.

The profiling of global endometrial gene expression revealed altered endometrial gene expression in most of the patients with COVID-19. Among the 24 initial endometrial biopsy samples, 5 out of the 14 samples from the COVID-19 group had low RNA integrity (DV 200 <30%) and were excluded from the RNA-seq analysis, whereas 1 out of the 10 samples from the control group was excluded because of low quality ( Supplemental Table 1 ). Thus, 18 endometria samples were sequenced as detailed in Figure 1 A. The endometrial biopsy samples from the patients with COVID-19 and controls had similar numbers of good-quality reads per sample ( Supplemental Fig. 1 , available online). The subsequent analysis of read quality in all the samples revealed 16,546 good-quality measured transcripts. Exploratory transcriptomic analyses of the 18 endometrial biopsy samples revealed that age, corrected menstrual phase, and RNA integrity did not affect endometrial gene expression. After removing 2 outliers identified using PCA ( Fig. 1 B), the samples were segregated into 2 groups based on gene expression related to COVID-19 ( Fig. 1 C). Two patients with COVID-19 (named COV2 and COV11), exceptionally clustered with the controls, were considered to have a “control-like” transcriptomic behavior and were removed from subsequent analysis because their endometria appeared transcriptomically unaffected by COVID-19 despite the other symptoms that they presented with. These findings led us to assume that COVID-19 variably affected this cohort (in terms of affected organs and manifested symptoms), and we estimated that 75% of our recruited symptomatic patients with COVID-19 had altered endometrial gene expression profiles. Our DEA revealed 235 DEGs (163 up- and 72 down-regulated) in the altered endometria of the patients with COVID-19 ( Supplemental Table 2 ). The top 10 up-regulated DEGs exhibited over fivefold increases (e.g., UTS2 exhibited a 16.126-fold increase), whereas the changes in the down-regulated genes ranged from −3- to −5.165-fold ( Supplemental Table 2 ). Notably, only 3 down-regulated DEGs passed our |FC| cutoff of ≥5, whereas 13 up-regulated genes met this criterion ( Fig. 2 ), suggesting that COVID-19-associated dysregulation of endometrial genes is mediated through increased rather than decreased expression. FIGURE 2 Results of the gene set enrichment analysis. A heatmap depicts 13 out of 235 significantly differentially expressed genes with a false discovery rate of 5. DEG = differentially expressed gene; FC = fold change. Results of the gene set enrichment analysis. A heatmap depicts 13 out of 235 significantly differentially expressed genes with a false discovery rate of 5. DEG = differentially expressed gene; FC = fold change. Next, to validate the expression of DEGs in the RNA-seq results, we analyzed the gene expression of 2 up-regulated ( UTS2 , FC = 8.871; IFI6 , FC = 3.740) ( Fig. 3 A) and 2 down-regulated ( CPZ , FC = −1.866; CDH3 , FC = −1.511) ( Fig. 3 B) DEGs using qRT-PCR. We additionally validated the expression of IFIH1 (FC = 1.601), involved in defense response to symbionts, positive regulation of type I interferon (IFN-I) production, and the COVID-19 KEGG pathway ( Fig. 3 A and Supplemental Fig. 2 , available online); BNIP3 (FC = 1.601), implicated in response to viral infections ( Fig. 3 A); and IRF4 (FC = −1.786), described in the regulation of anti-inflammatory interleukin 4 production ( Fig. 3 B). Comparison with the RNA-seq values is presented in Figure 3 C. FIGURE 3 RNA sequencing validation using quantitative reverse-transcription polymerase chain reaction. The graphs show gene expression, presented as mean fold changes ± standard deviations, between endometria of the coronavirus disease 2019 and control groups for 4 up-regulated ( A ) and 3 down-regulated ( B ) differentially expressed genes. ( C ) A comparison between the RNA sequencing and quantitative reverse-transcription polymerase chain reaction results is also shown. ∗ P <.050. COVID-19 = coronavirus disease 2019; DOWN = down-regulated; qRT-PCR = quantitative reverse-transcription polymerase chain reaction; RNAseq = RNA sequencing; UP = up-regulated. RNA sequencing validation using quantitative reverse-transcription polymerase chain reaction. The graphs show gene expression, presented as mean fold changes ± standard deviations, between endometria of the coronavirus disease 2019 and control groups for 4 up-regulated ( A ) and 3 down-regulated ( B ) differentially expressed genes. ( C ) A comparison between the RNA sequencing and quantitative reverse-transcription polymerase chain reaction results is also shown. ∗ P <.050. COVID-19 = coronavirus disease 2019; DOWN = down-regulated; qRT-PCR = quantitative reverse-transcription polymerase chain reaction; RNAseq = RNA sequencing; UP = up-regulated. The functional enrichment analysis was performed to identify underlying biologic functions associated with the significantly dysregulated DEGs in the endometria of the patients with COVID-19. We identified 161 enriched, dysregulated functions (36 up- and 125 down-regulated) in the endometria of the patients with COVID-19 ( Fig. 4 and Supplemental Fig. 3 , available online), involving a myriad of biologic processes (15 up- and 52 down-regulated), molecular functions (12 down-regulated), cellular components (12 up- and 7 down-regulated), and/or signaling pathways (9 up- and 54 down-regulated). Specifically, the endometrial transcriptomes obtained from the patients with COVID-19 were enriched in biologic processes associated with response to viral infections, positive regulation of IFN-I production, the COVID-19 KEGG pathway ( Supplemental Fig. 2 , available online), NET formation ( Fig. 4 ), and ribosome-guided messenger translation. Further, the DEGs were mainly enriched in the mitochondria and cilia ( Fig. 4 and Supplemental Fig. 3 ). Remarkably, the molecular functions related with the immune system (i.e., T helper 17 cell differentiation as well as T-cell and platelet activation), anti-inflammatory cytokine production (i.e., regulation of interleukin 4 production), and gene transcription (i.e., DNA-binding transcription activator activity and transcription regulatory region nucleic acid binding) were down-regulated ( Fig. 4 and Supplemental Fig. 3 ). FIGURE 4 Functional enrichment analysis. The figure shows the main functions that were differentially up- and down-regulated (false discovery rate<0.05) in the endometria of patients with COVID-19 (compared with those in the control group; COV11 and COV2 were considered “control-like” samples). Both the Kyoto Encyclopedia of Genes and Genomes pathways and Gene Ontology biologic processes are included. BMP = bone morphogenetic protein; COVID-19 = coronavirus disease 2019; GF = growth factor; IFN = interferon; IL = interleukin; NES = normalized enrichment score; P13k-Akt = phosphatidylinositol 3-kinase and protein kinase B; PCR = polymerase chain reaction; Ras = rat sarcoma virus; Th = T helper. Functional enrichment analysis. The figure shows the main functions that were differentially up- and down-regulated (false discovery rate<0.05) in the endometria of patients with COVID-19 (compared with those in the control group; COV11 and COV2 were considered “control-like” samples). Both the Kyoto Encyclopedia of Genes and Genomes pathways and Gene Ontology biologic processes are included. BMP = bone morphogenetic protein; COVID-19 = coronavirus disease 2019; GF = growth factor; IFN = interferon; IL = interleukin; NES = normalized enrichment score; P13k-Akt = phosphatidylinositol 3-kinase and protein kinase B; PCR = polymerase chain reaction; Ras = rat sarcoma virus; Th = T helper.

A total of 24 endometrial samples were prospectively collected for this study. Fourteen biopsy samples (COVID-19 group) were collected from patients hospitalized for COVID-19 at Hospital Universitari i Politècnic La Fe (Valencia, Spain) ( 17 ); these patients were symptomatic (patients with mild-to-severe symptoms were included, whereas asymptomatic or very severe cases were not recruited) and had a positive test result (cycle threshold<40) for SARS-CoV-2 infection, as indicated by quantitative reverse-transcription polymerase chain reaction (qRT-PCR) of nasopharyngeal swab samples, as previously described ( 17 ). The other 10 endometrial samples (control group) were derived from patients admitted to the same hospital who underwent hysteroscopic procedures for benign gynecologic disorders and tested negative for COVID-19 based on qRT-PCR of nasopharyngeal swab samples. Notably, the different gynecologic disorders in these control patients did not change the results between the groups because of normal endometrial tissue diagnoses in all the cases. The patients’ characteristics (including age and menstrual cycle stage at the time of endometrial biopsy sample collection, menstrual cycle length, and COVID-19 symptoms) are presented in Supplemental Table 1 (available online). Samples with good RNA quality (percentage of RNA fragments >200 nucleotides [DV 200 ] >30%) were analyzed using RNA-seq for the detection of changes in gene expression using an untargeted approach that compared both the groups (COVID-19 and control). Raw sequencing data were preprocessed, normalized, and functionally interpreted using bioinformatic pipelines that were used to identify disease-altered genes and their functions in the endometria. The study design is depicted in Figure 1 A. FIGURE 1 Study design. ( A ) The workflow indicates the size of each group throughout the study and the reasons for excluding specific samples. ( B ) An initial principal component analysis showed how the 2 experimental groups overlapped with the presence of 2 outliers. ( C ) After removal of these outliers, the samples were evidently segregated according to the 2 experimental groups: coronavirus disease 2019 and control. This principal component analysis also identified 2 samples from the coronavirus disease 2019 group that behaved similar to the ones from the control group. ∗Outliers ∗∗Control-like samples. COVID-19 = coronavirus disease 2019; DEA = differential expression analysis; GSEA = gene set enrichment analysis; PC = principal component; PCA = principal component analysis. Study design. ( A ) The workflow indicates the size of each group throughout the study and the reasons for excluding specific samples. ( B ) An initial principal component analysis showed how the 2 experimental groups overlapped with the presence of 2 outliers. ( C ) After removal of these outliers, the samples were evidently segregated according to the 2 experimental groups: coronavirus disease 2019 and control. This principal component analysis also identified 2 samples from the coronavirus disease 2019 group that behaved similar to the ones from the control group. ∗Outliers ∗∗Control-like samples. COVID-19 = coronavirus disease 2019; DEA = differential expression analysis; GSEA = gene set enrichment analysis; PC = principal component; PCA = principal component analysis. The research ethics committee of Hospital Universitari i Politècnic La Fe (Valencia, Spain) approved this study (registration number: 2020-268-1) on May 12, 2020 (addendum: November 18, 2020). All study participants received written and oral information on the characteristics of the study and signed informed consent forms (or consented orally in the presence of at least 1 witness in the case of patients with COVID-19). All endometrial biopsy samples (from both the COVID-19 and control groups) were collected via endometrial aspiration using an intrauterine cannula (Cornier Pipelle, ref. 4164, Gynetics, Lommel, Belgium) and preserved in RNAlater (ref. AM7021; Invitrogen, Vilnius, Lithuania) until cryopreservation. After the removal of RNAlater, the endometrial biopsy samples (median size, 3 mm 2 ) were frozen at −80°C until their analysis. Total RNA was extracted from the thawed endometrial biopsy samples using the miRNeasy Mini Kit (ref. 1038703; Qiagen, Hilden, Germany) according to RNeasy FFPE Handbook protocol instructions (Qiagen, Hilden, Germany). The amount of total RNA isolated from each tissue was quantified using the NanoDrop One spectrophotometer (ThermoFisher, Vilnius, Lithuania), and the integrity of the extracted RNA was assessed using the Fragment Analyzer System (Agilent, Santa Clara, CA). Nine endometrial samples from each group (9 from the COVID-19 group and 9 from the control group) with the highest DV 200 values were selected for RNA-seq analysis. Sequencing libraries were constructed using the TruSeq stranded messenger RNA protocol (Illumina, San Diego, CA), with poly-A selection for ribosomal messenger RNA depletion, followed by chemical fragmentation and complementary DNA synthesis, as per the manufacturer’s instructions. The libraries were sequenced in 2 runs using the NextSeq 500 platform (Illumina, San Diego, CA) at a single-end read length of 75 nucleotides using the High Output 75 cycles kit (Illumina, San Diego, CA). Raw data were trimmed using fastp (version 0.21) ( 19 ) and evaluated for read quality using the FastQC, version 0.11.5, tool ( 20 ). STAR algorithm ( 21 ) was employed for mapping data quality, using GRCh38 as the reference human genome, and obtaining raw counts. Genes with zero counts in all measured samples were excluded, whereas low-count genes were filtered using the NOISeq R-package ( 22 ) in R (version 4.0.5 [March 31, 2021]) ( 23 ). The remaining data were normalized using the Remove Unwanted Variation methodology ( 24 ) implemented in the RUVseq R-package ( 25 ). Briefly, this approach uses general linear model regression to remove unknown sources of data or variation from batch or other variable effects to improve the detection of differentially expressed genes (DEGs). We then corrected for the menstrual cycle phase to remove transcriptomic variability associated with the progression of the menstrual cycle ( 26 ). Finally, the principal component analysis (PCA) was performed to evaluate any remaining batch effects, look for outliers, and find transcriptomically related samples in our cohort. To identify genes in the endometrium that are most affected by COVID-19, the differential expression analysis (DEA) was performed using the limma R-package ( 27 ). DEGs that had an adjusted P value of <.05 and a false discovery rate of <0.05 were considered to have been significantly affected by COVID-19 and were subjected to a subsequent functional analysis using the Gene Ontology (version 2021-07-02, experimental and propagated annotation; to elucidate the affected biologic processes, molecular functions, and cellular components) ( 28 ) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases (release 99.0, August 1, 2021) ( 29 ). The genes were arranged in descending order according to their fold change (FC), which was obtained from DEA, and the gene set enrichment analysis was implemented in the clusterProfiler R-package ( 30 ). Only significantly enriched functions (false discover rate<0.05) were considered to be affected by COVID-19. All plots were generated using the ggplot2 package in R ( 31 ). To validate the DEGs identified using the RNA-seq analysis, the complementary DNA from the COVID-19 (n = 9) and control (n = 9) samples was analyzed using qRT-PCR with Power-Up SYBR Green (ref. A25742 ; ThermoFisher, Vilnius, Lithuania) on the StepOnePlus system (Applied Biosystems, Carlsbad, CA). Based on their high |FC | values ( Supplemental Table 2 , available online), we assessed UTS2 and IFI6 as up-regulated DEGs and CPZ and CDH3 as down-regulated DEGs. Further, based on their involvement in response to viral infections or inflammatory processes, we assessed IFIH1 , BNIP3, and IRF4 . The relative level of gene expression or FC was determined using the ΔΔCt method and normalized to the level of the expression of the housekeeping gene GAPDH . The primer sequences (Integrated DNA technologies, Leuven, Belgium) are listed in Supplemental Table 3 (available online). The qRT-PCR results are expressed as means ± standard deviations. We used the unpaired t test and Wilcoxon test for normally and nonnormally distributed data, respectively, to identify significant DEGs between the COVID-19 and control groups. GraphPad Prism 8.3 (San Diego, CA) was used for statistical analysis and visualization. In all the analyses, P <.05 was considered statistically significant.

In conclusion, our transcriptomic results revealed that human endometrial gene expression and, possibly, endometrial function are affected by COVID-19. The changes in the endometria of patients with COVID-19 seem to occur in the absence of direct viral infection, along with the down-regulation of other vital functions and host immune responses. Still, more studies are needed to elucidate whether these molecular changes also occur in patients with aggressive COVID-19 symptoms and/or asymptomatic women. Finally, evaluating whether these changes persist after disease remission or resolve in the following menstrual cycle(s) as the endometrium cyclically regenerates itself would be of high clinical interest.

Although our sample size was limited, this was the first study to analyze changes in global gene expression in the endometrium of women with COVID-19 during the first 10 months of the global pandemic (2020). It has been widely reported that COVID-19 produces heterogenic clinical phenotypes in terms of severity and complications (involving only pulmonary or additional renal, hepatic, and/or even cardiac pathophysiologies). The presentation of the disease also varies with age, with older patients having a higher prevalence of severe disease and risk of complications ( 32 , 33 ). Consistent with COVID-19 being a complex, multiorgan disease, this study highlighted that this disease can also impact the endometrium, at least at the level of gene expression. Indeed, we found that most (75%) of the patients with COVID-19 included in this study presented with an altered endometrial gene expression profile, revealing that the endometrium is commonly affected by this disease. In addition, using a nontargeted approach, we identified key components of the pathogenesis of COVID-19 ( Supplemental Fig. 2 ) among the most altered DEGs in the affected endometria. Overall, this study showcased that this systemic disease affects the normal endometrial gene expression profile in the absence of direct viral infections ( 17 ). By correcting for menstrual cycle biases, we overcame the limitation of unequal representation of the different menstrual phases in our COVID-19 and control groups and ensured that the gene dysregulations that we found were truly related to COVID-19 (empowering the results of our low-sample-size study). The latter was verified using PCA, which additionally corroborated that age did not result in a bias in our results. Other groups have observed some of the same genes and patterns during SARS-CoV-2–mediated disease. For example, NET formation is enhanced in various tissues of patients with COVID-19 because of severe inflammatory reactions ( 6 , 34 , 35 ) and is associated with collateral tissue injury during respiratory viral infections. In this study, we found that the endometrium of the women with COVID-19 exhibited remarkable up-regulation of genes associated with NET formation, suggesting that COVID-19 causes collateral damage to the endometrium as well. In terms of other host immune responses, we also found dysregulated IFN-I production in the endometria of the patients with COVID-19. This is consistent with previous studies that reported that IFN-I production is impaired in such patients ( 36 ). Interestingly, specific interferon-stimulated genes appear to be up-regulated in the peripheral blood of patients with COVID-19 ( 36 ). These implicated genes fall into 2 subgroups: those dysregulated in different viral infections (including SARS-CoV-2, severe acute respiratory syndrome coronavirus, or influenza) and those uniquely dysregulated in COVID-19. The second subgroup includes IFI2 , IFI6 , and IFI44 ( 36 ), which were also up-regulated in the endometria of our patients, indicating an altered inflammatory environment. Notably, both the IFI27 gene and the IFN-I signaling pathway are up-regulated in the endometria of women with recurrent pregnancy loss ( 37 ), and IFI6 and IFI44 are overexpressed in patients with implantation failure ( 38 ). Indeed, the entire reproductive system can be altered in response to an inflammatory environment, as has been reported for premature ovarian insufficiency ( 39 ). The other observed changes that we can associate with inflammatory processes include the remarkable 16.126-fold up-regulation of UTS2 in the endometria of patients with COVID-19. This gene encodes a peptide that acts as a vasoconstrictor and is involved in the pathogenesis of inflammatory diseases such as atherosclerosis and hypertension ( 40 ). Interestingly, a recent meta-analysis revealed that several preeclampsia-related pathways were affected in pregnant women with COVID-19, corresponding with an imbalance of vasoactive peptides resulting from the up-regulation of various genes, including UTS2 ( 41 ). Further, studies of other systemic diseases, such as lupus erythematosus ( 42 ), several rheumatic diseases ( 43 ), irritable bowel syndrome ( 44 ), and diabetes ( 45 ), revealed dysregulated function of the female reproductive system. Interestingly, these studies also support the fact that the endometrium is influenced by cytokine storms (triggered by viral infections) rather than by the virus itself. Changes in lipid metabolism are also consistent with SARS-CoV-2 infection. We found that the endometrial expression of genes involved in this process was altered in the patients with COVID-19. Despite not being directly linked to inflammation, this correlates with the differential lipid profile observed in the serum of patients with COVID-19, which, in turn, seems to be associated with the severity of the disease ( 46 , 47 ). We found the expression of genes with direct roles in endometrial function to be altered by COVID-19. In particular, AGR3 was among the top 10 up-regulated DEGs in the endometria of the patients with COVID-19, and the protein it codes for has been shown to be overexpressed in the endometrial ciliated luminal epithelium of women with recurrent implantation failure and recurrent pregnancy loss ( 48 ). Alterations in ciliated cells also occur in other reproductive disorders such as endometriosis ( 49 , 50 ). Further, these reports describe the alteration of ribosome and phosphorylation processes during an endometrial pathology ( 49 ), supporting our findings of these 2 mechanisms being up-regulated in our group of patients with COVID-19. The up-regulation of genes associated with reproductive disorders indicates that COVID-19 could negatively affect endometrial function, thereby potentially disrupting menstrual cycles and consequently compromising embryo implantation and fertility. However, to address this hypothesis, it would be necessary to investigate whether proinflammatory reactions (that can alter the normal function of the female reproductive system) ( 51 ) subside after the infection or persist for a short or even long term ( 52 ). Nevertheless, a more detailed dataset from a larger cohort that takes into account other factors that cause irregular menstrual cycles (i.e., baseline ovarian reserve and hormone measurements for the given cycle as well as follow-up of the menstrual cycle) is required to refine the correlation. Because alterations in the menstrual cycle have previously been described in other systemic disorders ( 42 , 43 , 44 , 45 ), it could be worthwhile to conduct an in silico analysis that compares our data with their transcriptomes to elucidate whether COVID-19 directly produces unique genetic signatures or indirectly induces molecular changes with the general inflammation it elicits. Although this study emphasized the up-regulation of genes and their pathways because of their high |FC| values, the down-regulated processes were also noteworthy. For example, smell and taste disorders are common in patients with COVID-19 ( 53 ), indicating the down-regulation of the taste transduction pathway, as shown in our analysis. Further, CPZ, a carboxypeptidase involved in the Wnt signaling pathway ( 54 ), was among the most down-regulated DEGs and pathways in the COVID-19–affected endometria. The Wnt pathway has well described roles in multiple endometrial processes and is critical for implantation ( 55 ). Two other down-regulated DEGs, CDH3 and IRF4 , are involved in cell adhesion during embryo implantation ( 56 ) and decidualization ( 57 ), respectively. Taken together, the down-regulated DEGs we identified in this study also suggest that COVID-19 could alter endometrial function. This study provides novel transcriptomic findings elucidating the potential impact of COVID-19 on the endometrium. Our bioinformatic approach for processing the global endometrial transcriptomes was concise and powerful, highlighting our findings despite the small sample size. Nevertheless, our small cohort (resulting from the risk of collecting high-quality samples from patients with the infection) implied several limitations, such as interpatient variability in terms of menstrual cycle parameters, fertility status, day of biopsy sample collection, and age, among others. By eliminating the effect of 2 of these confounding variables (menstrual phase and gynecologic disorders), our endometrial COVID-19 signature became more reliable. In addition, because of health measures and the risk of infection, we only collected biopsy samples from symptomatic, hospitalized patients (asymptomatic women were not recruited). Finally, regarding the 2 endometria with control-like expression profiles from the COVID-19 group, we carefully reviewed the medical records of COV11 and found that she might have had a false-positive result in her SARS-CoV-2 diagnostic test. Interestingly, the 2 patients with control-like expression profiles were the youngest in the cohort. Because of our small sample size, we cannot draw conclusions based on age at this point, especially considering that all the patients reported respiratory symptoms. However, age has been shown to be correlated with the severity of COVID-19 ( 58 , 59 ), and in this regard, it might be possible that the endometrium of younger patients is not affected to the same degree. Overall, additional research with a larger sample size would be needed to further validate the results of our RNA-seq and functional analyses and advance the understanding of the implications of COVID-19 on the human endometrium.