The present study was approved by the Ethical Committee of the Jiangxi Maternal and Child Health Hospital (Nanchang, China) and all participants provided written informed consent. All experiments were carried out in accordance with the guidelines provided by the Jiangxi Maternal and Child Health Hospital.

Samples of seminal plasma were collected at the Jiangxi Maternal and Child Health Hospital between October and December 2016 from three groups of men without leukocytospermia as follows: The control (C) group (22–36 years old), the obstructive azoospermia (OA) group (24–38 years old) and the non-obstructive azoospermia (NOA) group (25–37 years old). A total of 5 fertile donors with normal sperm parameters (seminal volume, 2–6 ml; sperm count, 0.6–1.5×10¹²/l; progressive motile spermatozoa ≥30% or progressive motile spermatozoa + motile spermatozoa ≥50%) (4) were recruited as the C group. A total of 12 patients with azoospermia who underwent surgical sperm retrieval via percutaneous epididymal sperm aspiration or testicular sperm extraction were recruited and were divided into the OA group (n=6) and the NOA group (n=6). In OA, sperm are produced normally inside the testicle, but are prevented from reaching the ejaculate due to a blockage or obstruction in the reproductive tubing. In non-obstructive azoospermia group, the tubes are not blocked; however, there is a spermatogenesis problem whereby there is a very low level of sperm production or a total lack of production (8). Patients with abnormal karyotypes and those who had previously suffered from an injury to the genitals were excluded.

Prior to sampling, the glans penis of each patient was washed using soap and water. Semen was obtained by masturbation, ejaculated into a sterile collection tube and incubated at 37°C for 30 min for liquefaction. Basic semen parameters and leukocytospermia were detected as previously described and the samples were stored at −80°C prior to DNA extraction (20).

Genomic DNA was extracted from each semen sample using a TIANamp Genomic DNA kit (Tiangen Biotech Co., Ltd., Beijing, China) and combined with bead beating as previously described (14). The concentration and quality of extracted genomic DNA were tested prior to sequencing using a spectrophotometer at 230 nm (A 260) and 260 nm (A 260) (NanoDrop; Thermo Fisher Scientific, Inc., Waltham, MA, USA). The extracted genomic DNA was used as a template to amplify the V4 region of the 16S rRNA genes of all samples using 515F/806R primers (515F, 5′-GTGCCAGCMGCCGCGGTAA-3′; 806R, 5′-GGACTACVSGGGTATCTAAT-3′). The PCR reaction conditions were as follows: 2 min initial denaturation at 94°C and 30 cycles of 30 sec denaturation at 94°C, 30 sec annealing at 52°C, 30 sec extension at 72°C, followed by a 5 min final extension at 72°C. The samples were maintained at 4°C when the PCR reaction finished. Each 25 µl reaction made up of 2.5 µl 10X Ex Taq buffer (Mg²⁺ free; Takara Biotechnology Co., Ltd., Dalian, China), 2 µl (2.5 mM) dNTP mixture, 1.5 µl Mg²⁺ (2.5 mM), 0.25 µl Ex-Taq DNA polymerase (Takara Biotechnology Co., Ltd.), 0.5 µl barcode primer V4F (10 µM), 0.5 µl barcode primer V4R252 (10 µM), 2 µl of fecal DNA and 15.75 µl milli-Q water. These PCR products were sequenced using an Illumina HiSeq 2000 platform (21).

Paired-end reads from the original DNA fragments were merged using FLASH (version 1.2.7; ccb.jhu.edu/software/FLASH/) when any of the reads overlapped that generated from the opposite end of the same DNA fragment. Paired-end reads were assigned to each sample according to the unique barcodes using Trimmomatic (version 0.33) (22). Sequence analysis was subsequently performed using the UPARSE software package (version 7.0.100) (23) with the UPARSE-OTU and Chao1, the Chao1 estimator and the Shannon index (24) with UPARSE-OTUref algorithms. In-house Perl scripts were used to analyze α (within samples) and β (among samples) diversity. Sequences with ≥97% similarity were assigned to the same OTUs. A representative sequence was selected for each OTU and the RDP classifier (version 2.2) (25); was used to annotate taxonomic information for each representative sequence. Cluster analysis was preceded by weighted UniFrac distance using QIIME software package (version 1.8.0) (26), partial least squares discriminate analysis (PLS-DA) was preceded using SIMCA-P software version 11.5 (Umetrics; Sartorius Stedim Biotech, Malmö, Sweden) and the species that were differentially abundant were characterized for their metabolic capacity using Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt), Version 1.0.0 (27,28). The statistical significance was set at P<0.05 for correction of multiple comparisons.

