From January 2015 to February 2021, IVF/ICSI cycles with live births at the Reproductive Medicine Center of Tongji Hospital were retrospectively identified. Inclusion criteria were shown below: a) cycles resulting in singleton births; b) fresh and FET cycles; c) cycles with gonadotropin releasing hormone (GnRH) stimulation protocols in fresh cycles; d) natural cycles and cycles with hormone replacement in FET cycles. Moreover, exclusion criteria were as follows: a) cycles with donor eggs; b) cycles with non-gonadotropin releasing hormone (GnRH) stimulation protocols; c) cycles with preimplantation genetic diagnosis; d) cycles with missing baseline and follow-up data. All information was retrieved from the electronic medical record system. This study was approved by the Medical Ethics Committee of Tongji Hospital, Tongji Medical College, affiliated Huazhong University of Science and Technology (TJ-IRB20220624).

As mentioned above, fresh cycles with GnRH stimulation protocols were included in our analysis. In light of ovarian reserve, female age, and other characteristics, the individualized COS protocol for each patient was determined. Generally, women with a normal ovarian response were provided with a long GnRH agonist (GnRH-a) regimen, while a flexible GnRH antagonist (GnRH-ant) protocol was offered to women with a diminished ovarian reserve or a known poor response (15). There are four COS protocols, including the GnRH-a ultralong protocol, the depot GnRH-a protocol, the GnRH-a long protocol, and the GnRH-ant protocol. All these protocols were described previously (16–18).

In short, in the GnRH ultralong protocol, women were injected subcutaneously with 3.75 mg long-acting GnRH-a (Decapeptyl; Ferring, SaintPrex, Switzerland) on days 1-3 of cycles to achieve pituitary suppression (16). The depot GnRH-a protocol was performed by subcutaneously administering 3.75 mg long-acting GnRH-a (Decapeptyl; Ferring, Saint-Prex, Switzerland) on Day 2 of the cycle. Moreover, in the GnRH-a long protocol, patients were provided with a daily subcutaneous injection of 0.1 mg GnRH-a (Decapeptyl; Ferring, Saint-Prex, Switzerland), initiating from the mid-luteal phase. In GnRH-a protocols, when complete pituitary desensitization was confirmed (a low serum estradiol (E₂) level of ≤30 pg/mL and a serum luteinizing hormone level of ≤ 2 IU/L), gonadotropins (Gn) were applied. While in the GnRH-ant protocols, the stimulation process starts with the administration of recombinant follicle stimulating hormone (r-FSH) (Gonal-F; Merck-Serono, Geneva, Switzerland) from Day 2 or 3 of the cycle. Then, the GnRH-ant cetrorelix acetate (Cetrotide; Merck-Serono) was injected at 0.25 mg/d to prevent premature ovulation when one of the criteria was met (E₂ levels >300 pg/ml, LH levels >10 IU/l, and the diameter of at least one follicle >14 mm).

When two leading follicles reached a mean diameter of 18 mm or three follicles reached a mean diameter of 17 mm, 250 µg human chorionic gonadotrophin (hCG) (Ovidrel; Merck-Serono, Geneva, Switzerland) was used to trigger ovulation. Oocytes were retrieved 34-36 h after triggering.

The processes of oocyte fertilization, embryo culture, and embryo transfer have been described previously (15, 16). In brief, oocytes were fertilized by IVF, ICSI, or rescue ICSI (RICSI). RICSI was performed when there was no second polar body (PB) in oocytes or <25% of the oocytes presented a second PB after insemination with IVF (19). Normal fertilization was defined as the appearance of two pronuclei (2PN). Moreover, a maximum of two embryos were transferred.

On day 3 (D3) after oocyte retrieval, all of the cleavage embryos were checked. Good-quality embryos on D3 comprised seven to eight blastomeres without multinucleation and had < 20% fragments (16). The quality of blastocysts (D5 or D6) was evaluated using the Gardner Blastocyst Scoring System (20). In short, according to the cell number and junction, the inner cell mass and trophectoderm of blastocysts were graded from A to C. Good-quality blastocysts were defined as AA, AB, BA, and BB (20).

Based on the instructions of the vitrified freezing/resuscitation solution (Japan Kato), the embryos were frozen and thawed. The survival of thawed embryos was observed after 2H. The surviving cleavage embryos reached or exceeded half of the blastomere, while the surviving blastocyst partially or fully dilated (21).

The endometrial preparation was described previously (21). In brief, in the natural cycles, the follicle, endometrium and ovulation were monitored during the menstrual cycle. When the peak of luteinizing hormone (LH) was observed, progesterone (P) was given. The cleavage embryos and blastocysts were thawed on D3 and D5 respectively after ovulation. As for the hormone replacement cycles, estradiol valerate (Bayer, Germany) was taken orally from D2 to D12 during the menstrual cycle. When the endometrium was ≥7 mm and the estrogen action time was ≥10 days, P was administrated. The cleavage embryos and blastocysts were thawed on D4 and D6 respectively after ovulation.

