A multilevel analysis identifies the different relationships between amino acids and the competence of oocytes matured individually or in groups

High-protein diets contribute to an increase in urea follicular concentrations associated with decreased fertility. Urea has been shown to interfere with the epidermal growth factor (EGF)/EGFR system, which has been shown to have a beneficial effect during in vitro maturation (IVM) of oocytes. Of note, the number of cumulus-oocyte complexes (COCs) in the maturation medium can change the maturation and the developmental competence of COCs. Therefore, it was hypothesized that, the presence of urea and EGF may have a differential effect on the depletion/appearance of AAs and competence of COCs matured individually (I-IVM system) or in groups (G-IVM system). In the G-IVM system, COCs increased consumption (depletion) of AAs compared with other groups in the presence of high-level urea (40 mg/dl) + EGF (10 ng/ml). In the I-IVM system, the non-cleaved COCs depleted more AAs than the cleaved COCs, in particular in the presence of urea. The combination of urea and EGF increased the depletion of AAs in the G-IVM system. However, the EGF abrogated the urea-induced depletion of AAs by the I-IVM COCs. The use of N-acetyl-l-cysteine as an EGFR inhibitor canceled urea-induced depletion of AAs. This shows the inhibiting effect of urea over the EGF/EGFR system. In the presence of urea + EGF, COCs had a lower degree of developmental competence than control in both I- and G-IVM systems. Arginine had the best predictive power to identify highly competent COCs in the G-IVM system, while glutamine was the best predictor of the cleavage in the I-IVM system. In conclusion, this multi-level study shows that COCs matured individually or in groups may have different association with AAs metabolism. These findings provide new insights into the relationships between AA metabolism and the subsequent developmental competence of COCs.


Intracellular urea content in bovine denuded oocytes (DOs).
The results revealed a concentration of 3.43 ± 1.4 mg/dl (mean ± SD) of urea in DOs cultured for 4 h with 40 mg/dl of urea (Fig. 1a).
Calculation of depletion or appearance of amino acids (AAs). Differences between the concentration of AA in the test medium (after 24 h of incubation) and that of the fresh medium were used to calculate the depletion or appearance of each AA 18,19 . Negative values indicate that AAs have disappeared from the maturation medium after 24-h incubation (depletion). Positive values indicate that AAs appeared at higher concentrations in the maturation medium after 24-h incubation (appearance). For the calculation of "total depletion", "total appearance" of AAs, the negative (depletion) or positive (appearance) values obtained for AAs were summed up 18,19 ( Supplementary Fig. 1, Phase 1). Next, the "total net balance" of all AAs was calculated by subtracting "total depletion" from "total appearance". In addition, the "total turnover" of all AAs was determined by the sum of "total depletion" and "total appearance". Therefore, there were 6 data (6 time points) for the analysis of "total depletion", "total appearance", "total net balance" and "total turnover" of AAs.

Univariate analysis to evaluate the impact of EGF and urea on the depletion/appearance of EAAs by COCs matured in groups.
Compared to COCs matured in the presence of a high urea level (40 mg/dl) plus EGF (10 ng/ml), higher amounts of Thr (p < 0.001), Phe (p < 0.01), Ile (p < 0.001), Val (p < 0.001), Leu (p < 0.001), Lys (p < 0.001), and Trp (p < 0.001) were released by COCs in the presence of 10 ng/ml EGF (Fig. 1b). Furthermore, Ile and Lys were significantly depleted by COCs matured in the presence of EGF and a high level of urea; however, these EAAs were produced by COCs in the presence of EGF (p < 0.001). COCs matured in the presence of EGF (10 ng/ml) consumed more Met than other experimental treatments (p < 0.01). It was also noted that the total appearance (p < 0.001), total turnover (p < 0.001) and total net balance (p < 0.001) of EAAs was higher in EGF-incubated COCs, but the presence of EGF plus a high level of urea during in vitro maturation (IVM) removed this effect from the EGF (Fig. 1b).

Univariate analysis to determine the effect of EGF and urea on the depletion/appearance of non-EAAs, semi-EAAs and all AAs by COCs matured in groups. Compared with COCs matured
during IVM at a high urea level (40 mg/dl) plus EGF (10 ng/ml), COCs matured in the presence of EGF (10 ng/ ml) produced more Ala (p < 0.001) and Gly (p < 0.001) (Fig. 2a,b). In comparison, COCs matured in the presence of EGF plus a high level of urea increased (p < 0.001) the consumption of Arg and Gln (Fig. 2b). In addition, COCs matured in a culture medium with a high level of urea plus EGF had a higher net balance of non-EAAs (p < 0.001) than the EGF-incubated COCs (Fig. 3a). COCs matured in the presence of a high urea level plus EGF had a higher total depletion (p < 0.001) and the lower total appearance (p < 0.001) of semi-EAAs and all AAs compared to the EGF group (Fig. 3b,c). COCs incubated with a high urea level plus EGF showed a significant decrease (p < 0.001) in the total net balance of semi-EAAs/all AAs and in the total turnover of all AAs (Fig. 3b,c).
Scientific RepoRtS | (2020) 10:16082 | https://doi.org/10.1038/s41598-020-73225-7 www.nature.com/scientificreports/ Developmental competence of COCs and mRNA expression of genes in day-7 blastocysts resulting from bovine COCs matured in groups. All parameters of developmental competence were reported on the basis of the initial number of COCs cultured in the maturation medium. The results showed a significant increase in 2-PN (p < 0.05) and the percentage of COCs reached the cleavage stage (2-16 cell stage, p < 0.001) by applying EGF to the maturation medium (Fig. 4a). Adding EGF and urea to the culture medium during IVM decreased the percentage of 2PN, cleavage, and blastocyst formation (p < 0.01), whereas the percentage of oocyte degeneration (by day 3 post-insemination) increased (p < 0.001) compared to the EGF group (Fig. 4a). The mRNA expression of BCL2 was up-regulated in day-7 blastocysts resulting from COCs matured in the presence of EGF compared to the control group (p < 0.05, Fig. 4b). However, the presence of high urea level and EGF during IVM resulted in BCL2 mRNA down-regulation in the resulting blastocysts (p < 0.001). In the presence of EGF and a high level of urea in the medium during IVM, the mRNA expressions of BAX, NANOG, OCT4, and DNMT1 were significantly up-regulated (p < 0.001) in the resulting blastocysts (Fig. 4b).
Linear and quadratic impact of urea on the depletion/appearance of AAs by COCs and subsequent developmental competence. Linear and quadratic analyses were conducted as a result of linear use of urea in this study (0, 20, and 40 mg/dl, Supplementary Table 1). With an increase in urea, there was a quadratic decrease in the appearance of Glu (p = 0.02), His (p < 0.01), Gly (p = 0.04), Cit (p = 0.02), Ala (p < 0.01), Tyr (p = 0.03), Lys (p < 0.05) and Orn (p < 0.0001) (Supplementary Table 1). (b) Using the group-IVM culture (10 COCs per droplet), the 24-h incubation of COCs with urea and EGF altered the depletion/appearance of essential amino acids (EAAs). Differences between the concentration of AA in the test medium (after 24 h of incubation) and the fresh medium were used to calculate the depletion or appearance of each AA. Negative and positive values indicate the "depletion" and "appearance" of the EAAs, respectively. The "total depletion", "total appearance", "total net balance" and "total turnover" of EAAs was calculated using data from the depletion/appearance of the EAAs grouped into six replicates (based on the 6 time points of IVM conducted, Supplementary Fig. 1, Phase 1). The negative (depletion) or positive (appearance) values obtained for EAAs were summed up for each time point. Boxes range from the 25th to the 75th percentile of the distribution of values for each group. The median is represented by the line inside the box. The whiskers represent the maximum and minimum data point inside 1.5 times the interquartile range. Dots denote values outside the range of adjacent values. The data was obtained from 17 measurements and the medians were analyzed using the Kruskal-Wallis test followed by the Dunn's multiple comparison test. Urea and EGF are shown as U and E, respectively. Cont: control group; EGF group; 10 ng/ml EGF; E + U20: 10 ng/ ml EGF + 20 mg/dl urea; and E + U40 group: 10 ng/ml EGF + 40 mg/dl urea group; Asterisks indicate significant differences between groups: *: p < 0.05; **: p < 0.01; ***: p < 0.001. Scientific RepoRtS | (2020) 10:16082 | https://doi.org/10.1038/s41598-020-73225-7 www.nature.com/scientificreports/ Total appearance, total net balance and total turnover of non-EAAs (p = 0.009) as well as total turnover of semi-EAAs (p < 0.05) decreased quadratically with an increase in urea levels. The appearance of Val (p = 0.006) and Trp (p = 0.03) and the total appearance (p < 0.001), total net balance (p = 0.004) and total turnover (p = 0.0006) of EAAs as well as 2-PN (p < 0.0001), cleavage (p < 0.0001) and the percentage of blastocyst formation (as a proportion of the initial cultured COCs) decreased linearly (p < 0.0001) with increasing urea levels. Total depletion of non-EAAs (p = 0.009) and the percentage of COCs degeneration by day 3 post-insemination increased linearly (p < 0.0001) with an increase in urea levels (Supplementary Table 1). Total appearance, total net balance and total turnover of all AAs and total appearance of semi-EAAs decreased linearly (p < 0.001, p < 0.0001, p < 0.05, and p < 0.01, respectively) and quadratically (p = 0.006, p = 0.008, p < 0.05, and p = 0.0009), with increasing urea levels. Arg depletion increased linearly (p < 0.01) and quadratically (p = 0.005) with an increase in urea levels (Supplementary Table 1).
In the urea group, the non-cleaved COCs consumed significantly more Ile, Lys, Cit, Glu, Asp, Arg, and Tyr and had a higher total depletion, net balance and turnover of all AAs (p < 0.05) compared to the cleaved COCs (Table 1).
In the EGF-urea (EU) group, Met, Ile, Lys, Glu, Asp, Arg, Gly, and Gln depletion was higher in the non-cleaved COCs than in the cleaved COCs (p < 0.05). Non-cleaved COCs displayed a higher total depletion and net balance of all AAs and produced a lower amount of Ala compared to cleaved COCs (p < 0.05, Table 1).
Univariate analysis for the evaluation of cleavage and depletion/appearance of AAs by COCs matured individually. In the individual-IVM system, urea abrogated the EGF-induced the percentage of COCs reached to the cleavage stage. It was found that the percentage of cleavage decreased from 66.7% (the EGF group) to 44.8% (the EGF + urea group, approximately 66%, Table 1).

