Identifying a causal link between prolactin signaling pathways and COVID-19 vaccine-induced menstrual changes

COVID-19 vaccines have been instrumental tools in the fight against SARS-CoV-2 helping to reduce disease severity and mortality. At the same time, just like any other therapeutic, COVID-19 vaccines were associated with adverse events. Women have reported menstrual cycle irregularity after receiving COVID-19 vaccines, and this led to renewed fears concerning COVID-19 vaccines and their effects on fertility. Herein we devised an informatics workflow to explore the causal drivers of menstrual cycle irregularity in response to vaccination with mRNA COVID-19 vaccine BNT162b2. Our methods relied on gene expression analysis in response to vaccination, followed by network biology analysis to derive testable hypotheses regarding the causal links between BNT162b2 and menstrual cycle irregularity. Five high-confidence transcription factors were identified as causal drivers of BNT162b2-induced menstrual irregularity, namely: IRF1, STAT1, RelA (p65 NF-kB subunit), STAT2 and IRF3. Furthermore, some biomarkers of menstrual irregularity, including TNF, IL6R, IL6ST, LIF, BIRC3, FGF2, ARHGDIB, RPS3, RHOU, MIF, were identified as topological genes and predicted as causal drivers of menstrual irregularity. Our network-based mechanism reconstruction results indicated that BNT162b2 exerted biological effects similar to those resulting from prolactin signaling. However, these effects were short-lived and didn’t raise concerns about long-term infertility issues. This approach can be applied to interrogate the functional links between drugs/vaccines and other side effects.

The occurrence of post-vaccine menstrual cycle disturbances has gone unnoticed during clinical trials phase of COVID-19 vaccines.Then, thousands of reports pointing to menstrual changes started to surface following worldwide vaccination campaigns 14,16,25 .Healthcare authorities dismissed these claims at the beginning and considered them unjustified.But such reports continued to appear from various countries around the world which led to increased vaccine hesitancy in young women [26][27][28][29] .
Recently, serious concerns have been raised about the effects of COVID-19 vaccines on menstruation [27][28][29][30][31] , and these fears keep escalating.Thousands of women reported post-COVID-19-vaccine menstrual changes to health care authorities around the world, and several published studies indicated an association between menstrual abnormalities and COVID-19 vaccines [32][33][34] .Women feared that menstrual changes suggest long-term adverse effects on fertility and pregnancy, which led to hesitation against vaccination among women.In fact, menstrual changes have been reported after receiving both mRNA and adenovirus vectored COVID-19 vaccines which led to the speculation that it is the nature of the immune response to vaccines, rather than vaccine components, that led to these adverse events 13,14,16,31,[35][36][37] .Furthermore, previous reports suggested that vaccine-associated period changes occur due to transient perturbations to the hypothalamic-pituitary-ovarian (HPO) axis [38][39][40] .In fact, changes in menstrual cycles have been reported for non-COVID-19 vaccines including the human papillomavirus and typhoid vaccines 41,42 .
Herein, we devised a workflow to assess menstrual adverse events in response to treatment with mRNA vaccine BNT162b2.Our results revealed a causal link implicating prolactin signaling and hormone-induced effects on the menstrual cycle and endometrium resulting from post-vaccine gene expression perturbances.Fortunately, gene expression perturbations were short-term and therefore are not expected to cause long-term menstrual irregularities.The approach devised and implemented herein can be applied to assess other vaccines and other vaccineinduced biological effects.

Systems biology findings
We undertook a systems biology approach to derive transcriptional signatures for COVID-19 mRNA vaccines relying on BNT162b2 transcriptomics data.Our workflow is shown in Fig. 1.Previously, we applied similar approaches to explore the network pharmacology of drugs and vaccines 43 , as well as investigating disease pathogenesis pathways for prioritizing biomarkers and drug targets 11,12 .Each study workflow is tweaked to suit the scientific questions we are asking as well as the types of data we have.In this study, the analysis of transcriptional raw data extracted from GSE169159 indicated that gene expression alterations on day 22 of receiving the first vaccine dose (i.e., the day after receiving the second vaccine dose) affected genes that are known biomarkers or drug targets for menstrual cycle disturbances.All details on deriving gene signatures from transcriptomics data of GSE169159 are described elsewhere 12 .Additionally, no significant DEGs were observed on day 28 after receiving the first vaccine dose.

Vaccine Gene Signatures
Two transcriptomic gene signatures (GS) for BNT162b2 vaccine were derived from gene expression profiling experiments in response to treatment with vaccine (GSE169159) 44 .The first gene signature (GS1) consisted of 1853 differentially expressed genes (DEGs), 884 upregulated and 969 downregulated DEGs that satisfied the prioritization criteria of a false discovery rate (FDR) ≤ 0.05, and a log2 fold change (log 2 FC) ≥ 2.00 or ≤ −2.00.The second gene signature (GS2) consisted of 108 DEGs, 74 upregulated and 18 downregulated DEGs that satisfied the prioritization criteria of an FDR ≤ 0.05, and a log 2 FC ≥ 5.00 or ≤ −5.00.
In order to get a better idea about the biological significance of the DEGs in response to treatment with BNT162b2, we applied a bioinformatics workflow relying on both downstream enrichment analysis, and upstream analysis for putative regulators responsible for causing the gene expression changes observed in transcriptomics data.

