Comprehensive Analysis of the Association Between Human Non-obstructive Azoospermia and Plasticisers via Single-Cell and Traditional RNA Sequencing Methods

A father’s lifetime experience is a major risk factor for a range of diseases in an individual. The influences of exposure can also be transmitted to offspring. Previous studies demonstrated that plasticisers can damage the male offspring reproductive system, but the link between mammalian research and human non-obstructive azoospermia remains underexplored. Here, we analysed reproduction-related genes from four publicly available single-cell RNA-Seq datasets and RNA-Seq datasets on GEO databases to investigate the correlation between human non-obstructive azoospermia and exposure to plasticisers during pregnancy. The R software was used in statistical analyses. A total of 9 co-upregulated genes and 1 co-downregulated gene were found. The Gene Ontology enrichment analyses were performed using the “clusterProfiler” package. Single-cell analyses were conducted to visualise the gene expression pattern in testis cell subgroups. Immunohistochemical images were used to evaluate the level of protein in testicular organs from The Human Protein Atlas. GSVA results provided further functional annotations. Three signature genes, i.e. COL1A1, CYP17A1 and KIF11, might serve as promising diagnostic biomarkers in non-obstructive azoospermia caused by plasticisers. Our results provided a potential new protocol to evaluate the feasibility of life or clinical intervention in patients with non-obstructive azoospermia. We believe that these observations will improve our understanding of the potential mechanisms of plasticiser contributions to human non-obstructive azoospermia and help identify potential targets for immunotherapy in the future.


Introduction
Given the increasing level of industrialisation and urbanisation, plastics, pesticides, plasticisers and synthetic drugs are widely used (Kasonga 2021). Although these materials have improved the certain quality of human life, increasing chemicals are released into the environment and cause potential harm to animals and humans through the food chain. Endocrine-disrupting chemicals (EDCs), which are exogenous substances that cause adverse health effects in undamaged organisms or changes in the endocrine function of the offspring of organisms, include artificial compounds, chemical pollutants, natural oestrogens, plants and fungi (Golshan and Alavi 2019;Li et al. 2021). Many EDCs have stimulating hormone activity, which can affect steroid metabolism and change the oestrogen/androgen balance in the body (Stewart 2020). After the twentieth century, plastic products develop rapidly due to their wide use, low price and Xu Zhang, Xiaohan Ren, Tongtong Zhang and Xiang Zhou have contributed equally to this work.
1 3 durable characteristics. Plasticisers are a kind of polymer material additive that can increase the flexibility of a polymer. Given their versatility, robustness and low production costs, plastics are used in a wide variety of applications. Plasticisers are mixed with polymers to increase the flexibility of plastics. However, plasticisers are not covalently bound to plastics and leach from products into the environment. Di-(2-ethylhexyl)-phthalate (DEHP) is a compound used as a plasticiser that can easily leach from the plastic into the environment. Mono-(2-ethylhexyl) phthalate (MEHP) is the metabolite of DEHP (Kuroda et al. 2020). Trioctyl trimellitate (TOTM), a new plasticiser, may be the preferred plasticiser instead of DEHP for medical polyvinyl chloride instruments due to its good plasticising property, migration resistance, high temperature resistance and low toxicity (Sheikh et al. 2016).
The male semen quality is closely related to male reproductive health and have been paid attention by people (Rahban and Nef 2020). In the past few decades, many studies reported a decline in human sperm quality. Wang et al. demonstrated that many semen parameters of Chinese young men show a decreasing trend within 15 years (Wang et al. 2019). Lidia Mínguez-Alarcón et al. found that sperm parameters, such as sperm concentration, total sperm count, motility and morphology, have significantly decreased amongst subfertile men from 2000 to 2017 (Minguez- Alarcon et al. 2018). The decline in sperm quality is affected by many factors, including tobacco and alcohol use, irregular daily routine and environmental toxicants. Recently, many studies focussed on the effects of plasticiser on male sperm quality due to the wide exposure to plasticiser through the food chain. The distribution of plasticiser or their metabolites can inevitably be detected in various organs of the human body (Sun et al. 2020). Animal studies confirmed that exposure to DEHP phthalate during pregnancy can lead to the deterioration of sperm parameters and imbalance of sex hormone secretion in male rat offspring (Lucas et al. 2012). One in 100 men is estimated to have azoospermia, which is the complete lack of sperm in the ejaculate. Currently, 20% of azoospermia cases remain idiopathic. The non-obstructive azoospermia (NOA) is mostly explained by congenital factors leading to spermatogenic failure, such as chromosome abnormalities and gene mutation 5 . However, knowledge on the link between NOA and exposure to plasticisers during pregnancy remains lacking. Knowledge on the causes of NOA is limited. High genetic heterogeneity due to the complexity of spermatogenesis and testicular function, lack of nonconsanguineous familial cases and confirmatory studies challenge this field. In clinics, patients with NOA often need to be treated with testicular and epididymal incision and sperm extraction under a microscope, but the success rate is relatively low. Therefore, increasing attention has been paid to the study of the pathogenesis of NOA and the further exploration of prevention and treatment measures. Further studies confirmed the associations of exposure to plasticisers during pregnancy with reproductive hormones and semen parameters (Dobrzynska and Tyrkiel 2019). The positive relationship between sperm concentration and the metabolite of DEHP is found (Barakat et al. 2019). Besides, with increased DEHP metabolites found in urine, the concentrations of serum FSH and LH increase, indicating impaired function of the testis and infertility (Hauser 2008). Iman Al-Saleh et al. confirmed that men who excrete high levels of MEHP relative to the other DEHP metabolites may be at risk of having low sperm concentration (Al-Saleh et al. 2019).
In the present study, two gene chips are downloaded from the NCBI-Gene Expression Omnibus (GEO) database (GSE41901 and GSE9210 datasets) to search for differentially expressed genes (DEGs) in the testes of offspring exposed to plasticiser during pregnancy and human testes with NOA. The possible mechanism of plasticiser on the male spermatogenic function is explored. The Gene Ontology (GO) functional annotation analysis is applied. Immunohistochemical (IHC) and single-cell analyses are used to compare the expression of these key genes in testicular cells. We hope to explore the link between NOA and exposure to plasticisers during pregnancy. Corresponding biomarkers are found to guide the diagnosis of male infertility caused by the plasticiser. The decline in human sperm level is an indisputable fact. Besides, male reproductive health should be protected.

