Co-expression of the SARS-CoV-2 entry molecules ACE2 and TMPRSS2 in human ovaries: Identification of cell types and trends with age.

This work was supported by grants from the 10.13039/501100001809 National Natural Science Foundation of China (No. 82001514, 81671394 and 81873824) and the Clinical Research Pilot Project of Tongji Hospital , 10.13039/501100003397 Huazhong University of Science and Technology (No. 2019CR205).

We first searched the online datasets of the GTEx project and then analyzed the expression level of ACE2 and TMPRSS2 in ovary. As lung is known to be the main target of SARS-CoV-2, we compared the ACE2 and TMPRSS2 expression between lung and ovary; there was no significant difference in the mRNA levels of ACE2 and TMPRSS2 between the ovaries and lungs ( Fig. 1 A, B). We accessed the scRNA-Seq data for macaque ovaries and human follicles in the GEO database ( Table 1 ) to further investigate the cell type-specific expression of ACE2 and TMPRSS2 . As M. fascicularis and humans share similar genetic and physiological features [ 24 ], the results could also be suggestive. In GSE130664 , eight clusters of cells in ovaries of M. fascicularis were observed ( Fig. 2 A). Overall, approximately 13.31% of all 2601 cells expressed ACE2 ; the highest proportion of these cells were oocytes (95.58%), and other cell types represented relatively lower proportions ( Table 2 ). However, the expression profile of TMPRSS2 indicated that 10.02% of all cell types expressed TMPRSS2, including 68.81% oocytes and other cell types ( Table 2 ). The violin plot ( Fig. 2 B), feature plot ( Fig. 2 C) and dot plot ( Fig. 2 D) further demonstrated that  ACE2  and  TMPRSS2  were mainly expressed in oocytes. The cluster of oocytes was confirmed by the co-expression of DDX4 and GDF9 ( Fig. 2 E), the canonical markers of oocytes, indicating a similar pattern as ACE2 . Another important cell in the ovary, granulosa cells with the marker of AMH and INHA ( Fig. 2 F), but did not express the ACE2 and TMPRSS2 ( Fig. 2 D). Similarly, the expression of ACE2 and TMPRSS2 was observed to be extremely low in the subset of macrophages and natural killer (NK) cells ( Fig. 2 D). Fig. 1 mRNA expression levels of ACE2 (A) and TMPRSS2 (B) in multiple organs; data obtained from the Genotype-Tissue Expression (GTEx) database (a human organ database); * p  < 0.05, ** p  < 0.001, *** p  = 0.0001, **** p  < 0.0001. The vertical coordinates represent the log transformed (log2 [x + 0.001]) relative expression level of ACE2 and TMPRSS2 . Fig. 1 Fig. 2 Gene expression profiles of distinct ovarian cell subpopulations based on ovarian scRNA-Seq data from dataset GSE130664 (macaque ovarian single cells). (A) T-distributed stochastic neighbor embedding (t-SNE) clustering of adult ovarian cells into eight clusters based on the scaled data. (B) Average expression level of ACE2 (upper) and TMPRSS2 (lower) in each cell cluster. The vertical coordinates represent the log-transformed [ log (x + 1)] relative expression level of ACE2 and TMPRSS2 . (C) Scatter plot of ACE2 and TMPRSS2 gene expression of each cell cluster. The ellipse indicates the cell cluster of oocytes that co-expressed ACE2 and TMPRSS2 . The number and gradient color of the legend bars show the log-transformed values of gene expression matrix raw data of each cell cluster. (D) Dot plot of the well-known markers including ACE2 and TMPRSS2 in each cell cluster. The color of the legend bar of average expression reflects the mean expression. The size of the dots corresponds to the percentage of cells expressing the listed genes in each cell cluster. (E) Scatter plots of cell clusters with canonical markers of oocytes and (F) granulosa cells. Fig. 2 Table 2 The distribution of ACE2 and TMPRSS2 in each type of ovarian cells in two scRNA-seq. Table 2 ACE2 TMPRSS2 GSE130664 (%) GSE107746 (%) GSE130664 (%) GSE107746 (%) Oocytes 95.58 86.42 68.81 16.05 Granulosa cells 2.55 67.14 3.19 10.00 meiotic prophase germ cells 0.31 – 0.00 – Stromal cells 0.58 – 0.58 – Smooth muscle cells 2.84 – 4.26 – Endothelial cells 2.33 – 2.33 – Macrophages 0.70 – 0.00 – Natural killer cells 0.00 – 0.24 – mRNA expression levels of ACE2 (A) and TMPRSS2 (B) in multiple organs; data obtained from the Genotype-Tissue Expression (GTEx) database (a human organ database); * p  < 0.05, ** p  < 0.001, *** p  = 0.0001, **** p  < 0.0001. The vertical coordinates represent the log transformed (log2 [x + 0.001]) relative expression level of ACE2 and TMPRSS2 . Gene expression profiles of distinct ovarian cell subpopulations based on ovarian scRNA-Seq data from dataset GSE130664 (macaque ovarian single cells). (A) T-distributed stochastic neighbor embedding (t-SNE) clustering of adult ovarian cells into eight clusters based on the scaled data. (B) Average expression level of ACE2 (upper) and TMPRSS2 (lower) in each cell cluster. The vertical coordinates represent the log-transformed [ log (x + 1)] relative expression level of ACE2 and TMPRSS2 . (C) Scatter plot of ACE2 and TMPRSS2 gene expression of each cell cluster. The ellipse indicates the cell cluster of oocytes that co-expressed ACE2 and TMPRSS2 . The number and gradient color of the legend bars show the log-transformed values of gene expression matrix raw data of each cell cluster. (D) Dot plot of the well-known markers including ACE2 and TMPRSS2 in each cell cluster. The color of the legend bar of average expression reflects the mean expression. The size of the dots corresponds to the percentage of cells expressing the listed genes in each cell cluster. (E) Scatter plots of cell clusters with canonical markers of oocytes and (F) granulosa cells. The distribution of ACE2 and TMPRSS2 in each type of ovarian cells in two scRNA-seq. Next, we downloaded transcriptomic sequencing data for human ovary follicles ( GSE107746 ) and divided 151 single cells into two sub-clusters according to the canonical markers of oocytes and granulosa cells ( Fig. 3 A). We found ACE2 expression in most oocytes (86.42%) and expression in some granulosa cells (67.14%) ( Fig. 3 B-D, Table 2 ). Moreover, ACE2 was highly expressed in some oocytes, which is similar to the expression of oocyte markers such as DDX4 and DAZL ( Fig. 3 E) Low expression of ACE2 was observed in human ovarian granulosa cells with the marker of STAR and AMH ( Fig. 3 ).However, 16.05% of all oocytes and 10.00% of granulosa cells were estimated to express TMPRSS2 , showing a significantly lower expression rate than ACE2 ( Fig. 3 D). Fig. 3 Gene expression profiles of oocytes and granulosa cells based on the ovarian scRNA-Seq data from dataset GSE107746 (human oocytes and granulosa cells). (A) T-distributed stochastic neighbor embedding (t-SNE) clustering of adult ovarian cells into two clusters—oocytes and granulosa cells. (B) Average expression level of ACE2 (upper) and scatter plot of ACE2 expression (lower) in each cell cluster based on the canonical cell markers of oocytes and granulosa cells. GCs: granulosa cells. The number and gradient color of the legend bars showed the log-transformed values of gene expression matrix raw data of each cell cluster. (C) Average expression level of TMPRSS2 (upper) and the scatter plot of TMPRSS2 expression (lower) of each cell cluster. (D) Dot plot of the well-known markers including ACE2 and TMPRSS2 in each cell cluster. The color of the legend bar of average expression reflects the mean expression. The size of the dots corresponds to the percentage of cells expressing the listed genes in each cell cluster. (E) Scatter plots of cell clusters with the markers of oocytes and (F) granulosa cells. Fig. 3 Gene expression profiles of oocytes and granulosa cells based on the ovarian scRNA-Seq data from dataset GSE107746 (human oocytes and granulosa cells). (A) T-distributed stochastic neighbor embedding (t-SNE) clustering of adult ovarian cells into two clusters—oocytes and granulosa cells. (B) Average expression level of ACE2 (upper) and scatter plot of ACE2 expression (lower) in each cell cluster based on the canonical cell markers of oocytes and granulosa cells. GCs: granulosa cells. The number and gradient color of the legend bars showed the log-transformed values of gene expression matrix raw data of each cell cluster. (C) Average expression level of TMPRSS2 (upper) and the scatter plot of TMPRSS2 expression (lower) of each cell cluster. (D) Dot plot of the well-known markers including ACE2 and TMPRSS2 in each cell cluster. The color of the legend bar of average expression reflects the mean expression. The size of the dots corresponds to the percentage of cells expressing the listed genes in each cell cluster. (E) Scatter plots of cell clusters with the markers of oocytes and (F) granulosa cells. In addition, we analyzed the expression profile of two other potential SARS-CoV-2 infection-related genes, BSG and CTSL (Fig. S1). In GSE130664 , BSG was mainly expressed in granulosa cells; However, CTSL was undetectable (Fig. S1—B). In GSE107764 , both BSG and CTSL were mainly expressed in parts of the granulosa cells (Fig. S1—D). Together, our analysis of the public datasets indicates that SARS-CoV-2 receptor proteins are expressed in M. fascicularis and human ovaries, especially in oocytes and partial granulosa cells. To further explore the difference in ACE2 and TMPRSS2 gene expression between ovaries of different ages, we compared ovarian single cells from old (18–20 yr) and young (4–5 yr) M. fascicularis (corresponding to human females aged ~70 yr and young ovaries ~20 yr, respectively) from dataset GSE130664 ; human ovarian cumulus cells of old (>40 yr) and young (<35 yr) females from dataset GSE81579 ; and selected data on human ovarian tissues of different ages (20–79 yr) retrieved from GTEx. Dividing each cell cluster into old and young groups showed that the percentage expression levels ( Fig. 4 A) and average expression levels of ACE2 and TMPRSS2 did not significantly differ between old and young ovarian cells ( Fig. 4 C). Fig. 4 B presents the visual outliers on a scatter plot with the significant DEGs displayed in an ellipse; however, ACE2 and TMPRSS2 expression were not included. Although ACE2 was identified as a DEG in the oocyte cluster, the average logFC value was 0.321, and the adjusted P -value was 1, indicating non-significant differential expression. The DEGs of other clusters between old and young groups did not include ACE2 or TMPRSS2 . The specific values of differential expression analysis are listed in Table 3 . Fig. 4 ACE2 or TMPRSS2 gene expression did not significantly differ between old and young ovaries. (A) Dot plot of each cell cluster in GSE130664 (macaque ovarian single cells) in old and young groups, showing the expression level and percentage of cells in a cluster expressing conserved markers (including ACE2 and TMPRSS2 ). (B) Comparative analysis of oocytes and granulosa cells in GSE130664 (macaque ovarian single cells) to visualize ACE2 and TMPRSS2 expression with the outlier scatter plot. The ellipse indicates the differentially expressed genes between old and young ovarian single cells. (C) Violin plots show the expression level and distribution of ACE2 and TMPRSS2 in all ovarian single cells from all seven cell types of each ovarian cell cluster in old and young groups in GSE130664 (macaque ovarian single cells). The vertical coordinates represent the log-transformed [ log (x + 1)] relative expression level of ACE2 . The average expression level of ACE2 and TMPRSS2 showed no significant differences between old and young ovarian single cells. (D) Volcano plot of aged and young human cumulus cells in dataset GSE81579 (bulk RNA-Seq of human cumulus cells), with ACE2 (black circle) included among the genes that showed non-significant differential expression (grey circles); however, TMPRSS2 was not detected. (E) Volcano plot of bulk RNA-Seq of aged and young ovarian tissues in Genotype-Tissue Expression (GTEx) project (bulk RNA-Seq of human ovarian tissues); neither ACE2 nor TMPRSS2 (black circles) was differentially expressed. (F) Volcano plot of oocytes with low and high ovarian reserve in dataset GSE87201 (bulk RNA-Seq of human oocytes), neither ACE2 nor TMPRSS2 (black circles) was differentially expressed. Fig. 4 Table 3 The expression status of ACE2 and TMPRSS2 between old and young oocytes/granulosa cells. Table 3 Gene Cell type Average expression value Percentage of expression (%) Scaled average expression TMPRSS2 Old granulosa cells 0.029124 1.536984 4 TMPRSS2 Young granulosa cells 0.014594 1.526718 3 TMPRSS2 Old oocytes 0.250906 64.78405 14 TMPRSS2 Young oocytes 0.473481 81.11888 20 ACE2 Old granulosa cells 0.025505 1.056676 3 ACE2 Young granulosa cells 0.012305 1.272265 3 ACE2 Old oocytes 3.003148 97.00997 16 ACE2 Young oocytes 4.520156 95.8042 20 ACE2 or TMPRSS2 gene expression did not significantly differ between old and young ovaries. (A) Dot plot of each cell cluster in GSE130664 (macaque ovarian single cells) in old and young groups, showing the expression level and percentage of cells in a cluster expressing conserved markers (including ACE2 and TMPRSS2 ). (B) Comparative analysis of oocytes and granulosa cells in GSE130664 (macaque ovarian single cells) to visualize ACE2 and TMPRSS2 expression with the outlier scatter plot. The ellipse indicates the differentially expressed genes between old and young ovarian single cells. (C) Violin plots show the expression level and distribution of ACE2 and TMPRSS2 in all ovarian single cells from all seven cell types of each ovarian cell cluster in old and young groups in GSE130664 (macaque ovarian single cells). The vertical coordinates represent the log-transformed [ log (x + 1)] relative expression level of ACE2 . The average expression level of ACE2 and TMPRSS2 showed no significant differences between old and young ovarian single cells. (D) Volcano plot of aged and young human cumulus cells in dataset GSE81579 (bulk RNA-Seq of human cumulus cells), with ACE2 (black circle) included among the genes that showed non-significant differential expression (grey circles); however, TMPRSS2 was not detected. (E) Volcano plot of bulk RNA-Seq of aged and young ovarian tissues in Genotype-Tissue Expression (GTEx) project (bulk RNA-Seq of human ovarian tissues); neither ACE2 nor TMPRSS2 (black circles) was differentially expressed. (F) Volcano plot of oocytes with low and high ovarian reserve in dataset GSE87201 (bulk RNA-Seq of human oocytes), neither ACE2 nor TMPRSS2 (black circles) was differentially expressed. The expression status of ACE2 and TMPRSS2 between old and young oocytes/granulosa cells. Moreover, we explored bulk RNA-Seq data of old and young human ovaries retrieved from the GEO and GTEx databases. The upregulated and downregulated gene expression patterns in dataset GSE81579 demonstrated that ACE2 was not included in the DEGs of cumulus cells between young and old females ( Fig. 4 D). The data of human ovarian tissue from GTEx was also analyzed; both ACE2 and TMPRSS2 were not differentially expressed during ovarian aging ( Fig. 4 E). In addition, the analysis of dataset GSE87201 revealed that both ACE2 and TMPRSS2 were not differentially expressed in oocytes between low and high ovarian reserve ( Fig. 4 F). From the results mentioned above, we concluded that there was no disparity in ACE2 or TMPRSS2 expression between ovaries of different ages or with different ovarian reserve conditions in all available datasets. The specific description of each dataset is listed in Table 1 . Next, we validated the results of transcriptome analysis by using immunohistochemistry in human ovaries. Immunoreactivity for ACE2 was found to be distributed in multiple regions of human ovary. As shown in Fig. 5 A, ACE2 was highly expressed in the ovarian cortex from sexual maturity to menopause, especially in oocytes and stromal cells. In addition, immunohistochemistry of human ovaries from 14 to 60 yr-old individuals revealed that TMPRSS2 was also expressed in the ovarian cortex of all ages, especially in oocytes ( Fig. 5 B). We found no difference in the expression of TMPRSS2 in ovaries of different ages, which was consistent with the results of RNA sequencing analysis. To confirm the specificity of the antibodies against ACE2 and TMPRSS2, we performed negative controls with PBS, and used testis sections in positive controls since were reported to overexpress ACE2. Fig. 5 Representative images of ACE2 (A) and TMPRSS2 (B) immunohistochemical staining in the human ovarian cortex of different ages (14 to 60 yr). The positive control was human testis. The negative control was carried out by omitting the primary antibody. Fig. 5 Representative images of ACE2 (A) and TMPRSS2 (B) immunohistochemical staining in the human ovarian cortex of different ages (14 to 60 yr). The positive control was human testis. The negative control was carried out by omitting the primary antibody. Next, we explored the expression of ACE2 and TMPRSS2 in human ovarian medulla of different ages. The results showed that ACE2 was highly expressed in the ovarian medulla from sexual maturity to menopause but not during puberty (Fig. S2). Approximately 80% of the ovarian medullary cells were positive for ACE2. Blood vessels and matrix cells showed intense immunoreactivity for ACE2 (Fig. S2A). The immunoreactivity data suggested that TMPRSS2 was highly expressed in human ovarian medullary blood vessels and stromal cells (Fig. S2B). In conclusion, ACE2 and TMPRSS2 are expressed in human ovarian medulla. As ovarian follicles play an essential role in ovarian function, we further explored the expression of ACE2 and TMPRSS2 in ovarian follicles of different stages. As shown in Fig. 6 , ACE2 (Fig. 6A) and TMPRSS2 (Fig. 6B) immunoreactivity was observed in all phases of follicle development. In primordial follicles, intense ACE2 and TMPRSS2 immunostaining appeared in the oocytes. Strong immunoreactivity for ACE2 and TMPRSS2 was also observed both in oocytes and the granulosa layer of primary follicles. In secondary follicles, the granulosa cells were negative for ACE2, whereas the theca-interstitial cells showed moderate staining for ACE2, and TMPRSS2 was highly expressed in some granulosa cells. In antral follicles, ACE2 was highly expressed in the outer granulosa cells, but the inner granulosa cells were moderate staining. Expression of TMPRSS2 was clearly localized in some granulosa cells and theca-interstitial cells. The cytoplasm of small and large luteal cells in corpora lutea showed intense staining for ACE2 and TMPRSS2. Taken together, ACE2 and TMPRSS2 were expressed in all stages of the follicles, which was consistent with the results of the RNA-Seq analysis. Fig. 6 Immunohistochemical staining indicated that the expression of ACE2 (A) and TMPRSS2 (B) occurred in all types of human follicles. The red arrows indicate follicles. The negative control was carried out by omitting the primary antibody. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 6 Immunohistochemical staining indicated that the expression of ACE2 (A) and TMPRSS2 (B) occurred in all types of human follicles. The red arrows indicate follicles. The negative control was carried out by omitting the primary antibody. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Furthermore, a double immunofluorescence study was performed to explore the relationship between ACE2 and TMPRSS2 expression. The results suggested co-expression of ACE2 and TMPRSS2 in the ovarian cortex, medulla, and follicles ( Fig. 7 , Fig. S3). Together, these results suggest the co-expression of ACE2 and TMPRSS2 in human ovaries. Fig. 7 Immunofluorescence staining showed the co-expression of ACE2 and TMPRSS2 in human ovaries. The yellow arrows indicate follicles. The blue channel is DAPI staining, the green channel is ACE2 expression, and the red channel is TMPRSS2 expression. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 7 Immunofluorescence staining showed the co-expression of ACE2 and TMPRSS2 in human ovaries. The yellow arrows indicate follicles. The blue channel is DAPI staining, the green channel is ACE2 expression, and the red channel is TMPRSS2 expression. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) We compared the features of ACE2/TMPRSS2-positive and ACE2/TMPRSS2-negative ovarian cells in dataset GSE130064 (macaque ovarian single cells) through a DEG functional analysis. The GO functional enrichment analysis indicated that eight virus-related biological process terms were enriched in ACE2-positive oocytes/granulosa cells, including viral transcription, viral gene expression, viral entry into host cells and viral life cycle ( Fig. 8 A, Table S2). The specific genes enriched in the GO terms were visualized in a heatmap-like functional plot ( Fig. 8 C), mainly focusing on the L ribosomal protein family. KEGG pathway analysis showed that a series of viral infection pathways (human papillomavirus infection, human T-cell leukemia virus 1 infection, human cytomegalovirus infection, and viral carcinogenesis) were enriched in ACE2-positive ovarian cells ( Fig. 8 B). Moreover, the focal adhesion pathway may be engaged in regulating virus entry and host response. Furthermore, the functional enrichment analysis of TMPRSS2-positive and TMPRSS2-negative ovarian cells showed that the enriched pathways were mainly involved in the female reproductive process, including oogenesis and oocyte development ( Fig. 8 D, Table S3). The significant KEGG pathways were enriched in oocyte meiosis, homologous recombination, and mismatch repair pathway, indicating the underlying affected meiotic prophase biological events (Fig. 8E). On the basis of the above analysis, we speculate that SARS-CoV-2 may potentially target ACE-positive ovarian cells and therefore influence oogenesis; however, this speculation still lacks direct evidence. Fig. 8 Differentially expressed gene (DEG) functional analysis of ACE2 and TMPRSS2 in macaque single ovarian cells ( GSE130664 ). (A) Gene ontology (GO) enrichment results (only biological process terms shown) of DEGs in ACE2-positive and ACE2-negative ovarian cells. (B) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of DEGs in ACE2 -positive and ACE2 -negative ovarian cells. (C) Visualization of the relationships of GO terms and corresponding mapped genes. (D) GO enrichment result (only BP terms shown) of DEGs in TMPRSS2 -positive and TMPRSS2 -negative ovarian cells. (E) KEGG pathway enrichment analysis of DEGs in TMPRSS2 -positive and TMPRSS2 -negative ovarian cells. BP: biological process. Fig. 8 Differentially expressed gene (DEG) functional analysis of ACE2 and TMPRSS2 in macaque single ovarian cells ( GSE130664 ). (A) Gene ontology (GO) enrichment results (only biological process terms shown) of DEGs in ACE2-positive and ACE2-negative ovarian cells. (B) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of DEGs in ACE2 -positive and ACE2 -negative ovarian cells. (C) Visualization of the relationships of GO terms and corresponding mapped genes. (D) GO enrichment result (only BP terms shown) of DEGs in TMPRSS2 -positive and TMPRSS2 -negative ovarian cells. (E) KEGG pathway enrichment analysis of DEGs in TMPRSS2 -positive and TMPRSS2 -negative ovarian cells. BP: biological process.

