{"paper_id":"aa8c12af-b695-48f1-821d-f633241a0b21","body_text":"Article https://doi.org/10.1038/s41467-023-41979-z\nA molecular staging model for accurately\ndating the endometrial biopsy\nW. T. Teh 1,2,3,13,J .C h u n g1,4,13,S .J .H o l d s w o r t h - C a r s o n1,2,5, J. F. Donoghue1,2,\nM. Healey1,2,H .C .R e e s2,6,S .B i t t i n g e r2,6,V .O b e r s7,C .S l o g g e t t4,8,\nR. Kendarsari9,10,J .N .F u n g11,S .M o r t l o c k9,G .W .M o n t g o m e r y9,\nJ. E. Girling 1,12 &P .A .W .R o g e r s 1,2\nNatural variability in menstrual cyclelength, coupled with rapid changes in\nendometrial gene expression, makes it difficult to accurately define and\ncompare different stages of the endometrial cycle. Here we develop and\nvalidate a method for precisely determining endometrial cycle stage based on\nglobal gene expression. Our ‘molecular staging model’ reveals significant and\nremarkably synchronised daily changes in expression for over 3400 endo-\nmetrial genes throughout the cycle, with the most dramatic changes occurring\nduring the secretory phase. Our study significantly extends existing data on\nthe endometrial transcriptome, and for thefirst time enables identification of\ndifferentially expressed endometrialgenes with increasing age and different\nethnicities. It also allows reinterpretation of all endometrial RNA-seq and array\ndata that has been published to date. Our molecular staging model will sig-\nnificantly advance understanding of endometrial-related disorders that affect\nnearly all women at some stage of their lives, such as heavy menstrual bleed-\ning, endometriosis, adenomyosis, and recurrent implantation failure.\nThe endometrium plays an essential role in embryo implantation,\nplacentation and success or otherwise of pregnancy in all mammals. A\nfundamental understanding of human endometrial biology underpins\nour knowledge of everyday physiological and pathological processes\nthat include uterine receptivity, pregnancy, menstruation, heavy\nmenstrual bleeding, recurrent implantation failure, endometriosis,\nadenomyosis, endometrial cancer and pelvic pain. Nearly all women\nduring their lifetime will see their gynaecologist for one or more\nendometrial-related health problems\n1. Despite this impact on quality\nof life for most women, the endometrium remains understudied rela-\ntive to the healthcare burden of endometrial-related disorders, with\nendometriosis\n2, menstrual problems 1 and contraceptive-related\nbleeding issues3 as prime examples.\nThere are two methodological challenges that have a profound\nimpact on endometrial research: the large normal variation in men-\nstrual cycle length, and the huge variability in gene expression across\nthe cycle. Compared to other tissues in the body, the endometrium\nundergoes dramatic cyclical changes in gene expression\n4,5. Daily and\nsometime hourly changes in expression are driven by increased cir-\nculating estrogen from the developing ovarian follicles during the\nproliferative phase of the cycle, then by progesterone from the corpus\nluteum following ovulation, and in non-conception cycles, during\nReceived: 15 July 2022\nAccepted: 26 September 2023\nCheck for updates\n1University of Melbourne Department of Obstetrics and Gynaecology, Melbourne, Victoria, Australia.2Royal Women’s Hospital, Melbourne, Victoria, Australia.\n3Melbourne IVF, Melbourne, Victoria, Australia.4Melbourne Bioinformatics, University of Melbourne, Melbourne, Victoria, Australia.5Julia Argyrou Endo-\nmetriosis Centre, Epworth HealthCare, Richmond, Victoria, Australia.6Royal Children’s Hospital, Melbourne, Victoria, Australia.7Melbourne Pathology,\nCollingwood, Victoria, Australia.8Microbiological Diagnostic Unit Public Health Laboratory, Department of Microbiology and Immunology, University of\nMelbourne at the Peter Doherty Institute, Melbourne, Victoria, Australia.9Institute for Molecular Bioscience, University of Queensland, St Lucia, Queensland,\nAustralia. 10Illumina Inc. 11 Biopolis Way, Singapore 138667, Singapore.11School of Biomedical Sciences, University of Queensland, St Lucia, Queensland,\nAustralia. 12Department of Anatomy, School of Biomedical Sciences, University of Otago, Dunedin, Aotearoa, New Zealand.13These authors contributed\nequally: W. T. Teh, J. Chung. e-mail: parogers@unimelb.edu.au\nNature Communications|         (2023) 14:6222 1\n1234567890():,;\n1234567890():,;\n\nmenstruation after the demise of the corpus luteum and the loss of\ncirculating progesterone. This rapidly changing gene expression\nwithin a highly variable length menstrual cycle has made accurate\ncomparisons between matched samples dif ficult at best, and often\nimpossible. As a consequence, studies linking endometrial gene\nexpression to various endometrial-related pathologies such asfibroid-\nrelated heavy menstrual bleeding, reduced endometrial receptivity for\nimplantation, and endometriosis, seldom replicate\n6–9.\nA critical variable in assessing differential endometrial gene\nexpression between samples is accurate menstrual cycle staging.\nThere is large variability between women in overall cycle length, as well\nas days of menstrual bleeding, and follicular and luteal phase lengths\n10.\nIn a study of over 30,000 women, only 12.4% had a 28-day cycle11.M o s t\nhad menstrual cycle lengths between 23 and 35 days, with a normal\ndistribution centred on day 28, and over half had cycles that varied by\n5 days or more from cycle to cycle. There was a 10-day spread of\nobserved ovulation days for a 28-day cycle, with the most common day\nof ovulation being day 15. Another large study of 612,613 ovulatory\ncycles reported a mean length of 29.3 days from 124,648 subjects\n12.\nThe mean follicular phase length was 16.9 days (95% CI: 10 –30) and\nmean luteal phase length was 12.4 days (95% CI: 7 –17). Part of the\nvariability in cycle length between women was due to age, with a\nconsistent shortening of the average cycle length by about 3 days from\n30 down to 27 days between ages 25 and 45\n12,13.\nMethods currently in use for estimating endometrial cycle stage\nhave limitations. Endocrine methods measuring the luteinising hor-\nmone (LH) surge or peripheral blood estrogen and progesterone are\nindirect and do not allow for variability over time in endometrial\nresponse. Ultrasound scans to detect follicle size and/or ovulation do\nnot provide an obligatory correlation with endometrial development.\nRecording the commencement of last menstrual period (LMP) gives an\naccuratefix on a major endometrial event, but as a singlefixed point in\nthe cycle is of limited use for accurately comparing different stages of\ncycles of variable length. Histopathology of the endometrium is the\nmost direct measure of endometrial stage and normalcy\n14, although\nthis is a subjective method with inherent inaccuracy even among\nexperts\n15. Although significant advances have been made using endo-\nmetrial gene expression to determine cycle stage, particularly in the\nmid-luteal phase around the time of embryo implantation 4,16,17,t h e s e\nmethods do not cover the whole cycle.\nA more precise method for normalising endometrial gene\nexpression across the menstrual cycle will provide a major contribu-\ntion to understanding endometrial function and provide foundational\ninformation to investigate the pathophysiology of common gynaeco-\nlogical conditions such as heavy menstrual bleeding, recurrent\nimplantation failure and endometriosis.\nTherefore, thefirst aim of this study was to develop and validate a\nnew method for accurately determining menstrual cycle stage based\non changing endometrial gene expression. The second aim was to\ndemonstrate the functional utility of the new method by using nor-\nmalised endometrial gene expression data to identify genes that\nchange expression most rapidly across the menstrual cycle, as well as\ninvestigate the effects of increasing age and ancestry on differential\ngene expression in the endometrium. We have previously demon-\nstrated strong genetic effects on endometrial gene expression with\nsome evidence for genetic regulation of gene expression in a men-\nstrual cycle stage-speci fic manner\n18,19. However, to date no-one has\nidentified differentially expressed endometrial genes between women\nof different ancestries, despite well-established differences in genetic\nmakeup.\nResults\nSubject Details\nThe median age of all subjects (study 1, 2 and Illumina HT-12 validation\nstudy) at time of endometrial biopsy was 33 years (range 18–49). Of the\ntotal of 358 subjects, 214 had con firmed endometriosis, 131 did not\nhave endometriosis and in 13 endometriosis status was unknown.\nSimilarly, 167 had had a prior clinical pregnancy, 183 had never been\npregnant, and pregnancy status information was unavailable for the\nremaining 8.\nThe average age at time of endometrial biopsy of the 236 sub-\njects in Study 1 from which the final molecular model was developed\nwas 31.1 years (range 18 –49). All these women provided endometrial\nbiopsies at the time of diagnostic laparoscopy for suspected endo-\nmetriosis, with the primary symptom for investigation being pelvic\npain. 168 (71%) had endometriosis and 68 (29%) did not. Of the\n236 subjects, 60 (25%) had had a successful prior pregnancy (Sup-\nplementary Table 1). Fertility intention information was available for\n136 of the women (Supplementary Table 2). Only 28 women repor-\nted problems conceiving (defined as trying for more than 12 months\nto conceive), of whom 10 went on to have successful pregnancies,\nand only 6 out of 96 subjects reported pregnancy loss due to mis-\ncarriage, although this number could have been higher due to\nmissing data. All subjects reported regular menstrual cycles and\nnormal endometrium as assessed by at least one experienced\npathologist.\nAnalysis 1: Development of the ‘molecular staging model’ to\nassign cycle stage for secretory stage samples only\nSplines were fitted to RNA-seq expression data for each of 20,067\ngenes from 96 endometrial samples where 2 or 3 independent\npathology reports agreed to within 2 post-ovulatory days (Fig. 1,\npanel 1). For each endometrial sample, an estimated post-ovulatory\nday (POD) was obtained using the day which minimised mean squared\nerror (MSE) between the observed expression and the expected\nexpression across all genes. Examples of MSE plots are shown in the\nFig. 1, panel 2. There was a strong correlation between the POD cycle\ntime calculated from the lowest MSE value and the average of the\npathology estimates (r = 0.9297) (Fig. 1, panel 3). To illustrate that\nlarger, less precise, units of time can be used to estimate cycle time\nusing the same method, an additional model was built using the\npathology-assigned 3 secretory cycle stages (i.e., early-, mid-, and late-\nsecretory). The cycle time estimated from the 3 stages model showed a\nstrong correlation to the cycle time estimated from the 14-day POD\nmodel (r = 0.9807) (Fig. 1, panel 4).\nAnalysis 2: Molecular staging model using 7 pathology stages for\nthe whole cycle with RNA-seq and array expression data\nIn Analysis 2 we modelled RNA-seq expression data from all 236 sam-\nples collected in Study 1. These samples had been classified by routine\npathology into 1 of 7 cycle stages. Because the majority of the pro-\nliferative phase samples were not assigned as early, mid or late by the\npathologist, we re-assigned all samples labelled as proliferative into\nearly, mid, and late byfitting a penalised cubic regression spline (k = 3)\nusing gene expression data from samples classified by the pathologists\nas menstrual, proliferative, and early secretory (Fig. 2a). Then a pro-\nliferative time point was estimated from the minimised MSE between\nthe observed expression and the expected expression across all genes\n(Fig. 2b). The proliferative samples were then split into equal sized\ngroups of early, mid, and late using this time point (Fig.2c). A penalised\ncyclic cubic regression spline (k = 8) wasfit for all 20,067 genes using\nthe 7 stages of the menstrual cycle, which included the re-assigned\nearly, mid and late proliferative samples (Fig. 2d). Each endometrial\nsample was then assigned a‘day’ or ‘model time’ using the time which\nminimised the MSE between the observed expression data for all genes\nand their corresponding gene models (Fig. 2e). ‘Model time’ is a rela-\ntive timepoint in the cycle and does not correspond to a real day.\nUnder the assumption that all 236 women were approximately uni-\nformly distributed across the menstrual cycle, the data were trans-\nformed so that the distance in time between each sample was identical\nArticle https://doi.org/10.1038/s41467-023-41979-z\nNature Communications|         (2023) 14:6222 2\n\n(Fig. 2f). This ranked all the samples in order from the start to the end\nof the cycle, removing the need for cycle stages or an idealised 28-day\ncycle. At this point the x-axis was changed to show the percentage of\nthe way through the menstrual cycle that each sample was. The new\ntime points were also compared to the pathology-derived cycle stages\nto get an approximation how the model time corresponds to stages in\nthe menstrual cycle (Fig.2g ) .G e n ec u r v e sw e r et h e nr efitted using the\nnewly derived cycle times for each sample with a penalized cyclic cubic\nregression spline (k = 30) (Fig. 2h). For visualisation purposes, nor-\nmalisation of gene expression for cycle stage was then derived by\nsubtracting the expected expression from the observed expression\n(i.e., calculating the residuals) and re-adding the mean (Fig.2i).\nValidation of the molecular staging model\nVarious validation studies were undertaken using the molecular sta-\nging model. As an initial check, data from Analysis 1 using POD to\ndevelop the secretory model was plotted against secretory stage data\nfrom the final molecular staging model generated in Analysis 2\n(Fig.3a). A second comparison confirmed that using only 3 cycle stages\n(early, mid and late secretory from only 1 pathologist) gave similar\nresults to having more frequent daily POD information from 2 or 3\nindependent pathologists (Fig.1, panel 4). To assess the repeatability\nof the molecular staging model method, Analyses 1 and 2 were repe-\nated using Illumina HT-12 data and the results compared for the 198\nsamples that had both RNA-seq and Illumina HT-12 data (Fig.3b). There\nwas a high level of correlation in cycle stage determination using data\nfrom the 2 different gene expression platforms, with slightly more\nvariation being seen in the mid-proliferative phase. Additionally, vali-\ndation using unsupervised methods with initial groupings based on a\nPCA plot (Supplementary Fig 1) and not using pathology dating\ninformation at all also yielded a high level of concordance between the\nmolecular staging model and the validation model (Supplementary\nFig 2). The correlation between the two models was 0.989 and the\nmean absolute difference between estimated sample cycle times was\n1.67%. Peripheral blood estradiol and progesterone levels were not\nused to help determine cycle stage and could therefore be considered\nas an independent variable. Estradiol ( N =1 5 9 ) a n d p r o g e s t e r o n e\n(N = 187) data were plotted against molecular staging model cycle\nstage and showed typical expected menstrual cycle distributions\n(Fig. 3c, d).\nReanalysis of published data\nThe molecular staging model was used to re-analyse 3 published\nendometrial gene expression datasets available on GEO (GSE65099,\nendometrial samples from GSE141549, and GSE180485). These\ndatasets were chosen because they contained endometrium from\nnatural cycles with attached estimates of cycle stage and plots. We\nfirst produced a principal component analysis (PCA) plot using our\nown RNA-seq dataset ( N = 266) with cycle stage as determined by\nthe molecular staging model (Fig. 4a). This PCA plot has a char-\nacteristic pattern with samples clustering according to cycle stage\nas determined using the molecular staging model, with no outliers.\nWe then generated PCA plots from each of the 3 published datasets.\nThe PCA plot using data from GSE141549\n20 is shown in Fig. 4b, with\nsamples labelled as per information in GEO as menstrual, pro-\nliferative, secretory and unknown. There is mixing of proliferative\nsamples with menstrual and secretory ones within the PCA plot\nwhen using the GEO assigned labels. The same PCA plot using data\nfrom GSE141549 but with cycle stage assigned by our molecular\nstaging model has minimal overlap between different cycle stages\n(Fig. 4c), demonstrating that the molecular staging model accu-\nrately aligns with PCA analyses of endometrial gene expression data\nRepeat steps 1 and 2 \nusing the 3 secretory \nstages instead of POD\nValidation 2\n●\n●\n●\n●●\n●\n●\n●\n●●\n●\n●●\n● ●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\nr =0.98075\n6\n7\n51 0\nCycle Time from POD Secretory Model\nCycle Time from 3 Stages Model\nStep 2:\nEstimate post-ovulation day\nof 96 samples\n1.0\n1.5\n2.0\n1 2 3 4 5 6 7 8 9 1 01 11 21 31 4\nPost−Ovulation Day\nMean Square Error\nEstimate\nModel\nPathology\nSample X210475\nP2\n0.5\n1.0\n1.5\n123456789 1 0 1 1 1 2 1 3 1 4\nPost−Ovulation Day\nMean Square Error\nEstimate\nModel\nPathology\nSample X210470\nP2P1\n1.0\n1.5\n2.0\n123456789 1 0 1 1 1 2 1 3 1 4\nPost−Ovulation Day\nMean Square Error\nEstimate\nModel\nPathology\nSample X210244\nP2 P1\nValidation\n●●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\nr = 0.9297\n5\n10\n51 0\nAverage Post−Ovulation Day from Pathology\nPost−Ovulation Day from Model (5−Fold CV)●●\n●\n●\n●\n●\n●\n●\n●\n●● ●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●● ●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\nRadj\n2 = 0.8107\n2\n4\n6\n8\n10\n51 0\nPost−Ovulation Day\nlog2 (CPM)\nENSG00000106688\n●\n●\n●\n●\n●\n● ●\n●\n●\n●● ●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●● ●\n●\n●●\n●\n●\n● ● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●●●●●●\n●\n●\n●\n●●\n●\n●●\n●\n●\nRadj\n2 =0.65534\n5\n6\n7\n51 0\nPost−Ovulation Day\nlog2 (CPM)\nENSG00000188313\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n● ●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●●●●●\n●●●\n●\nRadj\n2 =0.1728\n−4\n−2\n0\n2\n4\n51 0\nPost−Ovulation Day\nlog2 (CPM)\nENSG00000163464\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n● ●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n● ●\n●\nRadj\n2 = 0.0977\n−0.5\n0.0\n0.5\n1.0\n1.5\n51 0\nPost−Ovulation Day\nlog2 (CPM)\nENSG00000260852\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ● ●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●● ●●●●●\nRadj\n2 = 0.309\n2.8\n3.2\n3.6\n4.0\n4.4\n51 0\nPost−Ovulation Day\nlog2 (CPM)\nENSG00000103047\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●●\n●\n● ●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\nRadj\n2 = 0.482\n1\n2\n3\n4\n5\n51 0\nPost−Ovulation Day\nlog2 (CPM)\nENSG00000167984\nStep 1:\nModel gene expression of 20,067 genes\n1. 2. 3.\n4.\nP1\nFig. 1 | Analysis 1: Development of the molecular staging model to assign cycle\nstage for secretory stage samples only. Panel 1. Examples of regression spline\nfitting to expression data for individual genes from 96 endometrial samples taken\nbetween post-ovulatory days (POD) 1–14. Splines were fitted to a total of 20,067\ngenes. The error bands surrounding the regression lines represent the 95% con-\nfidence interval. Panel 2. Plots showing post-ovulatory time that gives lowest Mean\nSquared Error (MSE) using spline data for all 20,067 probes for 3 different\nendometrial samples (solid line). Dotted lines show POD estimates from 2 inde-\npendent evaluations by experienced pathologists. Panel 3. Correlation between\naverage POD from 2 or 3 independent pathology evaluations and the POD time at\nwhich the lowest MSE occurred. Panel 4. Correlation between the estimated POD\nfrom the POD secretory model and the estimated cycle time from the 3 stages\nsecretory model. Source data are provided as a Source Data file.\nArticle https://doi.org/10.1038/s41467-023-41979-z\nNature Communications|         (2023) 14:6222 3\n\nacross the menstrual cycle. In a similar fashion, but with a smaller\ndataset from GSE65099 21, samples reported as LH + 6 to LH + 10 do\nnot group in a consistent fashion by PCA (Fig. 4d). When the same\nsamples are assigned cycle stage times by the molecular staging\nmodel, the same PCA analysis shows consistent grouping according\nto cycle stage for all samples, with 2 outlying samples on the PCA\nplot being reassigned as proliferative and not secretory (Fig. 4e). We\nrepeated the PCA comparison approach using 36 RNA-seq endo-\nmetrial data sets deposited in GSE180485 (Supplementary Fig 3).\nThese samples are from a study called EndoTime to determine\nwhether the accuracy of timing of luteal phase endometrial biopsies\nbased on urinary ovulation testing could be improved by measuring\nthe expression of a small number of genes using a continuous, non-\ncategorical modelling approach\n17. We produced identical PCA plots\nwith the first showing samples labelled with percentage cycle times\nderived from our molecular model, and the second using the\nEndoTime model eLH+ estimates. While our percentage data shows\na steady progression through the cycle in concordance with the PCA\nplot, the eLH+ days 4-5 data points are widely spread (i.e., these\nsamples have signi ficantly different molecular pro files despite the\nEndoTime model calling them as similar), and several of the eLH+\ndays 8 –10 samples are grouped together suggesting that this ver-\nsion of the EndoTime model cannot reliably discriminate between\nthese post LH surge days.\nChanges in Endometrial Gene Expression with Increasing Age\nand Different Ancestries\nUsing our RNA-seq data ( n = 266, 20,067 genes analysed) with\nmenstrual cycle staging calculated using the molecular staging\nmodel, a total of 60 endometrial genes showed signi ficant changes\nin expression with increasing age. Examples of 2 signi ficant genes\nare shown in Fig. 5a. Re-running the age analysis using the original 7\ncycle stage pathology data instead of the staging from the mole-\ncular staging model reduced the number of age-related signi ficant\ngenes from 60 to 32, providing evidence that the molecular staging\nmodel provides a superior approach for identifying differentially\nexpressed genes. To further explore the effects of aging on endo-\nmetrial gene expression, an additional n = 87 Illumina microarray\nendometrial samples from GSE141549 were analysed and combined\nwith our RNA-seq differential expression results as part of a meta-\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●●●●●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n● ●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●● ●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●\n●\nRadj\n2 = 0.7396\n4\n5\n6\n7\n8\n9\n12345\nPathology Cycle Stage\nlog2 (CPM)\n●●\n●●\n●\n●\n●\n●\n●●\n●\n● ●● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●\n●\n● ●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●● ● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n● ●●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n● ●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n1\n2\n3\n4\n5\n12345\nPathology Cycle Stage\nModel Estimate of Cycle Stage\n●● ●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●●\n●\n●\n●\n●●\n● ●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n1\n2\n3\n4\n5\n12345\nReclassified Pathology Cycle Stage\nModel Estimate of Cycle Stage\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●● ●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●●● ●\n●●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●\n●●\n●●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n● ●\n●●●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●●\n●●●\n●●\n●\n●\n●\n●\n●\n●●●●●● ●\n●\n●●\n●\n●\n●●●●●●●●●●●●\n●●●●●●●●●●●●●●●●●●●\n●●\nRadj\n2 = 0.80580.0\n2.5\n5.0\n7.5\n10.0\n1234567\nReclassified Pathology Cycle Stage\nlog2 (CPM)\n●\n●\n●●\n● ●\n●\n●\n●●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n● ● ●●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●\n● ●\n●\n●\n●\n● ●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n● ●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n1\n2\n3\n4\n5\n6\n7\n1234567\nReclassified Pathology Cycle Stage\nModel Estimate of Cycle Stage ●\n●\n●●\n●●\n●\n●\n● ●\n●\n●\n●● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●●●● ●\n●\n●\n●\n●\n●●\n●\n●\n●\n● ●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●\n●●\n●\n●\n●\n● ●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n● ●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n1\n2\n3\n4\n5\n6\n7\n02 5 5 0 7 5 1 0 0\nModel Cycle Time\nModel Estimate of Cycle Stage\n● ●●● ●\n● ● ●● ●● ●●● ● ●● ●●● ●●● ●●● ●● ● ●●● ●● ●● ●●\n●●\n● ●● ●● ● ● ●\n●●● ●● ● ●●\n●●●● ●● ●\n● ●●\n● ● ● ●● ●●● ●●● ● ●● ● ●●●\n●● ●● ●●\n●●● ● ●● ●\n●●● ●● ●●●\n●\n●●● ●● ● ●●● ●● ●●● ●\n● ●\n●● ●●\n● ●● ●●\n● ● ●●\n●\n●●\n● ● ●● ●●● ●●●\n● ●●\n●● ●●\n●\n●●●\n●\n●● ●●\n●\n●\n●\n●\n●\n● ●● ●\n●\n●●● ●●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n● ●\n●\n●\n●●● ●\n●\n●\n●\n●● ●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●● ●●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n7\n6\n5\n4\n3\n2\n1\n0 1 02 03 04 05 06 07 08 09 0 1 0 0\nModel Cycle Time\nPathology Cycle Stage\nsecretory\n●\n●\n●\n● ●●\n● ●\n● ●\n●\n●\n●● ●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●●\n●\n●\n●● ●● ●\n●\n●\n● ●\n●\n● ●\n●\n●\n● ●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●● ●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●●\n●\n●\n●●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n● ●\n●●\n●\n●\n●\n●●\n●\n●\n●\n●●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●●\n●\n●\n●\n●\n● ●\n●\n●\n● ●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●●● ●\n●\n●\n● ●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●●●●●\n●●●●\n● ●●\n●●\n●●\n●\n●\n●\n●\n●\n● ●\n●\n●\nRadj\n2 = 0.8944\n5\n6\n7\n8\n9\n02 5 5 0 7 5 1 0 0\nModel Cycle Time\nlog2(CPM)\n●\n●\n●\n●\n●●\n●\n●● ●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●●\n●\n●\n●\n● ●● ●\n●\n●●\n●\n●\n●\n●\n●\n●\n● ●\n● ●\n●\n●\n●\n●\n●\n● ●\n● ●\n●\n●\n●\n●● ●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●● ●●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●●\n●\n●●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n● ●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n● ●\n●● ●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n5\n6\n7\n8\n9\n02 5 5 0 7 5 1 0 0\nModel Cycle Time\nCycle Normalised log2(CPM)\na\nd e f\nb c \ng ih\nmenstrual proliferative\nFig. 2 | Analysis 2: Development of the molecular staging model to assign cycle\nstage for the whole menstrual cycle.Figure 2a. Example of a spline curvefitted to\nmenstrual, combined proliferative and early secretory expression data. The error\nbands surrounding the regression line represent the 95% confidence interval.\nFigure 2b. Menstrual, proliferative and early secretory samples with estimated\nmolecular cycle stage calculated from minimum mean squared error data. Figure\n2c. Combined proliferative samples reassigned into early, mid and late proliferative\ngroups containing equal numbers. Figure 2d. Example of a spline curve fitted to\ndata from all 7 stages of the cycle using reclassified proliferative cycle stage\ninformation and pathology-derived menstrual and secretory staging. The error\nbands surrounding the regression line represent the 95% confidence interval.\nFigure 2e. Comparison of cycle time within the 7 cycle stages estimated from the\nmodel and pathology-derived cycle stage. Figure 2f. Under the assumption that 236\nwomen underwent surgery at random stages of the menstrual cycle, data from\n236 samples were transformed to be uniformly spaced along the x-axis on a\n0–100 scale. This transformation allows each sample to be identified as being a\npercentage of the way through the menstrual cycle. Figure 2g. Plot showing rela-\ntionship between pathology staging and the percentage of cycle from the mole-\ncular staging model. Menstrual is 0–8%, proliferative is 8–58% and secretory is\n58–100% of the molecular staging model cycle respectively. Figure 2h, i. Expression\ndata for ENSG00000187231 replotted using derived‘percentage’ cycle times and\nthen normalised across the menstrual cycle. The error bands surrounding the\nregression line represent the 95% confidence interval. Source data are provided as a\nSource Data file.\nArticle https://doi.org/10.1038/s41467-023-41979-z\nNature Communications|         (2023) 14:6222 4\n\nanalysis. Considering only genes in both datasets, this reduced the\nnumber of genes tested to 12,868, which still included 32 of the\n60 signi ficant genes from our original dataset. 65 signi ficant genes\nwere found in the GSE141549 data when analysed on its own. How-\never, when the 2 data sets were combined ( n = 353), 206 signi ficant\ngenes were identi fied across the whole menstrual cycle (Supple-\nmentary Table 3). We then split the samples into 3 cycle stages;\nmenstrual, proliferative and secretory (equivalent to 0 –8%, 8 –58%\nand 58 –100% of the molecular staging model cycle respectively)\nand analysed each stage of the cycle separately (Fig. 5b). Of note,\nnearly all (218/222 [98%]) of the genes showing signi ficant changes\nwith age were found in samples taken in the secretory phase of the\ncycle (Supplementary Data 1). A gene ontology enrichment analysis\nwas run using the 218 genes from secretory samples that changed\nsignificantly with age (Supplementary Data 2). The top biological\nprocesses enriched with upregulated genes were related to axo-\nnemes, cilia and microtubules while the top processes enriched with\ndownregulated genes were related to blood vessels, endothelial\ncells and angiogenesis.\nAncestry of subjects as defined by a previous study\n19 was used to\nlook for differences in endometrial gene expression using pairwise\ncomparisons of each ancestry group. In our Australian population the\nmajority of subjects were of European ancestry, however, despite small\nnumbers in other groups signi ficant differences in gene expression\nwere identified between the groups (Supplementary Fig 4, Supple-\nmentary Table 4 & Supplementary Data 3).\n0\n25\n50\n75\n100\nCycle Time\n●\n●\n●\n●\n●\n●\n●\n●\n● ●●\n●●●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●● ●\n●●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●● ●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●\n●\n● ●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n● ●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n● ●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●● ●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n−100\n0\n100\n−100 0 100\nPC1\nPC2\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n−60\n−30\n0\n30\nPC2\n−50 −25 0 25\nPC1\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n−50\n0\n50\nPC2\n−50 0 50 100\nPC1\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n−60\n−30\n0\n30\nPC2\n−50 −25 0 25\nPC1\n0\n25\n50\n75\n100\nCycle Time\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n−50\n0\n50\nPC2\n−50 0 50 100\nPC1\n0\n25\n50\n75\n100\nCycle Time\nPROLIFERATIVE\nMENSTRUAL\nSECRETO\nR\nY\na b c\nd e\nCycle Stage\n●\n●\n●\n●\nMenstruation\nProliferative\nSecretory\nUnknown\nLH+ Day\n●\n●\n●\n●\n●\nLH+6\nLH+7\nLH+8\nLH+9\nLH+10\nFig. 4 | Reanalysis of published data. Figure 4a. PCA plot using 266 endometrial\nsamples coloured to show the molecular staging model cycle stage. The PCA plot\nhas a characteristic pattern with the samples in cyclic pattern corresponding to\napproximate cycle stage order. Figure 4b. PCA plot using Illumina microarray data\nfrom GSE141549\n20, with samples identified as menstrual, proliferative, or secretory\nby the authors. Note significant mixing of proliferative samples with menstrual and\nsecretory ones. Figure 4c. The same PCA plot using data from GSE141549 with cycle\nstage assigned by the molecular staging model. Note minimal overlap between\ndifferent molecular staging model cycle stages. Figure 4d. PCA plot using RNA-seq\ndata from GSE65099\n21 with samples identified as LH + 6 to LH + 10 as reported in the\nstudy. Figure 4e. The same PCA plot using data from GSE65099 with cycle stage\nassigned by the molecular staging model. Note reassignment of 2 outlying samples\non the PCA plot as proliferative and not secretory. Source data are provided as a\nSource Data file.\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\nr = 0.9790\n5\n10Cycle Time from POD Secretory Model\n50 60 70 80 90 100\nModel Cycle Time\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\nr = 0.9744\n0\n25\n50\n75\n100Cycle Time from Microarray Model\n02 5 5 0 7 5 1 0 0\nCycle Time from RNA−Seq Model\n●●●\n●\n●● ●●● ●●●●● ●●●● ●●●●●● ● ●● ●● ● ●●● ●● ● ●● ●●● ●● ●● ●●● ●● ●●●● ● ●● ● ●● ●●●● ●● ●● ● ●● ●● ●●● ●\n●●●● ●\n●●●●●●\n●●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●●\n●\n●\n●\n●\n● ●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n● ●\n●\n●\n●●● ●\n●\n● ●\n●\n●●\n● ●●● ●● ● ●●\n●\n●●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n0\n5\n10\n15\n20\n25\n02 5 5 0 7 5 1 0 0\nModel Cycle Time\nP4 ng/mL\n●●●● ●\n●\n●●●\n●\n●●●\n●\n●\n●●● ●\n●\n●\n●\n●● ●\n●\n●\n●● ●\n●\n●●\n●\n● ●\n●● ●\n●\n● ●\n●\n●● ●\n●\n●\n●\n●\n●\n●\n●\n● ● ●●\n● ●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●● ●\n●\n●\n●\n● ●\n●\n●\n●\n● ●\n●\n●●●\n●●●●\n●\n●\n●\n●\n●\n●\n●●\n● ●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●●● ●●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n● ●\n●\n●●● ●●\n●\n●\n●\n● ●●\n0\n100\n200\n300\n02 5 5 0 7 5 1 0 0\nModel Cycle Time\nEstradiol pg/mL\na b c\nd\nFig. 3 | Validation of the NGS Molecular Staging Model. Figure 3a. For endo-\nmetrial secretory samples, POD predicted from the secretory molecular staging\nmodel was compared against the full molecular staging model. Note, POD 1 is\napproximately 58% of the way through the cycle. Figure 3b. Illumina HT-12 micro-\narray data were also available from 198 endometrial samples used to generate the\nmolecular staging model. A validation study was run comparing NGS vs Illumina\ndata. The molecular staging model results showed a strong correlation when\ncomparing the 2 different gene expression platforms. Figure 3c, d. Peripheral blood\nprogesterone (n = 187) and estrogen (n = 159) data plotted against the molecular\nstaging model cycle stage showed expected typical menstrual cycle-stage dis-\ntribution. Source data are provided as a Source Data file.\nArticle https://doi.org/10.1038/s41467-023-41979-z\nNature Communications|         (2023) 14:6222 5\n\nDifferential gene expression across the cycle using the mole-\ncular staging model\nTo investigate changing gene expression across the cycle, all samples\nwith RNA-seq data from Study 1 (n = 236) were ranked in chronological\norder from start to end of the molecular staging model cycle. A‘sliding\nwindow’ approach was then used to compare differential gene\nexpression (DGE) between samples 1–8 versus samples 9–16, followed\nby samples 2 –9v e r s u s1 0–17, then 3 –10 vs 11 –18 and so on for all\nsamples across the menstrual cycle. Group sizes were arbitrarily set at\n8 because this represents 3.4% of the 236 samples or approximately\n1 day assuming a mean cycle length of 28 days. Moderated t-tests were\nused to identify differentially expressed genes withP < 0.05 following\nmultiple testing correction, at each window. Using adjusted P values,\n488 unique genes significantly changed expression during menstrua-\ntion, 44 during the proliferative phase, and 2921 during the secretory\nphase. Peak times of rapid change in gene expression occurred during\nmenstruation (3% of the way through the cycle), late proliferative\n(51%), POD3 (66%), POD5 (71%), POD11 (94%) and POD13 (98%) (Fig.6a).\nExamples of 12 endometrial genes showing significant and very rapid\nchanges in expression across different stages of the menstrual cycle\nare provided in Fig. 6b.\nThe original Endometrial Receptivity Analysis (ERA) publication\nidentified 238 genes that show major changes in expression before,\nduring and after the time of embryo implantation at POD 3 –7\n22.O f\nthese 238 genes, 207 were recognised in our NGS data, and 70% of\nthese (145/207) changed expression significantly between cycle times\n66 ± 2% and 76 ± 2% (POD 3–7). Supplementary Fig 5 shows the 6 most\nsignificantly down-regulated genes and the 6 most signi ficantly up-\nregulated ERA genes that we identified.\nDiscussion\nWe have developed and validated a method for accurately determining\nendometrial cycle stage based on global gene expression. We did this\nby generating mathematically de fined curves fitted to RNA-seq\nexpression data from 236 endometrial samples for each of 20,067\ngenes. From these curves we then found the time of the cycle that\nminimised a loss function (mean squared error) which gives the bestfit\nfor all genes simultaneously for any individual endometrial sample. By\nplacing the 236 samples in chronological order from the start to end of\nthe menstrual cycle and converting the x-axis to percentage, we cre-\nated a method for both de fining how far through the cycle any given\nsample is, as well as being able to normalise gene expression for cycle\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n● ● ●\n● ●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●●● ●●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n● ●\n●\n● ●\n●\n●\n● ●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●● ●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●● ●\n● ●\n●\n●\n●\n●\n●\n●\n● ● ●●●●\n● ●\n●\n●\n●● ●\n● ●●●● ●●●●● ●●●●● ●●\n● ●●●●●●●●●●\npadj < 0.001\n4\n5\n6\n20 30 40 50\nAge\nNormalised log2(CPM)\nENSG00000180543\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n● ●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●●\n●\n●\n●\n●●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n● ●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●●\n●\n●\n●\n●● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n● ●●\n●● ●\n●\n●\n●\n●\n● ●\n●●●●●●●●\n●●●●●●●●●●\n●●\n●●●●\n●●●●●●●●\npadj < 0.001\n−4\n0\n4\n8\n12\n20 30 40 50\nAge\nNormalised log2(CPM)\nENSG00000129988\n●\n●\n●\n●●\n●● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●● ●● ●\npadj = 1\n●\n● ●\n●\n●\n●\n●\n●●\n●\n●●●●\n●\n●\n●\n●\n●\n●\n● ●●\n●\n●\n●\n●\n● ●\n●\n● ●\n●●\n●\n●\n●\n● ●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●● ●\n●\n●\n●\n●●\n●\n●\n● ●\n●\n● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●●\n●\n● ●\n●\n●\n●\n● ●● ●●● ●●●●●● ●●●●●●●\n● ●●● ●●\n●●●\n●●\n●●\n●\npadj = 1\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●\n● ●\n●\n●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●●●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●● ●\n●\n●●\n● ●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●\n● ●\n●\n●\n●\n●\n● ●\n●● ●\n●\n●\n●●●\n●\n●\n●\n●\n●\n● ●\n● ●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n● ●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●●●\n●● ● ●●●● ●\npadj = 0.010\nMenstrual Proliferative Secretory\n20 30 40 50 20 30 40 50 20 30 40 50\n3\n4\n5\nAge\nNormalised log2(CPM)\nENSG00000147113\n● ●\n●\n●●\n●\n● ●● ●● ●\n●\n● ●\n●\n●\n●\n● ●\n●● ●●●●●●● ●● ●●\npadj = 1\n●\n● ●● ●\n●\n● ●●\n●\n●●\n●\n● ●\n●\n●\n●\n● ●\n●\n●● ●\n● ● ●● ●\n●\n●\n● ●\n●●\n●\n●●\n●●●\n●●\n●\n● ●●\n● ●\n● ●\n●\n●\n●\n● ●\n● ●\n●\n● ●\n●●\n●\n●\n●\n●\n●\n● ●●\n●\n●\n●\n●\n●\n●●\n●\n● ●\n● ●\n●● ●\n●\n●\n●\n●\n● ●\n●\n●\n●\n● ●\n●\n●\n●●\n●\n●●\n● ●●● ●●\n●● ●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●●●●●●●●●● ●● ●● ●● ●●●●● ●●●● ●●●●●●●●●●●● ●●●● ●●●●\n●\n●\n●●●●●●●● ●●● ●●\npadj = 1\n●\n●\n●\n●\n●\n● ●●●\n●\n●\n●\n●\n●\n●\n●\n●\n● ●● ● ●\n● ●\n●\n●\n●\n●\n●●\n●\n●\n●●●\n●\n●●\n●\n●\n● ●\n●\n●●\n●\n●\n●\n●●\n●●\n●●\n●\n●\n●\n●\n●● ●\n●\n● ● ●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n● ● ● ●\n●\n●\n●\n●●\n●●\n● ●● ● ●●\n●\n●\n● ●\n●●●\n● ●●\n●\n●\n●\n●● ●\n●\n●\n●\n●\n●● ●● ●\n● ●●●●●●●● ●\n● ●●●●●●●●●● ●●●●●●● ●●●●●● ● ●●●●●●\n●\n●●●●●●●●●● ●●●● ●●●●●●●●●\npadj = 0.013\nMenstrual Proliferative Secretory\n20 30 40 50 20 30 40 50 20 30 40 50\n−8\n−4\n0\n4\nAge\nNormalised log2(CPM)\nENSG00000174175\na\nb\nFig. 5 | The impact of age on endometrial gene expression. Figure 5a. Two\nexamples from the 60 endometrial genes that change expression significantly with\nincreasing age. Statistical analysis used two-sided empirical Bayes moderated\nt-tests implemented in limma35 with corrections for multiple comparisons per-\nformed using the Benjamini-Hochberg method (ENSG00000180543: t = 6.025,\np =5 . 8 2 e - 9 ,p\nadj = 5.84e-5; ENSG00000164778: t =−5.886, p =1 . 2 9 e - 8 ,padj =8 . 2 2 e -\n5). Expression data were plotted following normalisation for changing expression\nacross the menstrual cycle (n = 266 RNA-seq samples, 87 samples from\nGSE14154920). Figure 5b. Expression data for 2 genes plotted separately for men-\nstrual, proliferative, and secretory samples (n = 266). Statistical analysis used two-\nsided empirical Bayes moderated t-tests implemented in limma with corrections for\nmultiple comparisons performed using the Benjamini-Hochberg method (secre-\ntory phase ENSG00000147113: t =−5.074, p =1 . 4 3 e - 6 ,padj = 0.010; secretory phase\nENSG00000174175: t =−4.429, p =2 . 0 9 e - 5 ,padj = 0.013). Source data are provided\nas a Source Data file.\nArticle https://doi.org/10.1038/s41467-023-41979-z\nNature Communications|         (2023) 14:6222 6\n\nstage so that gene expression differences between any 2 samples can\nbe compared.\nOur work has immediate translational application. Being able to\naccurately identify differentially expressed genes from relatively small\nnumbers of samples will substantially increase the prospect of identi-\nfying genes linked to different endometrial pathologies, thus helping\nto identify biomarkers for diagnosis and potential targets for therapy.\nSimilarly, detailed knowledge of normal gene expression pro files\nthroughout the endometrial cycle will assist in identifying ‘off-target’\neffects during clinical trials of new therapeutics. Accurate cycle staging\nis also critical for synchronisation of frozen-thaw embryo replacement\nduring in vitro fertilization treatment. Longer-term, improved under-\nstanding of endometrial disorders such as heavy menstrual bleeding,\nendometriosis and adenomyosis will lead to improved outcomes for\nmany women.\nOur molecular model for staging the endometrium creates many\nnew research opportunities in benign gynaecology. There is a com-\npelling argument for more research in thisfield given the considerable\nfinancial and quality of life burdens that endometrial-related disorders\nhave on patients and their families. The total direct and indirect costs\nof abnormal uterine bleeding in the UK and USA were reported more\nthan 10 years ago as greater than £1.2 billion and $37 billion,\nrespectively\n1. More recently in the USA, heavy menstrual bleeding\n(HMB) with uterine fibroids had a mean all-cause total cost of $16,762\nper patient, while HMB only was $11,135 compared to a control cohort\nat $6,691\n23. In 2011, in the UK, it was estimated that in England and\nWales approximately 1.5 million women were affected by HMB\nannually. HMB is the fourth most common reason for referral to\ngynaecology services in the UK and approximately 20% of 1.2 million\nreferrals to specialist gynaecology services concern women with HMB\nin the UK\n1.I na d d i t i o nt ot h e financial impact, women with\nmenstruation-related symptoms have lower scores on several domains\nof quality of life such as general health and physical, mental, social and\noccupational functioning during their periods\n24. A recent review of the\neconomic burden of endometriosis determined the overall direct\nmedical cost range to be from US$1459 to US$20,239 per patient per\nyear, and the indirect costs between US$4572 and US$14,079\n2.\nEndometrial cycle stage classi fication is currently based on his-\ntopathological criteriafirst published over 70 years ago14.T h eo n g o i n g\nneed for endometrial cycle staging is evidenced by the fact that this\n1950 paper has received over 3,500 citations25, the most for any study\nin obstetrics and gynaecology. Accuracy of pathology dating has well\nestablished limitations. It has been reported that accuracy in dating the\nsame POD between different pathologists using Noyes criteria was\n0\n250\n500\n750\n1000\n0 1 02 03 04 05 06 07 08 09 0 1 0 0\nModel Time\nNumber of Significant Genes\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●●\n●\n●●\n●\n●●\n●●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●●\n●●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●●●●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●●●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●●●\n●●●●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●\n−5\n0\n5\n0 2 04 06 08 0 1 0 0\nModel Cycle Time\nlog2 (CPM)\nleft−right determination factor 1\nLEFTY1\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●●\n●\n●●●●●●●●●\n●\n●●●●●●\n●\n●●●●●●●●●●●●●●●●●●●●●●\n●\n●●●\n●\n●●●●\n●\n●●●\n●\n●●●\n●\n●●●\n●\n●●●●\n●\n●●●\n●\n●●●●●●●●●●\n●\n●●●\n●●\n●\n●\n●●●\n●\n●●●●\n●●\n●●●●●●●●●●●●●●●●●●●●●●●\n●\n●●●●●●●●●●●●\n●\n●●●●\n●●\n●●●●●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n−5.0\n−2.5\n0.0\n2.5\n02 0 4 0 6 0 8 0 1 0 0\nModel Cycle Time\nlog2 (CPM)\nlong intergenic non−protein coding RNA 2274\nLINC02274\n●\n●●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●●●\n●●\n●●●\n●\n●\n●●\n●●●\n●●\n●●\n●●●●\n●\n●\n●●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●●\n●●\n●●\n●●\n●●\n●\n●●●●●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n2.5\n5.0\n7.5\n10.0\n02 0 4 0 6 0 8 0 1 0 0\nModel Cycle Time\nlog2 (CPM)\nDEPP1 autophagy regulator\nDEPP1\n●\n●\n●●\n●●\n●\n●\n●\n●\n●●\n●●\n●●●●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●●\n●●\n●●●\n●\n●\n●\n●\n●\n●\n●●●●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●●\n●●●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●●●\n●\n●●●\n●●●\n●\n●\n●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●\n−5.0\n−2.5\n0.0\n2.5\n5.0\n7.5\n02 0 4 0 6 0 8 0 1 0 0\nModel Cycle Time\nlog2 (CPM)\ncytochrome P450 family 3 subfamily A member 5\nCYP3A5\n●●●●●●●●●●●●●●\n●\n●●\n●\n●\n●●●●●●●●●●●●●●●●●●\n●\n●\n●●●●●●●●\n●\n●●●●\n●\n●\n●●\n●●●●●●\n●\n●●●●●●●●●●●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●●●●\n●\n●●●\n●\n●\n●\n●\n●\n●●●\n●\n●●●●●●\n●●\n●\n●\n●\n●●●●●●●●●●●●●●●●●●●●●●●●\n●\n●\n●\n●\n−5.0\n−2.5\n0.0\n2.5\n5.0\n02 0 4 0 6 0 8 0 1 0 0\nModel Cycle Time\nlog2 (CPM)\n3−hydroxy−3−methylglutaryl−CoA synthase 2\nHMGCS2\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●●\n●\n●\n●●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●●●\n●●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●●\n●●\n●●\n●\n●●\n●●●●\n●\n●●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●●\n●\n●\n●●●●●\n●\n●\n●●\n●●●\n●●●●\n●●\n●\n●\n●\n●\n●●\n●●\n●●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●●\n●●●\n●●\n●●\n●\n●\n●\n●●●●\n●●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n−5\n0\n5\n02 0 4 0 6 0 8 0 1 0 0\nModel Cycle Time\nlog2 (CPM)\nphospholipase A2 group IVF\nPLA2G4F\n●\n●●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●●●\n●●●\n●●●●●●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n−4\n−2\n0\n2\n4\n02 0 4 0 6 0 8 0 1 0 0\nModel Cycle Time\nlog2 (CPM)\ncoagulation factor II, thrombin\nF2\n●●●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●●\n●\n●●●\n●●\n●\n●●\n●●\n●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●●●●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n−2.5\n0.0\n2.5\n5.0\n02 0 4 0 6 0 8 0 1 0 0\nModel Cycle Time\nlog2 (CPM)\nsolute carrier family 24 member 4\nSLC24A4\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●\n●●●●\n●\n●●●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●●●\n●\n●●\n●\n●\n●●●\n●\n●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●●●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●\n−5.0\n−2.5\n0.0\n2.5\n5.0\n02 0 4 0 6 0 8 0 1 0 0\nModel Cycle Time\nlog2 (CPM)\nserpin family A member 4\nSERPINA4\n●\n●●●●\n●\n●\n●\n●\n●●\n●●\n●\n●●\n●\n●\n●\n●●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●●\n●●●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●\n●●\n●\n●\n●●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●●●\n●\n●\n●●●●●●●●\n●\n●●\n●\n●●●\n●\n●\n●\n●\n●\n●●\n●●●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●●\n●\n●●●●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●●●●\n●\n●\n●●\n●\n●●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●●●●\n●\n●\n4\n6\n8\n02 0 4 0 6 0 8 0 1 0 0\nModel Cycle Time\nlog2 (CPM)\nOCIA domain containing 2\nOCIAD2\n●\n●\n●\n●\n●\n●\n●\n●●●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●●●\n●●●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●●●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●●●\n●●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●●●●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n−5.0\n−2.5\n0.0\n2.5\n02 0 4 0 6 0 8 0 1 0 0\nModel Cycle Time\nlog2 (CPM)\nshisa family member 3\nSHISA3\n●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●\n●●\n●●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●●\n●●\n●\n●\n●●\n●\n●\n●●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●●●●●\n●\n●\n●\n●●\n●●\n●\n●\n●●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●\n●\n●\n●\n●\n●\n●\n●\n●●●●\n●\n●\n●\n●●\n●\n●\n●●●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●\n●●\n●\n●●●\n●\n●\n●●●●\n●\n●\n●\n●\n●\n−3\n0\n3\n6\n02 0 4 0 6 0 8 0 1 0 0\nModel Cycle Time\nlog2 (CPM)\nLHFPL tetraspan subfamily member 3\nLHFPL3\na\nb\nFig. 6 | Changing gene expression across the cycle. Figure 6a. Genes from RNA-\nseq analysis that significantly change expression (Padj <0.05) over 3.4% of the cycle\n(approximately equal to a 24-hour window) at different stages of the menstrual\ncycle. Statistical analysis used two-sided empirical Bayes moderated t-tests, tested\nrelative to a threshold fold change of 1.2, implemented in limma\n35 with corrections\nfor multiple comparisons performed using the Benjamini-Hochberg method. Using\nadjusted P values, 488 unique genes significantly change expression during\nmenstruation, 44 during the proliferative phase, and 2921 during the secretory\nphase. Peak times of rapid change in gene expression approximately correspond to\nmenstrual (3% of the way through the cycle), late proliferative (51%), POD3 (66%),\nPOD5 (71%), POD11 (94%) and POD13 (98%). Figure 6b. Examples of 12 genes that\nchange expression significantly at different times across the cycle. Source data are\nprovided as a Source Data file.\nArticle https://doi.org/10.1038/s41467-023-41979-z\nNature Communications|         (2023) 14:6222 7\n\npoor, ranging from 18 to 40%, rising to 60–68% for POD dating within\n1 day and to 65–81% for dating within 2 days15. The same study reported\nthat accuracy and interobserver reproducibility were unaffected by\nrefresher training, suggesting the limits of histological dating have\nbeen reached. By contrast, our model reliably identi fies signi ficant\ndifferences in gene expression from day to day through the cycle (see\nFig. 6)a n dc a na s s i g nat i m ep o i n to nas c a l em o r er efined than a single\nday. In addition, our model accurately determines endometrial stage\nregardless of cycle length and does not rely on measurement of\nexternal events such as changes in circulating hormones orfixed time\npoints such as ultrasound detection of ovulation or first day of men-\nstrual bleeding. How accurate these results are in terms of cycle stage\nis more difficult to ascertain. The comparison of NGS versus Illumina\ndata from the same samples (Fig. 3b) shows that at times of the cycle\nwhere gene expression is changing rapidly, such as menstruation and\nimplantation (see Fig.6a), there is very high concordance between the\n2 different methods. This agreement suggests that results at these\nstages of the cycle are accurate to within less than 3-4% or approxi-\nmately 1 day. At other times, such as during most of the proliferative\nphase, there is less agreement, suggesting an accuracy of 7–11% or 2-3\ndays. By contrast, pathology dating during the proliferative phase is\nlimited to early, mid or late. The use of percentage to describe distance\nthrough the cycle removes dif ficulties in comparing different length\ncycles. The human menstrual cycle is inherently variable in overall\nlength and in the relative lengths of the menstrual, proliferative and\nsecretory phases\n10–13. This variability makes accurate comparisons\nbetween endometrial samples challenging. For example, 10 days post\nlast menstrual period (LMP) can be proliferative endometrium in one\nwoman and secretory in another. Similarly, POD6 can be pre-receptive\nin one case and post-receptive in another. By de fining cycle stage in\nterms of the percentage of the way through the menstrual cycle that\nthe endometrium is (which still allows precise mapping to events such\nas menstruation and ovulation), it is possible to accurately synchronise\nall endometrial samples to their own well-de fined point in the cycle.\nFrom our data, menstrual is 0 –8% of the cycle, proliferative is 8 –58%\nand secretory is 58–100%.\nAnother advantage of the molecular staging model is that it\ngenerates mathematically de fined curves (splines) for expression\nlevels of all genes across the cycle. This allows gene expression from\ndifferent stages of the cycle to be readily normalised for comparison\nstudies. It also provides information on expression profiles for every\ngene and can be used across different gene array platforms. As part of\nthe validation studies undertaken for the molecular staging model, we\nused the same methods to derive a similar staging model using Illu-\nmina HT-12 data from 198 samples that also had RNA-seq data. This\ne x e r c i s ed e m o n s t r a t e dah i g hd e g r e eo fa g r e e m e n ti nt h et w om o d e l s’\ncycle stage between the RNA-seq and Illumina HT-12 data sets, con-\nfirming that the molecular staging model method gives robust out-\ncomes when used with different transcriptomics platforms.\nEndometrial expression data in GEO and other databases can now\nbe precisely staged and reinterp reted using this methodology. To\nsupport this claim, we have reinterpreted 3 published data sets and\ngenerated PCA plots to demonstrate that the molecular model pro-\nvides accurate timing for all samples in each data set. We show good\nagreement between the molecular staging model and most of the\npublished staging but identify several samples in each published data\nset that are misclassi fied. With relatively small sample sizes, a few\nmisclassified samples can severely reduce the ability to identify dif-\nferentially expressed genes. Based on these results, we suggest the\nneed for caution in interpreting endometrial gene expression studies\nwith small sample sizes and/or where the authors have not been able to\nincorporate accurate cycle staging into their analyses.\nThis day-to-day volatility in gene expression during the secretory\nphase provides a compelling explanation for the consistent lack of\nagreement between studies comparing gene expression where\ncomparison groups contain all stages of the secretory phase com-\nbined. Use of the molecular staging model to accurately assign cycle\nstage and account for these rapid changes in gene expression data will\nsignificantly enhance our ability to compare endometrial gene\nexpression between experimental groups and identify real differences\ndue to phenotype or treatment.\nThe molecular staging model reveals remarkably synchronised\nand dramatic changes in expression for many genes at different stages\nof the cycle. The model delivers highly accurate cycle staging because\nmost of the variance in endometrial gene expression is due to cycle\nstage, and not other factors. Over 50% of the gene models had an\nadjusted R-squared value greater than 0.19 and 30% were over 0.36,\nindicating that the proportion of variance in gene expression\nexplained by modelling cycle time is considerable. Based on this var-\niance, we have been able to identify 3453 genes that change expression\nsignificantly within a 3.4% window of the endometrial cycle (which\napproximately equates to a 24-hour period in a 28-day cycle). Of these,\n488 were during menstruation, 44 during the proliferative phase, and\n2921 during the secretory phase. These changes are remarkably con-\nsistent between samples and synchronised between the 3453 genes\nwithin each sample (see for example Fig.6 and Supplementary Fig 5). It\nis these cycle-driven changes that provide the statistical strength and\nremarkable chronological consistency of our modelling across the\ncycle. If very few genes changed expression across the cycle, most of\nthe expression curves would be relatively flat and there would be no\nobvious global minimum in the mean squared error plots. Instead, the\nhighly dynamic patterns of gene expression ensure that the mean\nsquared error plots give robust results for predicting a time in the\ncycle. It is also important to reiterate that the mean squared error\nsolution is simultaneous for all 20,067 genes, giving con fidence that\nsignificant changes in expression for individual genes are real.\nBecause most of the variance in endometrial gene expression is\nlinked to cycle stage, cycle stage determination using the molecular\nstaging model is unlikely to be in fluenced by other variables that\nimpact on gene expression during the menstrual cycle. This is because\neven if a pathological condition changed the expression 100 ’so f\ngenes, it will only minimally influence the accuracy of the model which\nis based on significant changes in the expression of several 1000’so f\ngenes. The overwhelming contribution of cycle stage variance to the\naccuracy of the molecular model provides the rationale for the gen-\neralisability of the model to endometrial samples taken from women\nwith a wide range of uterine phenotypes and pathologies during the\nmenstrual cycle.\nThe Endometrial Receptivity Analysis (ERA) is a clinical test used\nto determine if the endometrium is at the optimal stage of the men-\nstrual cycle for embryo implantation. The test was originally based on a\ngroup of 238 genes that change expression rapidly around the time of\nimplantation\n22, thus allowing the relative expression pro files to be\nequated with ‘pre-receptive’, ‘receptive’ and ‘post-receptive’ endome-\ntrium. We correctly predicted that using the molecular staging model\nto identify endometrial genes with rapidly changing expression pro-\nfiles would also identify many of the ERA genes. As well as identifying\n70% of the ERA genes we also found many hundreds of additional\ngenes that change expression signi ficantly around the time of\nimplantation, reinforcing how complex endometrial function is at this\ncritical stage of the menstrual cycle.\nUsing the molecular staging model to normalise gene expression\nacross the cycle we were able to identify, for the first time, numerous\nendometrial genes that change expression significantly with increasing\nage (Supplementary Table 3). Importantly, when we performed a meta-\nanalysis by adding a published dataset from another laboratory and\nanalysing it using our method, the number of signi ficant genes was\nincreased; a result that supports our original finding. An interesting\nfinding from this study was that when examining menstrual, pro-\nliferative and secretory stages separately, nearly all the genes that\nArticle https://doi.org/10.1038/s41467-023-41979-z\nNature Communications|         (2023) 14:6222 8\n\nshowed DGE with increasing age did so in endometrial samples taken\nin the secretory phase. This suggests that it is the endometrial\nresponse to progesterone that changes with increasing age, or at least\ngenes with biological actions related to secretory phase events such as\nimplantation, early pregnancy or the late secretory transition to\nmenstruation. One explanation for not identifying significant changes\nin gene expression during the menstrual phase is low subject numbers\n(total N = 28 with 0 genes signi ficant). But for the proliferative phase\nfrom N = 165 subjects we found only 4 signi ficant genes while for the\nsecretory phase withN = 160 subjects there were 218 significant genes.\nThese data highlight the fundamental difference between proliferative\nand secretory endometrium and bring in to question the biological\nrelevance of comparing these two very different tissues.\nThe most signi ficant up and down regulated GO pathways with\nincreasing age have strong biological plausibility (Supplementary\nData 2). Upregulated pathways include axoneme assembly, micro-\ntubule bundle formation and cilium assembly. Microtubules are pre-\nsent in nearly all cells and play a prominent role, along with ciliated\ncells, in well characterised changes to the endometrial luminal epi-\nthelium during the secretory phase. Although the functional relevance\nis not clear, the transition from luminal microvilli to pinopodes around\nthe expected time of embryo implantation is thought to be a key\nmorphological marker for endometrial receptivity\n26.M o r er e c e n te v i -\ndence has found an increased incidence of luminal epithelial micro-\ntubule abnormalities in the endometrium of women with\nadenomyosis, possibly due to increased inflammation\n27. Blood vessel\ndevelopment and regulation of endothelial cell proliferation were the\ntop 2 down regulated GO biological pathways. The endometrium is a\nsite of rapid cyclical angiogenesis, vascular development and break-\ndown, with spiral arterioles developing their classical shape and other\nendometrial vessels growing to match the increasing thickness of the\nendometrium\n28. Age-related alterations in endometrial angiogenic or\nvascular potential could be expected to have profound effects on\nendometrial function, and in particular fertility and subsequent\nplacentation.\nFuture endometrial DGE studies should account for age when\nperforming statistical modelling. It seems likely that use of the mole-\ncular staging model will identify other gynaecological conditions that\ninfluence endometrial DGE and hence need to be accounted for in\nsubsequent analyses.\nThere are limitations to the current study that need to be con-\nsidered. The 236 endometrial samples used to develop the menstrual\ncycle staging model in analysis 2 were nearly all taken from women\nundergoing diagnostic laparoscopy to investigate pelvic pain and\nsuspected endometriosis (IVF subjects with secretory stage endome-\ntrium that were used in analysis 1 were not included in thefinal model).\nOf the 236 women, 173 were diagnosed with endometriosis. Whether\nthis population can be considered as having normal endometrium is\nopen for debate. However, as reported in Supplementary Tables 1 & 2,\n60 of the 236 women had delivered 1 or more live births, with 38 of\nthese 60 also being diagnosed with endometriosis. The ability to carry\na normal pregnancy to a successful outcome is evidence that the\nendometrium was functionally‘normal’ in these women, and suggests\nthat if there are pathological differences in the endometrium of\nwomen with endometriosis, they are either not universal, or not severe\nenough to preclude normal endometrial function. Further, the 140\nwomen who had not been pregnant were on average 6.2 years younger\nthan those that had had a successful pregnancy (29.2 vs 35.4 years of\nage), suggesting that many had not yet tried to conceive and could well\nhave ‘normal’ endometrium. While it is possible that gynaecological\npathologies that affect endometrial gene expression may exist in some\nof our 236 subjects, the conclusion from this work, as discussed earlier,\nis that the predominant driver of changing endometrial gene expres-\nsion is menstrual cycle stage, and n ot phenotypic factors. If pheno-\ntypes that do influence gene expression are identified in the future, for\nexample age and ethnicity in the current work, then these can be\ncorrected for in the model and this will reduce variability and increase\nsensitivity of future studies.\nThe relative lengths of the menstrual, proliferative and secretory\nphases in the molecular staging model are determined by the dis-\ntribution of women at each stage of the menstrual cycle. We made the\nassumption that our population of 236 women were randomly dis-\ntributed across the cycle although it is possible that this was not the\ncase. Hence some phases in the molecular staging model could be\ncompressed or extended over time compared to the normal popula-\ntion. This does not alter the validity of the model for the relative\nposition of the samples in chronological order across the cycle but\ndoes mean that mapping of molecular staging model onto cycle stages\nmay update as larger studies are performed.\nIn conclusion, we have developed and validated a method for\naccurately determining endometrial cycle stage based on global gene\nexpression. Our ‘molecular staging model ’ reveals signi ficant and\nremarkably synchronised daily changes in expression for over 3400\nendometrial genes at different stages of the cycle, with most change\noccurring during the secretory phase. These major day-to-day differ-\nences in endometrial gene expression provide a compelling explana-\ntion for the failure of studies that lack accurate cycle staging to reach\nconsensus on genes of interest. Our study supports selected previous\nfindings and significantly extends existing data. Using the molecular\nstaging model to normalise expression data we can demonstrate sig-\nnificant changes in endometrial gene expression with increasing age.\nThe molecular staging model provides a wealth of new data on\nendometrial gene expression and establishes a method for investi-\ngating the role of the endometrium in critical biological events such as\nuterine receptivity for embryo implantation as well as gynaecological\npathologies such as endometriosis and menstrual disorders.\nMethods\nSubject recruitment and tissue collection\nThis work was approved by the Human Research Ethics Committee of\nthe Royal Women’s Hospital, Melbourne, Australia (Projects 11–24 and\n16–43 for Study 1) and Melbourne IVF (Project 13/17 for Study 2), and all\nsubjects gave written informed consent. Study subjects did not receive\ncompensation for participating. A total of 358 endometrial samples\nwere collected for this study, comprising 264 samples taken at the time\nof diagnostic surgery from women reporting pelvic pain with sus-\npected endometriosis ( ‘Study 1 ’) and 94 samples from individuals\nundergoing IVF (‘Study 2’).\nSome of the 358 subjects have had data published as part of\nprevious studies investigating genetic regulation of endometrial gene\ntranscription. Specifically, 123 Illumina Human HT-12 v4.0 samples\nwere included in a 2017 study\n18, 229 Illumina Human HT-12\nv4.0 samples in a 2018 study 29, and 169 & 206 RNA sequencing sam-\nples in studies from 202019,30.\nEndometrial tissue samples (collected by curette or Pipelle\nbiopsy) were obtained for gene expression analysis, along with blood\nsamples for DNA extraction and hormone assays, patient ques-\ntionnaires, past and present clinical histories, pathologyfindings and\nsurgical notes. All subjects were premenopausal and free from hor-\nmone treatment at the time of biopsy. Endometrial tissue samples\nwere split and either stored in RNA later (Life Technologies, Grand\nIsland, NY, USA) at 4\noC before being stored at -80 oC for total RNA\nextraction, or formalin fix e da n dp r o c e s s e dr o u t i n e l yf o rh i s t o l o g i c a l\nassessment.\nHistological dating of endometrium\nAll 358 endometrial samples were routinely evaluated14 by at least one\nexperienced pathologist and allocated to one of seven menstrual cycle\nstages. Menstrual cycle stage de finitions, assuming a standardised\n28 day cycle, and numbers of samples in each group were as follows:\nArticle https://doi.org/10.1038/s41467-023-41979-z\nNature Communications|         (2023) 14:6222 9\n\nStage 1 = menstrual (n =1 8 ,d a y s1–4), 2 = early proliferative (n = 5, days\n5–7), 3 = mid proliferative ( n =1 0 4 , d a y s 8–11), 4 = late proliferative\n(n =2 9 ,d a y s1 2–15, includes‘interval’), 5 = early secretory (n =6 4 ,d a y s\n16–19, or post ovulation days 2 –5), 6 = mid secretory ( n =7 6 , d a y s\n20–23, or post ovulation days 6 –9), 7 = late secretory ( n =4 0 , d a y s\n24–28, or post ovulation days 10 –14). Twenty-two biopsies taken by\nPipelle did not have adequate tissue for reliable pathology reporting. If\nthe pathology report crossed two definitions (e.g. mid/late secretory),\nthen the later cycle stage was used. The exception to this was in the\nsecretory phase where post-ovulatory day (POD) range crossed 2 cycle\nstages with the majority of the days in the earlier stage, i.e. POD 4–6=\nearly secretory. If the pathology report only recorded‘proliferative’ or\n‘secretory’, then mid-proliferative or mid-secretory was assigned,\nresulting in elevated numbers in these 2 cycle stages. A subset of\nsecretory stage samples (n = 164) underwent additional evaluation and\nwere assigned an individual post-ovulatory day by a further 1 or 2\npathologists working independently of each other and the previous\nassessment. Endometrial samples for which the pathologist reported\nany abnormalities or evidence of exogenous hormone had already\nbeen excluded from the study.\nRNA-seq and gene expression array data preparation\nOf the 358 endometrial samples, 290 had Illumina Human HT-12 v4.0\nperformed and 266 underwent RNA-seq (198 samples had both tech-\nniques performed). Total RNA was isolated from endometrial samples\nusing the Allprep DNA/RNA Mini Kit (Qiagen, CA) as per the manu-\nfacturer’si n s t r u c t i o n s\n29. Brie fly, RNA quality was checked using a\nBioanalyzer 2100 (Agilent Technologies, CA) and RNA concentration\nwas measured using a NanoDropND-6000 (Thermo Fisher Scientific,\nUSA). All samples were high quality with an RNA integrity number\ngreater than 8. Expression profiles in endometrial tissue were gener-\nated by hybridizing 750 ng of cRNA to Illumina Human HT-12 v4.0\nBeadchips.\nRNA samples were treated with Turbo DNA-free kit (Thermo\nFisher Scientific, USA) prior to RNA-seq library generation\n19. Stranded\nRNA-seq libraries were prepared using the Illumina TruSeq Stranded\nTotal RNA Gold protocol which includes ribosomal depletion\n(Illumina, USA).\nRNA-seq data processing\nRaw sequencing reads were quality checked using FastQC v0.11.7 and\nMultiQC v1.6. Low quality reads and contaminating HiSeq Illumina\nadapter sequences were trimmed using Trimmomatic v0.36\n31.T r i m -\nmed reads were aligned against the human reference genome\n(Ensembl Homo sapiens GRCh38 release 84) using HISAT2 v2.0.5\n32.\nTranscript assembly was performed using StringTie v1.3.1 and the\nEnsemblHomo sapiensGRCh38 release 91 reference annotation. Reads\nmapping to each known transcript were quanti fied in StringTie\n33 to\ngenerate transcript-, exon- and intron-level expression matrices for\neach individual. Raw gene count matrices were also produced using a\nPython script provided by StringTie.\nNormalisation of RNA-Seq and array expression values\nGenes expressed at a low level by RNA-seq, i.e. genes with counts per\nmillion (CPM) < 0.5 in > 80% of the samples, were removed. Raw gene\ncounts were normalized for composition bias and total raw reads\n(library size) using the Trimmed Mean of M (TMM) method in the\nedgeR R package\n34. Normalized counts were converted to CPM and\nlog2 transformed (log2-CPM). Batch effects from sequencing were\nremoved using the ComBat function from the sva R package (Leek\net al., 2020). To load and normalise the microarray data, the R\npackages limma\n35 and lumi36 were used. Background correction and\nrobust spline normalization (RSN) produced logged values of probe\nintensity and ComBat was used to remove microarray batch effects.\nFor probes to be included in the array analysis, annotation probe\nquality was required to be “Good” or “Perfect” and detection p-\nvalue < 0.05 in at least 20% of samples.\nGenotyping\nFor determination of ancestry, DNA samples from each of the 358\nindividuals were genotyped on HumanCoreExome or In finium Psy-\nchArray chips (Illumina, USA)\n19. Quality control (QC) was performed in\nPLINK29. Following QC, a total of 282,625 SNPs (hg19) were phased\nusing Shapelt V2 and taken forward to imputation using the haplotype\nreference consortium reference panel (version r1.1 2016) on the\nMichigan Imputation Server. SNPs with low imputation quality\n(R\n2 < 0.8), missing rate>5%, minor allele frequency (MAF) < 1×10–4,a n d\nHardy–Weinberg equilibrium<1×10−6 after imputation were removed.\nThe remaining SNP positions were lifted over to the Ensembl genome\nbuild 38 (GRCh38) using CrossMap v.0.2.8. SNPs failing to lift-over\nwere assigned to their new GRCh38 position manually based on\ndbSNP151 GRCh38 patch release 7 (GRCh38.p7), leaving 6,230,993\nSNPs for further analysis.\nHormone assays\nEstradiol and progesterone concentrations were measured in bloods\ntaken at the time of endometrial sampling. Some of these hormone\ndata have been published previously\n37. An additional 28 bloods were\nassayed for progesterone (Serum P was tested on the Roche Cobas\ne601 immunoanalyser, utilising electrochemiluminescence (ECLIA).\nThe lower limit of detection was 0.06 ng/mL. The inter-assay at a target\nmean of 1.4 ng/mL returned a CV% of 3.7. The intra-assay at a target\nmean of 1.5 ng/mL returned a CV% of 6.5). This gave a total of 187\nprogesterone results and 159 estradiol results that could be plotted\nagainst the molecular staging model cycle stage.\nAnalysis 1\nDevelopment of the‘molecular staging model’ to assign cycle stage for\nsecretory stage samples only: All secretory stage samples with RNA-seq\ndata from Study 1 and Study 2 where 2 or 3 independent pathology\nassessments agreed on the post-ovulatory day (POD) to within 2 days\n(n = 96 of a possible 180 secretory stage samples) were selected for\nanalysis 1. For each gene, batch-corrected expression was fitt oP O D\nusing a penalised cubic regression spline using 3 knots implemented\nwith the generalised additive model (gam) function from the mgcv R\npackage\n38. Each curve was used to obtain the expected expression\nvalue for each gene for any given day. For each sample, an estimated\nPOD was obtained using the day which minimised the mean squared\nerror (MSE) between the observed expression and the expected\nexpression across all genes.\nAlternatively, this procedure can be described minimisingd in the\nloss function:\nLdðÞ =\nX\ng2G\nyg /C0 f g dðÞ\n/C16/C17 2\nð1Þ\nwhere d is the POD constrained between 1 and 14 days, g is a gene in\ngene set G,y g is the observed expression of gene g,a n d fg(d) is the\nspline function that describes the expected expression of gene g for\nday d. Additionally, K-fold cross-validation where K = 5 was performed\nto ensure the model was not overfitting.\nTo illustrate that using larger, less precise, units of time can be\nused to estimate cycle time using the same method, an additional\nmodel was built using the pathology-assigned 3 secretory cycle stages\n(i.e. stages 5, 6, and 7 corresponding to early-, mid-, and late-secretory\nrespectively) instead of the pathology-assigned POD (i.e. 1 –14 days).\nUsing the RNA-seq batch-corrected expression data, each gene wasfit\nusing the same penalised cubic regression spline (k = 3) as a function of\nstage, and curves were generated for each gene. Cycle time was esti-\nmated for each sample by selecting the time point in the stage (from\nArticle https://doi.org/10.1038/s41467-023-41979-z\nNature Communications|         (2023) 14:6222 10\n\n4.5 to 7.5) that minimised the MSE between the observed expression\nand the expected expression across all genes. As validation, cycle time\nestimated from the 3 stages model was compared to the cycle time\nestimated from the 14-day POD model.\nAnalysis 2\nMolecular staging model using 7 pathology stages for the whole cycle\nwith RNA-seq and array expression data: The method for developing\nthe POD of cycle prediction model was replicated with some mod-\nifications using all samples from Study 1 ( N = 236 for NGS) classi fied\ninto 7 cycle stages by histopathology. Because the majority of the\nproliferative phase samples were not assigned as early, mid or late by\nthe pathologist, we re-assigned all proliferative samples into early, mid\nor late by fit t i n gap e n a l i s e dc u b i cr e g r e s s i o ns p l i n e( k = 3 )u s i n g\nmenstrual, proliferative, and early secretory gene expression data. A\nproliferative time point was estimated using the time-point which\nminimised the mean squared error between the observed expression\nand the expected expression across all genes. The proliferative sam-\nples were then split into equal groups of early, mid, and late using this\ntime point. This approach assumes that patients presenting for surgery\nin a public hospital system will approximate a uniform distribution\nacross the menstrual cycle, so the number of early, mid, and late\nproliferative samples will be approximately equal. Once the pro-\nliferative samples were assigned to early, mid and late stages, a pena-\nlised cyclic cubic regression spline (k = 8) was fit using the 7 stages.\nEach endometrial sample is then assigned a‘day’ from the model using\nthe time which minimises the mean squared error between the\nobserved expression data for all genes and their corresponding gene\nmodels. This is equivalent to minimising the loss function in Eq. ( 1),\nexcept now d can represent any timepoint in the cycle. This ‘day’ is a\nrelative timepoint in the cycle and does not correspond to a real day.\nContinuing with the assumption that all 236 women were approxi-\nmately uniformly distributed across the cycle, the data are trans-\nformed so that the distance in time between each sample is identical.\nThis process in effect ranks all the samples in order from the start to\nthe end of the cycle, and no longer relies on assigning days from an\nidealised 28-day cycle. The x-axis was then scored from 0–100 so that\nthe individual scores for each sample represented the percentage of\nthe way through the menstrual cycle that each sample was.\nThe gene curves were then re fitted using the new derived cycle\ntimes for each sample with a penalized cyclic cubic regression spline\n(k = 30). For visualisation purposes, normalisation of gene expression\nfor cycle stage was then derived by subtracting the expected expres-\nsion from the observed expression (i.e. calculating the residuals) and\nre-adding the mean.\nValidation of the molecular staging model\nSeveral studies were conducted t o assess how endometrial gene\nexpression data that were generated using the molecular staging\nmodel performed.\nThe first validation study as part of Analysis 1 was to compare\nresults from the secretory stage that had daily POD pathology, and\nhence substantially more accurate cycle stage information, with results\nthat only used 3 pathology stages across the secretory phase. The\nsimilar results from the 2 different pathology inputs validated the use\nof pathology data dividing the endometrium into 7 cycle stages to\ndevelop the molecular staging model across the whole menstrual cycle.\nThe second validation study was to repeat Analyses 1 and 2 using\nIllumina HT-12 data and then compare results for the 198 out of\n358 samples that had both RNA-seq and Illumina HT-12 data.\nThe third validation study was to repeat Analysis 2 using unsu-\npervised methods (i.e. without the use of pathology cycle stage\ninformation). K-means clustering (k = 5) was performed using thefirst\n2 dimensions after dimensionality reduction principal component\nanalysis (PCA) (Supplementary Fig 1). Clusters were ordered by their\ndistance between each other, and for each gene, expression wasfitted\nusing a penalised cyclic cubic regression spline (k = 6) using the\nordered cluster number as a proxy for cycle time. Identical to Analysis\n2, for each sample, an optimal cycle time was obtained by calculating\nthe time that minimised the MSE, and subsequently all times were\ntransformed to be uniformly distributed on a scale from 0 to 100.\nCycle time estimates from the validation model were compared to the\nAnalysis 2 cycle times, however, since time 0 for the validation model\nwas not expected to correspond to the start of menstruation, time was\noffset by the earliest menstrual sample (identified by Analysis 2) before\nassessing performance (Supplementary Fig 2).\nThe fourth validation study was to determine whether indepen-\ndent endocrine data supported the molecular staging model. For a\nsubset of samples where endocrine data were available, peripheral\nblood estradiol (n = 159) and progesterone (n = 187) levels were plotted\nagainst the cycle stage from the molecular staging model.\nThe molecular staging model was also used to re-analyse 3 pub-\nlished datasets available in GEO (GEO DataSets ID: GSE65099\n21,\nGSE14154920 and GSE18048517). Principal component analysis (PCA)\nplots from these data sets were replotted with the unaltered cycle\nstage from the original publication and our molecular staging model\ncycle stage for comparison.\nApplication of the molecular staging model\nThe molecular staging model for gene expression across the menstrual\ncycle was applied to 3 questions: Does (1) age or (2) ancestry have any\ninfluence on endometrial gene expression, and (3) at what stages of the\ncycle does gene expression change most rapidly? Differential expres-\nsion analysis was performed with the predicted cycle time as an\nadditional factor in the linear model, modelled as a regression spline.\nEmpirical Bayes moderated t-tests, implemented in the R limma\npackage, were used to assess if genes were differentially expressed.\nAge of patient at biopsy was analysed as a continuous variable\nusing all subjects with RNA-seq data ( N = 266). Since differential\nexpression effects due to age were identified, subsequent differential\ngene expression (DGE) analyses included age as a factor when fitting\nthe linear model. To further explore the effects of age on endometrial\ngene expression n = 87 subjects from GSE141549 were also analysed,\nand a meta-analysis was performed using the weighted Fisher ’s\nmethod for combining p-values implemented in the metapro R\npackage\n39, where weights were proportional to each study ’ss a m p l e\nsize. Ensembl ID’s from the current data were matched with the Illu-\nmina probe ID’s from GSE141549 resulting in 12,868 genes in common\nbetween the 2 data sets. If multiple probes matched to the same\nEnsembl ID, the probe with the greatest mean expression was used.\nAnalyses were run for the whole menstrual cycle, and separately for the\nmenstrual, proliferative and secretory phases. After multiple hypoth-\nesis correction, genes with false discovery rate (FDR) corrected\nP < 0.05 were considered to be differentially expressed. Gene ontology\nenrichment analysis was performed using the clusterPro filer R\npackage\n40.\nAncestry of subjects as defined by a previous study19 was used to\nlook for differences in endometrial gene expression using pairwise\ncomparisons of each ancestry group. Genetic ancestry of subjects was\nassigned using principal component analysis (PCA). Briefly, genotype\ndata from participants was merged with the 1000 Genomes P3v5\nreference data using markers common to both cohorts. Population\nclusters were determined using the first 5 PCs and were annotated\naccording to the five 1000 genomes super populations (European,\nEastern Asian, Admixed America, African, Southern Asian).\nT oi n v e s t i g a t ec h a n g i n gg e n ee x p r e s s i o na c r o s st h ec y c l ea l l\nsamples with RNA-seq data from Study 1 ( n = 236) were ranked in\nchronological order from start to end of the cycle. A ‘sliding window’\napproach was then used to compare DGE between samples 1–8v e r s u s\nsamples 9–16, followed by samples 2 –9 versus 10 –17, then 3 –10 vs\nArticle https://doi.org/10.1038/s41467-023-41979-z\nNature Communications|         (2023) 14:6222 11\n\n11–18 and so on for all samples across the menstrual cycle. Group sizes\nwere arbitrarily set at 8 because this represents 3.4% of the 236 samples\nor approximately 1 day assuming a mean cycle length of 28 days.\nModerated t-tests with a fold-change cut-off of 1.2, as implemented in\nlimma’s treat function, were used to identify differentially expressed\ngenes with P\nadj <0 . 0 5a te a c hw i n d o w .\nTo validate a subset of genes we identi fied as having rapidly\nchanging expression around the time of embryo implantation, we ran a\ncomparison with genes identified in the original Endometrial Recep-\ntivity Analysis (ERA) publication\n22. The ERA identi fied 238 genes that\nshow major changes in expression before, during and after the time of\nembryo implantation at POD 3–7. We identified all endometrial genes\nthat showed significant changes in expression within any given 24 hr\nperiod at the same time of the menstrual cycle to confirm that our list\ncontained a high proportion of the ERA genes.\nReporting summary\nFurther information on research design is available in the Nature\nPortfolio Reporting Summary linked to this article.\nData availability\nThe RNA-seq data generated in this study have been deposited in the\nGene Expression Omnibus database under accession number\nGSE234354 and the Illumina HT-12 microarray data under accession\nnumber GSE234368. Additional public datasets used are available on\nGEO with accession numbersGSE65099 (PMID: 26418742),GSE141549\n(PMID: 32859947) and GSE180485 (PMID: 35092277). Homo sapiens\nreference genome GRCh38 and genome annotation were obtained\nfrom Ensembl (release 91). Source data are provided with this paper.\nCode availability\nThe code for our model is available as an R package at https://github.\ncom/jessicachung/endest (DOI: 10.5281/zenodo.8321573) and an R\nShiny application is available at https://github.com/jessicachung/\nendspect. Analysis scripts for this manuscript can be found at\nhttps://github.com/jessicachung/endo_model_paper.E n d e s ti sp u b -\nlished pursuant to the terms located in the applicable repository at\ngithub.com which terms permit reproduction, publication and adap-\ntation of endest solely for non-commercial purposes. Publication of\nendest at github.com is not subject to the publication terms applicable\nto Nature publications.\nReferences\n1. Jain, V. et al. Uterine bleeding: how understanding endometrial\nphysiology underpins menstrual health.Nat. Rev. Endocrinol. 18,\n290–308 (2022).\n2. Darba, J. & Marsa, A. Economic Implications of Endometriosis: A\nReview. Pharmacoeconomics40, 1143–1158 (2022).\n3. Hoppes, E. et al. Global research and learning agenda for building\nevidence on contraceptive-induced menstrual changes for\nresearch, product development, policies, and programs.Gates\nOpen Res 6, 49 (2022).\n4. Ponnampalam, A. P. et al. Molecular classi fication of human endo-\nmetrial cycle stages by transcriptional profiling. Mol. Hum. Reprod.\n10,8 7 9–893 (2004).\n5. Kao, L. C. et al. Expression pro fil i n go fe n d o m e t r i u mf r o mw o m e n\nwith endometriosis reveals candidate genes for disease-based\nimplantation failure and infertility.Endocrinology144,\n2870–2881 (2003).\n6. Girling, J. E. et al. Differential Gene Expression in Menstrual Endo-\nmetrium From Women With Self-Reported Heavy Menstrual\nBleeding. Reprod. Sci. 24,2 8–46 (2017).\n7. Koot, Y. E. et al. An endometrial gene expression signature accu-\nrately predicts recurrent implantation failure after IVF.Sci. Rep. 6,\n19411 (2016).\n8. Aghajanova, L. et al. Stanniocalcin-1 expression in normal human\nendometrium and dysregulation in endometriosis.Fertil. Steril.106,\n681–691.e1 (2016).\n9. Aghajanova, L. et al. Effects of noncavity-distorting fibroids on\nendometrial gene expression and function.Biol. Reprod. 97,\n564–576 (2017).\n10. Najmabadi, S. et al. Menstrual bleeding, cycle length, and follicular\nand luteal phase lengths in women without known subfertility: A\npooled analysis of three cohorts.Paediatr. Perinat. Epidemiol.34,\n318–327 (2020).\n11. Soumpasis, I., Grace, B. & Johnson, S. Real-life insights on men-\nstrual cycles and ovulation using big data.Hum. Reprod. Open\n2020, hoaa011 (2020).\n12. Bull, J. R. et al. Real-world menstrual cycle characteristics of more\nthan 600,000 menstrual cycles.NPJ Digit Med 2,8 3( 2 0 1 9 ) .\n13. Tatsumi, T. et al. Age-Dependent and Seasonal Changes in Men-\nstrual Cycle Length and Body Temperature Based on Big Data.\nObstet Gynecol. 136, 666–674 (2020).\n1 4 . N o y e s ,R .W . ,H e r t i g ,A .T .&R o c k ,J .D a t i n gt h eE n d o m e t r i a lB i o p s y .\nFertil. Steril.\n1,3 –25 (1950).\n15. Duggan, M. A. et al. The accuracy and interobserver reproducibility\nof endometrial dating.Pathology 33,2 9 2–297 (2001).\n16. Ruiz-Alonso, M. et al. Endometrial Receptivity Analysis (ERA): data\nversus opinions.Hum. Reprod. Open 2021, hoab011 (2021).\n17. Lipecki, J. et al. EndoTime: non-categorical timing estimates for\nluteal endometrium.Hum. Reprod. 37,7 4 7–761 (2022).\n18. Fung, J. N. et al. The genetic regulation of transcription in human\nendometrial tissue.Hum. Reprod. 32,8 9 3–904 (2017).\n19. Mortlock, S. et al. Tissue speci fic regulation of transcription in\ne n d o m e t r i u ma n da s s o c i a t i o nw i t hd i s e a s e .Hum. Reprod. 35,\n377–393 (2020).\n2 0 . G a b r i e l ,M .e ta l .Ar e l a t i o n a ld a tabase to identify differentially\nexpressed genes in the endometrium and endometriosis lesions.\nSci. Data 7, 284 (2020).\n2 1 . L u c a s ,E .S .e ta l .L o s so fE n d o m e t rial Plasticity in Recurrent Preg-\nnancy Loss. Stem Cells 34,3 4 6–356 (2016).\n22. Diaz-Gimeno, P. et al. A genomic diagnostic tool for human endo-\nmetrial receptivity based on the transcriptomic signature.Fertil.\nSteril. 95,5 0–6 0 ,( 2 0 1 1 ) .6 0 . e 1 - 1 5 .\n23. Wang, A. et al. Health Care Costs and Treatment Patterns Asso-\nciated with Uterine Fibroids and Heavy Menstrual Bleeding: A\nClaims Analysis.J. Women’sH e a l t h( L a r c h m t . )31,8 5 6–863 (2022).\n24. Schoep, M. E. et al. Productivit y loss due to menstruation-related\nsymptoms: a nationwide cross-sectional survey among 32 748\nwomen. BMJ Open 9, e026186 (2019).\n25. Niederberger, C. et al. 25 historic papers: an ASRM 75th birthday\ngift from Fertility and Sterility.Fertil. Steril. 112,e 2–e27 (2019).\n26. Quinn, K. E. et al. Pinopodes: Recent advancements, current per-\nspectives, and future directions.Mol. Cell Endocrinol. 501,\n110644 (2020).\n27. Khan, K. N. et al. An axonemal alteration in apical endometria of\nhuman adenomyosis.Hum. Reprod. 36,1 5 7 4–1589 (2021).\n28. Rogers, P. A. et al. Endometrial angiogenesis, vascular maturation,\nand lymphangiogenesis.Reprod. Sci. 16,1 4 7\n–151 (2009).\n29. Fung, J. N. et al. Genetic regulation of disease risk and endometrial\ngene expression highlights potential target genes for endome-\ntriosis and polycystic ovarian syndrome.Sci. Rep. 8, 11424 (2018).\n30. Holdsworth-Carson, S. J. et al. Obesity does not alter endometrial\ngene expression in women with endometriosis.Reprod. Biomed.\nOnline 41, 113–118 (2020).\n31. Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trim-\nmer for Illumina sequence data. Bioinformatics30,\n2114–2120 (2014).\n32. Kim, D. et al. Graph-based genome alignment and genotyping with\nHISAT2 and HISAT-genotype.Nat. Biotechnol.37,9 0 7–915 (2019).\nArticle https://doi.org/10.1038/s41467-023-41979-z\nNature Communications|         (2023) 14:6222 12\n\n33. Kovaka, S. et al. Transcriptome assembly from long-read RNA-seq\nalignments with StringTie2.Genome Biol. 20,2 7 8( 2 0 1 9 ) .\n34. Robinson, M. D. & Oshlack, A. A sc aling normalization method for\ndifferential expressionanalysis of RNA-seq data. Genome Biol. 11,\nR25 (2010).\n35. Ritchie, M. E. et al. limma powers differential expression analyses\nfor RNA-sequencing and microarray studies.Nucleic Acids Res. 43,\ne47 (2015).\n36. Du, P., Kibbe, W. A. & Lin, S. M. lumi: a pipeline for processing\nIllumina microarray.Bioinformatics24,1 5 4 7–1548 (2008).\n37. Marla, S. et al. Genetic risk factors for endometriosis near estrogen\nreceptor 1 and coexpression of genes in this region in endometrium.\nMol. Hum. Reprod. 27, gaaa082 (2021).\n38. Wood, S. N., Pya, N. & Säfken, B. Smoothing Parameter and Model\nSelection for General Smooth Models.J. Am. Stat. Assoc. 111,\n1548–1563 (2016).\n39. Yoon, S. et al. Powerful p-value combination methods to detect\nincomplete association.Sci. Rep. 11,6 9 8 0( 2 0 2 1 ) .\n40. Yu, G. et al. clusterPro filer: an R package for comparing biological\nthemes among gene clusters.OMICS 16,2 8 4–287 (2012).\nAcknowledgements\nWe acknowledge the work of Research Nurses Ranita Charitra, Tracy\nMiddleton and Irene Bell who recruited all the participants in this study\nand Leonie Cann who provided laboratory support. Thank you to the\nsurgeons and theatre booking staff at the Royal Women’s Hospital and\nStaff at Melbourne IVF. We also acknowledge and thank the many\nwomen who provided tissues and clinical data that were used in this\nstudy. Research reported in this publication was supported in part by\nNational Health and Medical Research Council (NHMRC) project grants\nGNT1012245 and GNT1105321 (PAWR, GWM, JEG, SJH), NHMRC Fel-\nlowship GNT1078399 and Investigators Grant GNT1177194 (GWM),\nNHMRC Medical Postgraduate Scholarship No. 1055814 and RANZCOG\nFellows’ Clinical Research Scholarship (2013) (WTT), funds from Mel-\nbourne IVF, and the Medical Research Future Fund MRF1199715 (PAWR,\nMH, SJH, JFD).\nAuthor contributions\nWTT, JC and PAWR were responsible for the conception and design of\nthe study. PAWR, WTT, SJH, JFD, MH and JEG acquired the tissues along\nwith clinical and surgical data. HR, SB and VO provided expert patho-\nlogical assessment of endometrial biopsies. JC and CS provided bioin-\nformatics and statistical input. RK, JNF, SM and GWM undertook all the\ngene expression assays and initial QC of the expression data. All authors\nhave read, provided input and approved the final manuscript.\nCompeting interests\nThe method described in this study is the subject of a patentfiled at the\nAustralian Patent office (21/6/23 (PCT/AU2023/050559).), by the Uni-\nversity of Melbourne with WTT, JC and PAWR as inventors. The remaining\nauthors declare no competing interests.\nAdditional information\nSupplementary informationThe online version contains\nsupplementary material available at\nhttps://doi.org/10.1038/s41467-023-41979-z.\nCorrespondenceand requests for materials should be addressed to P.\nA. W. Rogers.\nPeer review informationNature Communicationsthanks Sascha Ott and\nthe other, anonymous, reviewer(s) for their contribution to the peer\nreview of this work. A peer review file is available.\nReprints and permissions informationis available at\nhttp://www.nature.com/reprints\nPublisher’s note Springer Nature remains neutral with regard to jur-\nisdictional claims in published maps and institutional affiliations.\nOpen Access This article is licensed under a Creative Commons\nAttribution 4.0 International License, which permits use, sharing,\nadaptation, distribution and reproduction in any medium or format, as\nlong as you give appropriate credit to the original author(s) and the\nsource, provide a link to the Creative Commons licence, and indicate if\nchanges were made. The images or other third party material in this\narticle are included in the article’s Creative Commons licence, unless\nindicated otherwise in a credit line to the material. If material is not\nincluded in the article’s Creative Commons licence and your intended\nuse is not permitted by statutory regulation or exceeds the permitted\nuse, you will need to obtain permission directly from the copyright\nholder. To view a copy of this licence, visithttp://creativecommons.org/\nlicenses/by/4.0/.\n© The Author(s) 2023\nArticle https://doi.org/10.1038/s41467-023-41979-z\nNature Communications|         (2023) 14:6222 13","source_license":"CC0","license_restricted":false}