Upstream regulation analysis
Upstream analysis was performed using the DEGs in GS1 (applying a threshold of log 2 FC) and GS2 using causal reasoning 45 .GS1 consisted of 1853 DEGs and therefore was trimmed by our upstream analysis algorithm to reduce the complexity of generated results.The algorithm automatically applied log 2 FC ≥ +2.62 or ≤ −2.62 threshold which reduced the number of DEGs to 1107 (743 upregulated and 612 downregulated); since the applied causal reasoning algorithm requires a query list of about thousand genes on average.This analysis resulted in the prediction of 826 activated and 480 inactivated upstream regulators including transcription factors, kinases, phosphatases, and microRNAs (Supplementary Table 1).All prioritized upstream regulators have prediction activities with p-values ≤ 0.05 and a calculation distance = 1-3.
To focus on upstream regulators of genes with the maximal differential expression in response to vaccination, we applied the same analysis to GS2 consisting of 92 DEGs (74 upregulated and 18 downregulated) that had log 2 FC ≥ +5.00 or ≤ −5.00.Our analysis resulted in predicting 625 activated proteins 389 inactivated (Supplementary Table 2).
All prioritized upstream regulators were scored based on the number of differentially expressed genes that can be reached via the shortest paths, and the correctness of the regulation.The activity prediction correctness is assessed based on the activation and inhibition edges along the paths and the expected and observed directionality of fold changes of the DEGs.It should be noted that the calculation distance is one of the most important parameters that can distinguish direct regulation effects from indirect effects.For example, a calculation distance of one means that the upstream regulator is one step away from the transcriptional event indicating that the regulation event in direct.

Filtering upstream regulators by molecular function and distance
The changes in gene expression that we observe in response to perturbagens (e.g., vaccines), often result from interactions between gene regulatory regions and regulatory proteins such as transcription factors, kinases, phosphatase, RNA molecules and others.Transcription factors (TFs) are considered topologically more important than DEGs especially for purposes of mechanism reconstruction; or at least complementary to DEGs for reconstructing the biological pathways and networks responsible for a phenotype of interest (e.g., menstrual irregularities in response to vaccination).Although transcription factors are not the only master regulators, but they are probably one of the easiest to validate experimentally using in vitro assays.
Hence, we sought filtering the predicted causal regulatory proteins by "molecular function" = "transcription factor" and "calculation distance" = 1.A calculation distance=1, indicates direct effects on transcription.This filtering step resulted in 13 transcription factors (IRF1, STAT1, STAT2, RELA, IRF9, SPI, NFKB1, IRF3, IRF7, BCL6, PRDM1, GATA3) for DEGs in GS1 with log 2 FC ≥ +2.62 or ≤ −2.62 (Supplementary Table 3).Applying the same filters on the causal regulatory proteins predicted for GS2 resulted in the prediction of five transcription factors (IRF1, STAT1, RELA, IRF3 and STAT 2), with a calculation distance = 1 (Supplementary Table 4).The overlapping TFs resulting for the previous two causal reasoning analyses are listed in Table 1.Such hits can be considered higher confidence transcription factors for causing the observed phenotype.
Two of the five higher confidence transcription factors (IRF1 and IRF3) belong to the interferon regulatory factors, and other two transcription factors (STAT1 and STAT2) belong to the signal transducers and transcription activators which mediate cellular responses to interferons.The fifth transcription factor, RelA, belongs to the Rel homology domain/immunoglobulin-like fold and is a regulator of NF-kB activity.As a validation of these findings, we reported examples of supporting evidence in the biomedical literature linking these five predicted higherconfidence transcription factors to menstrual cycle irregularities (Table 1).In fact, supporting studies cited in Table 1 brought our attention to significant interactions between these transcription factors and prolactin/PRL gene.Among all predicted transcription factors for GS1 and GS2, IRF1 had the smallest activity prediction p-values values (4.237E-16 for GS1, and 7.63E-06 for GS2).Therefore, IRF's causal reasoning network shown in Fig. 2a, b.This network serves as an example of the causal reasoning networks we relied in this work.Additional networks for STAT1, RELA, IRF3 and STAT 2 are provided in Supplementary Material (Supplementary Fig. 1-4).
Next, we checked overlaps between the 213-biomarker set and 3 gene sets of interest: 1) DEGs with log 2 FC ≥ +2.00 or ≤ −2.00, 2) causal hubs predicted using DEGs with log 2 FC ≥ +2.00 or ≤ −2.00, and 3) causal hubs predicted using DEGs with log 2 FC ≥ +5 or ≤ −5.Our results indicated that TNF was the only common gene between our biomarkers list and all other three gene lists.Nine genes (ARHGDIB, LIF, FGF2, MIF, IL6R, IL6ST, RHOU, BIRC3 and RPS3) were common between the biomarkers list and two gene sets.Finally, there were 33 common genes between menstrual cycle biomarkers and vaccine-induced DEGs and/or predicted causal proteins (Table 2), and 35 gene overlaps between prolactin biomarkers and vaccine-induced DEGs and/or predicted causal proteins (Table 3).All gene overlaps are shown in Fig. 3a, b.We could look at overlaps with more gene lists as a method of filtering higher confidence "topological" genes that may be driving the menstrual irregularities in response to COVID-19 vaccines (i.e., increased confidence due to biological relevance).Correct/total network predictions show the number of genes in the dataset predicted correctly over the total number of genes in the causal reasoning network.
c Calculation distance from the upstream regulatory key hubs and downstream genes.For example, distances of 2 and 3 identify key hubs that are distant key hubs, while a distance of 1 identify closest one-step away transcriptional factors.
d p-value of the predicted protein activity calculated using the polynomial test. e Evidence for an existing experimental link between the transcription factor and menstrual irregularity.f p-value of the predicted protein activity calculated using the polynomial test for GS1.g p-value of the predicted protein activity calculated using the polynomial test for GS2.*For GS1.**For GS2.