Data were collected on female age, female BMI, women smoking, duration of infertility (years), antral follicle counts (AFC), baseline FSH, type of infertility (primary or secondary infertility), and infertility factors (tubal, pelvic, endometriosis, uterine, cervical, ovulatory, diminished ovarian reserve (DOR), male, or unexplained factors) in fresh or FET cycles.

In fresh cycles, the following factors were included in the analysis: stimulation protocols (GnRH-a ultralong, GnRH-a depot, GnRH-a long, or GnRH-ant protocols), days of Gn, total dose of Gn, E₂ and P level on the triggering day, endometrial thickness (EMT) on the triggering day, No. of oocyte retrieved, type of fertilization (IVF, ICSI, or RICSI), No. of embryo transferred (double or single embryo transferred, DET or SET), type of embryo (cleavage embryo or blastocyst), and morphology of the embryos transferred.

In FET cycles, the factors as follows were included: type of endometrial preparation, such as natural cycles and hormone replacement cycles (HRT), EMT on progesterone administration, type of fertilization (IVF, ICSI, or RICSI), DET or SET, type of embryo (cleavage embryo or blastocyst), and morphology of the embryos transferred.

During pregnancies, gestational sacs, obstetric complications, such as placenta previa (PP), placenta abruption (PA), hypertensive of pregnancies (HDP), or premature rupture of membrane (PROM), neonatal gender, birth weight and gestational age were all collected. In our study, PTD was defined as delivery with gestational age reached 28 weeks but less than 37 weeks. All relevant information in the electronic medical record system was collected and checked by two authors (Z.L. and L.C.).

All selected risk factors were enrolled to establish the nomogram model using the R package “rms”. The predictive model was calculated using the R package “lrm” and displayed as an interactive nomogram by “regplot”. The multivariate logistic regression coefficient (β) generated a risk score with each factor. Each factor had its point (corresponding line is β(X-m) terms), and the total scores were calculated by adding all the points together. The total scores corresponded to the probability of PTD occurring.

In this study, continuous data with normal distribution are presented as the mean ± SD, while non-normal distribution data are displayed as median (IQR). Categorical data are presented as numbers (percentages). T tests or Mann‐Whitney U tests were utilized to compare the continuous data, while chi-square tests or Fisher’s exact tests were used to analyze the categorical data. Univariate values with P value < 0.1 were screened for the next step of multivariate regression analysis. In the multivariate logistic regression analysis, adjusted odds ratios (AORs) and 95% confidence intervals (CIs) were used to indicate the level of association between risk factors and PTD. P value < 0.05 was considered statistically significant. SPSS 26 software was used for statistical analysis.

All statistically significant risk factors in the multivariate logistic regression analysis were further used to build a nomogram predictive model. The discrimination calibration abilities of this model were evaluated by calculating receiver operating characteristic curves (ROC), area under the ROC curve (AUC), Harrell’s concordance index (C-index), and the calibration plots in the training and validation groups (14). The value of the C-index and AUC ranges from 0 to 1, with 0.5-0.7 indicating poor performance of this model, 0.7-0.8 demonstrating fair performance, and 0.8-0.9 suggesting good performance (22). R software (version 4.0.5) was used to construct the nomogram model and analyze the ROC, AUC, C-index, and calibration curve.

A total of 22809 cycles resulting in live births were identified ( Figure 1 ). Among these cycles, 11296 cycles were not available for analysis, as one of the following exclusion criteria was met: (a) multiple live births (4527cycles), (b) donor oocyte (1913 cycles), (c) non-GnRH stimulation protocols in fresh cycles (948 cycles), (d) cycles with preimplantation genetic diagnosis (629 cycles), and (e) missing baseline or follow-up data (3279 cycles). Finally, the number of IVF/ICSI available for analysis was 11513. Then, we randomly allocated 9208 cycles (80%) to the training group and 2305 cycles (20%) to the validation group.

Table 1 described the baseline characteristics of the IVF/ICSI cycles with singleton live births between the training group and validation group. No significant difference in most characteristics was observed between the two groups ( Table 1 ). Furthermore, Table 2 demonstrated a comparison of baseline characteristics between the PTD and non-PTD groups in the training group. There were 953 cycles (10.3%) in the PTD group and 8255 cycles (89.7%) in the non-PTD group ( Table 2 ). As shown in Table 2 , the median (IQR) of birth weight was 3050 (1750) and 2060 (1210) in non-PTD and PTD group respectively. The incidence of PTD was higher if there was one of the following factors: (a) female age ≥ 35, (b) female BMI > 24 kg/m², (c) smoking, (d) high baseline FSH, (e) secondary infertility, (f) uterine factors, (g) cervical factors, (h) ovulatory factors, (i) EMT ≤7 mm, (j) two embryo transferred, (k) blastocyst transfer, (l) FET, (m) non-high quality embryo transfer, (n) two gestational sacs, (o) obstetric complications (PP, PA, GDMHDP, and PROM), and (p) a male fetus ( Table 2 ).

