Candida pathogens induce protective mitochondria-associated type I interferon signalling and a damage-driven response in vaginal epithelial cells

Vaginal candidiasis is an extremely common disease predominantly caused by four phylogenetically diverse species: Candida albicans; Candidaglabrata; Candidaparapsilosis; and Candidatropicalis. Using a time course infection model of vaginal epithelial cells and dual RNA sequencing, we show that these species exhibit distinct pathogenicity patterns, which are defined by highly species-specific transcriptional profiles during infection of vaginal epithelial cells. In contrast, host cells exhibit a homogeneous response to all species at the early stages of infection, which is characterized by sublethal mitochondrial signalling inducing a protective type I interferon response. At the later stages, the transcriptional response of the host diverges in a species-dependent manner. This divergence is primarily driven by the extent of epithelial damage elicited by species-specific mechanisms, such as secretion of the toxin candidalysin by C. albicans. Our results uncover a dynamic, biphasic response of vaginal epithelial cells to Candida species, which is characterized by protective mitochondria-associated type I interferon signalling and a species-specific damage-driven response. Dual RNA sequencing and a time course infection model of vaginal epithelial cells with four Candida pathogens reveal a homogenous protective mitochondria-associated type I interferon host response at early time points, followed by distinct species-specific pathogenicity patterns.

V ulvovaginal candidiasis (VVC) is among the most common fungal infections, affecting 70-75% of women at least once in their lifetime 1 . VVC is characterized by acute inflammation of the vaginal mucosa due to the overgrowth of normally commensal Candida species [2][3][4] . Although Candida albicans is the predominant cause of VVC, the prevalence of species like Candida glabrata, Candida parapsilosis and Candida tropicalis has increased (reviewed in Makanjuola et al. 5 ). Despite the shared genus name, these species are phylogenetically diverse and often have non-pathogenic close relatives, indicating that their ability to infect humans has emerged independently 6 . How these diverse Candida species interact with host cells has rarely been addressed on a comparative basis. Improved knowledge of similarities and species-specific characteristics of infection processes is crucial to understand the pathogenesis, improve diagnostics and therapy of candidiasis 7 .
Research on host-fungi interactions has mainly focused on immune cells, which are considered crucial players in the defence against fungal infections 8,9 . However, epithelial cells play a fundamental role in shaping the host defence against fungi, which goes beyond their function as a physical barrier [10][11][12][13][14][15] .
Studies in infection biology often focus on either the pathogen or host response, yet microbial pathogenesis can be best interpreted in the framework of dynamic host-microbe interactions. Dual RNA sequencing (RNA-seq) enables the combined assessment of the transcriptional responses of host and pathogen 16,17 and provides insights into the interactions of fungal pathogens with different host cells [17][18][19][20][21][22][23] .
To elucidate general and species-specific interactions between vaginal epithelial cells and the four most prevalent VVC-causing Candida species, we applied dual RNA-seq and an in vitro infection model. Our experimental design allows pathogens to deploy their arsenal of pathogenicity factors without restriction by the immune system. Furthermore, it facilitates specific investigation of epithelial recognition and defence mechanisms, which constitute the first line of defence against infecting fungi.
Our results reveal that fungal transcriptomes show species-specific patterns during infection, probably reflecting the independently evolved pathogenic potential of Candida species. Vaginal epithelial cells display a biphasic response: an early protective type I interferon (IFN) response, which is mediated by sublethal