Pathway enrichment and interconnectivity analysis
We used pathway enrichment analysis to assess whether the identified causal regulators work collectively to affect certain biological pathways.Pathway enrichment results using five higher confidence upstream causal regulators as a query gene list (IRF1, STAT1, RELA, IRF3 and STAT 2), highlighted the prolactin signaling pathway as one of the significantly enriched pathways with an enrichment p-value of 1.60E-05 (Fig. 4a).Next, we included PRL in the query list (PRL, IRF1, STAT1, RELA, IRF3 and STAT 2) for pathway enrichment which led to the prioritization of prolactin singling pathway as the top enriched pathway with an enrichment p-value of 8.92E-07 (Fig. 4b).Similar analyses were performed on the 33 filtered biomarker genes/proteins (Fig. 4c), and biomarker proteins in addition to PRL, IRF1, STAT1, RELA, IRF3 and STAT 2 (Fig. 4d).
It is well known that changes in prolactin signaling can result in menstrual cycle irregularities [47][48][49][50][51] .Searching PubMed using terms "prolactin" AND "menstrual" returned 1950 results, while search terms "prolactin" AND "menstrual cycle" returned 1270 results, and search terms "prolactin" AND "menstrual irregularity" returned 32 results.This is a clear indication that the studied mRNA COVID-19 vaccine BNT162b2 has the potential to cause menstrual irregularities by inducing perturbations in the genes and/or proteins involved in prolactin signaling pathways (i.e., through affecting the activity of key transcription factors involved in this pathway).Prolactin signaling pathway affects a wide range of physiological processes ranging from reproduction and lactation to growth and development, from endocrinology and metabolism to brain and behavior, as well as immune regulation (Fig. 5a, b).
Prolactin is a polypeptide hormone encoded by the PRL gene and secreted by the anterior pituitary gland.It is known as a Fig. 2 Causal reasoning networks.a Causal reasoning network of highest confidence transcription factor IRF1 using DEGs in GS1.b Causal reasoning network of highest confidence transcription factor IRF1 using DEGs in GS2.Gene expression changes are shown in green and red sectors around each molecule.Increased expression value corresponds to the green sector which size increases clockwise around the molecule icon.Decreased expression value corresponds to the red sector which size increases counterclockwise.Supportive data panel contains over and under-expressed genes from the experimental data set which support a hypothesis that IRF1 is in a predicted predominant "active" state.Conflicting data panel contains over and under-expressed genes from the experimental data set which are discordant with the hypothesis that IRF1 is in predicted predominant "active" state.
growth regulator for many tissues, including cells of the immune system.It functions as a cytokine with immunomodulatory activities.It may also play a role in cell survival by suppressing apoptosis, and it is essential for lactation.Chemically, prolactin's structure is similar to those of the growth hormone and the placental lactogen hormone, which form together what is known as the "prolactin/growth hormone/placental lactogen" family, and they all originated from one ancestral gene.
Signaling pathway impact analysis (SPIA) SPIA was performed on the combined list consisting of the DEGs and the predicted causal hubs.Enrichment analyses performed on the combined list of DEGs, and key hubs have previously been shown to highlight more biologically relevant results 52 .The top five enriched pathway maps were: 1) immune response interferonalpha/beta (IFN-alpha/beta) signaling via JAK/STAT, 2) regulation of antiviral response by SARS-CoV-2, 3) immune response antiviral actions of interferons, 4) immune response induction of apoptosis and inhibition of proliferation mediated by interferon-gamma (IFN-gamma), and 5) immune response IFN-gamma in macrophages activation with impact p-values of 4.31E-25, 1.24E-11, 1.04E-10, 5.61E-10, and 1.59E-9 respectively.All predicted pathways are highlighting the role of interferons.Hence, SPIA results  All molecules were ranked based on their activity prediction p-values as well their overlap confidence score.An overlap confidence score of 3 indicates that a specific gene/protein is overlapping between the biomarker set and 3 other gene sets, while a score of 1 indicates that gene/protein is overlapping between the biomarker list and one other gene set.
revealed a prominent role of interferon signaling in the signaling pathways impacted by BNT162b2 vaccines.