When only included fresh cycles, we also found the difference in female age, female BMI, smoking, cervical factors, ovulatory factors, EMT, gestational sacs, obstetric complications, and neonatal gender between PTD and non-PTD groups ( Table S1 ). In addition to the patient’s characteristics, data related to ovarian stimulation in the fresh cycle were also presented in Table S1 , including stimulation protocols, average dose of Gn, E2 or P on hCG day, and No. of oocyte retrieved. The result showed that women with GnRH-agonist long protocol had lower risk to deliver a preterm baby, while women with GnRH-antagonist protocol had higher risk ( Table S1 ).

When only included FET cycles, there was also a higher risk for PTD in women with obesity, uterine factors, DET, non-high quality embryo transfer, two gestational sacs, obstetric complications, and neonatal gender between PTD and non-PTD groups ( Table S2 ). However, no difference was observed in the type of endometrial preparation (natural cycles and HRT cycles) between PTD and non-PTD groups ( Table S2 ).

To eliminate confounding factors, multivariate logistic regression analysis was used to analyze the risk factors for PTD, which were detailed in Table 3 . Factors (female age, female BMI, smoking, baseline FSH, type of infertility, uterine factors, cervical factor, ovulatory factor, male factor, EMT, No. of embryo transferred, blastocyst transfer, FET, high-quality embryo transfer, gestational sacs, PP, PA, GDM, HDP, PROM, and neonatal gender) with P<0.1 in Table 2 were entered into the multivariate logistic regression analysis. Thirteen factors remained statistically significant. Although P value in cervical factor was 0.053 (P>0.05), we still included this factor in further analysis as this factor had very important clinical significance in the occurrence of PTD.

The results showed that compared to women whose BMI was ≤ 24 kg/m2, overweight patients had a significantly higher risk of PTD after IVF/ICSI (OR 1.25, 95% CI: 1.06-1.47). Women who had a smoking habit were more likely to deliver a preterm infant (OR 4.97, 95% CI: 3.21-7.67). Moreover, women with uterine, cervical, and ovulatory infertility factors displayed a rising PTD ratio (OR 1.47, 95% CI: 1.23-1.76; OR 2.02, 95% CI: 0.99-4.11; OR 1.24, 95% CI: 1.04-1.47, respectively). In comparison to the SET group, women who had DET had a higher risk of PTD (OR 1.37, 95% CI: 1.17-1.61). Besides, women with blastocyst transfer (OR 1.56, 95% CI: 1.19-2.03) and FET (OR 1.56, 95% CI:1.20-2.03) had higher risk of PTD. During pregnancies, women who had two gestational sacs had a greater chance of delivering a preterm baby (OR 4.17; 95% CI: 3.46-5.03). In addition, women who had obstetric complications, such as PP (OR 3.06, 95% CI: 2.24-4.19), PA (OR 10.22, 95% CI: 6.35-16.46), HDP (OR 4.26, 95% CI: 3.21-5.66), and PROM (OR 20.32, 95% CI: 13.55-30.47), also had an increasing risk of PTD. Interesting, we found a male fetus (OR 1.182, 95% CI: 1.019-1.371) was more likely to be born prematurely ( Table 3 ).

Furthermore, the above prognostic factors were incorporated to establish a predictive model for predicting the probability of PTD incidence in IVF/ICSI singleton live birth ( Figure 2 ). The nomogram developed from the training group was shown in Figure 2 . The usage of the nomogram is illustrated with an assumptive woman with BMI > 24, smoking, cervical factor, DET, FET, two gestational sac, and have a male fetus. The total score added up to 4.92 for this patient, which represents approximately 0.792 of probability of PTD incidence ( Figure 2 ).

The performance of this nomogram was assessed by AUC, C-index, and calibration plots. In the training and validation groups, both calibration plots indicated good agreement between the predicted probability and actual observation of PTD ( Figures 3A, C ). In addition, the AUC for this predictive model was 0.774 and 0.770 in the training and validation groups, respectively ( Figures 3B, D ). The optimal threshold was -2.239 and -2.011 in training and validation groups, respectively ( Figures 3B, D ). At the above thresholds, the sensitivity was 77.4% and the specificity 6 was8% in training group, while in validation group, the sensitivity was 88.1% and the specificity was 57.4% ( Figures 3B, D ). The C-index of the predictive model was 0.774 (95% CI: 0.757-0.791) in the training group and 0.770 (95% CI: 0.731-0.809) in the validation group ( Table 4 ). These results suggested that the model exhibited fair performance.