Figure 2.
Using the group-IVM culture (10 COCs per droplet), the 24-h incubation of COCs with urea and EGF altered the depletion/appearance of (a) non-essential amino acids (non-EAAs) and (b) semi-EAAs. Differences between the concentration of AA in the test medium (after 24 h of incubation) and the fresh medium were used to calculate the depletion or appearance of each AA. Negative and positive values indicate the "depletion" and "appearance" of the AAs, respectively. Boxes range from the 25th to the 75th percentile of the distribution of values for each group. The median is represented by the line inside the box. The whiskers represent the maximum and minimum data point inside 1.5 times the interquartile range. Dots denote observations outside the range of adjacent values. The data was obtained from 17 measurements and the medians were analyzed using the Kruskal-Wallis test followed by the Dunn's multiple comparison test. Urea and EGF are shown as U and E, respectively. Cont: control group; EGF group: 10 ng/ml EGF; E + U20 group: 10 ng/ml EGF + 20 mg/dl urea; and E + U40 group: 10 ng/ml EGF + 40 mg/dl urea. Asterisks indicate significant differences between groups: *: p < 0.05; **: p < 0.01; ***: p < 0.001. Figure 3. Total depletion, total appearance, total turnover and total net balance of (a) non-essential amino acids (EAAs), (b) semi-EAAs, and (c) all AAs of bovine COCs. COCs were matured for 24 h in the group-IVM culture supplemented by EGF and urea. The "total depletion", "total appearance", "total net balance" and "total turnover" of EAAs was calculated using data from the depletion/appearance of the EAAs grouped into six replicates (based on the 6 time points of IVM conducted, Supplementary Fig. 1, Phase 1). Then, the negative (depletion) or positive (appearance) values obtained for EAAs were summed up for each time point. Therefore, there were 6 data from 6 time points for the analysis of "total depletion", "total appearance", "total net balance" and "total turnover" of AAs. The data was obtained from 17 measurements and the medians were analyzed using the Kruskal-Wallis test followed by the Dunn's multiple comparison test. Boxes range from the 25th to the 75th percentile of the distribution of values for each group. The median is represented by the line inside the box. The whiskers represent the maximum and minimum data point inside 1.5 times the interquartile range. Dots denote observations outside the range of adjacent values. Urea and EGF are shown as U and E, respectively. Cont: control group; EGF group: 10 ng/ml EGF; E + U20: 10 ng/ml EGF + 20 mg/dl urea; and E + U40 group: 10 ng/ml EGF + 40 mg/dl urea group. Asterisks indicate significant differences between groups: *: p < 0.05; **: p < 0.01; ***: p < 0.001. www.nature.com/scientificreports/ N-acetyl-l-cysteine completely abrogates urea-induced depletion of AAs. Data showed that the NAC completely abrogated the urea-induced depletion of AAs (p < 0.05, Fig. 4c).