Supporting evidence
We mined the VAERS database, PubMed using Abstract Sifter, and the Connectivity Map to gather supporting evidence for the prioritized hypotheses regarding vaccine-induced menstrual irregularities and prolactin or HPO mimicking effects.
VAERS database.We searched the VAERS database for vaccine adverse events that are relevant to menstrual irregularities.First, we extracted all COVID-19 vaccine adverse events and ranked them according to their proportions of the overall reported adverse events for COVID-19 vaccines.VAERS annotates menstrual irregularity under 'menstruation irregular', in addition to using more specific "symptom" terms, including heavy menstrual bleeding, dysmenorrhea, intermenstrual bleeding, amenorrhea, postmenopausal hemorrhage, premenstrual pain, menstrual discomfort, menstruation normal, premenstrual dysphoric disorder, menstrual cycle management, premenstrual headache, retrograde menstruation, premenstrual syndrome and menstrual disorder.Our analysis indicated that none of the menstruation disturbances listed above were among the most frequently reported adverse events for COVID-19 vaccines that we detailed before 12 .An updated adverse event report for COVID-19 vaccines is provided in Supplementary Table 5.However, some crosssectional studies reported high frequencies of these side effects.There could be many explanations for this discrepancy between adverse event databases and cross-sectional studies.Women tend to neglect seeking medical attention for what they perceive as a mild non-threatening short-term menstruation irregularity [53][54][55] .In fact, women participating in cross-sectional studies are unlikely to report changes to periods unless specifically asked 16 .
In fact, mining VAERS data for menstruation irregularities resulted in 35,386 adverse events that were not restricted to COVID-19 vaccines (Supplementary Table 6).The top five vaccines that had the highest share in these events were COVID-19 vaccine (26,714 events comprising 85.82%), human papillomavirus recombinant vaccine (1198 events comprising 3.85%), hepatitis B vaccine (1013 comprising 3.25%), trivalent influenza virus vaccine (581 events comprising 1.87%) and the zoster vaccine (566 events comprising 1.82%).It is noteworthy that all these vaccines are given later rather than the first few years of life, permitting adverse event reporting by menstruating women.
PubMed Abstract Sifter.To increase the confidence in the prioritized causal hits, we examined the relationship(s) between the prioritized genes and menstruation irregularities in more depth and complexity using the PubMed Abstract Sifter 56 .On the Landscape sheet we built queries that revealed the number of articles satisfying a variety of queries related to menstrual cycle, abnormal menstruation, vaccines, and prioritized upstream causal regulators.These results (Supplementary Table 7) revealed a sizable publication record for the relationship between menstruation and genes affected by COVID-19 vaccines (DEGs and/or causal hubs).The results of two queries consisting of the higher confidence list of prioritized causal transcription factors, and biomarkers causal hubs and DEGs are shown in Fig. 6a, b.
The Connectivity Map (CMap).The Connectivity Map analysis suggests that drugs capable of inducing transcriptomics effects opposite to those induced by mRNA vaccines could reverse vaccine side effects 57,58 .To identify small-molecule drugs that could prevent or reverse vaccine's side effects, we ranked all DEGs in response to vaccination by BNT162b2 according to their expression levels using log 2 FC values, to query the Connectivity Map database 59 .The CMap query gene signature consisted of the 50 most upregulated genes and the 50 most downregulated genes in response to vaccination with BNT162b2.Compound hits that produced opposite transcriptional signatures to the mRNA vaccine BNT162b2 are listed in Table 4.These compounds can reverse the transcriptomic signature of the vaccine, which will prevent or reduce side effects.In this study, we wanted to increase the confidence in the computational hypotheses derived from the enrichment and network analyses described earlier.