M.W., L.M., J.Z. and S.W. conceived and coordinated the study; M.W. and L.M. conducted all the major experiments and wrote the manuscript; L.X., S.Z., J.D. and W.Y. performed and analyzed the experiments; Q.Z. collected human ovarian tissues.

The experiments conducted in the study were ratified by the Human Research Ethics Committee of Tongji Hospital, Tongji Medical College, Huazhong University of Science and Technology. All patients provided written informed consent. Ovarian tissues were obtained from six women at Tongji Hospital for the following indications: adenomyosis, cervical adenocarcinoma cancer, endometrial cancer, and benign ovarian mass. No pathological abnormality was observed in the ovarian samples by the pathologist. The bulk RNA-Seq data from the Genotype-Tissue Expression (GTEx) database from 7845 samples of diverse normal human tissues were downloaded from UCSC Xena ( https://xenabrowser.net/datapages/ ). As GTEx provides the complete data regarding sample collection, RNA-Seq experiments, quality control, and gene expression normalization, such information is not reproduced here [ 17 , 18 ]. The downloaded data were analyzed using “ggstatsplot” package in R to visualize the expression level of ACE2 in multiple tissues. The Kruskal-Wallis test was used to analyze the difference of expression level between the lung and other organs. To identify the differentially expressed genes (DEGs) from old and young ovaries, we retrieved 178 ovarian samples (Table S1) collected from of different age stages (20–29 yr, 22 samples; 30–39 yr, 14 samples; 40–49 yr, 36 samples; 50–59 yr, 61 samples; and 60–79 yr, 45 samples) from GTEx and divided the samples into old (>40 yr) and young (<40 yr) groups; further analysis was performed using the DESeq2 package. Moreover, the DEGs in dataset GSE81579 were also analyzed using the “limma” package. DEGs were identified using the criteria |log 2 (fold change) | >1 and adjusted P _value <0.05. The results of the DEG analysis were displayed as volcano plots, with visualization of the genes ACE2 and TMPRSS2 . Two single-cell RNA sequencing (scRNA-Seq) datasets GSE130664 and GSE107746 were downloaded from the Gene Expression Omnibus (GEO) database ( https://www.ncbi.nlm.nih.gov/geo/ ) with accession numbers GSE130664 (species: Macaca fascicularis ) and GSE107746 (species: Homo sapiens ) ( Table 1 ). The scRNA-Seq dataset analysis followed the process described below. The unique molecular identifier (UMI) counts per gene per cell were processed using the “Seurat” package in R (version 3.1.4) [ 19 ]. After creating the Seurat object, further quality control was conducted. Cells were discarded if  5% of reads were mapped to the mitochondrial genome. Genes expressed in fewer than three cells were also removed. Then, the gene expression data were normalized by log transformation using the “NormalizeData” function (selection.method = “LogNormalize”). After extracting 2000 highly variable genes (HVGs) with the “FindVariableGenes” function, the data matrix was centered and scaled using the “ScaleData” function. For cell clustering and linear dimensional reduction, we performed principle component analysis (PCA) on the HVGs using the “RunPCA” function. To reduce the technical noise, the “JackStraw” procedure was used to determine the PCA dimensionality of the analysis through permutation calculation. Twenty principal components were selected using the “ScoreJackStraw” function for further analysis. The next step was performed using the “FindNeighbors” and “FindClusters” functions; cells were clustered based on a graph-based clustering approach with the resolution as 0.5. A non-linear dimensionality reduction technique was implemented by running the “RunTSNE” and “TSNEPlot” functions. T-distributed stochastic neighbor embedding (t-SNE) plots were thus produced by projecting these datasets onto a two-dimensional space. To find cluster biomarkers, the “FindMarkers” function was used, and the curated known cell markers according to the available publications [ 20 , 21 ] and CellMarker database [ 22 ] were utilized to annotate the cell clusters. Then, the ACE2/TMPRSS2 expression of distinct cell types was evaluated and visualized. To further analyze the DEGs between ACE2/TMPRSS2 -positive and ACE2/TMPRSS2 -negative cells, the function “FindMarkers” was used to generate a list of DEGs. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment pathway analysis were conducted using clusterProfiler [ 23 ]. Table 1 The sources of RNA-seq data from ovarian tissues used in this study. Table 1 GEO accession number Species Tissue or specific cells Sample Number Data Type Cell number Reference GSE130664 (scRNA-seq) Macaca fascicularis Ovarian single cells 8 Raw counts matrix 2601 [ 20 ] GSE107746 (scRNA-seq) Homo sapiens Oocytes and granulosa cells 7 Normalized expression matrix 151 [ 21 ] GSE81579 (bulk RNA-seq) Homo sapiens Cumulus cells 20 Raw counts matrix ⁎ [42] GSE87201 (bulk RNA-seq) Homo sapiens Oocytes 35 Normalized expression matrix ⁎ [43] GTEx portal (bulk RNA-seq) Homo sapines Ovarian tissue 178 Normalized expression matrix ⁎ [ 18 ] ⁎ Standard RNA-seq. The sources of RNA-seq data from ovarian tissues used in this study. Standard RNA-seq. The expression pattern of ACE2 and TMPRSS2 in human ovarian tissues was analyzed by immunohistochemistry. Ovarian tissues from human were formalin-fixed and paraffin-embedded. The 5-μm-thickness serial sections were processed by routine procedures. The sections were deparaffinized and rehydrated through gradient series of ethanol (100%, 95%, 85%, 75%). Then the slides were incubated in 3% H 2 O 2 for inactivation of endogenous peroxidase and were microwaved for 15 min in citrate buffer (pH = 6.0) for antigen retrieval. After cooling at room temperature, slides were blocked with 5% BSA without washing in phosphate-buffered saline (PBS) buffer and incubated with Rabbit-anti-ACE2 antibody (1:200, Proteintech, China) / TMPRSS2 (1:100, ABclonal, China) at 4 °C overnight. Subsequently, the slides were washed three times and incubated with horseradish peroxidase streptavidin solution (ZB-2404; ZSGB-Bio, Beijing, China). The positive signals were colored using DAB chromogenic agent (AR1022, Boster, Wuhan, China) and counterstained using hematoxylin. Negative controls were carried out by omitting the primary antibody. Ovarian tissues from humans aged 25–30 yr were formalin-fixed and paraffin-embedded. The 5 μm thickness serial sections were processed by routine procedures. The sections were deparaffinized and rehydrated through gradient series of ethanol (100, 95, 85, 75%). Then the slides were microwaved for 15 min in citrate buffer (pH = 6.0) for antigen retrieval. After cooling at room temperature, slides were blocked with 5% BSA without washing in phosphate-buffered saline (PBS) buffer and incubated with Rabbit-anti-ACE2 antibody (1:50, Proteintech, China) / Mouse-anti-TMPRSS2 (1:100, Santa Cruz Biotechnology, USA) at 4 °C overnight. Afterwards, the slides were washed three times and incubated with Fitc-conjugated Donkey Anti-Rabbit / Cy3-conjugated Donkey Anti-mouse (Servicebio, China). Then, sections were washed three times in PBS and mounted with Hoechst (Servicebio, China). Images were captured using an Axio Imager upright microscope with fluorescence microscope (Olympus, Germany). Negative controls were carried out by omitting the primary antibody.

The ongoing 2019–2020 COVID-19 pandemic has become a severe public health problem worldwide. As the COVID-19 pandemic continues, the health of more and more young women is threatened. The maintenance of the ovarian reproductive and endocrine function in females is related to the health of the individual and the next generation. Therefore, the crucial question that needs to be addressed is whether the SARS-CoV-2 infection would cause ovarian dysfunction. Recently, the pathological results from the systematic autopsy of COVID-19 victims in China demonstrated the presence of RNA and viral particles of SARS-CoV-2 in the ovary through PCR, transmission electron microscopy, and immunohistochemical staining [ 12 ]. However, it remains unclear how SARS-CoV-2 infects the ovary. In the present study, on the basis of scRNA-Seq data of human and non-human primate ovary and immunohistochemical experiments, we identified the co-expression of the two key entry molecules for SARS-CoV-2, ACE2 and TMPRSS2, in oocytes, luteal cells, blood vessels, and some matrix cells in human ovaries of different ages. Several types of viruses are reported to directly infect the ovary and cause decreased ovarian function, such as endogenous retroviruses (found in the human ovary) and human immunodeficiency virus (HIV; in human germ cells) [ 25 , 26 ]. In addition, hepatitis B virus, herpes simplex virus, and adenoviruses can also be detected in ovarian tissue, including oocytes, granulosa cells and other cells [ [27] , [28] , [29] ]. The mumps virus and HIV are the most recognized and extensively studied viruses that can potentially cause ovarian damage in humans. A multicenter study by Cejtin et al. [ 26 ] including 1139 HIV-seropositive and 292 HIV-seronegative women indicated that HIV-positive women were three times more likely to have protracted amenorrhea than HIV-negative women (10.5% vs. 5.5%). Moreover, the fertility of HIV-seropositive women who do not receive antiretroviral treatment may be lower than that of the seronegative control group [ 30 , 31 ]. Nonetheless, it remains unclear whether SARS-CoV-2 infection results in declined ovarian function. Another point of concern is whether the fertility of SARS-CoV-2-infected women will be adversely affected. While Rotshenker-Olshinka et al. reported that the COVID-19 infection did not seem to affect early first-trimester miscarriage rates in asymptomatic patients [ 32 ][], but there is still a lack of rigorous randomized controlled clinical studies. Therefore, clinicians should pay attention to ovarian function in female patients with SARS-CoV-2 infection. The major changes in the ovary due to viral infections include germ cell apoptosis, interstitial inflammation, leukocyte infiltration, and fibrosis [ 33 ]. Pathological studies have identified that SARS-CoV-2 infection lead to interstitial edema; monocyte, lymphocyte, and neutrophil infiltration; tunica intima inflammation, and cell degeneration in multiple organs [ 12 ]. In addition, some studies suggest that the severity and prognosis of COVID-19 is related to an excessive production of the pro-inflammatory cytokine cascades, also known as a cytokine storm [ 34 ]. Patients with COVID-19 with varying disease severity (mild to severe) had different levels of various cytokines and chemokines. Preliminary research indicates that plasma levels of interleukin (IL)-1β, IL-1RA, IL-7, IL-8, IL-10, IFN-ɣ, monocyte chemoattractant peptide-1, macrophage inflammatory protein (MIP)-1A, MIP-1B, granulocyte colony-stimulating factor, and tumor necrosis factor-alpha (TNF-α) were elevated in patients with COVID-19 [ 35 ]. Moreover, the plasma levels of IL-2, IL-6, IL-8, IL-10, and TNF-α were significantly higher in patients with severe infection than in those with non-severe infections [ 36 ]. A significant number of domestic and foreign studies have confirmed that inflammation mediates changes in the follicular function, which weakens the ability of oocytes to complete meiosis and undergo fertilization and embryo development up to the blastocyst stage and beyond [ [37] , [38] , [39] ]. Oocytes and ovarian cells are also susceptible to oxidative damage, which may deplete the ovarian reserve of primordial follicles and damage the surviving oocytes sufficiently to render them incapable of fertilization or normal development after fertilization; it follows that inflammatory factors may be a vital cause of POI [ 40 ]. Therefore, SARS-CoV-2 infection may affect ovarian function by directly binding to the ACE2 or TMPRSS2 receptor or indirectly through the cytokine storm. In conclusion, we demonstrated the co-expression of the SARS-CoV-2 receptors ACE2 and TMPRSS2 in human ovaries. Therefore, when evaluating the prognosis of COVID-19, ovarian damage and offspring health should be given significant attention. However, further clinical and fundamental research is still needed to monitor the long-term effects of the SARS-CoV-2 on the reproductive system in recovered female patients. The following are the supplementary data related to this article.Supplementary Fig. 1 The expression profile of BSG and CTSL in various datasets. (A) Feature plot of BSG expression in all ovarian cells in GSE130664 (macaque ovarian single cells). The legend bars represent a log-transformed gene expression value. (B) Dot plot of BSG , TMPRSS2 , and ACE2 expression for each major cell population in GSE130664 . (C) Feature plot of BSG and CTSL expression in all cells in GSE107746 (human oocytes and granulosa cells). (D) Dot plot of BSG, CTSL , TMPRSS2 , and ACE2 for the two subpopulations in GSE100746 . Supplementary Fig. 1 Supplementary Fig. 2 Representative images of ACE2 (A) and TMPRSS2 (B) immunohistochemical staining in the medulla of human ovaries of different ages. Supplementary Fig. 2 Supplementary Fig. 3 Quantification of immunofluorescence intensity of expression for both ACE2 and TMPRSS2. Supplementary Fig. 3 The expression of ACE2 and TMPRSS2 in human ovary medulla of different ages. Image 1 Immunofluorescence intensity for ACE2 and TMPRSS2. Image 2 The expression profile of BSG and CTSL in various datasets. (A) Feature plot of BSG expression in all ovarian cells in GSE130664 (macaque ovarian single cells). The legend bars represent a log-transformed gene expression value. (B) Dot plot of BSG , TMPRSS2 , and ACE2 expression for each major cell population in GSE130664 . (C) Feature plot of BSG and CTSL expression in all cells in GSE107746 (human oocytes and granulosa cells). (D) Dot plot of BSG, CTSL , TMPRSS2 , and ACE2 for the two subpopulations in GSE100746 . Representative images of ACE2 (A) and TMPRSS2 (B) immunohistochemical staining in the medulla of human ovaries of different ages. Quantification of immunofluorescence intensity of expression for both ACE2 and TMPRSS2. The expression of ACE2 and TMPRSS2 in human ovary medulla of different ages. Immunofluorescence intensity for ACE2 and TMPRSS2.

The coronavirus disease 2019 (COVID-19) caused by SARS-CoV-2 (previously known as 2019-nCoV) imposes a major threat to the world's public healthcare systems. The critically ill patients can rapidly progress to clinical complications, including acute respiratory distress syndrome (ARDS), septic shock, coagulation dysfunction and even multiple organ failure. Elderly patients and those with underlying health conditions are prone to have a relatively high mortality. Person-to-person transmission of COVID-19 has been repeatedly reported in both in-hospital and household clustering [ 1 ]. According to Zhou et al. [ 2 ], SARS-CoV-2 is a novel coronavirus belonging to the subgenus Sarbecovirus within the genus Betacoronavirus (lineage B) and shares 79.5% of the genetic sequence with SARS-CoV. The SARS-CoV-2 genome encodes four structural proteins; the most important of these—the spike (S) protein utilizes angiotensin-converting enzyme 2 (ACE2) as a cell receptor to facilitate the virus entry into the host cell and cell membrane fusion [ 3 , 4 ]. Recent studies suggest that the cellular entry of SARS-CoV-2 also depends on type II transmembrane serine protease (TMPRSS2) [ 5 , 6 ]. TMPRSS2 primes the viral spike protein, allowing for the potent binding of ACE2, and induces virus-cell fusion. The above findings indicate that the exploration of ACE2 and TMPRSS2 expression in human tissues can predict the potentially infected cells and their respective effects in COVID-19 patients. A series of bioinformatic analyses and experimental verifications have revealed that ACE2 and TMPRSS2 are expressed in the gastrointestinal tract, heart, lungs, and kidneys, demonstrating tissue-specific activity patterns [ [7] , [8] , [9] , [10] ]. In addition to causing lung damage such as ARDS, SARS-CoV-2 has been found to cause kidney, liver, and heart damage and gastrointestinal symptoms [ 11 ]. Most recently, Bian et al. [ 12 ] reported that SARS-CoV-2 RNA and viral particles were detected in multiple organs and tissues including the testis and ovary, which indicated that SARS-Cov-2 can reach the reproductive organs. Pathological analysis of the testes displayed varying degrees of spermatogenic cell reduction and damage. Using public datasets, Wang and colleagues [ 13 ] also provided evidence that the human testis are a potential target site of SARS-CoV-2 infection. The ovary is the core of the female reproductive system, and damage to it can cause infertility and premature ovary failure. One recent study did not detect significant expression of either ACE2 or TMPRSS2 in normal human ovaries by exploring a single-cell sequencing dataset [ 14 ], while another study reported that some oocytes in nonhuman primate ovarian tissue were found to express ACE2 and TMPRSS2 by analyzing transcriptomic and proteomic data [ 15 ]. However, such analyses have limitations; no verified experiments were conducted and they cannot be considered conclusive. Moreover, a detailed analysis of ACE2 and TMPRSS2 in the ovary from puberty to menopause was lacking; such an analysis is critically important to explore whether the ovaries of women of different ages will be affected by SARS-CoV-2. In the present study, we investigated the ACE2 and TMPRSS2 mRNA expression profiles in the human ovary at bulk and single-cell resolution. We also analyzed the co-expression of ACE2 and TMPRSS2 in human ovary sections of women of different ages. Our study revealed that approximately 80% of the ovarian cells were positive for ACE2 and TMPRSS2 ; they are predominantly enriched in oocytes and matrix cells of human ovaries. Additionally, there was no observed variation in ACE2 and TMPRSS2 expression in ovaries of different ages, consistent with the results of RNA sequencing analysis. In addition, we explored the expression of another receptor, basigin (BSG/CD147), and the cathepsins L (CTSL), which have residual cleavage activity toward the viral S protein [ 4 ] [ 16 ]; the results indicated that BSG and CTSL are expressed in somatic ovary cells. Our findings reveal that the receptors for SARS-CoV-2 are expressed in human ovaries and suggests that SARS-CoV-2 may potentially infect ovary.

The authors declare that they have no conflict of interest.