Scientific RepoRtS
Bivariate analysis to assess the association between AAs and the developmental competence of COCs matured individually or in groups. Using group-matured COC data, the Kendall correlation showed a correlation between 2-PN and Arg (R = 0.73, p < 0.05, Fig. 5a and Supplementary Fig. 2). It was also found that the percentage of COCs reached the cleavage stage (2-16 cell stage) was correlated with Arg (R = 0.58, p < 0.05). The Kendall correlation analysis did not find any significant correlation between the blastocyst formation and AAs. In addition, Arg (R = 0.40, p = 0.07), Leu (R = 0.34, p = 0.12) and Thr (R = 0.37, p = 0.09) tended to be positively associated with the blastocyst formation. Moreover, none of the AAs were associated with the The median is represented by the line inside the box. The whiskers represent the maximum and minimum data point inside 1.5 times the interquartile range. Dots denote observations outside the range of adjacent values. The data was obtained from 17 measurements and the medians were analyzed using the Kruskal-Wallis test followed by the Dunn's multiple comparison test. Asterisks indicate significant differences between groups: *: p < 0.05; **: p < 0.01; ***: p < 0.001. All parameters relating to developmental competence of the COCs were shown in relation to the initial number of COCs cultured in the maturation medium. www.nature.com/scientificreports/ blastocyst formation in each experimental treatment. The pairwise comparison heatmap (Fig. 5a) exhibited that oocyte degeneration (by day 3 post-insemination) was negatively correlated with Arg (R = − 0.62, p < 0.05). Table 1. Depletion/appearance of amino acids by COCs matured the individual-IVM culture. Differences between the concentration of AA in the test medium (after 24 h incubation) and the fresh medium were used to calculate the depletion or appearance of each AA. Negative and positive values indicate the "depletion (Dep)" and "appearance (App)" of AAs, respectively. The "total net balance (Net)" of all AAs was calculated by subtracting "total depletion" from "total appearance". In addition, the "total turnover (Turn)" of all AAs was calculated by the sum of "total depletion" and "total appearance". The data were obtained from 6 measurements. Asterisks indicate significant differences between cleaved vs. non-cleaved groups within each experimental treatment; data was normally distributed and analyzed using t-test with a Bonferroni correction (*: p < 0.05; **: p < 0.01; ***: p < 0.001). Different capital or small letters indicate significant differences among cleaved or non-cleaved groups across all experimental treatments, respectively (p < 0.05, one-way ANOVA followed by Tukey's multiple comparison test, n = 6). Data are means ± S.E.M. # indicates the percentage of cleaved or non-cleaved COCs in relation to the initial number of COCs cultured in the maturation medium. www.nature.com/scientificreports/ The Pearson correlation revealed that the percentage of cleaved COCs matured individually (relative to the initial cultured COCs) was correlated with Tyr (R = 0.71, p < 0.05) and Thr (R = 0.58, p < 0.05). Moreover, Glu (R = 0.57, p = 0.06), Cit (R = 0.45, p = 0.14), Val (R = 0.53, p = 0.08), Trp (R = 0.49, p = 0.10) and Phe (R = 0.55, p = 0.07) tended to be correlated with cleaved COCs (Fig. 5b and Supplementary Fig. 3).
Pearson analysis showed that the percentage of non-cleaved COCs matured individually (relative to the initial cultured COCs) was negatively correlated with the turnover of all AAs. Data showed a correlation between the percentage of non-cleaved COCs and the turnover of Asp (R = − 0.71, p < 0. The heatmap analysis supports the univariate analysis findings. The color intensity of the heatmap represents the depletion and appearance of the AAs and the developmental competence of the individual or group matured COCs ( Fig. 6a-c). Green represents a higher depletion (or lower appearance) of AAs or a lower percentage of other variables, i.e., oocyte competence. Red represents a higher appearance of AAs or higher levels of oocyte competence.
The heatmap analysis showed that, in the presence of EGF and a high level of urea, the turnover of all AAs (except Ser) by COCs (matured in groups) was higher than other experimental treatments. In addition, the heatmap analysis verified the results of the ANOVA test in such a way that oocyte degeneration increased in the presence of EGF and a high level of urea during IVM, but the percentage of COCs reached the cleavage, 2-PN and blastocyst stages decreased (Fig. 6a).
The heatmap analysis showed that when COCs were matured individually, the depletion of AAs by COCs was greater in the urea group (both cleaved and non cleaved COCs) than in other treatments whereas the noncleaved COCs depleted substantially more AAs than the cleaved COCs (Fig. 6b,c). In the case of individual-IVM culture, EGF abrogated the urea-increased depletion of AAs by COCs, but this effect was not achieved in the group-IVM culture (Fig. 6b,c).
A hierarchical cluster analysis, a multivariate approach, to classify the related variables. Using the Kendall and Pearson distance metric and a complete-linkage clustering method, cluster analysis was done to classify factors that are close together (Fig. 6a-c). Since the data from COCs matured in groups were not normally distributed, we used Kendall distance metric as a non-parametric test. The Pearson distance metric was used to identify the normal data obtained from individually matured COCs.
In the case of the group-IVM culture, the hierarchical cluster analysis revealed the similarity (lowest distance) between Arg and the percentage of 2PN and cleavage, since they were linked together (black box in Fig. 6a). Percentage of 2PN, cleavage, blastocyst formation, Arg, Trp, Val, Ser, and Gln all together formed a cluster (black box in Fig. 6a). Oocyte degeneration and His formed a distinct cluster together (red box in Fig. 6a).
In the case of the individual-IVM culture, the hierarchical dendrogram showed a correlation between the percentage of the cleaved COCs and Trp (black box in Fig. 6b). The hierarchical dendrogram showed the similarity between Orn, Ala and percentage of non-cleaved COCs (black box in Fig. 6c). The Anderson-Darling test showed that the group-IVM culture data was not normally distributed and that the individual-IVM culture data was normally distributed. Therefore, the Pearson correlation and Kendall rank correlation was used for normal and log-transformed data, respectively. Boxed squares mean a significant correlation (p < 0.05). The Bonferroni correction was used to test the significance of the correlations in the correlation matrix. Red and blue circles display negative and positive correlations, respectively. Scale bar on the right side of the pairwise heatmap shows the correlation between the factors. Deg percentage of oocyte degeneration, Cleav cleavage stage, Blast blastocyst stage, 2PN two-pronuclear zygote, Cl cleaved COCs, NC non-cleaved COCs.
Scientific RepoRtS | (2020) 10:16082 | https://doi.org/10.1038/s41598-020-73225-7 www.nature.com/scientificreports/ The principal component analysis, a multivariate analysis, to assess the association between AAs and the developmental competence of COCs. The principal component analysis (PCA) was performed in order to reduce the broad range of variables (high dimensional data) and to find a concurrent association between AAs and oocyte competence. Data of AAs depletion/appearance by COCs matured individually or in groups were analyzed using PAST software. The first three principal axis factors accounted for a sufficient amount of total variation (77.8%, 90.6%, and 96.4% found for data of COCs cleaved and non-cleaved in the individual-IVM culture or in group-IVM culture, respectively, Supplementary Table 2).
In the case of group-IVM culture, PCA identified a strong positive correlation between Trp, Arg, Gln, Ile, Val and the percentage of COCs reached the 2PN, cleavage and blastocyst stages (directional vectors < 45°, Fig. 6d). There was a strong negative correlation between Trp, Arg, Gln, Ile, Val and oocyte degeneration (directional vectors approaching 180°). The PCA showed a strong positive association between His and oocyte degeneration (Fig. 6d).
In the case of the individual-IVM culture, PCA identified a strong positive correlation between Orn, Trp, Ala, Glu, Asp, Phe, Tyr, Met and the percentage of cleaved COCs (directional vectors < 45°, Fig. 6e). The data showed a strong positive correlation between Orn, Ala, and the percentage of non-cleaved COCs (Fig. 6f). Moreover, PCA identified a strong negative correlation between Cit, Tyr, Arg, Ser, Gly, Lys,Asp, Glu, Leu and the percentage of non-cleaved COCs (Fig. 6f). Such observation of the PCA confirmed the results of the hierarchical cluster analysis (Fig. 6a-c).
Sub-network enrichment analysis describes the correlation between Arg and the incidence of high developmental competence. The nodes of the enriched sub-modules indicate AAs and parameters representing developmental competence (Fig. 7). The node size reflects the degree of connectivity of the node (i.e., AAs and parameters representing developmental competence) and the edges show the relationship between the two variables. The thicker edges reveal higher correlations between variables. Nodes with more links are close to each other.
In the case of the group-IVM culture, the amino acid-developmental competence interaction networks with the cutoff point > 80% showed that Arg was associated with the percentage of COCs reached the 2PN, cleavage www.nature.com/scientificreports/ and blastocyst stages (Fig. 7a). More specifically, when the cutoff point increased from 80 to 90%, only Arg had an association with the percentage of COCs reached the 2PN stage (Fig. 7b).
In the case of individual-IVM culture, the amino acid-cleavage interaction networks with the cutoff point > 82% showed that the cleaved COCs were associated with Trp, Met, His and Phe (Fig. 7c). In addition, a network analysis > 72% of the cutoff point showed an association between Trp, Orn, and non-cleaved COCs (Fig. 7d).
The receiver operator characteristic (ROC) curve analysis. To be more specific, the ROC curve study was performed to determine the predictive power of variables (i.e., AAs) and to calculate the optimal cutoff for the positive outcome prediction. In the case of group-IVM culture, AAs were divided into two separate groups based on IVF outcomes (i.e., high 2-PN percentage > 80% vs. low 2-PN percentage; high percentage of cleavage > 80% vs. low percentage of cleavage; high percentage of blastocyst formation > 30% vs. low percentage of blastocyst formation; and low percentage of oocyte degeneration < 15% vs. high percentage of oocyte degeneration). In the case of individual-IVM culture, the AAs were divided into two classes based on the incidence of cleavage as follows: the cleaved vs. the non-cleaved COCs. The AUC, P-value, optimal cutoff point, sensitivity, specificity, positive/negative predictive values and positive/negative likelihood ratios obtained for the significant AUC are shown in Table 2. The results showed that the positive likelihood ratio (LR) of all factors was much more than 1.0. However, the negative LR of all factors was less than 1.0, suggesting that there was an increased likelihood of a high percentage of developmental competence 26 (Table 2).
In the case of group-IVM culture, the ROC curve analysis for predicting the existence of a high 2PN percentage > 80% showed that Arg had the best predictive power with a significant area under the curve (AUC) of 1.00, a sensitivity of 100% and a specificity of 100%. Moreover, the depletion level of Arg > -0.60 pmol/h was the best optimal cutoff point to predict a high 2PN percentage ( Table 2, Supplementary Fig. 5a).
The AUC curve analysis was also used to determine the diagnostic value of AAs by COCs matured in groups for predicting the presence of low oocyte degeneration (i.e., < 15%). We observed that Arg had the highest www.nature.com/scientificreports/ predictive performance with a strong AUC of 0.91, a sensitivity of 85.7% and a specificity of 100%. The ROC analysis defined the depletion levels of Arg > − 0.60 pmol/h as the optimal cutoff for estimating the incidence of low oocyte degeneration ( Table 2, Supplementary Fig. 5b). Val (AUC = 1.00), Arg (AUC = 0.96), Ser (AUC = 0.93), Trp (AUC = 0.89), Thr (AUC = 0.89), Ile (AUC = 0.92) and the total turnover (AUC = 1.00), total appearance (AUC = 1.00) and total net balance (AUC = 0.89) of EAAs were the factors with the highest AUC for detecting a high percentage of cleavage > 80%. Moreover, the total depletion (AUC = 0.89) and net balance (AUC = 0.85) of all AAs and total depletion (AUC = 0.85) and net balance (AUC = 0.85) of semi-EAAs showed the best AUC for predicting a high percentage of cleavage > 80%. Furthermore, the optimal cutoff points for predicting a high percentage of cleavage were > 1.16, > -0.17, > -1.65, > 0.2, > 1.32, and > 0.32 pmol/h for Val, Arg, Ser, Trp, Thr and Ile, respectively ( Table 2 and Supplementary Fig. 5c). Table 2. Rank of amino acids turnover (depletion/appearance) based on their discriminating power (AUC) from ROC curve analysis. The ROC curve study was performed to determine the predictive power of factors (i.e., AAs) and to measure the optimal cutoff point for the positive outcome prediction. The AAs were divided into two separate groups based on the IVF outcomes of the group-IVM system (i.e., high 2PN percentage > 80%, high percentage of cleavage > 80%, high percentage of blastocyst formation > 30%; low percentage of oocyte degeneration < 15%) and the individual-IVM culture (i.e., incidence of cleavage vs. absence of cleavage). Dep depletion, Net net balance, App appearance, Turn turnover, Total all amino acids, Non non-essential amino acids, Semi semi-essential amino acids, EAA essential amino acids, 2PN twopronuclear zygote. www.nature.com/scientificreports/ Evaluation of the ROC curve for the presence of a high percentage of blastocyst formation (i.e., > 80%) showed that Arg, Leu, Gln, Ile, and Thr had the highest predictive ability performance with substantial AUC values of 0.84, 0.84, 0.84, 0.81, and 0.74, respectively. Moreover, the optimal cutoff for predicting a high percentage of blastocysts formation was estimated as > − 0.60, > 0.97, > − 8.23, > 0.57, and > 1.35 pmol/h for Arg, Leu, Gln, Ile, and Thr, respectively (Table 2, Supplementary Fig. 5d).
In silico molecular docking study reveals the molecular association of urea with DNA strains, EGF, EGFR and aquaporin (AQP) 3. It was evident from the ANOVA analysis that, in the presence of urea, oocyte degeneration increased and EGF-induced cleavage formation was abolished. These findings encouraged us to investigate whether urea interacts with DNA; this may help us to explain some of the increased percentage of oocyte degeneration found in this study. Next, we sought to find the feasibility of association between urea and EGF/EGFR; this was achieved to determine how urea would bind with these molecules. The potential interaction between urea and these molecules will allow us to better understand the interference of urea in the EGF-EGFR binding. In order to study changes in the energy binding of urea to DNA, EGF, EGRFR, and AQP3, with an increase in number of urea molecules, molecular docking was conducted 20, 12, 12 and 12 consecutive times, respectively. Indeed, each time the previous docked EGFR-urea, EGF-urea, or DNA-urea complex was used as a receptor for the next docking run under the same conditions. Molecular docking showed that free energy of binding was negative for all urea interactions with DNA, EGF, extracellular domain of EGFR, and AQP3 (Fig. 8a- Fig. 8e). It suggested, therefore, that urea had an increased affinity for DNA and EGF.
Next, in the presence of urea, molecular docking was performed again to compare the change in the free energy binding of EGF-EGFR complex. To this aim, the previous docked EGFR-urea or EGF-urea complex was used as a receptor for these docking processes. The total interaction energy(E tot ) was − 64.  (Fig. 8f). This finding suggested the strongest interactive impact of urea on the EGFR molecule.
The molecular docking analysis showed the interaction of urea with DNA bases (e.g., thymine and cytosine; Fig. 8g). Using LIGPLOT software, we found hydrogen bonds and non-covalent contacts between urea and some EGF residues and the extracellular domain of EGFR (Fig. 8h,i). The nitrogen atoms in urea molecules were mainly involved in hydrophobic interactions with the residues, such as Cys33, Gly39, Leu132 and Phe156. Furthermore, there are some hydrogen bonds between hydrogen atoms of urea molecules and oxygen atoms of residues (red arcs in Fig. 8h,i). In addition, the LIGPLOT showed the existence of hydrogen bonds between urea and some residues of EGF and EGFR, such as Cys31, Asn129, and Glu40 (green dashes; Fig. 8h,i). It was observed that there were several van der waals volume interactions between urea and AQP3, like Val43, Gly115 and Arg218 (Fig. 8k).

Discussion
In this study, the conditions of IVM modified the relationship between the metabolism of AAs and the developmental competence of COCs. Data showed that, in response to the presence of EGF and urea in the group-IVM culture, COCs released lower quantities of Thr, Phe, Val, Leu, Trp, Gly and Ala, consumed more Ile, Lys, Gln and Arg, and showed lower developmental competence. Moreover, the depletion of AAs was the highest in ureaincubated COCs with greater depletion in non-cleaved COCs than in cleaved COCs using an individual-IVM culture. Recently, we observed a greater depletion of AAs and a decrease in oocyte competence in the presence of high urea levels in the group-IVM system 17 . Similarly, Hemmings et al. reported that higher AA depletion is associated with a decrease in the developmental competence of oocytes in humans and bovines 18,19 . A growing body of evidence suggests that non-cleaved oocytes consume (deplete) more AAs for DNA repair and prevention of degeneration 18,27 . Greene et al. reported that the addition of Arg to the cell culture reduced the fragmentation of DNA in human endometrial cells 28 . It seems that COCs may increase the consumption of AAs to meet their needs in order to cope with the negative effects of urea as a cell-stressor 18,27 .
As mentioned above, more AAs were depleted by COCs in the presence of EGF and a high level of urea using the group-IVM system (Fig. 6a). In view of the fact that the EGF may inhibit AAs transport 23 , it is proposed that the increased depletion of AAs by COCs can result from the inhibiting action of urea on the EGF/EGFR system [9][10][11] . COCs were matured in the presence of NAC in order to test the potential interactive effect of urea on AAs metabolism via EGFR. It has been reported that NAC may inhibit the activation of EGFR 29 . Moreover, urea signaling has been shown to be sensitive to NAC 29 . Current data showed that NAC completely abrogated the urea-induced depletion of AAs by COCs. This meant that urea could use the EGFR system to increase the consumption of AAs by COCs. In addition, as shown by the present in silico molecular docking, both EGF and EGFR could be bound by urea, suggesting an interference role for urea in the EGF/EGFR system. Moreover, molecular docking showed that E tot for the EGF/EGFR-urea complex (− 62.96 kcal/mol) was the lowest compared to the EGF/EGFR complex (− 64.76 kcal/mol) and the EGF-urea complex/EGFR (− 69.49 kcal/mol). It seems that, Scientific RepoRtS | (2020) 10:16082 | https://doi.org/10.1038/s41598-020-73225-7 www.nature.com/scientificreports/ in the presence of urea, the extracellular domain of EGFR may have a lower affinity to the EGF. This indicated that urea interfered with the interaction between EGF and EGFR. Therefore, urea can reduce the inhibitory effect of EGF on AAs transport, resulting in an increase in the depletion of AAs by COCs when matured in the presence of urea and EGF. Next, an individual-IVM culture was used to see the potential differences in the metabolism of AAs based on the developmental competence of COCs. Compared to the group-IVM culture in which the depletion of AAs was the greatest in the presence of high-level urea + EGF, we found that the depletion of AAs was the highest in the high-level urea group in the individual-IVM culture. It has been shown that the culture conditions affect the metabolism of human oocytes 30 . EGF has been shown to increase the blastocyst formation in a single embryo culture system, although EGF did not affect the development of embryos in group culture 31 . The present results showed that the EGF abrogated urea-increased depletion of AAs in the individual-IVM system but synergistically increased depletion of AAs in the group-IVM culture with urea. Mutual interaction between oocyte and CCs may have contributed to different effects of EGF in group and individual IVP systems 32,33 . In other words, the CCs and oocyte maintain the supply of paracrine and autocrine factors, such as EGF, for each other 25 . For example, CCs have been shown to mediate the beneficial effects of EGF on bovine oocyte maturation 2-4 . Therefore, the Optimized urea was docked with these molecules by AutoDock 4. 2. Using Lamarkian Genetic Algorithm, molecular docking was performed to determine the binding free energy and the orientation of urea in its interaction with these molecules. The crystal form of the extracellular domain of EGFR (PDB ID: 3NJP), EGF (1jl9), and DNA (1 BNA) was taken from the Protein Data Bank. Computer modeling on AQP3 (ID: Q08DE6) was performed by providing amino acid sequences from the Universal Protein Resource (UniProt) database, followed by predicting a 3-D structure on the I-TASSER server. To assess the urea binding behavior, molecular docking was conducted at 20, 12, 12, and 12 sequential times for urea-DNA, urea-EGF, urea-EGFR, and urea-AQP3, respectively. Indeed, each time the previous docked urea-receptor complex was used as a receptor for the next docking run under the same conditions. These images (a-d) were produced by the VMD molecular visualization and analysis software (V. 1.9.2.). (e) Box plot to display the energy binding of urea to DNA, EGF, EGFR, and AQP3. (f) In order to investigate how the energy binding of EGF with EGFR shifts in the presence of urea, molecular docking was performed in four phases as follows: phase I: EGF/EGFR; phase II: (EGFurea complex)/EGFR; phase III: EGF/(EGFR-urea complex); and phase IV: (urea-EGF complex)/(EGFR-urea complex). The EGF-urea and EGFR-urea complexes obtained from the preceding step were used in these phases. In the panels (g) and (k), the ball represents the van der waals volume interactions. These images (g,k) were produced by the MGLTools software (V. 1.5.6.). Hydrogen bonds are shown in panels (h) and (i) by dashed green lines. Also, the distance between the receptor and the donor atoms is shown in the middle. Hydrophobic contacts are represented by arcs with spokes which radiate towards urea molecules. LIGPLOT software (V. www.nature.com/scientificreports/ cooperative relationship between CCs and oocytes is impaired in the individual-IVM system; which may lead to a different response to the presence of EGF in the maturation medium. In both individual-and group-IVM systems, urea abrogated the EGF-induced the percentage of COCs reached to the cleavage stage. It was found that the percentage of cleavage decreased from 90% found for the EGF group to 54% found for the EGF + urea group, approximately a decrease of 60% in the group-IVM culture (Fig. 4a). However, the percentage of cleavage decreased from 66.7% found for the EGF group to 44.8% found for the EGF + urea group, approximately a decrease of 66% (Table 1). In order to test the interference function of urea on the EGF/EGFR system, two in silico and in vitro experiments (with the use of an EGFR inhibitor, i.e., NAC) were carried out. The in silico data indicated an interference function for urea in binding of EGF to its receptor (EGFR). In addition, an in vitro experiment using NAC showed that urea can act through EGFR. Therefore, these data show that urea can interfere with the positive effects of the EGF/EGFR system in terms of oocyte developmental competence (i.e., the cleavage stage).

2.2.) was used to produce images and visualize hydrophilic interactions, hydrogen bonding, and hydrophobic contacts between urea and EGF (h) and EGFR (i). E, EGF
Data showed that the developmental competence of COCs matured in the group-IVM culture was improved upon maturation in the presence of the EGF. It was found that the percentage of cleavage increased from 70 to 90% (approximately 28%) in the presence of the EGF in the group-IVM culture. EGF has been reported to mediate oocyte maturation by phosphorylation of mitogen-activated protein kinase (MAPK) 34,35 . In fact, the suppression of EGFR signaling abolishes the acquisition of oocyte developmental competence 34,35 . Previous studies have reported that the MAPK pathway and the EGFR signaling are inhibited by urea and its derivatives 9-11,36 . Oyamada et al. reported that, based on the form of in vitro culture (IVC, i.e., individual-IVC or in group-IVC), the presence of EGF in the maturation medium changes differently the early embryonic development in bovines 33 . These findings again indicate that the nature of the IVM culture (individual system vs. group system) and the composition of the culture medium (i.e., the presence of EGF) may have an effect on the competence of the COCs.
As mentioned above, the presence of EGF in the maturation medium improved the subsequent fertilization and cleavage percentage of COCs but did not significantly improve the percentage of COCs reached the blastocyst stage compared to the control group. However, the percentage of blastocyst formation in the EGF group was numerically higher than in the control group (31.03% versus 28.47% for the EGF and the control group, respectively). Similarly, Merriman et al. observed that the fertilization was increased, while the subsequent developmental competence of mouse oocytes matured in the presence of EGF during IVM was impaired 37 . They claimed that the EGF would not be capable of supporting embryonic development 37 . On the other hand, Lonergan et al. reported that the presence of EGF during IVM significantly increased the percentage of cleavage and blastocyst formation in bovine COCs isolated from small follicles (2-6 mm) 2 ; however, it is known that COCs originating from small follicles (2-4 mm) 35 are generally non-responsive to the EGF 8,35,38 . In this study, the COCs were recovered from medium sized follicles (7-8 mm) and were expected to react more strongly to the EGF. In addition, in the Lonergan et al. study, COCs were matured in a maturation medium containing 10% FCS, which can contain a number of macromolecules, hormones and growth factors 2 . These unknown factors may have an effect on the consequent developmental competence of the COCs, while it was expected that small follicle-derived COCs would not respond to the EGF. In this study, heat-inactivated FBS was used to avoid unknown reactions. It is proposed that differences between the various reports may be due to the fact that conditions during IVM make the resulting embryos more reactive or non-responsive to the culture conditions used for embryo development 37 .
We have previously reported that high urea levels not only increased the percentage of oocyte degeneration, but also decreased the in vitro viability of CCs and oviduct epithelial cells 16,17,39 . In this study, the addition of urea to the maturation medium increased the percentage of oocyte degeneration by day 3 after insemination. Because bovine oocyte expresses aquaporin 3 (AQP3) 40,41 , small solutes, such as water and urea, can readily enter the oocyte. The present data showed that the intracellular content of urea in DOs (cultured for 4 h with 40 mg/dl of urea) was 3.43 mg/dl. In addition, using an in silico method, we found the binding affinity of urea to AQP3 (− 2.79 kcal/mol). This means that urea can enter oocytes using AQP3, providing the opportunity to interact directly with DNA. Moreover, we observed a strong interaction between the urea molecules and the DNA. Likewise, Qiu et al. also found a strong interaction of the hydrogen bond between the urea molecules and DNA bases 42 . Hence, it seems reasonable to assume that urea can bind directly to oocyte genetic material and cause oxidative damage to DNA 43 . This can, in part, explain the higher percentage of COCs degeneration.
In this study, the highest percentage of oocyte degeneration coincided with acute up-regulation of NANOG (16.7 folds), OCT4 (37 folds), DNMT1 (15.5 folds) and BAX (15.7 folds) in blastocysts derived from COCs matured in groups and in the presence of high-level urea and EGF. Previous studies have shown that bovine blastocysts with higher or lower expression of OCT4 and NANOG than normal expression have a poorer potential for further development 44,45 . In addition, the expression of BCL2 has been shown to be higher in bovine oocytes of good quality 46 . Over-expression of BAX is well known to be associated with accelerated apoptotic death 46 . Moreover, the BCL2 to BAX ratio can be used to estimate oocytes and embryos for apoptosis or survival 46 . Thus, the 24-h incubation of COCs with high urea levels and the EGF altered the normal pattern of gene expression in the resulting blastocysts. This was also part of the increased degeneration of COCs matured in the presence of urea and EGF.
More specifically, all data (regardless of the experimental treatments) were subjected to a bivariate and PCA analysis to assess the extent and direction of the relationship between AAs and oocyte competence. Indeed, the ANOVA tests cannot cope with complex treatment structures and multidimensional data, i.e., omics data 47 . The Kendall rank correlation identified a significant association between Arg and oocyte developmental competence or degeneration of COCs matured in groups. Interestingly, the bivariate analysis did not detect any significant correlation between the percentage of blastocyst formation and AAs. As regards the relationship between Trp, Val, Arg and the early stages of development (i.e., cleavage stage) but not the blastocyst stage, the advanced stages of reproduction (e.g., blastocyst development, implantation and pregnancy rates) must be recorded upon www.nature.com/scientificreports/ maturation of COCs. Indeed, embryos developing into the blastocyst stage undergo spontaneous selection during IVC 48 . It indicates that the stage of cleavage does not guarantee the development of embryos after a short-term culture (reflected at a high cleavage proportion) but may reduce the developmental competence of embryos after an extended culture (i.e., 6-7 days) 49 .
In the case of the individual-IVM system, the Pearson correlation identified a significant correlation between Tyr, Thr and the cleaved COCs; this result was different from that observed in the group-IVM system. Next, the PCA analysis was used to evaluate the complex interactions between AAs and the developmental competence of COCs matured individually or in groups 20 . In the case of the group-IVM system, the PCA analysis demonstrated a strong positive association between Arg, Gln, Ile, Trp, Val, and the percentage of COCs reached the 2-PN, cleavage and blastocyst stages (Fig. 6d). However, PCA detected a strong positive association between Ala, Orn, Trp, Glu, Asp, Phe and the proportion of cleaved COCs when COCs matured individually (Fig. 6e). This finding indicated that COCs maturation individually or in groups could change the metabolism of AAs by COCs. Moreover, multivariate analysis appears to outperform bivariate analysis in finding additional factors linked to oocyte developmental competence.
The ROC curve analysis was conducted to identify the best predictor for COCs with a high developmental competence in the group-IVM system and the incidence of cleavage in the individual-IVM system. The AUC analysis showed that, in the case of the group-IVM system, Arg was a biomarker with a good predictive ability to distinguish highly competent oocytes (i.e., high percentage of 2PN, cleavage and blastocyst formation > 80%) from those with low developmental competence (i.e., < 80%). Arginine displayed strong AUC (0.9, on average), sensitivity (87.4%, on average) and specificity (100%, on average) values. Similarly, cluster and network analysis also identified Arg and the percentage of 2PN, cleavage and blastocyst formation in the same cluster. In the network analysis, when an 80% cutoff point was set for network analysis, only Arg was linked to the proportion of COCs reached the 2PN and cleavage stages. Next, to increase the specificity, the cutoff point was increased from 80 to 90% and the results showed that only Arg was related to 2-PN (Fig. 7b). As mentioned above, the increase in Arg depletion coincided with a decrease in the competence of the COCs. Bódis et al. reported an increased level of Arg in follicular fluid, which adversely affected human embryos 50 . They suggested that increased activity of the L-Arg/NO system may have an adverse effect on IVF outcomes 50 . The present findings suggest that Arg depletion (i.e., > − 0.6 pmol/h) may be used as a biomarker during IVM to detect high-competent COCs.
In the case of the individual-IVM culture, the AUC analysis showed that Glu, Ala, Ile, Trp, Met, Orn, Phe, Tyr, and Arg were biomarkers with a good predictive ability to estimate the incidence of COC cleavage. Moreover, the AUC detected Orn and Ala as biomarkers with a good predictive ability to estimate the absence of COC cleavage. These findings were consistent with the results of the PCA, network analysis and the hierarchical dendrogram. The ROC analysis shows that the factors that predict the incidence of cleavage or identify high-competent COCs will also change with changes in the conditions of the IVM.
In conclusion, it can be deduced from these sets of analysis that the conditions of the IVM (i.e., the existence of interactive agents, such as urea and EGF, or the use of individual-IVM and group-IVM cultures) may somehow change the AAs depletion/appearance and developmental competence of the COCs. These findings may provide some insight into the role of urea in reduced fertility of high producing dairy cows fed high-protein diets, but further in vivo experiments are needed to confirm theses in vitro results.