DISCUSSION
This study describes the first attempt to provide a mechanistic insight for vaccine-induced menstrual cycle irregularities.Our approach combined the analysis of vaccine gene expression profiles with upstream predictions of causal regulatory proteins and RNAs, and downstream analysis of enriched biological pathways to provide a causal mechanistic insight for vaccineinduced menstrual irregularities.
Our analysis led to the prioritization of topologically significant genes, such as transcription factors and important enzymes (i.e., kinases) that were largely missed in the gene expression profiling experiment, and therefore were not among the prioritized DEGs.We used the 'Causal Reasoning' methodology to identify candidate proteins (i.e., hypotheses) in the network that can be reached through a pre-defined distance (i.e., maximum shortest path length) from the DEGs.Thus, this analysis was crucial for the reconstruction networks responsible for vaccine-induced menstrual irregularities.The top five transcription factors, listed from highest confidence to lower confidence based on their prediction p-values, were: IRF1, STAT1, RELA, IRF3 and STAT 2 (Table 1).All  were predicted to be activated, in response to BNT162b2 vaccination, based on the directionality of differential gene expression in GS1 and GS2.
IRF1 was ranked first as the highest confidence predicted activated transcription factor.To assess whether changes in IRF1 activation can affect menstrual cycle, we checked whether IRF1 is biomarker for "menstrual cycle irregularity" but we didn't find evidence to support that.Next, we reviewed the biomedical literature to search for possible links between IRF1 and the menstrual cycle.We used the PubMed's advanced search using query terms "IRF1" and "menstrual cycle" and found evidence that IRF1 is upregulated by prolactin during the secretory phase of the menstrual cycle 59,60 .Additionally, evidence pointed that IRF1 upregulation in the endometrium was linked to prolactin and is localized predominantly to the granular epithelial cells 59 .Network reconstruction using PLR in addition to seed nodes IRF1, STAT1, RELA, IRF3 and STAT led to the direct interactions network in Fig. 3a.Downstream enrichment analysis in biological pathways, highlighted the prolactin signaling pathway as the most significantly enriched pathway with the six network seeds mentioned above.Thus, upstream causal reasoning followed by downstream pathways analysis highlighted a putative role for prolactin signaling in modulating post-COVID-19-vaccine adverse events on the menstrual cycle.Prolactin is a multi-functional molecule; it is a transcription factor hormone, secreted from the pituitary glands, and it regulates diverse biological functions including female menstruation [61][62][63][64][65][66][67][68][69] .For example, high prolactin levels can interfere with the production of sex hormones including estrogen and progesterone which can further impact menstruation regulation [61][62][63][64][65][66][67][68][69] .In fact, women who experience menstrual cycle irregularities often have higher prolactin levels than others, a condition known as hyperprolactinemia 47 .Hyperprolactinemia, is the most prevalent endocrine dysfunction of the hypothalamicpituitary axis in young females, accompanied with ovulatory disorder and leading to menstrual irregularities 70,71 .High levels of prolactin in the body prevent the release of (luteinizing hormone) LH and follicle-stimulating hormone (FSH), leading to ovulation disturbances 62,65,69 .Symptoms of hyperprolactinemia include long or irregular cycles, anovulation, amenorrhea, oligomenorrhea, polycystic ovarian syndrome or amenorrhea [72][73][74][75][76][77][78] .In fact, hyperprolactinemia can be caused by some drugs, stress, and some Fig. 6 Screenshots from PubMed Abstract Sifter.a Landscape sheet of the PubMed Abstract Sifter showing relationships, in the form of article counts, between biological concepts highlighted in this study.The first column "id" lists the gene symbols of prioritized top five causal transcription factors.b Landscape sheet of the PubMed Abstract Sifter showing relationships, in the form of article counts, between causal biological concepts highlighted in this study.The first column "id" lists the gene symbols of prioritized causal genes and vaccine-induced DEGs that are known as also biomarkers for menstrual cycle according to the CDDI database 46 .
Although we perceive menstrual changes as adverse events, prolactin-mimicking effects of vaccine are not necessarily a negative consequence of vaccines.Recently, prolactin has been suggested as a promising immunomodulator for the treatment of COVID-19 85 .However, we caution that prolactin mimicking effects may worsen auto-immune disease symptoms in patients suffering from systemic lupus erythematosus (SLE), multiple sclerosis, rheumatoid arthritis, psoriatic arthritis, and AIDS.Caution should be also practiced in patients undergoing organ transplantation.Elevated PRL levels have been reported in the previous conditions 86,87 .Furthermore, a recent study showed that prolactin hormones in addition to FSH and LH of healthy vaccinated males were higher than non-vaccinated males or COVID-19 male patients, indicating that changes in prolactin signaling are not limited to females 88 .Prolactin levels were 27.86 ± 4.35 ng/L in vaccinated males, 5.35 ± 1.59 ng/L in non-vaccinated males, and 16.65 ± 6.15 ng/L in COVID-19 male patients 88 .
To identify biomolecules that are implicated in menstrual changes, or the pathological processes that underlie the observed vaccine-induced menstrual symptoms, we filtered all predicted causal molecules and DEGs based on their overlaps with "menstruation irregularity"/"menstruation abnormality" diagnostic and prognostic biomarkers found in CDDI (Fig. 2a).We had four gene lists: 1) all DEGs in GS1 and GS2, 2) causal hubs for DEGs predicted for GS1, 3) causal hubs predicted for GS2, and 4) known diagnostic biomarkers for menstrual irregularity.TNF was identified as a high-confidence hit, i.e., a causal protein leading to the observed changes in gene expression and the predicted prolactin mimicking effects of the vaccine.TNF was an overlapping gene among four gene lists: 1) a DEG (log 2 FC = 3.07), 2) a causal key hub considering DEGs in both GS1, 3) a causal hub considering DEGs in GS2, and 4) a diagnostic biomarker for menstrual irregularity.Furthermore, our causal reasoning results predicted TNF-alpha activation in response to vaccination with BNT162b2.It should be noted that all causal predictions (Supplementary Tables 1 & 2) are based on experimental gene expression data.
Finally, SPIA results combined the enrichment results of DEGs with the actual amount of perturbation which highlighted the role of interferons on the signaling pathways influenced by BNT162b2.In fact, mRNA and vector-based COVID-19 vaccines result in the formation of neutralizing antibodies and activation of immune cells via the release of pro-inflammatory markers like cytokines and interferons 89 .There is evidence indicating that the treatment of multiple sclerosis with beta interferons causes menstrual irregularities associated with increased levels of luteinizing hormone (LH) and/or hyperprolactinemia 90 .Furthermore, the upregulation of interferon-gamma perturbs calcium signaling pathways which can in turn impact hormonal balance 12 .
But what is the relationship between prolactin signaling, TNFalpha activation and interferons?In fact, TNF-alpha activates the human prolactin gene promoter via NF-κB signaling 91 .TNF-alpha activation also stimulates the hypothalamic-pituitary-adrenal axis while suppressing the hypothalamic-pituitary-thyroid and gonadal axes, and growth hormone release 92 .Menstrual bleeding (menses) is known to be regulated by hypothalamic and pituitary hormones, and even the slightest changes in hormone levels, e.g., in response to medication or stress, can result in menstrual cycle abnormalities 93 .There is evidence that TNF-alpha and interleukin 1 beta (IL-1B), both are upregulated DEGs in this analysis, exert significant inhibitory effects on the GnRH-LH system in rats 94 , which may be the case in humans too.Moreover, the occurrence of reproductive disorders in poultry is highly correlated with the HPO axis and neuro-endocrine-immune network molecules, such as TNF-alpha and interferon-gamma (IFN-γ, IFNG) 95 .Thus, integrating enrichment and causal reasoning results with SPIA findings uncovered causal relationships between BNT162b2-induced menstrual changes and all the following pathways: prolactin signaling pathways, TNF-alpha activation, interferons the hypothalamic-pituitary-gonadal/ovarian/testicular axis.These results agree with previous studies suggesting that stabilizing the hypothalamic-pituitary-ovarian (HPO) axis with combined hormonal contraception reduces the chance of experiencing vaccine-associated menstrual changes 38,96 .
Different lines of supporting evidence increased the confidence in the derived causal hypothesis implicating menstrual changes with prolactin signaling, TNF-alpha and the HPO axis.First, VAERS data showed that post COVID-19 menstrual changes occurred in response all COVID-19 vaccines included in the databases including mRNA and vector-based vaccines and were not tied to the vaccine platform.The menstrual changes reported in VAERS included wide range of symptoms and were not limited to the length of menstrual cycle or menses period.Secondly, PubMed Abstract Sifter results highlighted 299,927 articles linking DEG TNF to any menstrual symptoms and 141 articles linking TNF to abnormal menses.Other high-confidence causal DEGs were progesterone receptor (PGR) with 45,560 and 796 articles linking it to any menstrual symptoms or abnormal menses subsequently, IL-1B with 56,099 and 43 articles linking it to any menstrual symptoms or abnormal menses subsequently.Finally, chemogenomics evidence from the CMap highlighted significant links to the HPO axis per results shown is Table 4.
It should be noted that the transcriptomics perturbations in response to treatment with BNT162b2 diminished on day 28 after receiving the second vaccine dose of BNT162b2.This suggests that vaccine effects on gonadal hormones, for females and males, and the predicted prolactin-mimicking effects, TNF-alpha In conclusion, our integrative computational network biology approach revealed that BNT162b2 can induce transcriptomics changes which may induce menstrual cycle changes by several mechanisms including prolactin-mimicking effects resulting from changes in interferon signaling and associated hormonal imbalance particularly in the HPO axis.This remains a high-confidence biological hypothesis supported by different lines of computational evidence derived from transcriptomics studies, causal reasoning analysis, downstream pathway enrichment results, and additional supporting evidence from vaccine adverse event databases (e.g., VAERS) and the biomedical literature.Further experimental validation is warranted to assess whether postvaccine prolactin-mimicking effects are due to increased levels of prolactin or due to other networking biology events mimicking prolactinemia.These effects may not be restricted to COVID-19 vaccines and should be assessed for other vaccines as well.
This study sheds the light on post-vaccine menstrual irregularity by revealing short-term post-COVID-19 vaccine prolactin mimicking effects resulting from the transcriptomics irregularities induced by COVID-19 vaccines.Most women associate menstruation irregularities with infertility which is one of the leading causes of vaccine hesitancy among females 97 .By providing a mechanistic insight into post-vaccine menstrual irregularities, this study is promised to correct misinformation about the relationship between vaccine-induced period irregularities and infertility.Thus, it is expected to decrease vaccine hesitancy.
This study establishes a causal relationship between COVID-19 vaccine and menstruation regulation by highlighting perturbed gene expression or dysregulated transcription of known diagnostic or prognostic biomarkers for menstruation and menstruation irregularities.Additionally, top scoring key hubs provide valuable hypotheses explaining gene expression and can be explored further in laboratory tests.
This study explores the causal links between COVID-19 vaccines and menstruation regulation based on an integrative bioinformatics approach that analyzed vaccine-induced transcriptomics irregularities.Integrating COVID-19 vaccine transcriptomics data with menstruation biomarkers, reinforced the selection of biologically relevant hypotheses from an overwhelming number of statistically significant hypotheses by increasing the confidence in computational hypotheses predicted by several methods.The fact that our computational hypotheses were supported by multiple lines of evidence is considered a major strength for this study.In fact, our integrative informatics workflow has several advantages over relying solely on conventional enrichment analyses for identifying the biological mechanisms that underlie vaccine side effects.Our approach integrates hypotheses derived independently from pathway and network enrichments, causal reasoning, SPIA, and the CMap to prioritize high confidence computational hypotheses predicted independently by various computational approaches and using different data types.The CMap, for example, is considered a unique chemogenomics database capable of connecting genes, drugs, and diseases based on genes expression similarities between polypharmacologic drugs and studied vaccines.This permits the prediction of vaccine side effects as well as underling causal mechanisms based on gene expression similarities with well-studied drugs.Finally, mining VAERS and PubMed for adverse event reports and vaccine-relevant data, serves as a validation step for the computationally-derived hypotheses.Thus, computational hypotheses prioritized using our integrative informatics approach are inherently high-confidence hypotheses with potentially improved clinical outcome.
Conversely, the applied methodologies or public databases have a few limitations that should be pointed out.First, reports from VAERS may not be conclusive or sufficient to establish causal relationships adverse events and specific vaccines.Due to the voluntary nature of VAERS reporting system, the information provided about an adverse event can be imperfect, imprecise, coincidental, or unconfirmed, limiting the scientific use of such reports 11,12 .Secondly, bioinformatics techniques relying on gene expression, pathway over-representation and network biology have some limitations and biases that we reviewed previously elsewhere 43 .Herein, the main limitation for the generalizability of the bioinformatics results to other COVID-19 vaccines, was the reliance on transcriptional data for the mRNA COVID-19 vaccine BNT162b2, which was the only publicly available COVID-19 transcriptomics data in humans at the time of conducting this research.As a result, our bioinformatics results apply directly to BNT162b2 or and may be extended to other COVID-19 mRNA vaccines (e.g., mRNA-1273) since COVID-19 mRNA vaccines share common features of the nature, strength, and timing of the immune responses as well as similar vaccine compositions 7,8,12,18,44 .The dosing regimens of vaccines may affect the results as well 30,89 .Our integrative workflow can be used to assess the safety of other vaccines using their transcriptional signatures in vaccinated individuals.

