microRNA-based signatures obtained from endometrial fluid identify implantative endometrium.

STUDY QUESTION: Is it possible to use free and extracellular vesicle-associated microRNAs (miRNAs) from human endometrial fluid (EF) samples as non-invasive biomarkers for implantative endometrium? SUMMARY ANSWER: The free and extracellular vesicle-associated miRNAs can be used to detect implantative endometrium in a non-invasive manner. WHAT IS KNOWN ALREADY: miRNAs and extracellular vesicles (EVs) from EF have been described as mediators of the embryo-endometrium crosstalk. Therefore, the analysis of miRNA from this fluid could become a non-invasive technique for recognizing implantative endometrium. This analysis could potentially help improve the implantation rates in ART. STUDY DESIGN, SIZE, DURATION: In this prospective study, we first optimized different protocols for EVs and miRNA analyses using the EF of a setup cohort (n = 72). Then, we examined differentially expressed miRNAs in the EF of women with successful embryo implantation (discovery cohort n = 15/validation cohort n = 30) in comparison with those for whom the implantation had failed (discovery cohort n = 15/validation cohort n = 30). Successful embryo implantation was considered when pregnancy was confirmed by vaginal ultrasound showing a gestational sac 4 weeks after embryo transfer (ET). PARTICIPANTS/MATERIALS, SETTING,
METHODS: The EF of the setup cohort was obtained before starting fertility treatment during the natural cycle, 16-21 days after the beginning of menstruation. For the discovery and validation cohorts, the EF was collected from women undergoing frozen ET on Day 5, and the samples were collected immediately before ET. In this study, we compared five different methods; two of them based on direct extraction of RNA and the other three with an EV enrichment step before the RNA extraction. Small RNA sequencing was performed to determine the most efficient method and find a predictive model differentiating between implantative and non-implantative endometrium. The models were confirmed using quantitative PCR in two sets of samples (discovery and validation cohorts) with different implantation outcomes. MAIN RESULTS AND THE ROLE OF CHANCE: The protocols using EV enrichment detected more miRNAs than the methods based on direct RNA extraction. The two most efficient protocols (using polymer-based precipitation (PBP): PBP-M and PBP-N) were used to obtain two predictive models (based on three miRNAs) allowing us to distinguish between an implantative and non-implantative endometrium. The first Model 1 (PBP-M) (discovery: AUC = 0.93; P-value = 0.003; validation: AUC = 0.69; P-value = 0.019) used hsa-miR-200b-3p, hsa-miR-24-3p and hsa-miR-148b-3p. Model 2 (PBP-N) (discovery: AUC = 0.92; P-value = 0.0002; validation: AUC = 0.78; P-value = 0.0002) used hsa-miR-200b-3p, hsa-miR-24-3p and hsa-miR-99b-5p. Functional analysis of these miRNAs showed strong association with key implantation processes such as in utero embryonic development or transforming growth factor-beta signaling. LARGE SCALE DATA: The FASTQ data are available in the GEO database (access number GSE178917). LIMITATIONS, REASONS FOR CAUTION: One important factor to consider is the inherent variability among the women involved in the trial and among the transferred embryos. The embryos were pre-selected based on morphology, but neither genetic nor molecular studies were conducted, which would have improved the accuracy of our tests. In addition, a limitation in miRNA library construction is the low amount of input RNA. WIDER IMPLICATIONS OF THE
FINDINGS: We describe new non-invasive protocols to analyze miRNAs from small volumes of EF. These protocols could be implemented in clinical practice to assess the status of the endometrium before attempting ET. Such evaluation could help to avoid the loss of embryos transferred to a non-implantative endometrium. STUDY FUNDING/COMPETING INTEREST(S): J.I.-P. was supported by a predoctoral grant from the Basque Government (PRE_2017_0204). This study was partially funded by the Grant for Fertility Innovation (GFI, 2011) from Merck (Darmstadt, Germany). It was also supported by the Spanish Ministry of Economy and Competitiveness MINECO within the National Plan RTI2018-094969-B-I00, the European Union's Horizon 2020 research and innovation program (860303), the Severo Ochoa Centre of Excellence Innovative Research Grant (SEV-2016-0644) and the Instituto de Salud Carlos III (PI20/01131). The funding entities did not play any role in the study design, collection, analysis and interpretation of data, writing of the report or the decision to submit the article for publication. The authors declare no competing interests.

Introduction

Increasing embryo implantation rates is one of the greatest challenges in ART, as only 35% of embryo transfers (ETs) result in a clinical pregnancy (Matorras ; De Geyter ). Despite numerous studies focused on improving implantation rates, a reliable method of determining the competence of the endometrium, fundamental for successful implantation, is still lacking (Strowitzki ; Craciunas ). Currently, the endometrial biopsy is used to establish whether the endometrium is ready for ET (Casper, 2020). This is an invasive methodology, and the ET is not performed in the same cycle in which the sample is taken as it can have detrimental effects on implantation (van der Gaast ). If the biopsy shows that the endometrium is receptive, the results will be extrapolated to the next cycle. This assumption is not realistic, since the endometrial cycle is a dynamic process involving many factors affecting the receptivity of the endometrium. The analysis of endometrial fluid (EF) obtained in a non-invasive manner, without biopsy, is a promising alternative (van der Gaast ). It has been demonstrated that the aspiration of EF immediately before the ET does not affect the implantation. Moreover, the prompt analysis of EF composition might allow the ET in the same cycle (van der Gaast ; Azkargorta ; Matorras , 2020). The EF can be obtained several times during the cycle and its analysis could reveal whether the endometrium is ready for implantation or therapeutic intervention is necessary for a successful procedure.

The EF is a complex biological fluid that can modulate endometrial homeostasis and receptivity, it can sustain the preimplantation embryo and initiate the implantation process and it plays an important role in the embryo–endometrium communication (Ng ; Vilella ; Bhusane ; Nguyen ). microRNAs (miRNAs) are small non-coding RNA sequences (18–22 nucleotides) that are important regulators of genes at the post-transcriptional level (Bhaskaran and Mohan, 2014). They are essential during early embryonic development since they regulate cell proliferation and differentiation (Bhaskaran and Mohan, 2014). Some of these miRNAs have been associated with the extracellular vesicles (EVs), also present in the fluid obtained from the uterine cavity (Vilella ). EVs are widely known mediators of intercellular communication, transmitting information from one cell to a multitude of other cells and locations (Han ). Moreover, analyses of miRNA content of endometrium-derived EVs show that they are taken up by the embryos, modifying their transcriptomic and adhesive phenotypes (Ng ; Vilella ; Greening ; Balaguer ; Marinaro ). For example, the EV-associated hsa-miR-30d is internalized by mouse trophoectoderm and increases the embryo adhesion via upregulation of adhesive molecules (Vilella ).

One of the main challenges in ART is finding non-invasive tools for detecting the best time to perform the ET. Here, we developed a reproducible, sensitive, low-invasive method to comprehensively examine the miRNA landscape of the EF. First, we optimized the EF sample collection technique. Then, we established a robust method for analyzing vesicular and non-vesicular miRNAs from EF obtained in clinical settings, where sample size is limited and no sophisticated equipment is available. Finally, we applied these methods to a set of EF samples from women with different implantation outcomes. Our aim was to define a miRNA signature to identify the competence of the endometrium. If we could determine the state of the endometrium, it would then be possible to change the ET strategy when the results show an unfavorable implantative pattern. Thus, the implantation rates could potentially be improved and the loss of embryos minimized by avoiding their transfer to non-implantative endometrium.

Materials and methods

Ethical approval

Ethical approval for the study was obtained from the Cruces University Hospital Ethics Committee and Institutional Review Board (CEIC 11/45) and all the participants gave written consent for their participation.

Study population