Methodology
Animal studies and data statement. All reagents were obtained from Sigma Aldrich (St. Louis, MO, USA) unless otherwise stated. Animal experiments were conducted in accordance with the Guiding Principles for the Care and Use of Research Animals Promulgated by the Isfahan University of Technology, Iran. The protocol and methods were approved by the Committee on the Ethics of Animal Experiments of the Isfahan University of Technology (No. 390132). The datasets are available from the corresponding author upon request.
In vitro maturation and treatments. From November 2018 to February 2019, bovine ovaries were sourced from the local slaughterhouse and transported to the laboratory within 2 h at about 30 °C using a thermo box containing 0.9% NaCl, 0.1% penicillin, and streptomycin. Antral follicles (7-8 mm in diameter) were isolated from ovaries containing an active corpus luteum and kept at 37 °C in HEPES-TCM-199 (tissue culture medium; Invitrogen, Carlsbad, CA, USA).
In the group-IVM culture (group: 10 COCs per 50-μl droplet), a total of 960 COCs derived from 7-8 mm follicles (from 150 slaughtered cows) were randomized into four experimental groups. All treatments were repeated six times (6 independent experiments). At each replication, COCs from different cows were pooled and randomly assigned into treatment groups. On average, a total of 24 wells as the experimental units (24 droplets with 10 COCs per each) were assigned to each experimental treatment ( Supplementary Fig. 1, Phase 1). Only COCs with at least five cell layers were selected under a stereomicroscope (Olympus, Tokyo, Japan) and transferred to HEPES-TCM-199 supplemented by 50 mg/ml of kanamycin and 50 mg/ml of heparin. Intact COCs were then washed twice in HEPES-TCM199 and randomly placed in groups of 10 in 50-µl droplet of maturation medium containing bicarbonate-buffered TCM199 supplemented by 10% (w/v) of heat-inactivated fetal bovine serum (FBS, BioWhittaker, Walkersville, MD, U.S.A.) and 20 IU/ml follicle-stimulating hormone (FSH) in 35 mm cell culture dishes (Falcon brand, BD Biosciences, San Jose, CA, U.S.A.) for 24 h at 38.5 °C under 5.5% CO 2 , 20% O 2 , balanced N 2, and maximum humidity under mineral oil. The experimental treatments were as follows: (1) maturation medium without urea and EGF (control group, n = 250 COCs); (2) maturation medium with 10 ng/ ml EGF (E group, n = 230 COCs); (3) maturation medium with 10 ng/ml EGF and 20 mg/dl urea (EU20 group, n = 240 COCs); and (4) maturation medium with 10 ng/ml EGF and 40 mg/dl urea (EU40 group, n = 240 COCs). www.nature.com/scientificreports/ Following IVM, the 24-h spent maturation medium (24 samples from 24 wells per treatment) was kept at -70 until HPLC analysis.
In the individual-IVM culture (individual: 1 COC per 10-μl droplet), a total of 366 COCs were randomized to four experimental groups. Intact COCs were washed twice in HEPES-TCM199 and placed randomly individually in 10-µl droplet of the maturation medium overlaid with paraffin liquid. The individual-IVM culture was conducted under the same conditions as those described above for the group-IVM culture. The experimental treatments used for individual-IVM culture were as follows: (1) maturation medium without urea and EGF (control group, n = 90 COCs); (2) maturation medium with 10 ng/ml EGF (E group, n = 87 COCs); (3) maturation medium with 40 mg/dl urea (n = 93 COCs); and (4) maturation medium with 10 ng/ml EGF and 40 mg/dl urea (EU group, n = 96 COCs). Following IVM, the 24-h spent maturation medium (12 samples per treatment) was randomly selected from each treatment on the basis of the cleavage status (cleaved and non-cleaved COCs) and kept at − 70 until HPLC analysis.
NAC was used to test the possible effects of urea on AAs via EGFR system. NAC has been shown to inhibit EGFR activation 29 . Urea signaling has been shown to be sensitive to NAC 29 . Using the group-IVM culture (group: 10 COCs per 50-μl droplet), a total of 360 COCs were randomized into three experimental groups as follows: (1) maturation medium with no supplementation (control group, n = 60 COCs); (2) maturation medium with 40 mg/dl urea (n = 60 COCs); and (3) maturation medium with 40 mg/dl urea + 1 mM NAC (n = 60 COCs). The IVM was conducted under the same conditions as those described above for the group-IVM culture. After IVM, the 24-h spent maturation medium (6 samples per treatment) was randomly selected from each treatment and kept at -70 until HPLC analysis.
In the present study, COCs, but not DOs, MII-oocytes or in vitro-synchronized oocytes, were used for the assessment of AAs depletion/appearance. In addition, COCs were specifically drawn from 7-8 mm follicles to maintain the same maturation stage as possible and to minimize variability in the amount of cumulus layers. Since there is a cross-talk between oocyte and its accompanying CCs, which plays an important role in metabolism, in AAs uptake by oocytes, and in the subsequent developmental competence of oocytes 51 , COCs were not denuded in this study.
Based on the concentrations of urea in the blood and follicular fluid of dairy cows fed a low-or high-protein diet, 20 mg/dl urea (equivalent to 9.3 mg/dl BUN) and 40 mg/dl urea (equivalent to 18.7 mg/dl BUN) were used in this study 13,52,53 .Therefore, these levels of urea simulated the BUN values found in healthy dairy cows fed low-or high-protein diets, respectively 53,54 .In fact, the concentration of urea in follicular fluid is similar to that in plasma as it quickly spreads to body fluids, such as follicular fluid 13 . In this study, the EGF concentration (10 ng/ml) was based on the minimum effective dose reported in previous reports 2 . Follicular concentrations of EGF have been reported to vary from 2 to 15 ng/ml 2 .
Amino acid and urea determination. AAs were measured using HPLC, before HPLC analysis seven samples (seven out of 24 samples of 24-h IVM-wells per treatment) that were contaminated with oil were discarded; ultimately 17 droplets or wells (i.e., experimental units) per treatment were selected for the AAs analysis by HPLC (Supplementary Fig. 1, Phase 1). The determination of AAs in the 24-h spent media (n = 17) and fresh maturation media (n = 6) was conducted using the HPLC method (HyperClone ODS C 18, 250 mm × 4.6 mm, 5 μm, Agilent 1100, Agilent Technologies, Waldbronn, Germany) based on our previous work 17 . HPLC was used at a flow rate of 1.2 ml/min with fluorescence detection after pre-column derivatization (with o-phthaldialdehyde, OPA). The fluorescent signal was 450 nm when it was excited at 348 nm. As the reference (blank), maturation medium without COCs, EGF and urea supplementation were cultured under the same conditions as the test samples. After 24-h incubation, all blanks and test media were immediately kept at -70 °C and analyzed using the same protocol.
For the determination of AA concentrations, the peak areas of the samples and the standard AA mixture were compared. Using peak signals from internal standards (25 mM sarcosine and 25 mM norvaline both in 0.1 M HCl as internal standards), all peak signals were normalized. O-phthaldialdehyde in the presence of the thiol component (OPA method) used in this study has been shown to have a low response to cystine and none to proline 55 . The assay of asparagine was unlikely due to conversion to aspartic acid during the preparation of the sample 56 .
The accuracy of the AA assessment was checked by calculating the mean recovery of the biological medium with a certain AAs volume close to that of the AAs in the maturation medium. Accuracy was defined as the agreement between the measured value and the actual AAs value in the reference media. The average recovery for all AAs observed was 105.8%. The accuracy ranged from 91.1% to 109.1%, excluding aspartic acid (116.4%). Therefore, the limit for all AAs was set at between 90 and 117%. In order to test the precision of the method, the repeatability was calculated by injecting the same samples of fresh maturation medium (four times) in a row and calculating the relative standard deviation (RSD) for each AA. The RSD values obtained were within an acceptable range of 1.15 to 3.59%. The average statistical power of this experiment was adequate at 88.9 ± 10.0%, ranging from 68.0 to 100%.
In order to determine intracellular levels of urea, DOs in groups of 25 were cultured in the culture medium (50 µl droplet) for 4 h with 40 mg/dl urea. The harvested DOs were washed and centrifuged at 300 g for 5 min Scientific RepoRtS | (2020) 10:16082 | https://doi.org/10.1038/s41598-020-73225-7 www.nature.com/scientificreports/ and re-suspended in the culture medium (1 ml). The resultant mixture was mechanically lysed with ultrasonic homogenization for 1 min, and then urea was measured on the Hitachi 917 auto-analyzer (Roche Diagnostics, Mannheim, Germany) using an enzymatic method.