Data Set Download
Two transcription profile datasets (i.e. GSE41901 and GSE9210) were obtained from GEO databases (http:// www. ncbi. nlm. nih. gov/geo/), an open functional genome dataset. The platform for GSE41901 was GPL7294 (Agilent-014879 Whole Rat Genome Microarray 4 × 44 K G4131F). The microarray data of GSE9210 were based on the GPL887 platform (Agilent-012097 Human 1A Microarray (V2) G4110B). The GSE41901 dataset included 20 foetal testis samples and rats exposed to DEHP, MEHP and TOTM in utero at doses of 500 mg/kg administered through daily oral gavage to pregnant dams between gestational days 12 and 19. The RNA of their foetal testis was extracted. The GSE9210 dataset was based on genomic gene expression in testes of 47 patients with NOA and 11 patients with obstructive azoospermia (OA). Patients diagnosed with azoospermia underwent two semen examinations, and two centrifugal semen samples confirm the absence of sperm in the ejaculate. All patients underwent testicular sperm extraction (TESE) for diagnosis. OA is defined as motile sperm collected from microsurgical epididymal sperm aspiration, or a large number of mature sperm are found form testicular sperm extraction. NOA is defined as having no sperm from microsurgical epididymal sperm aspiration and testicular sperm extraction. In addition, azoospermia patients with varicocele, ejaculatory dysfunction and endocrinopathy were excluded (Okada et al. 2008). In addition, in order to verify the specificity of the screened genes as diagnostic biomarkers, two gene chips (GSE7755, GSE108194) were downloaded from GEO databases, which is correlated with cryptorchidism and human oligospermia. The GSE7755 dataset included 34 normal testis samples and 20 testis samples with cryptorchidism. The GSE108194 dataset was based on genomic gene expression of 18 patients in peripheral blood with at least two semen routine examinations showing sperm density lower than 15 million/ml and genomic gene expression of 6 normal people in peripheral blood.

Processing of Dataets and Screening of DEGs
We first downloaded the original comment file and TXT matrix file of the platform and used the original comment file of the platform to change the probe ID in the matrix to gene symbol ID through the Perl script. We divided the data of GSE41901 into two groups, i.e. control and exposure groups, to investigate their expression. The GSE9210 dataset was also divided into two groups (i.e. patients with NOA and OA). The R software (version 4.0.2, https:// www.r-proje ct. org/) was used in data mining and statistical analyses, and the LIMMA package was used to mine the mRNA of abnormal expression levels. Through the LIMMA analysis, genes with absolute logFC > 0.8 and P value < 0.05 were considered to be DEGs.

GO and Pathway Enrichment Analyses
The "clusterProfiler" package was used in R to perform the GO functional annotation pathway enrichment analysis and explore the biological and molecular functions of DEGs. The biological process (BP), molecular function (MF) and cellular component (CC) in the GO analysis were performed. P < 0.05 was selected as threshold in the analysis. The Venn diagram was used to identify co-downregulated and coupregulated genes in patients with NOA and rat offspring with plasticiser exposure.

Assessment of DEGs at the Single-Cell Transcriptional Level
Two adult testis tissues of GSE124263 from the GEO dataset were used to visualise the expression patterns of co-downregulated and co-upregulated genes in testis cell subgroups and examine the correlation between co-expression genes and the normal testis tissue at the single-cell level. The selected cell samples contained different cell populations, including Sertoli cell, spermatogonium, Leydig cell, endothelial cell, myoid cell, primary spermatocyte, macrophage, elongated spermatid, differentiating spermatogonia and sperm. The Seurat R package was used to perform the single-cell analysis. The gene features in each cluster were found by "scanpy". The cell-cell interactions were calculated using the "SingleCellSignalR" package in R.

IHC Detection in Testicular Tissue Sections
IHC images were used to evaluate the level of co-downregulated and co-upregulated genes and proteins in testicular organs from The Human Protein Atlas website (https:// www. prote inatl as. org/) and validate the different expression levels of co-expression genes in testis cells.

Screening of DEGs
We used the RNA-seq data from the GEO database to compare the differential expression of two transcription profiles. After data download and preprocessing by using the Perl script, we manually grouped the GSE9210 dataset to search for DEGs between DEHP-, MEHP-and TOTM-exposure and control groups. Compared with those in the control group, 244 upregulated and 513 downregulated genes were found in DEHP-exposure groups, 377 upregulated and 1230 downregulated genes were found in MEHP-exposure groups and 413 upregulated and 1476 downregulated genes were found in TOTM-exposure groups. Similarly, we then compared the DEGs in the testis of patients with NOA and OA (control group), and results showed 1217 downregulated and 483 upregulated genes in the testis of patients with NOA compared with those in the testis of patients with OA. The standard was P < 0.05 and absolute log2FC > 0.8 (Fig. 1).

Amalgamation of DEGs
We observed the expression of mRNA in testis and tried to further identify the genes that were upregulated or downregulated in the testis of male offspring after exposure to DEHP, MEHP and TOTM during pregnancy. A total of 227 co-upregulated and 91 co-downregulated genes were found. We then converted the genes expressed in Rattus norvegicus into human homologous genes to figure out the relationship between exposure to EDCs during pregnancy and azoospermia in offspring. Finally, nine genes, including TKT, NME2, CTSB, RPL10, FTL, COL1A1, CYP17A1, NDN and IGF2, were upregulated in male rat offspring exposed to plasticisers during pregnancy and human patients with NOA. Besides, only one gene, KIF11, was downregulated in male rat offspring after exposure to plasticisers during pregnancy and human offspring with azoospermia (Fig. 5).

DEG Investigation at Molecular and Functional Levels
GO BP, CC and MF pathway analyses were carried out for the functional analysis of DEGs. The analysis of DEGs in male rat offspring exposed to DEHO, MEHP and TOTM during pregnancy showed their function. In BP (Fig. 2), DEGs were enriched in cholesterol metabolic or biosynthetic process, sterol biosynthetic or metabolic process and cytoplasmic translation. In CC (Fig. 3), DEGs were enriched in the cytosolic ribosome, ribosomal subunit, ribosome, cytosolic large ribosomal subunit and cytoplasmic side of the rough endoplasmic reticulum membrane. In MF (Fig. 4), DEGs were enriched in the structural constituent of ribosome and rRNA binding. The results of DEG functional analysis between OA and NOA showed that DEGs were enriched cellular processes involved in reproduction in multicellular organism, germ cell development, spermatid development and differentiation and meiotic cell cycle and fertilisation in BP. In CC, DEGs were enriched in acrosomal vesicle, motile cilium and collagen-containing extracellular matrix. In MF, DEGs were enriched in tubulin binding, nucleocytoplasmic carrier activity and microtubule motor activity (Fig. 5).

Single-Cell Analysis of Co-upregulated and Co-downregulated Genes
The tissues of testis contain several cell populations, such as Sertoli cell, spermatogonium, Leydig cell, endothelial cell, myoid cell, primary spermatocyte, macrophage, elongated spermatid, differentiating spermatogonia and sperm. To investigate the expression levels of co-upregulated and co-downregulated genes in human testicular tissues, we performed a single-cell analysis of the GSE124263 dataset. The UMAP method was used to perform nonlinear dimensional reduction. Cells were clustered in accordance with their cell type. The microsatellite status was also shown by the UMAP plot. For all co-upregulated genes, NME2, TKT, CTSB, FTL and RPL10 showed no significant difference in gene expression amongst different clusters. NBN and COL1A1 showed an enrichment in Sertoli, Leydig and myoid cells, whereas We also investigated the KIF11, a co-downregulated gene in the exposure group and patients with OA, and cell distribution patterns by the single-cell analysis. Interestingly, KIF11 was enriched in spermatogonium and primary spermatocyte (Figs. 6 and 7).

IHC Detection of NBN, COL1A1, CYP17A1 and KIF11 in Testis Tissue Sections
In HPA databases, we evaluated the mRNA and protein expression levels of NBN, COL1A1, CYP17A1 and KIF11 in the testis tissue. CYP17A1 was stably expressed in the testis tissue. Elevated staining for CYP17A1 was observed in the Leydig cell, whereas no staining was found in the seminiferous cell. Low staining for COL1A1 was observed in the Leydig cell, whereas moderate staining was observed in the seminiferous cell. IHC results revealed that NBN was highly expressed in the seminiferous cell and moderately expressed in the Leydig cell, and this finding was not consistent with the results of our single-cell analysis. KIF11, the only co-downregulated gene, was observed to be overexpressed in the seminiferous cell and had low expression in the Leydig cell (Fig. 8).

Biological Relevance of COL1A1, CYP17A1 and KIF11
We performed the GSVA of COL1A1, CYP17A1 and KIF11 to figure out the GO BP, CC and MF pathways in their overexpressed cells. GSVA showed enrichment of COL1A1, CYP17A1 and KIF11 in their overexpressed cells. In spermatogonia, the GO pathway analysis showed that cells with low expression of KIF11 were significantly enriched in cyclin binding, apoptotic process, toxic substance binding, biological adhesion, DNA replication and cell cycle compared with cells with high expression of KIF11. In the primary spermatocyte, the GO pathway analysis showed that cells with low expression of KIF11 were significantly enriched in cyclin binding, cell cycle, toxic substance binding, identical protein binding and nuclear outer membrane of the endoplasmic reticulum membrane network compared with cells with high expression of KIF11. The COL1A1 gene was upregulated in Sertoli and Leydig cells. In Sertoli cells, specific GO categories, such as exon-exon junction complex, Fig. 2 Top 4 terms from the GO BP enrichment analysis. a DEHP-exposure group. b MEHP-exposure group. c TOTM-exposure group. d NOA group. The size of the dot represents counts, and the colour of the dot represents P value chaperone binding, RNA binding, cyclin binding and microtubule cytoskeleton, were significantly related to cells with high expression of COL1A1 than to cells with low expression of COL1A1. In Leydig cells, GO results showed that these cells with high expression of COL1A1 could participate in sequence-specific DNA binding, biological adhesion and protein dimerisation activity. Besides, in Sertoli cells, cells with high expression of CYP17A1 were enriched in kinase binding, misfolded protein binding and microtubule cytoskeleton compared with cells with low expression of CYP17A1 (Fig. 9).

Exploration of the Role of COL1A1, CYP17A1 and KIF11 in Cryptorchidism and Human Oligospermia
In order to verify the specificity of COL1A1, CYP17A1 and KIF11 as the diagnostic biomarkers for non-obstructive azoospermia caused by plasticisers, two transcription profile datasets (GSE7755, GSE108194) were obtained from GEO databases, which is correlated with cryptorchidism and human oligospermia. Then, the differential expressed analysis of GSE108194 dataset and GSE7755 dataset were performed with absolute logFC > 0.5 and P value < 0.05. We obtained a total of 527 co-upregulated and 340 co-downregulated genes in the GSE108194 dataset. In the GSE7755 dataset, a total of 68 co-upregulated and 119 co-downregulated genes were observed. The Venn diagram demonstrated that COL1A1, CYP17A1 and KIF11 are not involved in the DEGs of GSE108194 dataset and GSE7755 dataset (Fig. 10). The results revealed that COL1A1, CYP17A1 and KIF11 are not correlated with other reproductive disorders, including cryptorchidism and human oligospermia.

Discussion
Infertility affects about 1 in 6 couples attempting pregnancy, and the man is responsible in approximately half of the cases. Given that the pathophysiology underlying azoospermia is not elucidated, most male infertility is diagnosed as idiopathic. With the recent azoospermia pandemic, male reproductive health, as one of the important reasons that cannot be ignored, has attracted increasing attention. At present,  the aetiology of most NOA remains unclear. The continuous accumulation of EDCs in the human body undoubtedly causes irreversible harm to the male reproductive system (Cui et al. 2020). Plasticiser, a polymer auxiliary that is widely used in industrial production, can be found in food packaging, cosmetics and medical equipment. Previous studies confirmed that the damaging effect of plasticisers on the male reproductive system is intergenerational. Exposure to DEHP during lactation or pregnancy can damage the development of the male reproductive system of F1 offspring. A recent study evaluated sperm quality associations with urinary phthalate metabolites amongst 660 Chinese adult men. Results show that elevated urinary phthalate metabolites may affect semen quality.
The present study was based on the integrated analysis of single-cell RNA-Seq datasets and RNA-Seq datasets from GEO, which included testis samples from patients with NOA and OA and healthy adults and rats. We aimed to associate NOA with plasticiser exposure during pregnancy. After searching DEGs, we identified hundreds of genes whose expression was significantly altered in patients with NOA and rats with plasticiser exposure during pregnancy. Then, we first filtered nine co-expressed genes, including TKT, NME2, CTSB, RPL10, FTL, COL1A1, CYP17A1, NDN and KIF11. Then, we performed single-cell analysis on nine co-expressed genes and finally found three specific genes.
We found that KIF11 was significantly downregulated in patients with NOA and the plasticiser exposure group. Miki et al. found that the KIF11 protein is expressed in spermatocytes (Hara-Yokoyama et al. 2019). Besides, the immunolocalisation of KIF11 suggests the low expression in Sertoli cells. We analysed hundreds of thousands of cells in testis tissue and examined the overexpression of KIF11 in spermatogonium and primary spermatocyte. Therefore, we speculated that KIF11 played a key role in spermatogenesis. Besides, the downregulation of KIF11 could be one of the key factors that link NOA with plasticisers. A large proportion of KIF11 was correlated with the NOA led by plasticisers, which improved our understanding of NOA occurrence and helped in the development of therapeutic targets.
Testis tissues are mixtures of different compartments, e.g. Sertoli, Leydig, endothelial and myoid cells, which provide a good growth environment for sperm production and maturation. CYP17A1 is a steroidogenesis-related gene (Molaie et al. 2019). Testosterone is synthesised by Leydig cells, and the gene responsible for testosterone synthesis and cholesterol transport in testis is CYP17A1. CYP17A1 is one of the important genes in the microenvironment of sperm production (Yu et al. 2020). Sperm capacitation/acrosome reaction and testosterone biosynthesis occur in Leydig cells. We performed single-cell analysis and IHC detection and found that CYP17A1 was upregulated in patients with Fig. 6 Overview of the single cells from testis tissue samples. a The cell types identified by marker genes. UMAP plot of the marker genes. b-j Expression of co-expression genes for Sertoli cell, sper-matogonium, Leydig cell, endothelial cell, myoid cell, primary spermatocyte, macrophage, elongated spermatid, differentiating spermatogonia and sperm NOA and rats with plasticiser exposure. Our results showed that plasticisers might also damage sperm by affecting the function of Leydig cells to reduce testosterone production. Plasticisers may directly impair spermatogenesis and sperm maturation. Further study confirmed that promoter polymorphism may signify a genetic risk factor for male infertility (Rahali et al. 2020).
COL1A1 is a key gene implicated in epididymal immunoregulation, inflammation and fibrosis (Wijayarathna et al. 2020). The downregulation of COL1A1 in testis tissue often indicates the damage of immune system, inflammation or fibrosis (Morgan et al. 2010). By searching DEGs, the upregulation of COL1A1 was found in patients with NOA and F1 male rat offspring with plasticiser exposure. Our findings also suggested that the destruction of sperm quality in male rat offspring was associated with inflammation, fibrosis or immune damage. The elucidation of the mechanisms of immune damage, inflammation and fibrosis might improve our understanding of NOA and the reproductive damage of plasticisers.
After the enrichment analysis of DEGs, we found that plasticiser was involved in a diversity of functions, such as sterol biosynthetic and metabolic processes, cholesterol biosynthesis, metabolic process and ribosome biogenesis, and participated in a series of physiological and pathological processes, including RNA binding, coenzyme binding and amide binding. Besides, the GO pathway analysis was conducted on the DEGs of NOA, and results showed that the functions and pathway genes activated were predominantly related to spermatid development and differentiation, fertilisation and germ cell development. Figure 11 shows the pathway related to metabolism. A direct comparison of cells with high and low expression levels of CYP17A1, COL1A1 and KIF11 was performed on the basis of the single-cell sequencing. As shown in Fig. 11, apoptotic process, cell cycle, toxic substance binding and identical protein binding were the enriched pathways. Recent studies showed that cell cycle and apoptosis regulation are associated with NOA. The percentage of apoptotic cells differentiated from patients with NOA was higher than those from normal ones. Besides, the DEGs of patients with NOA were involved in cell cycle and regulation of apoptosis. These defects might be related to the genetic causes of male infertility. A previous study confirmed that exposure to plasticisers during pregnancy leads to the activation of peroxisome proliferator-activated receptors, increased fatty acid oxidation and reduction in the ability to cope with the augmented oxidative stress, resulting in reproductive organ malformations, reproductive defects and decreased fertility. In summary, exposure to plasticisers during pregnancy can affect the reproductive system of F1 male rat offspring. Besides, the long-term exposure to plasticisers during pregnancy is likely to be an important factor in male NOA. We explored DEGs and revealed the critical genes that might regulate cholesterol metabolism, transcription and protein translation in testis cells. In addition, single-cell RNA-seq had relatively good performance in distinguishing the molecular characteristics in each cell type. With the help of single-cell analysis, we aimed to figure out the expression of all co-downregulated and co-upregulated DEGs in different testis cells. Combined with IHC results, COL1A1, CYP17A1 and KIF11 were found to have a strong correlation in exploring the link between exposure to plasticisers Fig. 9 GSVA analysis in human testis tissue samples. a GSVA analysis of the GO pathways between COL1A1 lowly expressing Leydig cells and COL1A1 highly expressing Leydig cells. b GSVA analysis of the GO pathways between COL1A1 lowly expressing sertoli cells and COL1A1 highly expressing sertoli cells. c GSVA analysis of the GO pathways between CYP17A1 lowly expressing sertoli cells and CYP17A1 highly expressing sertoli cells. d GSVA analysis of the GO pathways between KIF11 lowly expressing primary spermatocyte and KIF11 highly expressing primary spermatocyte. e GSVA analysis of the GO pathways between KIF11 lowly expressing spermatogonia and KIF11 highly expressing spermatogonia during pregnancy and human NOA. In addition, we performed the validation analysis based on the dataset of cryptorchidism and human oligospermia. The results showed that COL1A1, CYP17A1 and KIF11 are not correlated with other reproductive disorders, such as cryptorchidism and human oligospermia. Furthermore, we speculated that plasticisers might induce the impairment of the structure of seminiferous tubules and spermatogenesis, which might be attributed to the imbalance of testosterone. The damage of immune system, inflammation or fibrosis might also be involved in this impairment. Given the prevalence of plasticisers in our atmosphere and environment and the impairment of male reproductive system, we suggest their replacement. With increasing efforts to develop viable replacement compounds, rigorous leaching, toxicity and impact assessment studies are needed before alternative plasticisers can be adopted as viable replacements.

Data Availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Conflict of interest
We declare that we have no financial and personal relationships with other people or organisations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the research of "Comprehensive analysis of the association between human non-obstructive azoospermia and plasticizers: via single-cell RNA sequencing and traditional RNA sequencing".
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Fig. 11
Top terms from the GO BP, CC and MF enrichment analysis. A DEHP-exposure group. B MEHP-exposure group. C TOTM-exposure group. D NOA group. The size of the dot represents counts, and the colour of the dot represents P value