The population under study consisted of a cohort of 162 women who attended the Human Reproduction Unit of Cruces University Hospital (Basque Country, Spain) from January 2018 to February 2021. For the setup and optimization of the techniques, the samples were collected before starting the fertility treatment. The samples were collected during the natural cycle, 16–21 days after the beginning of menstruation. To test the selected method, the samples were collected just before Day-5 frozen ETs, a practice which is performed increasingly often (Matorras ). Out of 162 women (Supplementary Fig. S1), 72 participated in the setup, 30 in the discovery of the predicted models and 60 in the validation of the models. Forty-five women became pregnant and were included in the implantative endometrium group. The other 45, who did not achieve pregnancy, were included in the non-implantative endometrium group. The endometrium was considered implantative when pregnancy was confirmed by vaginal ultrasound showing a gestational sac 4 weeks after ET. Cases with a positive β-hCG test where a gestational sac was not seen on vaginal ultrasound (biochemical miscarriages) were not included in the study.

The inclusion criteria in the setup study were: age between 18 and 37 years; cycle duration between 27 and 29 days; absence of ovulatory disorders, myomas, endometriosis, polyps, uterine scars or hydrosalpinges; normal uterine and ovarian ultrasound; serum anti-Müllerian hormone > 0.4 ng/ml; and no history of gynecological infections, immune disorders or gynecological surgery. The inclusion criteria for the discovery and validation cohorts also included: frozen ET on Day 5 (good quality embryos; Types A and B of the Spanish Society for the Study of Reproductive Biology (ASEBIR) classification (ASEBIR, 2015) and transfer of 1–2 embryos derived from the oocytes of the same subject.

The management of endometrial preparation was always carried out using the same protocol. A vaginal ultrasound was performed on Day 1 or 2 to confirm ovarian quiescence (absence of follicles >10 mm). An artificial cycle was started on Day 2 by administering 6 mg of estradiol daily (Progynova, Bayer, Barcelona, Spain). The development of the endometrium was monitored using serial vaginal ultrasounds. When the endometrium became 7-mm thick, the transfer day was scheduled. Vaginal progesterone at a dose of 400 mg/12 hr (Utrogestan, SEID, Barcelona, Spain) was started the next morning, and the ET was performed on the 5th day of progesterone administration. If pregnancy was achieved, the estradiol and progesterone treatment was maintained until the 12th week of gestation.

Embryo vitrification was performed on Day 4 or 5 using a Cryotop device (Kitazato BioPharma Co., Shizuoka, Japan). The embryos were cryopreserved and warmed using the Kitazato vitrification/warming kit (Kitazato BioPharma Co.), according to the manufacturer's instructions. Frozen Day-4 embryos were thawed and cultured for 24 hr before the ET and Day-5 blastocysts for 2 hr before the ET.

Sample collection and storage

The EF was aspirated with a catheter used for ET (Frydman, Instrumentos Médicos Estériles SA, Spain) connected to a 10-ml syringe under abdominal ultrasound guidance. Sample extraction was performed by gently applying a negative pressure with the syringe. The aspiration was interrupted at the internal cervical os to prevent contamination with cervical mucus. Special care was taken to avoid touching the uterine fundus or injuring the cervix and minimize sample contamination with blood and endometrial tissue. In cases with excessive vaginal secretions, the vagina was cleaned with saline solution before aspiration. Aspirate volumes ranged from 5 to 50 µl. After aspiration, the 10-ml syringe was replaced with a 2-ml syringe containing 1.5 ml of 1× Dulbecco’s PBS (DPBS) (Gibco, Thermo Fisher Scientific, # 14190250, MA, USA) to expel the EF. The aspirates were mixed with the 1× DPBS and expelled into a cryogenic tube (5–50 µl of EF + 1500 µl of 1× DPBS). The mixed samples were centrifuged to remove contaminants at 2500g for 5 min at room temperature, and the supernatants were then kept frozen at −80°C until processed. The dilution of the supernatants was 1:30, with a final volume between 400 and 1300 µl.

EV enrichment methods

Size-exclusion chromatography

A Poly-Prep chromatography column (BioRad, # 731-1550, Hercules, USA) was filled with 2.5 ml of Sepharose CL-2B cross-linked resin (Sigma, # CL2B300-100ML) and left packing overnight at 4°C. The column was then washed twice with 2.5 ml of 1× DPBS. Three aliquots of 400 µl from the setup cohort sample pool were used. Each aliquot was applied to the column, and then 4 ml of 1× DPBS was added. The size-exclusion chromatography (SEC) separated the sample into 12 fractions (F1–F12); the EVs were eluted mainly in F3 but also in F4 and F5 fractions, as described by Prieto-Fernández . F1 to F10 had a final volume of 200 µl, and F11 and F12 of 1 ml. The 12 fractions of one aliquot were each used for RNA extraction with mirVana™ PARIS™ Kit (Thermo Fisher Scientific, # AM1556). The RNA obtained from fractions F3 and F4 was further analyzed by small RNA-sequencing (RNA-Seq). The 12 fractions of the other two aliquots were characterized using western blot (WB).

Polymer-based precipitation method

Since there was no published protocol for using the Invitrogen Total Exosome Isolation Reagent with the EF, we compared the Total Exosome Isolation Reagent for the cell culture media (Invitrogen by Thermo Fisher Scientific, # 4478359) with Total Exosome Isolation Reagent for other body fluids (Invitrogen by Thermo Fisher Scientific, # 4484453). Although both worked well with the EF, we used the # 4478359 because of its better cost-effectiveness ratio. The optimized protocol was as follows: centrifuge the EF supernatants at 3000g for 30 min at 4°C; transfer the supernatants to a fresh tube and add an equal volume of the Total Exosome Isolation Reagent (1:1); stir the mixture by vortexing until there is a homogeneous solution and incubate the sample for 30 min at room temperature; after the incubation, centrifuge the samples at 10 000g for 1 hr at 4°C; aspirate the supernatant by pipetting and discard it; the EVs are contained in the pellet, which may not be visible at the bottom of the tube; and finally, add 100 µl of 1× DPBS to resuspend the pellet.

Ultracentrifugation

Ultracentrifugation (UC) was carried out in a single step (100 000g for 75 min at 4°C) using a Beckman-Coulter TLA 120.2 rotor. The EV pellets were resuspended in 100 µl of 1× DPBS.

RNA extraction methods

We used two RNA isolation methods; we followed the manufacturer's instructions for the mirVana™ PARIS™ Kit (Thermo Fisher Scientific, # AM1556) (DCT-M) and Norgen Plasma/Serum RNA Purification kit (DCT-N). Two different Norgen kits were used as needed; the midi kit (Norgen Biotek Corp., # 56100, Ontario, Canada) or the mini kit (Norgen Biotek Corp., # 55000). The RNA was eluted in nuclease-free water (Ambion, # AM9930 by Thermo Fisher Scientific).

cDNA synthesis and TaqMan miRNA assay

Following the manufacturer's recommendations, cDNA was synthesized from 2 µl of RNA using the TaqMan Advanced miRNA cDNA Synthesis kit (Applied Biosystems, # A28007, by Thermo Fisher Scientific). The TaqMan reactions used were the TaqMan Fast Advanced Master Mix (Thermo Fisher Scientific, # 4444557) and TaqMan Advance miRNA assays (Thermo Fisher Scientific, # A25576). The quantitative PCR was performed using a Viia7 or QS6 system, and the data were analyzed using the QuantStudio Real-Time PCR System version 1.3 (Applied Biosystems, by Thermo Fisher Scientific). The expression profiles of seven EV-associated miRNAs were used as reference (Thermo Fisher Scientific): hsa-let-7-5p (478579_mir), hsa-miR-17-5p (478447_mir), hsa-miR-200c-3p (478351_mir), hsa-miR-30c-5p (478008_mir), hsa-miR-30d-5p (478606_mir), hsa-miR-451a (478107_mir) and hsa-miR-92a-3p (477827_mir) (Supplementary Table SI). These miRNAs have been reported as secreted by endometrial epithelial cell lines (Ng ), found in the EF aspirates (Vilella ; Campoy ) and secreted in endometrial exosomes associated with early embryo implantation (Vilella ; Balaguer ). Two other miRNAs were selected after the small RNA-Seq analysis using the setup pool cohort sample. These were the hsa-miR-21-5p (477975_mir) and hsa-miR-155-5p (483064_mir) miRNAs, which were among the most and least abundant miRNAs in the pool, respectively (Supplementary Table SI). In addition, two exogenous miRNAs were used as internal controls. To examine the efficiency of the RNA extraction, 4 µl of cel-miR-39 (478293_mir, Thermo Fisher Scientific) of a 0.1 nM stock were added to the sample before each RNA extraction procedure. To test the differences between the cDNA synthesis reactions, 0.2 µl of ath-miR-159a (478411_mir, Thermo Fisher Scientific) of a 0.001 nM stock was added at the beginning of each cDNA synthesis reaction.

Comparing the miRNA extraction methods

Five different methods were compared to define a simple and effective strategy for detecting vesicular and non-vesicular miRNAs in small volumes of EF (Fig. 1). Two of these involved direct extraction using different RNA extraction kits, DCT-N (Norgen kit, # 56100) and DCT-M (mirVana PARIS kit, # AM1556). The other three required enrichment of EVs before RNA extraction. In one case, the enrichment was carried out by UC followed by RNA extraction with mirVana PARIS kit (UC-M). In the remaining two cases, the enrichment was carried out using the polymer-based precipitation (PBP) method and the RNA was extracted using the Norgen (# 55000) (PBP-N) or mirVana PARIS kit (PBP-M). In parallel, SEC was performed to characterize the protein and miRNA content of the EF (Fig. 1). The volumes of recovered EF samples varied depending on many factors, such as the operator collecting the sample and the EF volume or viscosity. In general, the volumes ranged from 400 µl to 1.3 ml. Therefore, we optimized the protocols to be used with the minimum volume available (400 µl) in all cases. All the tests were performed in triplicate.

Figure 1.

Workflow summarizing the different methods used to analyze microRNAs from the endometrial fluid of patients undergoing ART. We compared five different methods, two of which used the direct extraction of RNA from the endometrial fluid (EF) (DCT-N and DCT-M). The other three included the extracellular vesicle (EV) enrichment (UC-M, PBP-N and PBP-M) before RNA extraction. In parallel, we carried out a size-exclusion chromatography (SEC-M) to characterize the proteins and miRNAs in the EF. The samples came from the setup pool cohort, and each experiment was performed in triplicate, using sample aliquots of 400 µl. DCT-N: direct RNA extraction with Norgen Plasma/Serum RNA purification kit. DCT-M: direct extraction of RNA with mirVana PARIS kit. UC-M: EV enrichment by ultracentrifugation and RNA extraction using mirVana PARIS kit. PBP-N: EV enrichment with a polymer-based precipitation method and RNA extraction with Norgen Plasma/Serum RNA purification kit. PBP-M: EV enrichment using the polymer-based precipitation method and RNA extraction with mirVana PARIS kit. SEC-M: EV enrichment with SEC and RNA extraction with mirVana PARIS kit. miRNAs, microRNAs; PBP, polymer-based precipitation.

Technical reproducibility experiment

A technical reproducibility experiment was conducted using the PBP-M and PBP-N protocols. Two operators (J.I.-P. and M.C.-G.) performed the tests independently. The samples used in these experiments came from the setup pool cohort, and each of the operators tested 10 aliquots using each method. Quantitative PCR (qPCR) was used to examine reproducibility; nine miRNAs were analyzed (seven reference miRNAs and two miRNAs obtained from the small RNA-Seq analysis) (Supplementary Table SI).

Dithiothreitol treatment assay

Two aliquots from the setup pool cohort were used to perform the experiment. One of the aliquots was treated with a 1.4% dithiothreitol (DTT) solution in a 1:1 ratio (Miller ; Wang ) and the other served as control; 1× DPBS (1:1) was added. The samples were mixed by vortexing and incubated at room temperature for 15 min. Then, 1× DPBS was added until the EF samples were diluted to the ratio of 1:8, and the samples were centrifuged at 3000g for 15 min at 4°C. The supernatants were recovered, and 400-µl aliquots were taken. The EV enrichment was conducted using 400-µl aliquots, following the PBP method, and the pellet was resuspended in 100 µl of 1× DPBS. From this volume, 15 µl was reserved for WB, 5 µl for cryo-electron microscopy, 5 µl for nanoparticle-tracking analysis (NTA), and the rest of the suspension was used for RNA analysis with Norgen (# 55000). The isolated RNA was eluted in 100 µl of nuclease-free water. Two microliters of the eluate was used for the subsequent cDNA synthesis, and the rest was stored at −80°C.

RNase protection assay

The samples used in this step came from the setup pool cohort, and each aliquot tested had a final volume of 400 µl. All the samples were first EV-enriched using the PBP method (described above) and the EVs were resuspended in 200 µl of DPBS. In the RNase protection assay, four different procedures were compared. Samples were treated according to following protocols: RNase A (Sigma-Aldrich, # 10109142001, MA, USA) (RNase); Proteinase K (Sigma-Aldrich, # 03115879001) + RNase (PRT-K); Triton X-100 (Sigma-Aldrich, # T8787) + RNase (TX-100); and TX-100 + Proteinase K + RNase (TX + PRT). An untreated sample was used as a control. Samples were treated with TX-100 to a final concentration of 0.1%. Proteinase K (0.05 mg/ml concentration) was added and the mixture was incubated for 10 min at 37°C. The reaction was stopped by adding 5 mM of phenylmethylsulfonyl fluoride (Sigma-Aldrich, # 10837091001) and heating at 90°C for 5 min. The samples were finally treated with 0.1 mg/ml RNase A (RNase) for 20 min at 37°C. The control samples were kept at 4°C until RNA extraction. Before extraction, β-mercaptoethanol was used to inhibit RNases, as described by Norgen (# 55000). The RNA was eluted in 50 µl of nuclease-free water. Two microliters were used for the subsequent cDNA synthesis, and the rest was stored at −80°C. All the analyses were performed in triplicate with two technical duplicates, ending with six TaqMan qPCR replicates.

Analysis and quantification of the EF protein content

WB analysis

A sample of 15 µl was mixed with 5 µl of NuPAGE LDS Sample Buffer 4× (Invitrogen # NP0007, by Thermo Fisher Scientific). The fractions obtained by SEC were concentrated using 99.5% acetone (Panreac Applichem, # 161007, Darmstadt, Germany) and resuspended in 20 µl of 1× LDS sample buffer. They were heated for 5 min at 37°C, 10 min at 65°C and 15 min at 95°C and centrifuged for 10 min at 13 000g. Each protein preparation was loaded and separated under non-reducing conditions in 4–12% Bis–Tris precast gels (Invitrogen, # NP0336BOX, by Thermo Fisher Scientific) in MOPS SDS Running Buffer 1× (Invitrogen, # NP0001, by Thermo Fisher Scientific). Precision Plus Protein Dual Color Standard (BioRad, # 161-0374) was used as a marker for protein molecular weights. The proteins were transferred to an Immobilon-P Transfer membrane (Merck Millipore, # IPVH00010, MA, USA) in NuPAGE Transfer Buffer 1× (Invitrogen, # NP0006-1 by Thermo Fisher Scientific) for 1 hr at 100 V. The blocking was performed using 5% Blotting-Grade Blocker (BioRad, # 170-6404) and 0.2% Tween-20 (Sigma-Aldrich, # P2287, MA, USA) diluted in 1× DPBS, for 1 hr. Primary antibodies were incubated overnight and the membranes were washed three times for 10 min with 1× DPBS. Incubation with the secondary horse-radish peroxidase-conjugated antibody (1:6000) was performed at room temperature for 30 min. The chemiluminescence was detected using Pierce ECL Plus Western Blotting Substrate (Thermo Fisher Scientific, # 32132). The bands were visualized on high-performance films (GE Healthcare, # 28906844, IL, USA) employing the AGFA Curix-60 automatic processor (Agfa, Cologne, Germany). The primary antibodies used in this study were mouse anti-CD63 (1:500; clone H5C6 from Developmental Studies Hybridoma Bank, IA, USA), mouse anti-CD9 (1:500; clone 209306, R&D Systems, Minneapolis, MN, USA), mouse anti-CD81 (1:500, Clone JS-81, 555675, BD, NJ, USA), mouse anti-CD133 (1:500 clone W6B3C1, Miltenyi Biotec, North Rhine-Westphalia, Germany), mouse anti-Rab8 (1:1000; Clone 4, 610844, BD, NJ, USA), mouse anti-Flotillin-1 (1:500; Clone 18 610820, BD, NJ, USA), mouse anti-HSP90 (1:500; 610418, BD, NJ, USA) and rabbit anti-Limp II (1:500; ab16522, Abcam, Cambridge, UK). The intensity of the bands was quantified by densitometry using ImageJ software v. 1.52a (ImageJ software, MD, USA).

Coomassie blue staining

SimplyBlue™ SafeStain from Invitrogen (Cat. # LC6060, Thermo Fisher Scientific) was used following the manufacturer's recommendations. The intensity of the bands was quantified by densitometry using ImageJ software (v. 1.52a).

Spectrophotometer

Spectrophotometric measurements were performed using a NanoDrop™ One Microvolume UV–Vis Spectrophotometer (Thermo Fisher Scientific) in the wavelength range of 230–576 nm. Concentrations of RNA and proteins were obtained after measuring the absorbance of 1 µl of the sample.

Nanoparticle-tracking analysis

The size distribution of the EV preparations was analyzed by measuring the rate of Brownian motion using a NanoSight LM10 system (NanoSight, Amesbury, UK), equipped with fast video capture and particle-tracking software. NTA acquisition settings were the same for all samples, and each video was analyzed to obtain the mean and mode of vesicle size and estimate the particle concentration (Dragovic ).

Cryo-electron microscopy

EV preparations were directly adsorbed onto glow-discharged holey carbon grids (Quantifoil, Großlöbichau, Germany). The grids were blotted at 95% humidity and rapidly plunged into liquid ethane with the aid of Vitrobot (Maastricht Instruments BV, Maastricht, The Netherlands). Vitrified samples were imaged at liquid-nitrogen temperature using a JEM-2200FS/CR transmission cryo-electron microscope (JEOL, Tokyio, Japan) equipped with a field emission gun and operated at an acceleration voltage of 200 kV.

Real-time qPCR assay

The relative expression levels of the miRNAs obtained for the setup pool cohort were normalized to ath-miR-159 expression and calculated using the 2−ΔCt (Ct miRNA−Ct ath-miR-159a) method. The relative expression levels of the discovered and validated miRNAs were normalized to internal controls; the differences between the groups were calculated employing the 2−ΔCt (Ct miRNA−Ct mean internal controls) equation. Subsequently, the fold changes were obtained using the 2−ΔΔCt method (Rao ). Endogenous controls were selected from the reference miRNAs (Supplementary Table SI) using the NormFinder software (MOMA, Aarhus, Denmark). The NormFinder is an algorithm using a model-based approach to calculate the stability of a reference transcript; the calculation is based on the intergroup and intragroup variations. The stability score is a weighted measure of these two parameters, and the most stable reference transcript is the one with the smallest stability value (Andersen ).

Only the samples for which we could find the internal controls with fewer than 30 Ct cycles (in qPCR) were used in the regression study of the discovery cohort (Supplementary Table SII) and to validate the models in the validation cohort (Supplementary Table SIII).

Correlation analysis

The corrplot package (Wei ) of the R 3.6.2 program was used to analyze the correlations between the proteins (2019-12-12, R Foundation for Statistical Computing, Vienna, Austria).

Statistical analysis

GraphPad Prism v.8.0 (GraphPad Software, California, USA) was employed to analyze the data. The statistical significance of the experiments carried out with the setup pool cohort was determined using paired Student’s t-tests. For the results obtained for the discovery and validation cohorts, unpaired Student’s t-tests with Welch's correction were employed. Statistical differences were considered significant at a P-value smaller than 0.05 (two-sided). Sample sizes and P-values are all shown in the figures and figure captions.

Small RNA-Seq

The quantity and quality of the RNA were evaluated using Agilent RNA 6000 Pico Chips (Agilent Technologies, Cat. # 5067-1513, CA, USA). Sequencing libraries were prepared following the protocol included with the NEXTflex™ Small RNA-Seq Kit v3 (©Bioo Scientific Corp., Cat. # 5132-06, protocol V19.01, Austin, TX, USA). Briefly, the total RNA from each sample was incubated for 2 min at 70°C. Then, a 3′ 4 N adenylated adapter (adapter dilution 1/4) and ligase enzyme were added, and ligation was carried out by incubation overnight at 20°C. After removing the excessive 3′ adapter, 5’ adapter was added with the ligase enzyme and the mixture was incubated at 20°C for 1 hr. The ligation product was used for reverse transcription with the M-MuLV reverse transcriptase in a thermocycler for 30 min at 42°C and 10 min at 90°C. Next, the enrichment of the cDNA was performed using PCR cycling: 2 min at 95°C; 20–27 cycles of 20 s at 95°C, 30 s at 60°C and 15 s at 72°C, with the final elongation of 2 min at 72°C and a pause at 4°C. The PCR products were resolved on 8% Novex TBE polyacrylamide gels (Cat. # EC6215BOX, Thermo Fisher Scientific), and a band between 150 and 400 bp was cut out. Small RNAs were extracted from the polyacrylamide gel using an adapted protocol in which the DNA from gel slices was dissolved in ddH2O overnight at room temperature. Afterwards, the libraries were visualized employing an Agilent 2100 Bioanalyzer with an Agilent High Sensitivity DNA kit (Agilent Technologies, Cat. # 5067-4626) and quantified using a Qubit dsDNA HS DNA Kit (Thermo Fisher Scientific, Cat. # Q32854). The amount of cDNA in each library that was sent for sequencing was 10 nM. Sequencing was carried out in pools of isomolar libraries and all of them were sequenced in a HiSeq2500 (Illumina Inc) to achieve at least 10 million 50-nt single-reads per sample.

Alignment

The FASTQs were trimmed for the adapters following the recommendations of the NEXTflex™ Small RNA-Seq Kit manufacturers. We used the Bowtie program (Langmead ) to align the reads against the human genome (GRCh38), with a mismatch of 0 to avoid false positives. We chose miRBase v22 to quantify the mature miRNAs, employing the Partek Flow application (version 7.0).

Small RNA-Seq data analysis

We performed differential abundance analyses to identify miRNAs associated with different implantation outcomes. To avoid rare molecules, following the Trimmed Mean of M-values (TMM) normalization, we retained miRNAs with counts per million > 1, non-zero counts in at least 15 individuals, and at most 10 zero counts in each of the two subgroups, i.e., the successful (n = 15) and unsuccessful implantations (n = 15). Differential expression was then assessed employing the edgeR (Robinson ) using the SARTools R package (Varet ). The program fits a log-linear model for each miRNA that uses a group (implantative versus non-implantative) as the factor of contrast. Applying the edgeR default parameters for normalization and shrinkage, this gives a fold change estimate that corresponds to the mean expression level in the implantative samples divided by the mean expression level in the non-implantative group. For further analysis, we selected the miRNAs with logFC > 1.5 or logFC < −1.5 and the adjusted P-value < 0.05. The Benjamini–Hochberg procedure was used to calculate the false discovery rate for each comparison and obtain the adjusted P-values (Supplementary Tables SIV and SV).

Regression study

A subset of miRNAs was used to generate two linear regression models with k-fold cross-validations, one for each miRNA extraction protocol assessed. Samples were randomly divided into training and testing datasets (80–20%). Three miRNAs were used per modeling process. The hsa-miR-24-3p, hsa-miR-200b-3p and hsa-miR-148b-3p were selected for PBP-M and hsa-miR-24-3p, hsa-miR-200b-3p and hsa-miR-99b-5p for PBP-N. The resulting model reproducibility was further tested by bootstrap correction with 500 replications. The analysis was performed using R v4.0.0 software (R Development Core Team; http://cran.r-project.org) with ROCR (Sing ) and caTools packages.

Functional analysis of the miRNAs

The target genes of the validated miRNAs were obtained from the TarBase database, v7.0. The biological processes in which these miRNAs are involved were analyzed using the Kyoto encyclopedia of genes and genomes (KEGG) and gene ontology (GO) in terms of biological process categories employing Diana-miRPath tools v3.0 (Vlachos ). The Fisher's exact test and false discovery rate correction were performed to select enriched KEGG pathways and GO processes. We selected only the pathways and processes with P-values < 0.05. The results for the KEGG were merged by ‘pathway union’ and the results for GO by ‘category union’.

Results

Optimization of EF sample preparation

The treatment of the samples with 1.4% DTT (Miller ; Wang ) was useful for degrading the mucus pellet formed after centrifugation; most mucus disappeared, as shown in Supplementary Fig. S2A. The NTA and cryo-electron microscopy analyses (Supplementary Fig. S2B and C) revealed heterogeneous EV populations with diameters between 100 and 800 nm under both experimental conditions (with and without DTT). In the untreated samples, the average concentration was 1.9 × 109 ± 7.8 × 107 particles/ml, with a mean size of 291.5 ± 0.1 nm and mode 196.4 ± 3.2 nm. In the DTT-treated samples, we detected more particles (mean 2.7 × 109 ± 5.9 × 107 particles/ml), with larger mean size (mean 313.6 ± 2.5 nm and mode 248.5 ± 8.1 nm). The WB showed different patterns of vesicular markers in the two conditions (Supplementary Fig. S2D). The intensity of the Rab8 marker in the DTT-treated samples was stronger than in the untreated samples, in agreement with the number of particles detected in the NTA. In contrast, the intensities of the Limp II, CD133 and CD63 markers were stronger in the untreated samples. The seven reference miRNAs (Supplementary Table SI) were detected in all the samples for both conditions. In the DTT-treated group, the levels of the following miRNAs were significantly reduced compared to the untreated group: hsa-let-7e-5p, hsa-miR-17-5p, hsa-miR-200c-3p, hsa-miR-30c-5p and hsa-miR-451a (Supplementary Fig. S2E).

Characterization of the miRNAs in EF

The SEC method separated the EF into several fractions. The results of WB analysis of the fractions demonstrate that it was possible to detect exosomal markers in small-volume EF samples (Fig. 2A). The CD63 and CD81 markers were detected in the F3 and to a lesser extent in F4 and F5 fractions. Rab8 was also mainly detectable in F3–F5 fractions. Immunoglobulins were also found in fractions F6 to F11. The study of the distribution of the seven reference miRNAs (Supplementary Table SI) showed that the relative quantity of miRNAs increased in fractions F3 to F11 (Fig. 2B). In F3, which corresponds to the vesicular fraction, the most abundant miRNAs were hsa-miR-451a and hsa-miR-92a-3p, and the least abundant were hsa-miR-30d-5p and hsa-miR-200c-3p. This trend was maintained in the rest of the fractions, except for F7, where hsa-miR-200c-3p was the second most abundant miRNA.

Figure 2.

Characterization of the microRNAs (miRNAs) in the endometrial fluid of patients undergoing ART. (A) Western blot shows different EV markers (CD63, CD81 and RAB8) and soluble proteins (Igs) in the fractions of size-exclusion chromatography (SEC). The fractions obtained by SEC were numbered from F1 to F12. (B) Distribution of the seven reference miRNAs among the fractions of the SEC. Normalized relative quantification was used to detect the miRNAs in the fractions. To perform experiments A and B, a 400-µl sample aliquot from the setup pool cohort was added onto the column. The number of replicates for each fraction was six and the data show the mean with SEM. (C) RNase protection assay. Sample analysis to examine the association of miRNAs with proteins and EVs. The graphs show the Ct values of the reference miRNAs evaluated using the qPCR. The number of replicates for each condition was six and the data show the mean with SEM. The number of replicates in which each miRNA was detected is shown at the bottom of each column. Each aliquot (400 µl) came from the setup pool cohort. Statistical significance was determined using the paired Student’s t-test analysis. *,$,&,#P < 0.05; **,$$,&&,##P < 0.01; ***,$$$,&&&,###P < 0.001. * versus Control, $ versus RNase, & versus TX-100, # versus PRT-K. Control: control sample without treatment. RNase: samples treated with RNase. TX-100: samples treated first with Triton-X 100 (TX-100) followed by RNase treatment. PRT-K: samples treated first with proteinase K and then with RNase. TX-PRT: samples treated first with TX-100, then with proteinase K and finally with RNase. EVs, extracellular vesicles; qPCR, quantitative PCR.

The protective effect of the EVs on the miRNAs was confirmed by the RNase assay. In the samples treated with RNase or Triton X-100 (TX-100) with RNase, only the hsa-miR-30c-5p was significantly degraded compared to the control (Fig. 2C). In the samples treated with proteinase K (PRT-K) and RNase, there was a significant decrease in the levels of all the analyzed miRNAs (although all the miRNAs were detectable in all the replicates). The miRNAs were further degraded when the samples were treated with TX-100, PRT-K and RNase (TX + PRT). In this case, we could only detect hsa-miR-200c-3p and hsa-miR-92a-3p in all the replicates. Under these conditions, hsa-miR-451a, hsa-miR-17-5p and hsa-miR-30d-5p were detected in four of six replicas, hsa-let-7-5p in two of six replicas and hsa-miR-30c-5p was undetectable. In the TX + PRT treatment, the detection of all miRNAs was significantly reduced compared to the previous combinations.

Identification of an efficient method to perform a comprehensive analysis of miRNAs from EF

Although we detected all the reference miRNAs (Supplementary Table SI) in all the extraction replicates obtained using the DCT-M, DCT-N, PBP-M, PBP-N and UC-M protocols, differences in the abundance of each miRNA were found among them (Fig. 3A). The miRNA analysis showed that the methods employing the EV enrichment step with PBP (PBP-N and PBP-M) performed better than the others, with PBP-N being the most efficient method. The protocols using direct extraction (DCT-N and DCT-M) and the UC-M method obtained significantly smaller signals of our reference miRNAs than the PBP-N procedure (Fig. 3A). Small RNA-Seq analysis of the RNA obtained using PBP-N and PBP-M detected 151 and 251 unique miRNAs, respectively. Of these, 145 miRNAs were shared between the two methods (Fig. 3B). Moreover, the small RNA-Seq analysis of fractions F3 and F4 detected 204 and 149 unique miRNAs, respectively. The samples obtained using PBP methods and SEC fractions F3 and F4 shared a large number of miRNAs (Fig. 3B). The percentage of alignments with the human genome for PBP-M, PBP-N, F3 and F4 were 74 ± 5.9%, 68.4 ± 1.8%, 54 ± 6.8% and 67.7 ± 5.1%, respectively.

Figure 3.

Optimization of different methods for analyzing the miRNAs in endometrial fluid of patients undergoing ART. (A) Results for the seven reference miRNAs analyzed by quantitative PCR for each of the compared techniques. Normalized relative quantification revealed that the most efficient method was the PBP-N, while the UC-M method was the least efficient. Statistical significance was determined using paired t-test analysis. The number of replicates for each case was 12 and the data show the mean with SEM. * versus PBP-N; $ versus PBP-M; & versus DCT-N; # versus DCT-M. (B) The Venn diagram shows the number of unique miRNAs detected using small RNA-Seq for each method and the number of miRNAs common among them. The number of unique miRNAs detected by each technique was 251 for PBP-M, 151 for PBP-N, 204 for SEC F3 and 149 for SEC F4. The samples (400 µl) for experiments A and B came from the setup pool cohort, and each experiment was performed in triplicate. (C) A technical reproducibility experiment was conducted to compare the performance of PBP-M and PBP-N methods. The graphs show Ct values for each miRNA, each operator (a, JIP; b, MCG) and method (PBP-M or PBP-N). Box plots show the median, maximum and minimum values and all the points. The 400-µl samples came from the setup pool cohort. Each operator analyzed 20 aliquots, 10 by employing the PBP-M and 10 by PBP-N. Statistical significance was assessed using the paired t-test analysis of the total results obtained with PBP-N and PBP-M. *,$,&,#P < 0.05; **,$$,&&,##P < 0.01; ***,$$$,&&&,###P < 0.001; ****,,&&&&,####P < 0.0001. DCT-M: direct extraction of RNA with mirVana PARIS kit. DCT-N: direct RNA extraction with Norgen Plasma/Serum RNA purification kit. PBP-M: extracellular vesicle enrichment using the polymer-based precipitation method and RNA extraction with mirVana PARIS kit. PBP-N: extracellular vesicle enrichment with a polymer-based precipitation method and RNA extraction with Norgen Plasma/Serum RNA purification kit. SEC-M: extracellular vesicle enrichment with size-exclusion chromatography and RNA extraction with mirVana PARIS kit. UC-M: extracellular vesicle enrichment by ultracentrifugation before RNA extraction with mirVana PARIS kit. miRNAs, microRNAs; PBP, polymer-based precipitation; RNA-Seq, RNA-sequencing.

The technical reproducibility experiment showed that using the PBP-N method, the nine miRNAs were detected by both operators (Supplementary Table SI). However, using the PBP-M protocol, no operator could detect hsa-miR-155-5p and hsa-miR-21-5p in all the replicates (Supplementary Table SI). PBP-M showed bigger differences between coefficients of variation than PBP-N, both for the different aliquots and the operators performing the experiment (Supplementary Table SI). In addition, the normalized relative quantification of the miRNAs showed significant differences between the techniques, with the PBP-N method being the most efficient, as judged by the detection of all the miRNAs with lower Ct values (Fig. 3C).

Performance of the selected methods in a set of samples with different implantation outcomes

PBP-M and PBP-N, the most efficient of the tested protocols, were chosen for implementation in the discovery cohort. There were no significant differences between the clinical characteristics of the women for whom the implantation was successful and those for whom the procedure failed (Table I). We observed large variability in the total amounts of protein in the samples (Supplementary Fig. S3A and Supplementary Table SII). The WB analysis revealed that the expression of albumin and Igs was higher in the samples with higher total protein concentration, with correlation coefficients (r) of 0.65 (P-value = 0.0066) and 0.86 (P-value = 0.0001), respectively (Supplementary Fig. S3B). The vesicular markers, such as flotillin-1 (r = 0.76, P-value = 8.7 × 10−7), Rab8 (r = 0.89, P-value = 1.4 × 10−12) and HSP90 (r = 0.82, P-value = 2.3 × 10−7), were also detected in most of the samples and positively correlated with total protein content. However, the exosomal tetraspanins CD63 and CD9 were detectable in some samples but not correlated with the total protein quantity (Supplementary Fig. S3C). Both the Coomassie blue and WB results were analyzed by densitometry of non-saturated films (data summarized in Supplementary Table SII).

Table I

Main characteristics of the study population of women undergoing ART.

	Discovery cohort		Validation cohort
	Implantative endometrium	Non-implantative endometrium	Implantative endometrium	Non-implantative endometrium
	(n = 15)	(n = 15)	(n = 30)	(n = 30)
Woman’s age at transfer (years)	36.7 ± 2.6	36.3 ± 1.8	36.6 ± 2.3	36 ± 3.4
Woman’s age at cryopreservation (years)	35.5 ± 2.4	35.1 ± 1.7	35.8 ± 2.5	35.1 ± 3.4
BMI (kg/m2)	26.2 ± 4.4	25.3 ± 1.9	24 ± 5.2	23.9 ± 4.7
Smokers (%)	26.7	25	25	14.3
Primary infertility (%)	78.6	73.3	53.6	64.3
Previous insemination failure (%)	40	41.6	17.3	25
Male factor (%)	37.5	26.7	14.8	39.3
Tubal factor (%)	6.7	0	28.6	17.9
Estradiol on the day of hCG (pg/ml)	3580.1 ± 1820.2	4101.4 ± 1154.6	3807.9 ± 2560.3	3937.1 ± 1445.5
Oocytes obtained	13.5 ± 5.2	15.6 ± 7.5	14.1 ± 6.9	13.3 ± 6.7
Metaphase II oocytes	11.73 ± 5.1	13.2 ± 5.7	12.3 ± 5.9	11.2 ± 6.1
Fertilized oocytes	7.2 ± 2.4	7.9 ± 3.8	8 ± 4.8	7.5 ± 4.5
Frozen embryos	3 ± 1.7	3.9 ± 3.9	4 ± 3	3.6 ± 2.5
Embryos transferred	1.3 ± 0.5	1.3 ± 0.5	1.4 ± 0.5	1.1 ± 0.4
Twins (%)	6.7	NA	0	NA

Statistical significance was assessed using the unpaired Student’s t-test. There were no significant differences between the clinical characteristics of the women in the discovery and validation cohort (implantative versus non-implantative) or discovery versus validation cohort. Data are expressed as mean ± SD unless specified otherwise. NA, not applicable. The endometrium was considered implantative when pregnancy was confirmed by vaginal ultrasound showing a gestational sac 4 weeks after embryo transfer.

The post-alignment quality assurance/quality control of the small RNA-Seq analysis showed great variability among the samples (Supplementary Figs S4 and S5 and Supplementary Table SII). After the TMM normalization of the small RNA-Seq results, we considered the 341/910 and 231/845 unique miRNAs for further differential abundance analysis for the PBP-N and PBP-M datasets, respectively (Fig. 4A). Statistical analysis applied to the PBP-M data detected 13 miRNAs suitable for further validation by qPCR in the same samples. In the case of PBP-N, five miRNAs were deemed suitable for validation (fold changes and P-values are shown in Supplementary Tables SIV and SV). The NormFinder algorithm identified the hsa-miR-200c-3p (stability for PBP-M, 0.31 and for PBP-N, 0.24) and hsa-miR-92a-3p (stability for PBP-M, 0.27 and for PBP-N, 0.1) as the most suitable pair of normalizer miRNAs (Fig. 4B) (Supplementary Table SI).

Figure 4.

Performance of the selected methods (PBP-M and PBP-N) in the discovery cohort. These experiments were conducted using the discovery cohort samples (n = 30): 15 samples from women with successful implantation and 15 from women with implantation failure. The endometrium was considered implantative when pregnancy was confirmed by vaginal ultrasound showing a gestational sac 4 weeks after embryo transfer. (A) The Venn diagram shows the number of miRNAs common for the PBP-M and PBP-M (n = 230) detected by small RNA-Seq (after TMM normalization) and the number of miRNAs unique for each of the techniques (PBP-M, 1 and PBP-N, 111). (B) Correlation analyses were carried out to determine the suitability of the selected internal controls. The graphs show the results for each miRNA in each sample (n = 30), obtained using the small RNA-Seq and qPCR. The miRNAs selected as internal normalizers were hsa-miR-200c-3p and hsa-miR-92a-3p. (C) Receiver operating characteristic curves for the three miRNA-based predictive models tested by qPCR in the discovery cohort. The performance of the original model is shown in black, and the mean performance of the bootstrap output is shown in red. The shading indicates the extent of the standard deviation. The AUC is shown in the respective colors in the lower right-hand corner of the curves. PBP-N: extracellular vesicle enrichment using polymer-based precipitation and RNA extraction with Norgen Plasma/Serum RNA purification kit. PBP-M: extracellular vesicle enrichment using polymer-based precipitation method and RNA extraction with mirVana PARIS kit. miRNAs, microRNAs; PBP, polymer-based precipitation; qPCR, quantitative PCR; RNA-Seq, RNA-sequencing; TMM, Trimmed Mean of M-values.

The predictive model was chosen by conducting a regression study using bootstrapping correction of the normalized qPCR data (dCt) for the differentially expressed miRNAs in the discovery group. For the PBP-M method, Model 1 (AUC = 0.93; P-value = 3.3 × 10−3), based on three miRNAs (hsa-miR-24-3p, hsa-miR-200b-3p and hsa-miR-148b-3p), was robust against the data that fitted the given characteristics. The results for PBP-N also highlighted a predictive model based on three miRNAs, Model 2 (hsa-miR-24-3p, hsa-miR-200b-3p and hsa-miR-99b-5p) (AUC = 0.92; P-value = 2.3 × 10−4) (Fig. 4C). Both models were significantly predictive compared to random chance results (AUC = 0.5).

The regulatory functions of these miRNAs were broken down into biological process categories by analyzing their predicted target genes. The KEGG pathways showed many different enriched pathways, including adherens junction proteins, transforming growth factor (TGF)-beta signaling, fatty acid biosynthesis and fatty acid metabolism, all of which turned out to be significantly enriched (Supplementary Fig. S6). The GO analysis also showed TGF-beta signaling pathway enrichment (Supplementary Fig. S7). The in utero embryonic development, immune system processes, endosome and vesicle-mediated transport processes were also enriched. Overall, most detected pathways were closely related to embryo implantation and endometrial decidualization.

Validation of the models in an independent cohort of samples with different implantation outcomes

The performance of the two predictive models (Model 1 and Model 2) was validated in an independent cohort. There were no significant differences between the clinical characteristics of the women in the validation cohort (implantative versus non-implantative) or validation versus discovery cohort (Table I) or between the total RNA and protein amounts in the different groups (Supplementary Table SIII). The results showed that our two models based on miRNA signature remained predictive of the implantation status, especially Model 2.

Model 1 (hsa-miR-200b-3p, hsa-miR-24-3p and hsa-miR-148b-3p) has an accuracy of 0.68 (95% CI (0.54, 0.8)) and an AUC of 0.69 (95% CI (0.55, 0.86)) (Fig. 5). The test showed a statistically significant difference (P-value: 0.019) in AUC compared to random chance (AUC = 0.5). The receiver operating characteristic (ROC) curve analysis identified hsa-miR-148b-3p as the mostly likely variable to differentiate between implantative and non-implantative endometrium. A significant differential expression of this miRNA was seen in group comparisons (P-value 0.02); it was upregulated in the non-implantative group.

Figure 5.

Receiver operating characteristic curves of the validated models and box plots for the most significant miRNAs in each model. The predictive models designed using the qPCR results for the PBP-M and PBP-N methods applied to the discovery cohort were validated in a new group, the validation cohort (n = 60; 30 subjects in the implantative subgroup and 30 in the non-implantative subgroup). The endometrium was considered implantative when pregnancy was confirmed by vaginal ultrasound showing a gestational sac 4 weeks after ET. The analyses were carried out with those samples that pass the quality control; amplification of the reference miRNAs (hsa-miR-200c-3p and hsa-miR-92a-3p) less than 30 Cts. For PBP-M were: n = 28 in the implantative group and n = 25 in the non-implantative group. For PBP-N were: n = 27 in the implantative group and n = 25 in the non-implantative group. (A) Receiver operating characteristic curves for the qPCR data obtained for the two predictive models in the validation cohort. Model 1 PBP-M (hsa-miR-200b-3p, hsa-miR-24-3p and hsa-miR-148b-3p) had an AUC of 0.69 (95% CI, 0.55–0.86) and Model 2 PBP-N (hsa-miR-200b-3p, hsa-miR-24-3p and hsa-miR-99b-5p) had an AUC of 0.78 (95% CI, 0.6–0.89). (B) Box plots showing the most likely miRNAs to differentiate between the non-implantative and implantative endometrium. The microRNA levels in the EF at the time of embryo transfer. The horizontal line in the middle of the box plot represents the median, while the horizontal limits of the boxes represent the first and third quartiles. The levels of significance were assessed using unpaired t-tests with Welch's correction. ΔCt is inversely correlated with the amount of miRNA in the samples. PBP-N: enrichment with polymer-based precipitation method and RNA extraction with Norgen Plasma/Serum RNA purification kit. PBP-M: extracellular vesicle enrichment with a polymer-based precipitation method and RNA extraction with mirVana PARIS kit. Ct, cycle threshold; ET, embryo transfer; miRNAs, microRNAs; PBP, polymer-based precipitation; qPCR, quantitative PCR.

We established the cutoff point as 2.34 dCt (dCt = Cthsa-miR-148b-3p − Ctinternal controls) based on Youden's J statistic with sensitivity of 0.56 and specificity of 0.86. Values above the cutoff point would indicate an implantative endometrium (dCt > 2.34) (negative predictive value = 0.69) and values below, a non-implantative endometrium (dCt < 2.34) (positive predictive value = 0. 78).

Model 2 (hsa-miR-200b-3p, hsa-miR-24-3p and hsa-miR-99b-5p) had an accuracy of 0.77 (95% CI (0.63, 0.88)) and an AUC of 0.78 (95% CI (0.6, 0.89)) (Fig. 5). The test also showed a statistically significant difference (P-value: 0.0002) in AUC compared to random chance results (AUC = 0.5). The ROC analysis identified hsa-miR-99b-5p as the most able variable to distinguish the non-implantative from implantative endometrium. A significant difference in the expression of this miRNA (≥1.5-fold) was found between the groups at P-value 0.0004 (the miRNA upregulated in the implantative group). We established the cutoff point as 2.81 dCt (dCt = Cthsa-miR-99b-5p − Ctinternal controls) based on Youden's J statistic with 0.6 sensitivity and 0.93 specificity. Values above the cutoff point would indicate a non-implantative endometrium (dCt > 2.81) (positive predictive value = 0.88) and values below, implantative endometrium (dCt < 2.81) (negative predictive value = 0. 71) (Fig. 5).

Discussion

The endometrium undergoes a series of changes during the ovarian cycle until the endometrial glands achieve maximal secretory activity 6 days after ovulation, which is necessary to create an optimal uterine microenvironment for embryo implantation (Gellersen and Brosens, 2014). It has been widely believed that the embryo has the most important role in human implantation. However, in many cases, even when the maternal conditions are apparently optimal and the transferred embryo is chromosomally normal, implantation does not occur (Cozzolino ). The term ‘implantative endometrium’ was coined to signify the endometrium in which implantation succeeds in the same cycle as the EF aspiration (Matorras ). Thus, the aim of this study was to investigate predictive markers of the implantation failure. We have found differences in the miRNA patterns between implantative and non-implantative cycles, and these results have led us to define two predictive models of implantative endometrium.

In this study, we optimized the EF preparation starting from a small volume of sample. As described previously for sputum (Miller ; Wang ), pre-treatment of the samples with DTT increased the number of particles released from the mucus. However, it did not improve the relative quantification of the reference miRNAs in such samples (Supplementary Fig. S2). This finding and the report that DTT treatment might modify RNA–protein interactions (Zaman ) prompted us to collect and process the EF samples without DTT. The RNase experiment suggested that most of the miRNAs found in the EF were not only protein-associated but also protein-associated within the EVs, which protects them from degradation. Some authors have described the circulating Ago2-miRNA complexes in human plasma, which suggests that Ago2 protein might play an important role in the stability of secreted miRNA (Arroyo ; Groot and Lee, 2020). This protein has been identified within exosomes and has also been shown to protect the miRNAs within EVs from RNase degradation (Li ; Groot and Lee, 2020). However, other authors have not found the Ago2 in classical exosomes and they believe that there is no evidence that exosomes, or any other type of small EV, contain other major components of the miRNA biogenesis machinery (Jeppesen ). It is not clear which proteins are protecting these miRNAs from degradation; nevertheless, our results suggest a strong association between miRNAs and proteins and indicate the existence of miRNA-protein complexes within the EVs.

Once we had characterized the EF, we compared different vesicle enrichment and RNA extraction methods. A few studies have compared various EV-enriching protocols for EF samples (Campoy ; Li ), but none has evaluated different techniques for RNA extraction. A study published by Li reports that, for EV isolation, UC is superior to PBP. However, our findings do not agree with their results, and the differences could be due to the fact that different protocols are used. To isolate EVs from EF they (Li ) have used a 1:2 ratio of PBP and an overnight incubation at 4°C, while we used 1:1 ratio of PBP and incubation at room temperature for 30 min. The UC methods were also different; they (Li ) used a two-step UC procedure and resuspended the EVs in 35 µl, while we only performed a single step UC and resuspended the EVs in 100 µl.

In our search for efficient methods for a comprehensive analysis of miRNAs from EF, we found that the PBP-based EV enrichment techniques (PBP-M and PBP-N) increased the efficiency of miRNAs detection in EF samples (Fig. 3A). Moreover, our small RNA-Seq analysis revealed that different populations of miRNAs could be obtained depending on the RNA extraction method (Fig. 3B). Similar findings have been reported by other authors who have made such comparisons, showing different results with different RNA extraction kits (El-Khoury ; Wright ). Therefore, changing the RNA extraction protocols may lead to different results, complicating the comparisons between studies. Given the importance of selecting a robust methodology, we conducted a technical reproducibility experiment to compare the PBP-M and PBP-N methods. The results demonstrated that the PBP-N improves the qPCR amplification results and lowers the coefficients of variation (Fig. 3C).

Finally, we ended up with two predictive models, Model 1 using the PBP-M method and Model 2 based on PBP-N. The results obtained with Model 1 (PBP-M) showed differences between the discovery and validation cohorts. Specifically, the qPCR showed that hsa-miR-148-3p was significantly downregulated in the validation implantative endometrium group. This was consistent with the data obtained by the qPCR for the discovery group. However, these results did not agree with the data obtained by small RNA-Seq for the discovery cohort. In addition, hsa-miR-200b-3p did not follow the same trend in the discovery and validation cohorts. The only miRNA that followed the same trend in the different analyses and cohorts was hsa-miR-24-3p (downregulated in the implantative endometrium groups). In Model 2 (PBP-N), the results obtained for hsa-miR-24-3p were not consistent between different analyses; the qPCR results indicated that it was significantly downregulated in the implantative endometrium of the discovery cohort while, in contrast, it appeared upregulated in the small RNA-Seq analysis of that subgroup and in the qPCR results for the validation cohort. However, the results for the hsa-miR-200b-3p and hsa-miR-99b-5p were consistent, following the same trend in all the analyses, and these two miRNAs were upregulated in the implantative endometrium group (Supplementary Table SV). To be precise, the ROC analysis identified hsa-miR-99b-5p as the most likely variable to differentiate the non-implantative and implantative endometrium, and the results suggest that the dCt values above the cutoff point indicate a non-implantative endometrium (dCt > 2.81).

Overall, the PBP-N method delivered the best results in this optimization procedure, starting with the limited amounts of EF obtained in a clinical setting. In addition, it proved to be the most efficient and reproducible, with a simplified protocol and its application in two independent cohorts has shown consistent results. Currently, only 35% of the implantation attempts are successful (De Geyter ) and, as our test in the validation cohort achieved high specificity (Model 2, 0.93), we believe that we may be able to improve the success rate by using our predictive model PBP-N.

One limitation of our approach is the inherent variability of the women involved in the trial and the embryos transferred. When the pregnancy is achieved after ET, it can be said that the endometrium was implantative. However, when the implantation fails, it cannot be assumed that the problem was associated with the endometrium only, because the fault might lie with the embryo, the endometrium or both. In our study, we tried to minimize the issue of individual variability as all transfers were performed during an artificial cycle under the same hormone supplement protocol. Another limitation was the alignment rate of the sequenced data against the human genome, which showed great variability between the samples, with highly variable alignment percentages (Supplementary Table SII). We believe that these differences could be related to the different biological content of the samples, because the EF composition can vary as it may contain genetic material from the uterine microbiomes (Agostinis ). To overcome the problem of the biological variability effect on small RNA-Seq data, we used the TMM normalization at the time of selecting the differentially expressed miRNAs. In the case of qPCR data, we selected two endogenous miRNA controls (hsa-miR-200c-3p and hsa-miR-92a-3p) with the help of the NormFinder algorithm (Fig. 4B). These miRNAs were chosen because they were equally expressed in both groups (implantative and non-implantative endometrium groups) and positively correlated with the amount of protein in the samples. Moreover, they were detected at high levels by small RNA-Seq and had low Ct values in the qPCR. However, the normalization for analyzing the small RNA-Seq data and selecting an endogenous miRNA suitable for normalizing the qPCR data are critical points that might generate great uncertainty. Since there is no standardized method to normalize the RNA-Seq data, each group selects the method that they consider most appropriate. This is also the case with internal controls as the expression of such a control could also vary depending on the kit used for RNA extraction. The internal controls used should not be generalized, and each should be adapted to the experiment performed. The relatively low negative predictive value observed here in the models (Model 1 = 0.69, Model 2 = 0.71) is another limitation to be taken into account. We hypothesize that this was caused by the fact that we have only evaluated the endometrium, and this means that the test will fail in cases where the implantation failure is caused by the embryo. We selected the embryos on the basis of their morphology; however, we did not use genetic and molecular data to improve the selection. Obtaining such data would certainly improve our test accuracy and the AUC.

Functional analysis of the validated miRNAs showed strong associations with key processes involved in implantations. Thus, some of the pathways targeted by the differentially expressed miRNAs were related to adherens junctions, necessary for the initiation of implantation as they are required for cell attachment, adhesion and recognition (Buck ). These miRNAs were also associated with the TGF-beta signaling pathway, essential for decidualization of the endometrial stromal cells (Jones ). Furthermore, interactions with immune system processes (Lee ), vesicle-mediated transport and in utero embryonic development (Kurian and Modi, 2019), with key roles during implantation, were also found. In addition, the available hsa-miR-148b-3p data suggest that its activity may depend on tissue and cell types and is mainly involved in the regulation of cell progression (Dai ). Some studies have found that this miRNA inhibits malignant tumor progression (Wang ; Li ). Its overexpression has also been associated with osteogenesis (Mollazadeh ) and cancer cell progression (Dai ). Furthermore, this miRNA has been selected as a reference in a study designed to identify candidate miRNA markers of endometriosis, as its mean Ct values did not differ significantly between the women with endometriosis and the control group. Here, we observed that the expression of hsa-miR-148b-3p was upregulated in the EF samples of patients with implantation failure. This leads us to believe that high concentrations of this miRNA could be inhibiting processes related to embryo implantation. The role of hsa-miR-99b-5p in the EF could be related to its function as a pathway regulator, contributing to natural killer (NK) cell activation and effector function (Petty ). The NK cells have been described as the major leukocytes in the endometrium. They accumulate extensively around spiral arterioles in the mid-secretory-phase endometrium and early-pregnancy decidua in accordance with increasing levels of ovarian-derived estrogen and progesterone (Quenby and Farquharson, 2006). These findings suggest that the NK cells have a crucial role in implantation and decidualization. However, a meta-analysis using 22 studies examining the uterine NK-cell percentages in infertile versus fertile women showed no significant differences between the groups (Seshadri and Sunkara, 2014). Taken together with our results, these studies could indicate that in women with a low concentration of hsa-miR-99b-5p in the endometrium, the activation of uterine NK cells is suboptimal despite a normal cell count and, as a consequence, implantation does not occur.

In summary, this study introduces new protocols to analyze the miRNAs in very small volumes (5–50 µl) of EF collected just before Day-5 frozen ETs, which could be implemented in clinical practice. These new methods could be employed to assess endometrial competence using miRNA-based non-invasive tools. Hence, the professionals in assisted reproduction centers could use hsa-miR-99b-5p (employing the PBP-N detection method) to predict the endometrial status. This could potentially help to improve the implantation rates for women undergoing ART. It could be possible to change the ET strategy when the results showed an unfavorable implantative pattern and thus increase the implantation rates. Using this method may also reduce the loss of embryos, so common after ET to a potentially non-implantative endometrium.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.