Species-specific transcriptional responses to epithelial cells.
Subsequently, we investigated whether these differential pathogenicity patterns were reflected in the transcriptional responses of both fungal and epithelial cells (see the experimental set-up in Extended Data Fig. 1 and cross-mapping analysis between human and yeast sequencing reads in Supplementary Files 1-4). The transcriptional dynamics of each Candida species throughout the infection were analysed (Fig. 2). All species induced rapid transcriptional responses following infection with increasing numbers of differentially expressed genes over the course of infection ( processes and stress response pathways. C. parapsilosis upregulated, among others, genes related to iron transport, ribosome assembly and translation. In contrast, C. tropicalis differentially expressed genes were mainly related to RNA processing, ribosome biogenesis and ergosterol biosynthetic processes. At later stages, similar functional enrichments were observed across the species. The GO terms oxidation-reduction process, fatty acid beta-oxidation, iron homeostasis and acetate catabolism were enriched in at least three species throughout the infection. When comparing differentially expressed genes across species, a remarkably distinct pattern was observed for each pathogen (Extended Data Fig. 2). Analysis of the distribution of species-specific (without orthologues in the other species), partially shared (with orthologues in one or two of the other species) and fully shared (1-to-1 orthologous genes in all species) differentially expressed genes (Extended Data Fig. 2a) revealed that species-specific and partially shared genes constitute a substantial proportion (31-72%). Moreover, species-specific genes are more likely to be differentially expressed than fully shared genes (chi-squared test P < 0.05, except for C. tropicalis). Even orthologous genes present in all four species showed species-specific differential expression (Extended Data Fig. 2b) C . a lb ic a n s C . g la b r a ta C . p a r a p s il o s is C . tr o p ic a li s C . a lb ic a n s C . g la b r a ta C . p a r a p s il o s is C . tr o p ic a li s C . a lb ic a n s C . g la b r a ta C . on orthologous gene expression showed species-specific clusters (Extended Data Fig. 2c). Gene coexpression analysis was used as an independent approach to investigate commonalities and differences of fungal transcriptional responses. By constructing host-pathogen interaction coexpression networks, highly interconnected gene clusters (modules) were defined and their biological functions were inferred by GO term analysis. In each fungal infection scenario, we detected numerous modules of coexpressed genes (22-28 modules; Supplementary File 10).
Based on fully shared genes, we then assessed whether the fungal genes in the coexpressed modules were conserved across species. Distinct modules were observed for each Candida species with few shared genes. On average, only 5% of orthologous genes were shared between any modules of different species (Extended Data Fig. 3). Therefore, the genes in the coexpression modules functionally showed a large species specificity (Supplementary File 10) with few exceptions. The modules with the highest similarity, that is, module3 in C. glabrata and module3 in C. tropicalis, were both enriched for genes associated with DNA replication. Interestingly, we observed modules related to adhesion in C. albicans (module11) and C. tropicalis (module13), respectively, possibly related to shared virulence features of these two species.
Infection-specific differentially expressed genes of Candida species. Comparisons of C. albicans gene expression during infection of oral epithelium or vascular endothelium, and growth in the tissue culture medium, revealed that only a fraction of genes were specifically expressed during interactions with host cells 20 . This indicates that most of the genes induced during interaction with the host are also required for growth in culture media.
To investigate whether such a phenomenon also occurs during interaction with vaginal cells, controls of fungal cells grown in culture medium only were investigated (Extended Data Fig. 1 and Supplementary Table 1). A subtraction of the differentially expressed genes in medium from those expressed on epithelial cells revealed a large overlap between the differentially expressed genes during infection and in culture medium ( Fig. 4 and Supplementary Files 5-8). Infection-specific genes are mostly species-specific (Extended Data Fig. 4) and GO enrichment analysis identified different functional enrichments depending on the species (Supplementary File 9).
These were characterized by genes involved in mitochondrial electron transport and ATP synthesis for C. glabrata, ergosterol biosynthesis, sulphite and manganese ion transport for C. parapsilosis and ribosome biogenesis and ribosomal RNA processing for C. tropicalis. No functional enrichments were identified in the upregulated, infection-specific C. albicans genes, yet the downregulated genes showed enrichment in three GO terms related to white-opaque phenotypic switching. These downregulated genes include WOR1, a master regulator inducing the less virulent opaque state 24 .
In summary, distinct transcriptional patterns for each yeast species during infection were observed, suggesting highly species-specific strategies to cope with epithelial cells.
Epithelial transcriptomic responses to Candida species. To shed light on how the Candida species-specific pathogenicity and transcriptional patterns influence the host response, epithelial transcriptome responses to infections were analysed (Fig. 3a,b). The epithelial cell transcriptome dynamics showed a bias towards upregulation of genes at the initial stages of infection (Fig. 3a), which is consistent with previous findings 25 .
When compared against the total number of differentially expressed genes (Fig. 3c, top), the proportion of common differentially expressed genes induced by infection with any of the four species decreased throughout the time course-8.8% of common differentially expressed genes at 3 h post-infection, and 7.6 and 6.4% at 12 h and 24 h post-infection, respectively. A similar pattern is observed when comparing common differentially expressed genes to the differentially expressed genes induced specifically by each fungal species (Fig. 3c, bottom). The larger fraction of shared differentially expressed genes at the early time points suggests that the response to the different yeast species is more conserved at the early infection stages while increased species specificity is observed at the later stages (Fig. 3b,c).
PCA analysis of the gene expression of epithelial cells revealed a similar pattern (Fig. 3d). The tight clustering at 3 h post-infection indicates that epithelial cells exhibit a uniform transcriptional response to the four Candida species at the early stages of infection, in contrast to the fungal transcriptional profiles. However, the epithelial transcriptomes diverge from 12 h post-infection onwards; at 24 h post-infection we observed three distinct clusters of responses to the different species. The transcriptional response of epithelial cells to C. tropicalis and C. glabrata showed high similarity, being different from the responses to C. albicans or C. parapsilosis. Functional GO term enrichment analysis revealed a similar trend: at the early time point, GO terms associated with mitochondrial processes are enriched for all species and of type I interferon (IFN) responses for all species except C. tropicalis. At the later stages, species-specific terms appeared (Fig. 3e).
Based on these results, we decided to unravel the basis of the two observed key phenomena: (1) the uniform early transcriptional response related to mitochondria and type I IFN signalling; and (2) the divergence of the host transcriptome response at the later stages of infection.
Uniform early responses to Candida infections. The upregulation of genes associated with the respiratory electron transport chain (Fig. 3e) in epithelial cells at the early time points with all Candida species included induction of all mitochondrial genes (Fig. 4a). The response was dependent on viable fungi being in direct contact with the epithelial cells (RNA-seq, Extended Data Fig. 6 and Supplementary Files 5-8; quantitative PCR with reverse The numbers on the plots indicate the number of differentially expressed genes (log 2 (fold change) >1.5, P adj < 0.01, upregulated (red), downregulated (blue)). b, Venn diagrams showing similarities and differences of human differentially expressed genes in response to Candida species. c, Proportion between shared differentially expressed genes and the total number of differentially expressed genes (top) and proportion between shared differentially expressed genes and the differentially expressed genes induced exclusively by each fungal species (bottom). d, PCA biplot of all analysed human samples. The labels of the data points correspond to sample identifiers, where 'reseq' indicates that the sample was sequenced more than once (Supplementary Table 1 for further details). e, GO term enrichment analysis for upregulated genes (category 'biological process') of the host at different time points. The x axis indicates the infecting Candida species. Only significant (P adj < 0.05) GO enrichments are shown. Differentially expressed gene analyses were done using DESeq2 v.1.26.0; comparisons against time point 0 were done using a two-sided Wald test. GO enrichment analysis was done using clusterProfiler v.3.14.3, which was used to perform the hypergeometric test. Adjustments of P values for differentially expressed gene and GO enrichment were done by Benjamini-Hochberg procedure.    Host mitochondria have recently been identified as hubs of the innate immune responses 26,27 . In particular, mitochondrial signalling is known to activate type I IFN signalling pathways 28 . We observed enrichment of GO terms (Fig. 3e) and upregulation of interferon-stimulated genes (ISGs) associated with the type I IFN response 29 on exposure of epithelial cells to the four Candida species (Fig. 4a). Since type I IFN responses are implicated in antiviral host defence, additional metagenomic analyses were performed and no viral contamination was detected.
The connection between ISG expression and mitochondrial functions of the host was characterized during early Candida infection. Morphological changes of the mitochondrial network in epithelial cells were observed, changing from reticular (uninfected) to fragmented (infected), and an accumulation of mitochondria around the nucleus (Fig. 4b). Some mitochondria in infected epithelial cells lost their integrity and changed their shape but not in uninfected cells (Fig. 4c,d). Additionally, endoplasmic reticulum regions surrounded these altered mitochondria, suggesting mitophagy of damaged mitochondria. Interestingly, mitochondria were also localized frequently around the invading hyphae of C. albicans (Fig. 4c).
Mitochondrial membrane potential (ΔΨm), a key indicator of mitochondrial health, indicated depolarization in epithelial cells infected with any of the four species (Fig. 4e). This change in ΔΨm is associated with the production of mitochondrial reactive oxygen species (mtROS), which are critical players in the regulation of immune signalling pathways 30 . Epithelial mtROS levels were increased on infection with all four species (Fig. 4f). Finally, release of mitochondrial DNA (mtDNA) into the cytosol was observed during infection with all Candida species (Fig. 4g). The release of mtROS and mtDNA was not detected with killed C. albicans cells or when contact was restricted using a transwell system (Fig. 4e,f). This supports our notion that viable Candida cells in direct contact with epithelial cells induce mitochondrial dysfunctions at both transcriptional and biochemical levels (Fig. 4h).
During infection with bacteria or viruses, host mitochondria can release mtDNA into the cytosol acting as a damage-associated molecular pattern (DAMP) that activates immune pathways 31 . Cytosolic mtDNA can bind the DNA sensor cyclic GMP-AMP synthase and promote stimulator of interferon genes (STING)-IRF3-dependent signalling to induce a type I IFN response 28 . In line with this, depletion of epithelial mtDNA ( Fig. 5a) prevented upregulation of ISGs after Candida infections (Fig. 5b). In addition, transfection of uninfected epithelial cells with amplified mtDNA induced ISG expression (Fig. 5c). Transfection of uninfected epithelial cells with total DNA from epithelial cells only induced ISGs when the transfected DNA contained mtDNA (Fig. 5d), which supports the role of mtDNA in the induction of ISG expression (Fig. 5e).
Although mitochondrial dysfunction is a hallmark of cellular apoptosis, no apoptotic or necrotic epithelial cells were observed during the early stages of infections (Extended Data Fig. 5). This was expected since mitochondrial dysfunctions were only transiently observed. Later during infection, we observed necrotic cell death but no increase in apoptosis compared to the uninfected control (Extended Data Fig. 5). The A-431 cell line lacks functional p53, an important apoptosis inducer 32 . However, we observed similar mitochondrial depolarization in primary vaginal cells on Candida infections, while apoptosis levels did not differ between infected and uninfected cells (Extended Data Fig. 5). Additionally, treatment of A-431 cells with the apoptosis inducer staurosporine excluded that the observed mitochondrial signalling and induction of ISG expression were related to apoptosis (Extended Data Fig. 5).
The function of many ISGs is poorly characterized but induction is associated with protection against viral infections 33 . To gain insights into their role during Candida infection, selected ISGs (IFI6, MX2, CMPK2) were silenced in epithelial cells before infection with C. albicans. IFI6 was selected since it was previously observed to be induced by C. albicans 34 , MX2 was among the most highly upregulated common genes (this study) and CMPK2 encodes a protein with mitochondrial localization 35 . The level of epithelial damage was increased once these ISGs were silenced (Fig. 6a). While stimulation of cells with IFN-β (0.1 ng ml −1 ) before infection resulted in reduced damage (Fig. 6b); blocking IFN-α/β receptor (IFNAR) signalling led to increased damage (Fig. 6b). These data illustrate that type I IFN signalling increases epithelial resistance to Candida infection.
Neutrophil recruitment and activation is a hallmark of vaginal candidiasis 36 . Pro-inflammatory mediators, such as IL-6 or IL-1β, which are characteristic for these events, were not produced during Candida infection (Fig. 6c), while IL-1α and IL-8 levels increased (Fig. 6d), which correlates with the level of damage. IL-8 levels also increased on IFNAR blocking ( Fig. 6d), suggesting a relationship between epithelial type I IFN signalling and the initiation of pro-inflammatory responses, which can drive immunopathology in VVC.
When exposing neutrophils to the culture supernatants of epithelial cells infected with Candida species, IL-8 release was observed (Fig. 6e). IL-8 levels further increased when the IFNAR receptor was blocked in the same setting (Fig. 6e). This suggests that type I IFN signalling in epithelial cells restricts pro-inflammatory responses and subsequent neutrophil activation (Fig. 6f).

Damage-driven transcriptional responses.
While the initial phases of infection showed a conserved epithelial response, this response separated into different trajectories at later stages. Considering that host cell damage is a major determinant of pathogenicity [37][38][39] , we hypothesized that Candida species-specific differences in the damaging potential ( Fig. 1  . h, Schematic model of the events associated with mitochondrial dysfunctions resulting in mtROS production and mtDNA release. All data are derived from n = 3 independent experiments, unless indicated otherwise (d). Representative microscopy images (b-c) were taken from n = 3 biological replicates and similar results were observed. All values are presented as the mean ± s.d. relative to the uninfected (−) control (dotted lines on (d-g)). Statistical significance is indicated as *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001, ****P ≤ 0.0001 (Kruskal-Wallis test with two-sided Dunn's multiple comparison (d) and one-way ANOVA with Dunnett's multiple comparisons test (e-g). Credit: graphics in h adapted from Servier under a Creative Commons licence CC By 3.0. Data Fig. 6), reflected by the pattern observed in the PCA plot (Fig. 3d), drive the different transcriptional responses during late infection.

and Extended
Epithelial cell damage during C. albicans infection is mediated by the cytolytic toxin candidalysin 40,41 . Deletion of the ECE1 gene, which encodes candidalysin, renders C. albicans almost unable to To determine whether candidalysin-driven epithelial cell damage might dictate the transcriptional response of epithelial cells, the transcriptional response on interaction with the C. albicans ece1Δ/Δ mutant was investigated. The epithelial transcriptional response to the candidalysin-deficient mutant was notably similar to the response to the non-damaging species C. parapsilosis (Extended Data Figs. 6 and 7). This confirms a pivotal role for host cell damage as the major driver of epithelial transcriptional responses to Candida infections. GO term enrichment analysis of 774 genes, specifically upregulated on exposure to damaging wild-type C. albicans, showed no significant enrichment for any process. However, previous studies demonstrated that candidalysin induces c-Fos and mitogen-activated protein kinase-driven release of the pro-inflammatory cytokines IL-1α and IL-1β and the chemokine IL-8 in vaginal epithelial cells 12,25 . Manual inspection revealed similar responses including upregulation of HBEGF, CXCL1, CXCL2, IL1A, IL1B, CXCL8 and CSF2 and genes associated with the 'danger' response pathway FOS, JUN and DUSP1 (refs. 12,25 ). This confirms that epithelial damage and pro-inflammatory signals that drive neutrophil recruitment induced by C. albicans depend almost exclusively on candidalysin.

Discussion
In this study, we dissected the interaction of the main four Candida species that cause VVC with human vaginal epithelial cells. Large-scale dual transcriptomic analysis of human and fungal cells during the course of infection revealed common and species-specific Candida pathogenicity patterns. We observed a biphasic host response to Candida species, which is characterized by an early common mitochondria-induced type I IFN signalling and diverged responses at later stages depending on the species-specific capacities to inflict damage to the vaginal epithelial cells. It has previously been hypothesized that phylogenetically diverse Candida species independently acquired their ability to colonize and infect humans and thus are expected to use distinct sets of pathogenicity mechanisms 6,23 . Our study empirically supports this hypothesis by showing that the main VVC pathogens express species-specific transcriptional responses and pathogenicity patterns on contact with vaginal epithelial cells, even when only orthologous genes were considered.
Similarly, the epithelial transcriptional responses at the late stages of infection were specific depending on the Candida species. These diverse patterns paralleled the varying damaging capacities of the four Candida species. We confirmed fungal-induced damage as the major driver of epithelial responses by infecting epithelial cells with the non-damaging candidalysin-deficient C. albicans ece1Δ/Δ mutant. The epithelial transcriptional response to this mutant did not resemble the response to wild-type C. albicans, but instead was similar to the response to the non-damaging species C. parapsilosis at the late stages of infection. This finding confirms the crucial role of candidalysin during interaction of C. albicans with vaginal epithelial cells leading to DAMP release that can catalyse immunopathology during vaginal infections 25 .
In contrast, the epithelial response towards the different Candida species was highly uniform at the early stages. This initial response was driven by common epithelial processes rather than by convergent activities (such as virulence programmes) of the tested Candida species. For example, independent of their viability, DOCK8 expression was upregulated after infection with any of the four Candida species (Supplementary Files 5-8). Although DOCK8 has multiple   signalling functions, it was suggested that it promoted immune responses to diverse external stimuli 42 . Several studies associated DOCK8 with mucocutaneous candidiasis due to impaired T H 17 differentiation [43][44][45] . Thus, it is tempting to speculate that DOCK8 may regulate the recognition of Candida species by epithelial cells.
Mitochondria-encoded genes and genes associated with the type I IFN response pathway were uniformly upregulated. While type I IFN responses are associated with viral infections 46 , type I IFN responses have been observed in peripheral blood mononuclear cells infected with C. albicans 34 . Additionally, type I IFN responses were recently shown to dysregulate host iron homeostasis and enhance C. glabrata infection 47 ; the type I IFN-inducing RIG-I helicase MDA5 has been associated with systemic and chronic mucocutaneous candidiasis 48 . Finally, IFNαR1 signalling is crucial for efficient host defence against systemic candidiasis in mice 49 .
Our data show that type I IFN signalling, induced by vaginal epithelial cells in response to Candida species, increases epithelial resistance to infection and dampens pro-inflammatory responses.
Such an immune response may be relevant in host niches colonized by commensal microbes that need to be tolerated without induction of inflammation. For example, intestinal epithelial cells regulate the durability and specificity of immune responses and guide the immune system to differentiate between commensal and pathogenic microbiota via expression of type I IFN and ISGs [50][51][52] . Since Candida species are commensals of the vaginal mucosa 53 , the epithelial type I IFN pathway may maintain the threshold between commensalism and pathogenicity and regulate antifungal immunity. Our results show a protective role for ISGs and IFNAR signalling in increasing epithelial resistance to Candida-induced damage while reducing potentially detrimental pro-inflammatory responses. Supporting this, Li et al. 54 showed that administration of human IFNα-2b decreased the inflammation and vaginal epithelial damage in a rat VVC model. The combined effects, immunomodulation and epithelial antifungal resistance, may be crucial to restrict Candida species to commensalism and avoid inflammation-driven pathology. Moreover, this highlights the type I IFN response as a potential target for host-directed therapy aimed at improving epithelial resistance and preventing immunopathology. Such a therapy could benefit VVC patients that fail to recover with antifungal treatment alone, a phenomenon often observed in women with recurrent VVC 55 .
At the early infection stage, genes encoded by mtDNA, in particular genes coding for the respiratory electron transport chain, were upregulated. Apart from the well-established roles in metabolism and energy production, mitochondria are central hubs in innate immunity 26,27 . Mitochondrial dysfunction, resulting in mtROS and mtDNA release into the cytosol, can act as a DAMP and activate various signalling pathways [26][27][28][56][57][58] , including induction of cytokine production 59 and type I IFN responses 60 . Intriguingly, altered mitochondrial function at the early stages of infection were observed in different host cell types on infections with various bacterial pathogens 61-64 , including Chlamydia trachomatis 65 , Chlamydia pneumoniae 66 , Listeria monocytogenes 67 and the parasite Toxoplasma gondii 68 . However, these mechanisms have not yet been observed during fungal infections. We observed that mitochondria in Candida-infected vaginal epithelial cells changed shape and lost integrity, had decreased membrane potential and released mtROS and mtDNA. The release of mtDNA was observed to act as a DAMP that activates the type I IFN pathway during Candida infections of the vaginal epithelium. This activation may potentially occur through the STING pathway, as shown previously for Streptococcus pneumoniae 69 , at the level of post-translational modifications 70 .
To maintain epithelial integrity and mount an effective epithelial host defence while preventing detrimental inflammatory responses, such a mitochondrial response must occur on a sublethal level 59 . Accordingly, no changes in apoptosis were observed over the course of infection. The non-lethal mitochondrial dysfunction was also independent of necrosis, which was only observed at the later stages. Likewise, we observed consistent activation of the type I IFN pathway, which is suppressed by apoptotic caspases 71 . Consequently, induction of apoptosis abrogated expression of ISGs. Similar studies with L. monocytogenes and T. gondii showed that mitochondrial dysfunction was uncoupled from the apoptotic pathway 67,68 . During infection of epithelial cells with diverse microbes, the mitochondrial apoptosis apparatus can be activated at a low level, which is insufficient to induce apoptosis 59,72 . This phenomenon has been termed limited mitochondrial outer membrane permeabilization (MOMP), or 'minority MOMP' , and induces pro-inflammatory cytokine production via STING.
Viruses, bacteria and parasites can induce minority MOMP, thereby contributing to cytokine release during infection 59 . We propose that this mechanism plays a significant role in epithelial sensing of Candida species, induction of epithelial antifungal immunity and modulation of immune responses via type I IFN signalling.
It remains to be determined how Candida species initiate the mitochondria-induced epithelial type I IFN response. We observed mitochondrial signalling at stages when C. albicans had not yet invaded or damaged epithelial cells, whereas all other Candida species failed to invade epithelial cells. Therefore, we propose that the induction of mitochondrial signalling may rely on sensing of pathogen-associated molecular patterns (PAMPs).
In summary, we identified species-specific pathogenicity patterns of Candida species infecting vaginal epithelial cells, which are reflected at the transcriptional level during the course of infection. In contrast, vaginal epithelial cells exhibit a conserved response at early stages, but a diverse, damage-driven response at later stages. The conserved response was characterized by non-lethal mitochondrial signalling, which induced a type I IFN response that protects against Candida-induced damage and modulates pro-inflammatory responses. This acts as a common pathway of host-pathogen interactions between vaginal epithelial cells and Candida pathogens.

Methods
Fungal strains and culture conditions. C. albicans SC5314 (ref. 73 ), Candida glabrata ATCC 2001 (obtained from ATCC), C. tropicalis DSM 4959 (obtained from the German Collection of Microorganisms and Cell Cultures), C. parapsilosis 73-037 (ref. 74 ) and C. albicans ece1Δ/Δ 40 were used in this study. For all experiments, single colonies were picked from yeast extract peptone dextrose (YPD) agar plates and grown overnight in liquid YPD medium in an orbital shaker at 180 r.p.m. at 30 °C (C. albicans, C. tropicalis and C. parapsilosis) or 37 °C ( C. glabrata). Yeast cells were then collected by centrifugation (20,000g, 1 min), washed twice with PBS and adjusted to 2 × 10 6 yeast cells per ml −1 .

In vitro vaginal epithelial infection model. To mimic the vaginal epithelium,
A-431 epithelial cells (ACC 91) were used. These cells are derived from a vulva epidermoid carcinoma and routinely used to model the vaginal mucosa 75,76 . A-431 cells were authenticated using short tandem repeat analysis (DNA fingerprinting) and routinely tested for the absence of Mycoplasma contamination. Epithelial cells were cultivated in Roswell Park Memorial Institure (RPMI) 1640 medium (Thermo Fisher Scientific) supplemented with 10% fetal calf serum (FCS; Bio&SELL) in a humidified incubator at 37 °C and 5% CO 2 . For infection, epithelial cells were seeded in 6-well plates (3 × 10 5 cells per well) and cultured for 2 d. On the day of infection, the medium in each well was replaced with 1.5 ml of RPMI 1640 without FCS and incubated for 30 min to allow cells to adjust to the change of medium. In subsequent bioinformatics analyses, we considered the control samples at 30 min in this medium as the 0 h time point. Epithelial cells were subsequently infected with Candida cells (1.5 ml of 2 × 10 6 yeast ml −1 in RPMI 1640 without FCS) and incubated at 37 °C and 5% CO 2 . Samples for RNA isolation were collected at different time points: 3, 12 and 24 h post-infection. More specifically, the well content was removed and replaced with 500 μl of RLT buffer (QIAGEN), containing 1% β-mercaptoethanol (Roth). Cells were detached using a cell scraper (<3 min), immediately shock-frozen in liquid nitrogen and stored at −80 °C until further use (see 'RNA isolation and pooling'). As controls, Candida cells alone and epithelial cells alone were incubated for 30 min (0 h control: C0) and 24 h (24 h control: C24) and samples for RNA isolation were collected as described above.
RNA isolation and pooling. Collected samples were defrosted on ice and centrifuged for 10 min (20,000g, 4 °C). The supernatant was transferred to a new microcentrifuge tube and used to isolate human RNA (RNeasy Mini Kit; QIAGEN), according to the manufacturer's instructions. Fungal RNA was isolated from the pellet using a freezing-thawing method, as described previously 77 . Both human and fungal RNA concentrations were quantified using a NanoDrop 1000 Spectrophotometer (Thermo Fisher Scientific) and RNA quality was assessed with an Agilent 2100 Bioanalyzer (Agilent Technologies). To subsequently achieve sufficient sequencing depth of both counterparts for a robust differential gene expression analysis 78 , corresponding fungal and human RNA samples were pooled in a 2:3 ratio by weight for further library preparation and sequencing. Before pooling, we first checked whether this strategy would result in ambiguous read mappings between the host and pathogen data after sequencing and data analysis. To assess the rates of cross-mapping, that is, reads originated from human but mapped to fungi and vice versa, which can bias expression level quantifications, we used CROSSMAPPER v.1.1.0 (ref. 79 ), which simulates reads from multiple reference genomes/transcriptomes, maps the data back to the concatenated reference sequences and reports the rates of cross-mapping. We used the 'RNA' mode of CROSSMAPPER and simulated and back-mapped 20 and 40 million 2 × 50 and 2 × 75 reads for each fungal species and human, respectively. In all cases, the pooling and sequencing strategy resulted in virtually no cross-mapping between human and yeast data ( Supplementary Files 1-4).
Growth curves. Candida cells were adjusted to 10 6 yeast ml −1 either in YPD or RPMI medium. Growth was monitored in 96-well-plates by measuring the absorbance at 600 nm every 30 min for 24 h at 37 °C in a microplate reader (Tecan M-Plex). Before each measurement, plates underwent 10-s orbital shaking followed by 10-s waiting time. The OD 600 values were converted into log 2 and the generation time was calculated from the slope of the exponential growth phase. The experiment was repeated five times.
Ultraviolet killing of Candida. Candida cells from overnight cultures were collected by centrifugation, washed twice with PBS and adjusted to approximately 5 × 10 7 yeast ml −1 in PBS. The suspension was transferred to a Petri dish as a thin liquid layer (10 ml) and exposed to 4 doses of 100-120 mJ per cm 2 in an ultraviolet (UV) CROSSLINKER (CL 508S; Uvitec). The efficiency of UV killing was evaluated by plating 50 µl of the sample onto YPD agar and incubated for 48 h at 30 °C.
Adhesion assay. Epithelial cells were infected with Candida yeast cells as described above and incubated for 1 h. Non-adherent Candida cells were removed by rinsing with PBS. Subsequently, epithelial cells with adhered Candida were fixed with Roti Histofix 4% (Roth). Adherent Candida cells were stained with Alexa Fluor 647 conjugate of succinylated concanavalin A (Invitrogen) and visualized using a fluorescence microscope (Leica DM5500B; Leica DFC360 FX). Pictures of each sample were taken until a total of 100 adherent cells were counted. Adhesion was calculated based on the average of Candida cells counted in each picture with a defined area. This number was expressed as a percentage of adhered cells versus inoculated cells 77 .
Invasion assay and hyphal length. Epithelial cells were infected with Candida cells as described above and incubated for 3 h. Non-adherent Candida cells were removed by rinsing with PBS and samples were fixed with Roti Histofix 4%. Extracellular, non-invasive fungal components were stained by concanavalin A. After rinsing with PBS, epithelial cells were permeabilized in 0.5% Triton X-100 for 10 min. Next, fungal cells were stained with Calcofluor white (Sigma-Aldrich) and visualized by fluorescence microscopy. The total hyphal length was noted as well as the percentage of invasive hyphae (only Calcofluor white-stained), counted from at least 100 hyphae per strain for each biological replicate.
Epithelial damage assay. Epithelial cells were infected with Candida cells as described above and incubated for 24 h. Release of the cytoplasmic enzyme lactate dehydrogenase (LDH) was measured as a marker for necrotic epithelial damage 80 using a Cytotoxicity Detection Kit (Roche) according to the manufacturer's instructions. The background LDH value of uninfected epithelial cells (low control) was subtracted and the corrected LDH release was expressed as the percentage of high (full lysis) control (maximum LDH release induced by the addition of 0.25% Triton X-100 to uninfected epithelial cells for 5 min) unless otherwise stated. For the protection effect experiments, 0.1 ng ml −1 of IFN-β (InvivoGen) or neutralizing anti-human IFNAR2 antibody (4 ng ml −1 ; PBL InterferonSource) were added to epithelial cells 3 h before infection.
Transwell assay. Epithelial cells were seeded in 24-well plates (1 × 10 5 cells per well) in RPMI 1640 with FCS and incubated for 2 d at 37 °C and 5% CO 2 . After medium exchange with 750 µl of RPMI 1640 without FCS, transwell inserts (polycarbonate membrane inserts with 0.4 µm pore size; Corning), loaded with 250 µl of Candida suspension (4 × 10 6 yeast ml −1 ), were placed in the wells. After 3 h of incubation, the inserts were discarded and human RNA samples were collected and isolated as described above.
Measurement of epithelial cell mtDNA release. The release of mtDNA in response to infection was measured using the protocol of Bronner and O'Riordan 81 with some modifications. Briefly, epithelial cells were seeded in 6-well plates and infected as described above. After 6 h of infection, the medium and non-adherent Candida cells were removed. After the addition of 200 µl of the cell membrane detergent IGEPAL CA-630 (1%, NP-40; Sigma-Aldrich), cells were loosened by scraping. Lysates were incubated on ice for 15 min and centrifuged (12,000g, 15 min, 4 °C). The supernatant was used to isolate human mtDNA from the cytosolic fraction using the DNeasy Blood & Tissue Kit (QIAGEN), according to the manufacturer's instructions. Finally, cytosolic human mtDNA was measured by qPCR using 18S rRNA as a reference 81 . Results are presented relative to an uninfected control. Tunicamycin (10 µM; Sigma-Aldrich) was used as a positive control, as an endoplasmic reticulum stress inducer that leads to mitochondrial dysfunction 82 . The same procedure was carried out on yeast cells only to exclude that fungal cells would also release mtDNA following this protocol. The lysis step did not cause any lysis of yeast cells and no DNA was detected after the isolation procedure, confirming that the DNA obtained from the infected epithelial cells originated from epithelial cells only. RNA interference assay. The RNA interference assay was used to silence the expression of selected ISGs (IFI6, MX2 and CMPK2). Small interfering RNA, control siRNA (siRNA-A; catalogue no. sc-37007), siRNA transfection reagent and siRNA transfection medium were purchased from Santa Cruz Biotechnology. Epithelial cells were seeded in a 6-well plate and transfected with 1 µg of siRNA according to the manufacturer's instructions. After 48 h, cells were infected with C. albicans; 3 h post-infection, RNA was isolated and silencing of selected genes was confirmed using RT-qPCR. LDH release was measured 24 h post-infection. mtDNA depletion assay. Epithelial cells were seeded in 6-well plates and incubated for 2 d. Once confluent, 200 μM of 2′,3′-dideoxycytidine (Jena Bioscience) was added to the medium and epithelial cells were incubated for 6 d. mtDNA depletion was confirmed by quantifying mtDNA using qPCR, as described above.

Transfection of epithelial cells.
MtDNA was PCR-amplified from the entire human mitochondrial genome in 17 overlapping fragments as described previously 83 . Epithelial cells were transfected with amplified mtDNA fragments (Fig. 5c) or total DNA isolated from epithelial cells with and without their mtDNA depleted (Fig. 5d). Total DNA was isolated using the DNeasy Blood & Tissue Kit. A total of 2 μg ml −1 of DNA was used to transfect epithelial cells using the UltraCruz Transfection Reagent (Santa Cruz Biotechnology) according to the manufacturer's instructions. After 6 h, RNA samples were collected and the expression of ISGs was quantified by RT-qPCR.
Apoptosis induction. Epithelial cells were seeded in 96-well plates and infected with Candida as described above with the addition of staurosporine (1.2 µM) simultaneously with infection. After 6 h, RNA samples were collected and the expression of ISGs was quantified by RT-qPCR. Results were compared to infected cells incubated in the medium without staurosporine.

Collection of epithelial cell supernatants. Epithelial cells were infected with
Candida cells as described above and incubated for 24 h in the presence or absence of neutralizing anti-human IFNAR2 antibody. Supernatants were collected and stored at −80 °C until use (see Methods sections 'Cytokine release' and 'Neutrophil cytokine production'). The supernatants of only Candida cells grown in the absence of epithelial cells were included as control.
Cytokine release. Epithelial cells were infected with Candida cells as described above and incubated for 24 h. The release of IL-6, IL-8, IL-1α and IL-1β was measured with commercially available human enzyme-linked immunosorbent assay kits (IL-6, IL-8, IL-1β, Invitrogen; IL-1α, R&D Systems) according to the manufacturers' instructions.

Blood donors.
Human peripheral blood was collected from healthy volunteers (n = 3) with ethical approval and after obtaining written informed consent. This study was conducted according to the principles expressed in the Declaration of Helsinki version 2008. The blood donation protocol and use of blood for this study were approved by the institutional ethics committee of Jena University Hospital (permission no. 2207-01/08).
Neutrophil cytokine production. Primary human neutrophils were isolated from blood using a previously published protocol 84 and seeded in a 24-well plate (5 × 10 5 cells ml −1 ). Neutrophils were exposed for 24 h to the supernatants of epithelial cells that had been infected with each of the Candida species (24 h post-infection) to determine whether pro-inflammatory mediators released by epithelial cells played a role in neutrophil stimulation. The control supernatants of Candida cells alone were included to ensure that the neutrophils responded to the secretions of epithelial cells rather than the fungus (blue bars on Fig. 6e). After incubation, cytokine release was measured using an enzyme-linked immunosorbent assay as described above.
Fluorescence microscopy. Epithelial cells were seeded in µ-Slide 8 Well (IBIDI) and infected with Candida cells as described above. At 1 h post-infection, epithelial cells were stained with 100 nM of MitoTracker Deep Red FM for 15 min at 37 °C and washed and fixed with Roti Histofix 4%. Fluorescence imaging was done with the Cell Observer microscope (Carl Zeiss) with fluorescence settings at 644 and 665 nm. CCCP was used as a positive control (100 μM). Image acquisition was done in a fully blinded manner to avoid potential bias.
Transmission electron microscopy and imaging. Cells were fixed by adding glutaraldehyde (2.5% (v/v) final) to the growth medium. After 1 h, the cell layer was gently scraped off the surface, collected as a pellet by centrifuging at 600g and washed three times with PBS. After fixation in osmium tetroxide (1% (w/v) in distilled water) for 1 h, dehydration in ascending ethanol series with post-staining in uranyl acetate was performed. Afterwards, samples were embedded in epoxy resin (Araldite) and sectioned ultrathin (60 nm) using an ultramicrotome (Leica Ultracut E; Leica Biosystems). After mounting on filmed copper grids and post-staining with lead citrate, the sections were studied in a transmission electron microscope (EM 902 A; ZEISS) at 80 kV. Images were acquired with a 1k FastScan CCD camera (TVIPS). The mitochondrial aspect ratio (the ratio of length/width) was measured using ImageJ version 1.43h by analysing at least 80 mitochondria for each condition 85 in 1 biological replicate. Irregular structures were excluded from analysis. All transmission electron microscopy (TEM) analyses were conducted in a fully blinded manner to avoid potential bias in image acquisition and analyses.

Primary vaginal cells.
Primary human vaginal epithelial cells (catalogue no. PCS-480-010) were obtained from ATCC and cultured in vaginal epithelial cell basal medium (catalogue no. PCS-480-030; ATCC), supplemented with components from the Vaginal Epithelial Cell Growth Kit (catalogue no. PCS-480-040; ATCC). Cells were not authenticated but they were routinely checked for the absence of Mycoplasma contamination. Apoptosis/necrosis and mitochondrial membrane potential assays were performed as described for A-431 cells.
RNA-seq library preparation and sequencing. Library preparation for RNA-seq was performed with the TruSeq Stranded mRNA Sample Prep Kit v2 (catalogue no. RS-122-2101/2; Illumina) according to the manufacturer's instructions unless otherwise stated. One microgram of total RNA was used for poly(A)-mRNA selection using streptavidin-coated magnetic beads. Samples were then fragmented to approximately 300 base pairs (bp); subsequently, cDNA was synthesized using reverse transcriptase (SuperScript II; Invitrogen) and random primers. The second strand of the cDNA incorporated deoxyuridine triphosphate in place of deoxythymidine triphosphate. Double-stranded DNA was further used for library preparation. It was subjected to A-tailing and ligation of the barcoded TruSeq adaptors. All purification steps were done using AMPure XP beads (Agencourt). Library amplification was performed by PCR on size-selected fragments using the primer cocktail supplied in the kit. Final libraries were analysed using the Agilent DNA 1000 chip (Agilent) to estimate the quantity and check fragment size distribution; they were then quantified by qPCR using the KAPA Library Quantification Kit (Kapa Biosystems) before amplification with Illumina's cBot. To avoid potential batch effects, all samples were randomly distributed on the sequencing flow cells.
Libraries were sequenced with 2 × 50 (n = 21), 2 × 75 (n = 70) and 2 × 150 (n = 1) read lengths on the Illumina HiSeq 2500 system (2 × 50 bp) and HiSeq 3000 (the rest) at the Genomics Unit of the Centre for Genomic Regulation, Barcelona, Spain. Samples that contained mixed fungal and human RNA were sequenced for (on average) almost equal to 65 million reads (Supplementary Table  1 and Extended Data Fig. 1) to achieve sufficient sequencing depth for robust downstream analysis 78 .
For read mapping and quantification, we used the splice junction-sensitive read mapper STAR v.2.5.2b 89 using the basic two-pass mode and default parameters. For samples comprising either fungal or human RNA, reads were mapped to the corresponding reference genomes. In the case of pooled samples containing RNA from both host and pathogen, data were mapped to concatenated human and corresponding yeast reference genomes. For human data, we used the primary genome assembly GRCh38 and genome annotations from the Ensembl database release 89 (last accessed on 8 August 2017 (ref. 90 )). Reference genomes and genome annotations for C. albicans SC5314 (assembly 22), C. glabrata CBS138 and C. parapsilosis CDC317 were obtained from the Candida Genome Database (CGD, last accessed on 17 August 2017 (ref. 91 )). From the phased reference genome and annotations of C. albicans, we selected haplotype A to perform further analysis to avoid substantial rates of ambiguously mapped reads. The reference sequence and annotations for C. tropicalis were obtained from the RefSeq database (last accessed on 9 August 2017 (ref. 92 )). The genes missing from the RefSeq genome annotations were manually added from the Candida Gene Order Browser 93 . GFF genome annotation files were converted to GTF format using the gffread utility v.0.9.8 (ref. 94 ). We used Centrifuge v.1.0.4 (ref. 95 ) to test the presence of viral contamination in our dataset, by remapping the reads that mapped either to the human or fungal reference genomes to the whole National Center for Biotechnology Information nucleotide database (downloaded on 23 of March 2018). No traces of contamination were observed.
Differential gene expression analysis was performed using the Bioconductor package DESeq2 v.1.26.0 (ref. 96 ) using the read counts obtained with STAR mapping. For human samples and each fungal species, we compared time point 0 with the other time points throughout the course of infection by Wald test using the contrast option of DESeq2. To detect any statistically significant changes of expression throughout the course of infection, we also used a likelihood ratio test of DESeq2, by dropping the 'time' component of the formula design. Genes with a log 2 (fold change) > 1.5 and adjusted P (P adj ) < 0.01 were considered differentially expressed unless otherwise stated. To account for possible batch effects in the experiments involving C. albicans ece1Δ/Δ and non-viable fungal cells, we applied the RUVg function of the Bioconductor package RUVseq v.1.20 (ref. 97 ) using non-differentially expressed genes (base mean > 10 and P adj > 0.05 obtained by likelihood ratio test in DESeq2) across all samples and time points as negative controls. Since the optimal parameters for the batch effect removal algorithm are not defined a priori, we employed a strategy of incremental increase of k values (k = 1,2,…,n), until we observed disruption of the PCA clustering of original data from the first batch of sequencing. To perform differential expression analysis, the obtained matrix of batch effect coefficients was further supplied to the design formula of the DESeq2 object, which was subsequently run using original count data. To plot 'batch-free' PCA plots, we used batch-corrected counts retrieved from the RUV package.
The list of 1-to-1 orthologues between the four fungal species was obtained from the CGD. For interspecies gene expression comparisons, the raw read counts for each fungal species were normalized by gene length and library size. GO term enrichment analysis was performed using clusterProfiler v.3.14.3 (ref. 98 ). GO enrichment plots were produced with the dotplot function using showCategory set to 10 (for human data) and 8 (for fungal data) for better plot readability (the full list of GO enrichments is available in Supplementary File 9). Adjustment of P values was done using the Benjamini-Hochberg procedure. GO information for fungal species was obtained from the CGD, while for human data we used the Bioconductor package Genome wide annotation for Human (org.Hs.eg.db) v.3.10.0 in R 99 .
We assessed the patterns of host-pathogen gene coexpression across the infections using the weighted correlation network analysis approach implemented in WGCNA v.1.69 (ref. 100 ). For each infection, we combined fungal and corresponding human data at all available time points of infection, excluding the data for the C. albicans ece1Δ/Δ mutant. As recommended by the package developers, we selected genes that had 10 or more counts in more than 90% of samples for downstream analysis. As expression levels, we used variance-stabilized read count data obtained using the vst function of DESeq2. Before the actual network construction, we first selected the β power values using the pickSoftThreshold function implying an unsigned network. The minimum β values reaching 80% of scale-free network topology, namely 12, 20, 7 and 22 for infections with C. albicans, C. glabrata, C. parapsilosis and C. tropicalis, respectively, were used for downstream analysis. After network construction, we inferred modules (that is, highly interconnected clusters of genes) in the WGCNA networks using 1-topology overlap matrix values at the 0.25 hclust tree cut-off and identified eigengenes (that is, the first principal component of each module). For each identified module, we performed GO term enrichment analysis of fungal and host genes using clusterProfiler (Supplementary File 10), selecting the top 3 enrichments with the lowest P adj values. Then, the fungal gene content of each module of a given fungal species was compared against those of all modules of the other three species taking into account 1-to-1 orthology information (Extended Data Fig. 3). This analysis was done for modules that contained at least one fungal gene. Similarity between fungal gene contents of two given modules was defined as the intersection of the fungal gene lists of these modules divided by the union of these gene lists.
All custom calculations and visualizations were performed in R v.3.6.1 using various packages (all packages and their versions are available at our GitHub page https://github.com/Gabaldonlab/Host-pathogen_interactions).

Statistics and reproducibility.
Experiments were performed in biological triplicates (n = 3) with 3 different donors (neutrophil cytokine release) or 3 independent experiments. The growth curve experiments (Fig. 1b) were performed five times to ensure reproducibility. Only the mitochondrial aspect ratio was calculated based on 1 biological replicate ( Fig. 4d) but multiple mitochondria were measured for each condition (n > 80). All microscopy findings were reliably reproduced. Data were analysed using Prism 8 (GraphPad Software). Values are presented as the mean ± s.d. All the ratio data were log-transformed as indicated before statistical analysis in Prism and compared to 0 (uninfected/non-transfected/ non-treated control) using a one-way analysis of variance (ANOVA) with Tukey's multiple comparisons test (Figs. 1b and 6b,d) or Dunnett's multiple comparisons test (Figs. 4e-g, 5c and 6a and Extended Data Fig. 5c,d), Kruskall-Wallis test with two-sided Dunn's multiple comparison test (Fig. 4d), two-tailed one-sample t-test (Fig. 5a) or two-way ANOVA and Sidak's multiple comparisons test (Fig.  5b,d). Statistical significance is indicated in the figures as *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001 and ****P ≤ 0.0001. The exact P values are provided in the Source data.
Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The data supporting the findings of this study are available within the paper and its Supplementary Information. All relevant data, including further image and processed data are available by request from the corresponding authors, with the restriction of data that would compromise the confidentiality of blood donors. Raw sequencing data have been deposited in the Sequence Read Archive under accession nos. SRR10279972-SRR10280067. Mapped data from the four Candida species can be mined and browsed at Candidamine (http://candidamine.org/ candidamine/begin.do); the gene read counts from all samples can be found in our GitHub page https://github.com/Gabaldonlab/Host-pathogen_interactions along with the data analysis scripts for results reproducibility. Publicly available datasets/ databases used in the study can be accessed at: Ensembl (https://www.ensembl.org/ index.html); RefSeq (https://www.ncbi.nlm.nih.gov/refseq/); CGOB (http://cgob.ucd.ie/); NCBI FTP site (https://www.ncbi.nlm.nih.gov/home/ download/); CGD (http://www.candidagenome.org/); and Genome wide annotation for Human (https://bioconductor.org/packages/release/data/ annotation/html/org.Hs.eg.db.html). Source data are provided with this paper.

code availability
All transcriptome data analysis results, including figures, extended data and supplementary materials are fully reproducible using the scripts provided at our GitHub page https://github.com/Gabaldonlab/Host-pathogen_interactions.

Data exclusions
If no data were excluded from the analyses, state so OR if data were excluded, provide the exact number of exclusions and the rationale behind them, indicating whether exclusion criteria were pre-established.

Non-participation
State how many participants dropped out/declined participation and the reason(s) given OR provide response rate OR state that no participants dropped out/declined participation.

Randomization
If participants were not allocated into experimental groups, state so OR describe how participants were allocated to groups, and if allocation was not random, describe how covariates were controlled.

Ecological, evolutionary & environmental sciences study design
All studies must disclose on these points even when the disclosure is negative.

Study description
Briefly describe the study. For

Sampling strategy
Note the sampling procedure. Describe the statistical methods that were used to predetermine sample size OR if no sample-size calculation was performed, describe how sample sizes were chosen and provide a rationale for why these sample sizes are sufficient.

Data collection
Describe the data collection procedure, including who recorded the data and how.
Timing and spatial scale Indicate the start and stop dates of data collection, noting the frequency and periodicity of sampling and providing a rationale for these choices. If there is a gap between collection periods, state the dates for each sample cohort. Specify the spatial scale from which the data are taken

Data exclusions
If no data were excluded from the analyses, state so OR if data were excluded, describe the exclusions and the rationale behind them, indicating whether exclusion criteria were pre-established.

Reproducibility
Describe the measures taken to verify the reproducibility of experimental findings. For each experiment, note whether any attempts to repeat the experiment failed OR state that all attempts to repeat the experiment were successful.

Randomization
Describe how samples/organisms/participants were allocated into groups. If allocation was not random, describe how covariates were controlled. If this is not relevant to your study, explain why.

Blinding
Describe the extent of blinding used during data acquisition and analysis. If blinding was not possible, describe why OR explain why blinding was not relevant to your study.
Did the study involve field work?

Yes No
Field work, collection and transport

Field conditions
Describe the study conditions for field work, providing relevant parameters (e.g. temperature, rainfall).

Location
State the location of the sampling or experiment, providing relevant parameters (e.g. latitude and longitude, elevation, water depth).

Access and import/export
Describe the efforts you have made to access habitats and to collect and import/export your samples in a responsible manner and in compliance with local, national and international laws, noting any permits that were obtained (give the name of the issuing authority, the date of issue, and any identifying information).

Disturbance
Describe any disturbance caused by the study and how it was minimized.
Reporting for specific materials, systems and methods

Specimen deposition
Indicate where the specimens have been deposited to permit free access by other researchers.

Dating methods
If new dates are provided, describe how they were obtained (e.g. collection, storage, sample pretreatment and measurement), where they were obtained (i.e. lab name), the calibration program and the protocol for quality assurance OR state that no new dates are provided.
Tick this box to confirm that the raw and calibrated dates are available in the paper or in Supplementary Information.

Animals and other organisms
Policy information about studies involving animals; ARRIVE guidelines recommended for reporting animal research

Laboratory animals
For laboratory animals, report species, strain, sex and age OR state that the study did not involve laboratory animals.

Wild animals
Provide details on animals observed in or captured in the field; report species, sex and age where possible. Describe how animals were caught and transported and what happened to captive animals after the study (if killed, explain why and describe method; if released, say where and when) OR state that the study did not involve wild animals.

Field-collected samples
For laboratory work with field-collected samples, describe all relevant parameters such as housing, maintenance, temperature, photoperiod and end-of-experiment protocol OR state that the study did not involve samples collected from the field.