Systems biology informatics workflow
We have developed a network biology workflow to identify causal links between COVID-19 Vaccines and menstruation irregularities.This workflow (Fig. 1) incorporates three major components: (1) a module for mining and prioritizing gene signatures representative of a condition or a biological state; (2) a causal reasoning network module to identify upstream regulators of gene expression and (3) a pathway enrichment module to understand the biological processes regulated by DEGs and predicted causal regulators of gene expression.The resulting hypotheses are then evaluated based on existing evidence in vaccine reporting system databases and the biomedical literature.

Vaccine-induced differential gene expression
We searched the gene expression omnibus (GEO) [98][99][100][101] for transcriptional studies performed in response to treatment COVID-19 vaccines and we were able to identify one whole transcriptomics RNA-seq dataset (GSE169159) for COVID-19 vaccines in response to treatment with BNT162b2 at different time points.Our transcriptomics data analysis of GSE169159 raw data indicated that gene expression alterations from baseline were more prominent on day 22, which is the day after receiving the vaccine second dose.None of the genes analyzed at other time points (e.g., day 7, day 21, day, day 22.23, day 28) passed the applied thresholds for the prioritization of DEGs in this study (i.e., log 2 fold change (log 2 FC) of +2 or -2, and false discovery rate (FDR) 0.05.Therefore, we relied on differential gene expressions on day 22 for all our bioinformatics analyses.

Upstream transcriptomics analysis
Causal Reasoning 45,105 was used to identify upstream regulators (transcription factors, RNA molecules, kinases, phosphatases, and others proteins) that could cause/explain the observed postvaccine gene expression changes in transcriptomics experiments.We relied on MetaBase 106,107 as an interactions database, and the causal reasoning algorithm implemented in Clarivate's Key Pathway Advisor 103,108 .This method relies on a directed network which is annotated with activation and inhibition edges as well as biological mechanisms (transcription regulation).The significance of the predictions made by a particular hypothesis is assessed using a binomial test and calculating p-values as probabilities to get k successes in n predictions using binomial tests with pvalue = 0.50 according to the following equation: Here, k is the sum of correct predictions and n is the sum of correct and incorrect predictions.
Finally, p-values are assigned in the score matrix and hypotheses above the p-value threshold are filtered out of the score matrix.

Downstream pathway analysis
Pathway enrichment analyses were conducted in Cytoscape version 3.9.1 109 and MetaCore TM45 to interpret the consequences of vaccine-induced differential gene expression on biological processes.The significance of the enrichment was determined by calculating hypergeometric p-values 110 .All terms from the ontology were ranked based on their calculated p values.Ontology terms with p-values less than the p-value threshold 0.05 are defined as statistically significant and therefore relevant to the studied list of genes.All terms from the ontology were ranked according to their calculated p-values.
Signaling pathway impact analysis (SPIA) SPIA 111,112 was performed to identify the impact of the DEGs on the activity of the enriched pathway.This method aids in the identification of the most biologically relevant pathways and candidate causal genes.Herein, we identified perturbed pathways in response to vaccination by performing the enrichment analysis on the union gene list consisting of the experimentally derived DEGs in response to vaccination with BNT162b2, and the list of key hubs (e.g., activated, or inhibited proteins) using causal reasoning.
Vaccine adverse events database Raw data files were downloaded in comma-separated value (CSV) files from the CDC website 106,107 .CDC WONDER online search tool was used to mine VAERS data by vaccine type and symptoms 108 .The COVID-19 vaccines included in the databases were: BNT162b2, mRNA-1273 's and Janssen's 113,114 .

PubMed Abstract Sifter
The advanced literature retrieval tool PubMed Abstract Sifter was used to explore relationships between the biological concepts and molecular concepts that play roles in this research area.The steps in using the Abstract Sifter are described in the user guide.The tool and the user guide are available from the US EPA and downloadable from this webpage: https://comptox.epa.gov/dashboard/downloads 56,[115][116][117] .

The Connectivity Map (CMap)
The Connectivity Map analysis suggests that drugs capable of inducing transcriptomics effects opposite to those induced by mRNA vaccines could reverse vaccine side effects 57,58 .To identify small-molecule drugs that could prevent or reverse vaccine's side effects, we ranked all DEGs in response to vaccination by BNT162b2 according to their expression levels using log 2 FC values, to query the Connectivity Map database 59 .In fact, our transcriptomics data analysis of GSE169159 indicated that gene expression alterations from baseline were more prominent on day 22, which is the day after receiving the vaccine second dose.None of the genes analyzed at other time points passed the applied thresholds for identifying DEGs (i.e., log2 fold change (log2FC) of +2 or -2, and false discovery rate (FDR) 0.05.Therefore, we relied on differential gene expressions on day 22 for all our bioinformatics analyses.

Fig. 1
Fig. 1 Informatics systems biology workflow.A devised workflow for studying the mechanism(s) underlying the biological effects of vaccines.

Fig. 3
Fig. 3 Overlapping DEGs, causal upstream hubs and biomarkers.a Venn diagram showing overlaps between DEGs, predicted causal upstream regulatory hubs using DEGs in GS1 and GS2, and menstrual irregularity biomarkers.b Venn diagram showing overlaps between DEGs, predicted causal upstream regulatory hubs using DEGs in GS1 and GS2, and prolactin signaling biomarkers.

Fig. 4
Fig. 4 Interconnectivity between prioritized high confidence transcription factors.a Direct interactions network of five higher confidence causal transcription factors.b Direct interactions network of five higher confidence causal transcription factors in addition to prolactin (PRL).c Direct interactions network of 33 causal upstream regulators that are known biomarkers for menstrual disturbances.d Direct interactions network of 33 causal upstream regulators that are known biomarkers for menstrual disturbances, in addition to prioritized 5 topological genes and PRL.Thick edges correspond to confidence score !0.70 (i.e., high confidence score), while the thin edge corresponds to a confidence level 0.50 (i.e., low confidence score).Nodes are color-coded using a split pie chart coloring scheme indicating pathway/gene set contribution to each node from the top 5 most enriched pathways/gene lists.FDR values represent he significance of the predicted pathway.Generated based on STRING data on 27 September 2022.

Fig. 5
Fig. 5 Prolactin signaling pathway.a Prolactin signaling pathway map.A node (or object) on the map could be a gene, protein or chemical compound.Query genes from experimental data which intersect with pathway protein or chemical compound.Query genes from experimental data which intersect with pathway objects are denoted by thermometers.Thermometer 1 represents causal transcription factors.Thermometer 2 represents DEGs in response to treatment with vaccine, applying thresholds of log 2 FC ≥ 2.00 or ≤ −2.00, and FDR ≤ 0.05.b Biological processes involved in prolactin signaling pathway.The % refers to the percentage of network objects in the pathway map.The p-value is the process prediction p-value.

Table 1 .
Top commonly predicted upstream regulatory hubs at a distance of 1 for DEGs in GS1 and GS2.
a Predicted activity of the key hub by causal reasoning is denoted byif the hub is inhibited, and denoted by + if the hub is activated.b

Table 2 .
List of menstrual irregularity biomarkers which overlap with BNT162b2 vaccine-induced DEGs and/or predicted causal proteins.

Table 3 .
List of prolactin biomarkers which overlap with BNT162b2 vaccine-induced DEGs and/or predicted causal proteins.

Table 4 .
Small-molecule drugs and chemical compounds that regulate gene expression in an opposite manner to BNT162b2.