Calculation of amino acid depletion/appearance. Differences between the concentration of AA in
the test medium (after 24 h of incubation) and that of fresh medium were used to calculate the depletion or appearance of each AA 18,19 . The concentration of AAs in fresh maturation medium was based on concentrations present in the tissue culture medium (TCM)-199 and previous studies 59 . In recent studies, the depletion/ appearance of AAs in bovine and human MII oocytes was measured during the final 6 h of IVM with a decreased concentration of AA (one-sixth MEM and five-sixth Earle's Balanced Salt Solution) 18,19 . In previous studies by Hemmings et al., IVM was divided into two different phases, including phase I in which the maturation media had a normal AA concentration and phase II in which the AA concentration was reduced at the final 6 h of IVM 18,19 . They measured the AAs concentration of spent media of the last 6 h of IVM. On the other hand, we did not transfer COCs to a fresh maturation medium with a reduced concentration of AA (e.g., for the last 6 h of IVM) and measured the concentration of AA in a 24-h spent medium. It has been shown that, in the mammalian cell culture model, sudden nutrient reductions, such as AAs, can quickly stimulate autophagy within minutes 60 . In addition, it has been confirmed that MII oocytes have completed the final maturation processes, with reduced synthetic activity and a quiet stage 22 . As a consequence, we used the 24-h spent medium for AA determination to cover all phases of COCs metabolism. In addition, in order to reduce possible adverse reactions, such as autophagy 60 , we did not transfer COCs to the reduced AA concentration media during IVM (e.g., for the last 6 h of IVM).
Data were calculated in pmol/h and adjusted on the basis of the total number of COCs (n = 10) within each droplet, according to previous reports 18,19 . In addition, the concentrations of AAs in each 24-h spent medium were normalized to the number of CCs derived from 10 COCs inseminated in each droplet. It was performed to reduce the effect of variability (in the amount of CCs attached to COCs) on the depletion/appearance of AAs. To estimate the number of CCs per each droplet (consisting 10 COCs), the number of fertilizing sperm cells (2 × 10 6 / ml) was subtracted from the total cell numbers (including sperm and CCs) counted 8 h after IVF. Cumulus cells were completely removed from the oocytes through gentle vortex agitation. Approximately (2.6 ± 0.2) × 10 5 cells, varying from (1.9 ± 0.2) × 10 5 to (2.7 ± 0.8) × 10 5 , were counted per each droplet. The total number of cells (dead and live) was evaluated using Trypan blue staining.
In vitro fertilization and embryo culture. After IVM (using a group culture), COCs from two droplets of the same treatment were pooled in order to achieve a well (50-μl droplet) of 18 to 20 COCs. Therefore, COCs were allocated to 12 wells (50-μl droplet) and to four groups based on their own IVM treatments (Supplementary Fig. 1, Phase 2). Thus, the IVF and IVC media were not treated with any agents (e.g., EGF and urea). Next, COCs matured under different experimental treatments during IVM were inseminated with sperm cells at a final concentration of 2 × 10 6 /ml in the IVF medium (50-μl droplets) and categorized as follows: (1) fertilizing sperm + COCs matured in the absence of urea and EGF (the control group); (2) fertilizing sperm + COCs matured in the presence of 10 ng/ml EGF; (3) fertilizing sperm + COCs matured in the presence of 10 ng/ml EGF and 20 mg/dl urea; and (4) fertilizing sperm + COCs matured in the presence of 10 ng/ml EGF and 40 mg/ dl urea.
In the case of individual-IVM culture, following IVM, COCs were individually transferred to IVF medium (20-μl droplet), inseminated with sperm cells at a final concentration of 2 × 10 6 /ml and categorized as follows: (1) fertilizing sperm + COCs matured without urea and EGF (the control group); (2) fertilizing sperm + COCs matured in the presence of 10 ng/ml EGF; (3) fertilizing sperm + COCs matured in the presence of 40 mg/dl urea; and (4) fertilizing sperm + COCs matured in the presence of 10 ng/ml EGF and 40 mg/dl. The insemination of COCs was performed in the IVF medium, supplemented by BSA in 35 mm cell culture dishes (Falcon brand, BD Biosciences, San Jose, CA, USA) at 38.5 °C under 5% CO 2 in air. According to a previous study 61 , sperm cells were capacitated in a 96-well untreated polystyrene microplate (Corning Incorporated, New York, U.S.A.) using a modified Tyrod's albumin, lactate, and pyruvate medium (Sp-TALP). In brief, after 20 min of swim-down, sperm cells (5 × 10 6 sperm/ml) were suspended in Sp-TALP supplemented with 10 mg/ ml heparin, and incubated for 18 h.
Zygotes were recovered 8 h after insemination and CCs were manually separated by repeated pipetting in HEPES-buffered MEM (BioWhittaker, Walkersville, MD, U.S.A.) containing 80 IU/ml of bovine hyaluronidase using narrow-bore glass pipettes. The resulting embryos were transferred to 35 mm cell culture dishes containing synthetic oviduct fluid medium (50-μl droplet), supplemented by 30 µl/ml of EAAs solution (50 × ; 11130-051; Gibco), 10 µl/ml of non-EAAs solution (100 × ; 11140; Gibco), and 4 mg/ml BSA in an incubator at 38.5 °C, 5% O 2 , 6% CO 2 , and 90% N 2 . The proportion of COCs which did not reach the 2-PN stage was recorded using a stereomicroscope (Olympus, Tokyo, Japan) at 18 h post-insemination. The percentage of COCs reached the cleavage stage (2-16 cell stage, day 3 post-insemination) or blastocyst stage (day 7 post-insemination) was reported. All developmental parameters, including 2PN, cleavage, and blastocyst percentage, were calculated in relation to the initial number of COCs cultured in the maturation medium. Eventually, a total of nine pools were derived from 12 IVC droplets per treatment ( Supplementary Fig. 1, Phase 3). Each pool consisted of 7 to 8 blastocysts, which produced from COCs that matured under different treatments (i.e., urea and EGF) and kept at − 70 until qRT-PCR. These embryos had been cultured in IVC media based on initial experimental treatment used for COCs and had not been treated with any agents (e.g., EGF and urea). www.nature.com/scientificreports/ RNA extraction, reverse transcription, quantitative analysis of transcripts using qRT-PCR. RNA extraction and reverse transcription were conducted as reported earlier 16 . In short, total RNA was extracted from day-7 blastocysts using RNeasy Micro kit (Qiagen, Mississauga, Ontario, Canada) and then treated with DNase I (Ambion, Streetsville, Ontario, Canada) according to the manufacturer's manual. To determine the quality and quantity of RNA, the WPA Biowave spectrophotometer (Cambridge, United Kingdom) was employed. Reverse transcription was done for 10 min at 25 °C, for 1 h at 42 °C, and for 10 min at 70 °C. Using 1 µl of cDNA (50 ng), 5 µl of the SYBR Green/0.2 µl of ROX qPCR Master Mix (2X, Fermentas, Germany) and 1 µl of forward and reverse primers (5 pM) adjusted to a total volume of 10 ml using nuclease-free water, qRT-PCR was carried out. The reference gene was ACTB. Each replicate of qRT-PCR was repeated three times in order to minimize the technical errors. The primer sequences that were used for qRT-PCR were: 5′-CCT TCT  TTG AGT TCG GAG -3′,  Despite normal morphology, the condition of IVM has been shown to have a negative effect on the level of gene expression in the human preimplanted embryos 62 . Therefore, we assumed that different IVM conditions in this study could affect the gene expression in the resulting blastocysts. Warzych et al. reported that supplementation of the maturation medium with a macromolecule (PVP40) was associated with a higher apoptotic index in bovine blastocysts 63 . Moreover, because urea as a toxin was used during IVM in this study, we decided to choose BCL2 and BAX transcripts to determine whether the IVM environment could affect the major apoptosis regulators in the resulting blastocysts. In addition, EGFR inhibition increases the rate of apoptosis and changes the expression of specific BCL2 family members in favor of apoptosis 64 . EGFR has been shown to be clearly associated with the transcription factors OCT4 and NANOG, which are essential for determining the pluripotent status of embryonic stem cells 65 . EGFR inhibition has been reported to reduce the expression of OCT4 and NANOG, while the activation of EGFR acutely increases the activity of DNA methyltransferase (DNMT) 65,66 . Hence, these transcription factors were selected in this study to be evaluated in the resulting blastocysts.
Statistics. The Anderson-Darling test (EasyFit software, version 5.6, MathWave Technologies, Spokane, WA, USA) showed that the data from the group-IVM culture was normally distributed in each experimental treatment, but the data was non-normally distributed regardless of the experimental treatments. In the case of individual-IVM culture, data was normally distributed in each experimental treatment or regardless of experimental treatments. The medians of depletion/appearance of AAs and developmental competence of COCs were analyzed using the Kruskal-Wallis test followed by the Dunn's multiple comparison test. In the case of gene expression and use of NAC, the one-way ANOVA followed by Tukey's multiple comparisons test was used to evaluate the results. Differences at p < 0.05 were considered statistically significant. Orthogonal polynomial contrasts were used to determine the linear or quadratic effect of urea using SAS software (SAS Institute Inc., Cary, NC). For ANOVA analysis study, we also calculated q-value vs. p-value to address multiple testing problems. Multiple testing problems usually exist in metabolomic-based studies 67 . Indeed, a p-value of 0.05 is a false positive rate when one test is done, but with a large number of tests on the data, one needs to compensate with false positives. There are various approaches to correct for multiple testing. The use of false discovery rates (i.e., q-values) is most common. The q-value represents the false discovery rate; the lower q-value (q < 0.10) implies high confidence in the test 67 . In this analysis, we observed that the false discovery rate of various comparisons with significant p values was small (q < 0.08, Supplementary Fig. 6 and Supplementary Fig. 7). According to p-value vs .q-value plot, in these results, we found that q-value 0.05 corresponded p-value of about 0.03 ( Supplementary  Fig. 7a). For this analysis, threshold comparisons of q-values of less than 0.05 yield 113 significant differences between comparisons (Supplementary Fig. 7b,c). It indicates that only 6 of the 113 comparisons listed significant are likely to be false positives.
The bivariate analysis was performed using Kendall's tau or Pearson correlation (for the data derived from the group-and individual-IVM systems, respectively) to determine the associations between each parameter relevant to the developmental competence and AAs depletion/appearance. Next, multivariate methods, such as the Hierarchical Clustering Analysis (HCA) and the Principal Component Analysis (PCA), were used to consider the effects of inter-correlated variables (i.e., AAs) simultaneously on IVF outcomes. Bonferroni correction was used to test the significance of the correlations in the correlation matrix or in the t-test analysis.
Because, the data derived from the group-IVM culture were not normally distributed, these non-parametric data were transformed using the logarithm function. These log-transformed data were then applied to HCA and PCA and analyzed using the Kendall distance metric. However, Pearson metric was used in the case of normal data derived from the individual-IVM culture.
Using the "Heatmapper" web tool (https ://www2.heatm apper .ca/expre ssion /), an HCA analysis with Kendall or Pearson distance metrics (for data derived from group-and individual-IVM systems, respectively) and a complete linkage method was used to produce clustering patterns of AAs and developmental competence. The PCA analysis was used to consider the effects of all AAs on oocyte competence at the same time, irrespective of experimental treatments, to reduce the multiple data dimensions to two dimensions and to produce biplots.
Network analysis was performed to identify central nodes (i.e., AAs). Rho and Pearson similarity indexes were used for data derived from group-and individual-IVM systems, respectively. The Fruchterman-Reingold Scientific RepoRtS | (2020) 10:16082 | https://doi.org/10.1038/s41598-020-73225-7 www.nature.com/scientificreports/ algorithm was used as a force-directed layout algorithm. To this aim, the data was mapped to the amino aciddevelopmental competence interaction network. For the Kendall rank correlation, ANOVA, PCA, network mapping, and log transformation, the PAST program (accessible at: https ://folk.uio.no/ohamm er/past) was used.
Using the easyROC web-tool (https ://www.bioso ft.hacet tepe.edu.tr/easyR OC/), a Receiver Operating Characteristic (ROC) curve analysis was performed to identify the best predictive factors for the detection of COCs with high developmental competence. The optimal cutoff for high IVF outcome prediction was determined by maximizing the Youden index. In order to deal with potential biases in the results produced by assuming that any relationship between AAs and high oocyte competence could be linear, AAs were divided into two separate groups based on IVF outcomes derived from the group-IVM system , i.e., high percentage of 2PN (> 80%) vs. absence of high 2PN percentage (< 80%), high percentage of cleavage (> 80%) vs. absence of high percentage of cleavage (< 80%), high percentage of blastocyst formation (> 30%) vs. absence of high percentage of blastocyst formation (< 80%), and low percentage of oocyte degeneration (< 15%) vs. absence of low percentage of oocyte degeneration (> 15%). Also, AAs were divided into two separate groups based on the incidence of cleavage (incidence of cleavage vs. absence of cleavage) in the case of an individual-IVM system. Finally, the area under the curve (AUC) was determined by the Mann-Whitney statistic (a non-parametric unbiased likelihood estimator) to detect COCs with a high percentage of 2PN, cleavage and blastocyst formation, as well as low percentage of oocyte degeneration. In the case of an individual-IVM culture with normal data, the AUC was determined by the Asymptotic statistic (a parametric unbiased likelihood estimator) for the detection of the competent COCs.
In silico molecular docking. The initial structure of urea was constructed by GaussView 5.0, which is a graphical software widely used with Gaussian 68 . The optimized structure of urea was obtained using GAMESS-US package 69 at the B3LYP/6-31 + G (d) level through density functional theory (DFT). Molecular docking to predict urea interactions with DNA, EGF, EGFR, and AQP3. Molecular docking procedure is widely recognized as one of the important useful tools for computational biology and chemistry 70,71 . This approach can be used to study the binding free energy and ligand binding activity to receptors 72 . Using AutoDock 4.2 tools, the optimized urea was docked with the crystal structure of DNA, EGF, EGFR and AQP3 to predict the binding energy of urea with these molecules 73 . Urea (without any rotatable bond) and other molecules were respectively regarded as ligand and receptors. Molecular dockings were performed using a Lamarkian Genetic Algorithm to predict binding free energy and the direction of urea in its interaction with DNA, EGF, EGFR and AQP3. The grid sizes of "82 × 84 × 126", "82 × 84 × 126", "126 × 126 × 126", and "126 × 40 × 40" were used for the urea-DNA, urea-EGF, urea-EGFR, and urea-AQP3 complexes, respectively. Also, grid spacing of 0.39, 0.45, 1, and 0.375 Å was used for the urea-DNA, urea-EGF, urea-EGFR, and urea-AQP3 complexes, respectively. The crystal form of the extracellular domain of EGFR (PDB ID: 3NJP), EGF (1jl9), and DNA (1 BNA) was taken from the Protein Data Bank at a resolution of 3.30, 3.90, 1.0 Å, respectively. Computer modeling on AQP3 (ID: Q08DE6) was done by using AA sequences from the Universal Protein Resource (UniProt) database (https ://www.unipr ot.org), accompanied by 3-D structure predictions in the I-TASSER repository (https ://zhang lab.ccmb.med.umich .edu/I-TASSE R). The blind docking was completed to discover the suitable binding sites for the urea molecule's interaction with the receptors. The blind docking scans the entire receptor surface and calculates the positions of the highest binding affinity. The focus docking was performed after the blind docking had been completed. The total number of docking runs for all urea-DNA, urea-EGF, urea-EGFR, and urea-AQP3 docking procedures was set at 100.
In order to study the urea binding behavior to different receptors, the docking procedure was performed for the urea-DNA, urea-EGF, urea-EGFR, and urea-AQP3 complexes at 20, 12, 12, and 12 sequential times, respectively. To this aim, the previous urea-receptor complex was used as a receptor for the next molecular docking (under the same conditions). For example, the urea-EGFR complex with n urea molecules was used as the initial template for obtaining the urea/EGFR complex with n + 1 urea molecules. It should be noted that when the next urea molecule was added to the complex, the positions of the urea molecules in the urea/receptor complex obtained from the preceding docking were rigid.
Next, to investigate how EGF's binding energy with EGFR changes in the presence of urea, the molecular docking was performed as follows with four phases: phase I: EGF/EGFR; phase II: (EGF-urea complex)/EGFR; phase III: EGF/(EGFR-urea complex); and phase IV: (urea-EGF complex)/(EGFR-urea complex). The EGFurea and EGFR-urea complexes obtained from the preceding step were used as receptor in these phases. Using Lamarkian Genetic Algorithm, molecular dockings were performed to predict binding free energy of these complexes. The grid sizes of "126 × 104 × 94", "76 × 126 × 92", "105 × 126 × 82" and "82 × 96 × 68" were used for phases I, II, III and IV, respectively. Also, grid spacing of 0.669, 0.686, 0.664 and 1 Å was used for phases I, II, III and IV, respectively. LIGPLOT software was used to identify van der Waals interactions between urea and residues 74 . LIGPLOT is a bioinformatics computer program that produces 2-D ligand-receptor interaction representation 74  www.nature.com/scientificreports/