To compare microbes in the OA and NOA groups, the V4 hypervariable region of bacterial genomic DNA was sequenced using high-throughput sequencing and the effective tags with ≥97% similarity were considered as one OTU. In total, 675,539 usable sequences and 5,867 OTUs were obtained from the samples, with a mean of 345 OTUs per group (data not shown). Furthermore, the Chao1 and Shannon indices were nearly saturated and the rarefaction curve of every sample entered the plateau phase, which indicated that the sequencing data was reasonable (OTU number was sufficient for further analysis) and species composition was highly uniform in each group (data not shown).

At the phylum level, data from the 20 microorganisms that were most abundant in each group was analyzed. Firmicutes, Proteobacteria, Bacteroidetes and Actinobacteria were the predominant phyla, accounting for 98.14% in the C group, 98.26% in the OA group and 90.96% in the NOA group (Fig. 1A).

Compared with the C group, azoospermia reduced the OTU number of Cyanobacteria, Acidobacteria, Gemmatimonadetes, Planctomycetes, Chloroflexi, Crenarchaeota, Armatimonadetes, Elusimicrobia, Nitrospirae, Euryarchaeota, Spirochaetes, Chlorobi, Synergistetes, Chlamydiae and Verrucomicrobia, whereas the number of Bacteroidetes was increased (P<0.05; Fig. 1B).

Lactobacillus, Prevotella, Proteus, Pseudomonas, Veillonella, Corynebacterium, Rhodococcus, Staphylococcus and Bacillus were the dominant genera detected in the C (55.02%), OA (62.54%) and NOA (59.74%) groups (Fig. 2A). The distribution map indicated that Alicyclobacillus, Amaricoccus, Anaeromyxobacter, Aquicella, Arsenicicoccus, Azospirillum, Chitinimonas, Chlorobaculum, Coprococcus, Desulfovibrio, Dokdonella, Gallionella, Geobacter, Helicobacter, Idiomarina, Kaistia and Kribbella were significantly reduced in the C group compared with the OA and NOA groups (P<0.05; Fig. 2B).

There were 1,093, 925 and 840 OTUs detected in the C, OA and NOA groups, respectively (Fig. 3A). A total of 398 common OTUs were identified in the C, OA and NOA groups, of which 27 belonged to the Lactobacillus phylum (Fig. 3A). In addition, the partial least squares discriminant analysis (PLS-DA) indicated that patients in the OA and NOA groups scattered on the left, whereas the C group gathered together on the right of the coordinate axis (Fig. 3B). No overlapping of the C group with the OA or NOA groups was observed; however, some samples in OA and NOA groups overlapped (Fig. 3B), indicting the microbial diversity between the OA or NOA groups had a higher degree of similarity than the C group.

Significantly more Sneathia and Lysobacter sp. were present in the NOA group compared with the C and OA groups (P<0.05; Fig. 4). The number of Solibacillus, Campyiobacter, Campyiobacteraceae and Plesiomonas OTUs in the OA group was significantly higher compared with the NOA and C groups (P<0.05; Fig. 4). The number of bacteria in red-labeled genera was significantly higher in the C group compared with the OA and NOA groups (P<0.05; Fig. 4).

Correlation analyses were performed for various Kyoto Encyclopedia of Genes and Genomes pathways using the PICRUSt software. The association between seminal bacteria and health was investigated (Fig. 5A). The results indicated that patients in the OA and NOA groups were at an increased risk of having metabolic, infectious and immune diseases (Fig. 5A). Furthermore, the results indicate that azoospermia enhances the protein expression of genes associated with nucleotide metabolism, metabolism of cofactors and vitamins, glycan biosynthesis and metabolism, enzyme families and energy metabolism; however, the expression of genes associated with xenobiotic biodegradation and metabolism, metabolism of terpenoids and polyketides, lipid metabolism and amino acid metabolism was downregulated (Fig. 5B).