Vincent Estrade
email: vincent.estrade@gmail.com

Baudouin Denis De Senneville

Paul Meria
email: paul.meria@aphp.fr

Christophe Almeras
email: c_almeras@yahoo.fr

Franck Bladou
email: franck.bladou@chu-bordeaux.fr

Jean-Christophe Bernhard
email: jean-christophe.bernhard@chu-bordeaux.fr

Gregoire Robert
email: gregoire.robert@chu-bordeaux.fr

Olivier Traxer
email: olivier.traxer@aphp.fr

Michel Daudon
email: michel.daudon@aphp.fr

Toward improved endoscopic examination of urinary stones: a concordance study between endoscopic digital pictures vs. microscopy

Keywords: Morpho-constitutional analysis of urinary stones, flexible ureteroscopy training, endoscopic diagnosis, aetiological lithiasis, LASER fragmentation of stones 1. ESR is Endoscopic Stone Recognition, 2. LASER

is Light Amplification by Stimulated Emission of Radiation, 3. FTIR is Fourier Transform InfraRed spectroscopy. 4. CCD sensor is Charge-Coupled Device sensor ABSTRACT Objective: To improve endoscopic recognition of the most frequently encountered kidney stone morphologies for a better etiological approach in lithiasis by urologists.

Materials and methods: An expert urologist intra-operatively and prospectively (between June 2015 and June 2018) examined the surface, the section and the nucleus of all encountered kidney stones. Fragmented stones were subsequently analysed by a biologist based on both microscopic morphological (i.e. binocular magnifying glass) and infrared (i.e. FTIR) examinations (microscopists were blinded to the endoscopic data). Morphological criteria were collected and classified for the endoscopic and microscopic studies. The Wilcoxon-Mann-Whitney test was carried out to detect differences between the endoscopic and microscopic diagnoses. A diagnosis for a given urinary stone was considered "confirmed" for a non-statistically significant difference.

Results: A total of 399 urinary stones were included in this study: 51.4% of the stones exhibited only one morphological type while 48.6% were mixed stones (41% had at least two morphologies and 7.6% had three morphologies). The overall matching rate was 81.6%.

Diagnostics were confirmed for the following morphologies: whewellite (Ia or Ib), weddellite (IIa or IIb), uric acid (IIIa or IIIb), carbapatite-struvite association (IVb), brushite (IVd).

Conclusions: Our preliminary study demonstrates the feasibility of using endoscopic morphology for the most frequently encountered urinary stones and didactic boards of confirmed endoscopic images are provided. The current study constitutes the first step toward endoscopic stone recognition, which is essential in lithiasis. We provide didactic boards of confirmed endoscopic images which paves the way for automatic computer-aided in-situ recognition.

INTRODUCTION

A morpho-constitutional examination of urinary stones plays an essential role in the aetiological diagnosis [START_REF] Daudon | Urinary Calculi: Review of classification methods and correlations with etiology[END_REF][START_REF] Daudon | Recurrence rates of urinary calculi according to stone composition and morphology[END_REF][START_REF] Cloutier | Kidney stone analysis: "Give me your stone, I will tell you who you are![END_REF][START_REF] Estrade | Why should urologist recognize urinary stone and how? The basis of endoscopic recognition[END_REF][START_REF] Daudon | Comprehensive morpho-constitutional analysis of urinary stone improves etiological diagnosis and therapeutic strategy of nephrolithiasis[END_REF]. The international morpho-constitutional classification of urinary stones includes seven groups denoted by roman numerals (i.e. "I" to "VI") (Table 1). Each group is associated with a specific crystalline type: I = whewellite, II = weddellite, III = uric acid and urates, IV = calcium and non-calcium phosphates, and V = cystine (group VI other stones).

Each group is then divided into several subgroups to differentiate morphologies and aetiologies for a given crystalline type. Five morphological subtypes are encountered in whewellite group I (differentiated by subscripts in the Latin alphabet: Ia to Ie), each one being associated with a specific aetiology. Practically: Ia = excessive concentration of oxalate in the urine induced by diet or default of diuresis, Ic = inherited primary hyperoxaluria, Ie = enteric hyperoxaluria induced by inflammatory bowel disease or by-pass. Interestingly, the most recent lithogenic events (in chronological order) are located on the surface of the stone, whereas less recent events are observable on a section of the stone. The nucleus of the stone, which is the oldest part, corresponds to the initial lithogenic context. Urinary stones are mixed (i.e. include at least two morphologies) in almost half of cases [START_REF] Daudon | Comprehensive morpho-constitutional analysis of urinary stone improves etiological diagnosis and therapeutic strategy of nephrolithiasis[END_REF].

Daudon and Cloutier emphasised the importance of studying urinary stone morphologies [START_REF] Daudon | Urinary Calculi: Review of classification methods and correlations with etiology[END_REF][START_REF] Cloutier | Kidney stone analysis: "Give me your stone, I will tell you who you are![END_REF][START_REF] Daudon | Comprehensive morpho-constitutional analysis of urinary stone improves etiological diagnosis and therapeutic strategy of nephrolithiasis[END_REF] for an aetiological diagnosis of stone disease. A complete examination of the entire stone includes a visual morphological examination of the stone surface, the stone section and the nucleus, as well as a spectrophotometric infrared recognition (FTIR) analysis of crushed stone fragments. However, it is now well established that modern endoscopic treatment of urinary stones relies on LASER fragmentation [START_REF] Emiliani | Optimal Settings for the Noncontact Holmium:YAG Stone Fragmentation Popcorn Technique[END_REF][START_REF] Doizi | Dusting technique for lithotripsy: what does it mean?[END_REF][START_REF] Keller | Thulium fiber laser: ready to dust all urinary stone composition types?[END_REF][START_REF] Doizi | Evolution of the stone treatment over 30 years in a French academic institution[END_REF][START_REF] Carpentier | Management of Adult's Renal and Ureteral Stones. Update of the Lithiasis Committee of the French Association of Urology (CLAFU). General Considerations[END_REF][START_REF] Chabannes | Management of adult's renal and ureteral stones. Update of the Lithiasis Committee of the French Association of Urology (CLAFU). General considerations[END_REF][START_REF] Estrade | Place of the flexible ureterorenoscopy first choice for the treatment of kidney stones. Survey results practice committee of the AFU lithiasis completed in 2011[END_REF][START_REF] Raynal | In-hospital stays for urolithiasis: Analysis of French national data[END_REF]. Fragmentation, whether achieved with "popcorn" [START_REF] Emiliani | Optimal Settings for the Noncontact Holmium:YAG Stone Fragmentation Popcorn Technique[END_REF] or "dusting" modes [START_REF] Doizi | Dusting technique for lithotripsy: what does it mean?[END_REF], destroys the morphology of the targeted stone [START_REF] Keller | Thulium fiber laser: ready to dust all urinary stone composition types?[END_REF]. The morphological examination, which is the first essential diagnostic step, is impossible to redo once the stone is destroyed. Moreover, Keller et al. recently showed the impact of LASERbased dusting on changes in stone composition with significant changes in the infrared spectra (particularly for weddellite, carbapatite, struvite and brushite) [START_REF] Keller | Thulium fiber laser: ready to dust all urinary stone composition types?[END_REF]. Consequently, a FTIR examination of the stone powder itself does not provide sufficient information to determine correctly the lithogenic stage [START_REF] Daudon | Urinary Calculi: Review of classification methods and correlations with etiology[END_REF][START_REF] Daudon | Recurrence rates of urinary calculi according to stone composition and morphology[END_REF][START_REF] Cloutier | Kidney stone analysis: "Give me your stone, I will tell you who you are![END_REF][START_REF] Estrade | Why should urologist recognize urinary stone and how? The basis of endoscopic recognition[END_REF][START_REF] Daudon | Comprehensive morpho-constitutional analysis of urinary stone improves etiological diagnosis and therapeutic strategy of nephrolithiasis[END_REF]. This finding reinforces the need to observe the morphology of the stone before LASER-induced destruction to preserve an aetiological approach.

The current study constitutes the first step toward endoscopic stone recognition (ESR). We aimed to improve the recognition of the most frequently encountered kidney stone morphologies. Thus, an endoscopic-based examination of unfragmented stones was compared with microscopic observations of laser-fragmented stones, and a concordance study was conducted.

MATERIAL & METHOD

The study adhered to all local regulations and data protection agency recommendations (the National Commission on Informatics and Liberty (CNIL) dictates). Patients have been informed

for the use of their data anonymously.

Endoscopic study

An expert urologist (VE, 20 years of experience) intra-operatively and prospectively examined stones between June 2015 and June 2018 using a flexible digital ureterorenoscope (Olympus URF-V CCD sensor). We recall that visual cues that the surgeon has to consider for the estimation of a stone type are summarized in Table 1. The examination included a visual observation of the stone surface first, before LASER fragmentation, then visual observation of the section and the nucleus and after LASER stone section. LASER-based stone section in two parts was performed using the following LASER parameters: frequency = 5 Hz, energy = 1.2-1.4 J, power = 6-7 W, pulse length = short, fibre diameter = 230 or 270 µm [START_REF] Traxer | Flexible ureteroscopy with Holmium laser: technical aspects[END_REF][START_REF] Traxer | Flexible ureteroscopy with Holmium laser: the tools[END_REF]. That way, it was possible to fragment all types of pure and mixed stones. Once fragmented, the fragments were removed with a Nitinol basket (Zerotip 1.9) through the ureteral access sheat Coloplast Retrace 12/14 or 10/12. An additional fragmentation session was carried out when needed. Macro-fragments including at least one representative copy of the stone surface, the stone section and the nucleus were removed using a nitinol basket for subsequent examination in a dedicated laboratory.

Microscopic study

The fragmented stones were subsequently analysed by a biologist (MD, 40 years of experience) based on both morphological (i.e. binocular magnifying glass) and infrared (i.e. SPIR) analyses. Similar to the above-mentioned endoscopic analysis, the examination included the surface, the section and the nucleus of each stone.

Observations

The urinary stones were classified according to microscopic morphological instructions given by Daudon et al. in [START_REF] Daudon | Urinary Calculi: Review of classification methods and correlations with etiology[END_REF]. We collected the following eight observations for all urinary stones in the endoscopic and microscopic studies.

 Three observations about surface morphology referred to as "majority" (for the most visible morphology on the surface of the stone), "secondary" (for potentially secondary visible morphology) and "minority" (for another potentially visible morphology).

 Three observations about the inner structure: "central" (for the morphology visible in the centre of the stone), "peripheral" (for the morphology visible on the periphery of the stone) and "intermediate" (for the morphology in between).

 Two observations about the nucleus: "majority" (for the most visible morphology in the nucleus of the stone) and "minority" (for another potentially visible morphology).

Concordance between the endoscopic-and microscopic-based observations

All data were collected into a single Excel spreadsheet that was retrospectively analysed using Matlab software (MathWorks, Inc., Natick, MA, USA).

We selected a subset of urinary stones for which "Ia" or "Ib" morphologies were encountered in any microscopic-based observation. Among these, urinary stones for which the "Ia" or "Ib" morphologies were encountered in any endoscopic-based observation were listed and referred to as a "good match". The total number of urinary stones with a "good match" was calculated and recorded. The Wilcoxon-Mann-Whitney test was carried out to determine whether the differences between endoscopic and microscopic diagnoses were significant. A significance threshold of p = 0.05 was used. A diagnosis for a given urinary stone was considered "confirmed" for a non-significant difference between the endoscopic and microscopic examinations (p > 0.05). In such a case, the area under the ROC curve (AUROC), sensitivity, specificity, positive predictive value (PPV) and negative predictive value (NPV) were calculated (considering the presence of a given morphology in the microscopic study as a "positive case").

RESULTS

Demographic characteristics

In total, 399 urinary stones were included in this retrospective study. Among them, 48.6% were mixed stones (41% had at least two morphologies and 7.6% had three morphologies).

Our cohort was composed of 63% calcium oxalate stones (group I whewellite + group II weddellite), 20% phosphate stones (group IV: carbapatite, brushite and struvite), 15% uric acid and urate stones (group III), 0.2% cystine (group V) and 1.8% of other stones (group VI).

Typical images obtained for each of the six urinary stone groups are shown in Figures 1 to 9 respectively. Endoscopic-and microscopic-based observations of the stone surface and sections are reported, together with corresponding aetiological conditions.

Stone morphologies analysed

We identified 16 morphologies distributed over the following groups:

-For whewellite group I: Ia, Ib, Id and Ie -For weddellite group II: IIa and IIb -For the uric acid and urate group III: IIIa and IIIb -For the calcium and non-calcium phosphate group IV: IVa1 and IVa2 (carbapatite), IVb (carbapatite and struvite), IVc (struvite), IVd (brushite).

-For the cystine group V: Va -For the group VI aggregating other stones (protein matrices +/-whewellite): VIa and VIb For this study, we grouped the following morphologies that had similar aetiologies: Ia and Ib, IIa and IIb, IIIa and IIIb, Va and Vb and VIa and VIb.

Concordance between the endoscopic-and microscopic-based observations

Table 2 summarises the concordance between our endoscopic-and microscopic-based observations.

Morphologies validated by the microscopic and endoscopic examinations

Diagnoses were confirmed for the following morphologies: whewellite (Ia or Ib: concordance = 85%, n = 205 ; Id: concordance = 92%, n = 12 ; Ie: concordance = 80%, n = 5), weddellite (IIa or IIb: concordance = 85%, n = 178), uric acid (IIIa or IIIb: concordance = 91% n = 64), carbapatite-struvite association (IVb: concordance = 50%, n = 10), brushite (IVd: concordance = 65%, n = 23). Several other pure stone morphologies depicted good matching, although just below the significance threshold of 0.05, such as carbapatite (IVa1, concordance = 81%, n = 176). Other pure stone morphologies depicted excellent matching but insufficient number of cases, such as cystine (Va: concordance = 100%, n = 1).

DISCUSSION

Our study consisted of endoscopic recognition of the morphological elements constituting urinary stones before LASER stone destruction. ESR allows the morphological identification of an entire stone and is thus essential in lithiasis. While the main objective of ESR is to identify a stone type, it also provides examinations of the anatomy of the excretory pathway as well as the renal papillae. Flexible ureteroscopy is thus a great diagnostic and therapeutic candidate in lithiasis. We offer didactic boards to help with endoscopic recognition comprised of reference images and descriptions (surface, section and nucleus), and the lithogenic mechanisms and aetiologies associated with each morphology (Figures 123456).

It must be reported that a learning curve is needed to acquire the ESR skill which may limit its translation to practical use. In the current study, the urologist (VE) had, first of all, to learn the classification of the different types of stone surface, section and nuclei using the microscopic images provided by Daudon and al. [START_REF] Daudon | Urinary Calculi: Review of classification methods and correlations with etiology[END_REF][START_REF] Cloutier | Kidney stone analysis: "Give me your stone, I will tell you who you are![END_REF][START_REF] Daudon | Comprehensive morpho-constitutional analysis of urinary stone improves etiological diagnosis and therapeutic strategy of nephrolithiasis[END_REF]. Subsequently, he had to acquire the ability to recognize such types of stone surface, section and nuclei on endoscopic images, based on the learned microscopic images and associated descriptions. Along this line, one goal of our concordance study is to reduce possible subjectivity in ESR and urologist bias by history. To make easier the ESR for urologist, we propose didactic boards of confirmed endoscopic images for the most frequently encountered urinary stones on a daily practice. This follows the work of Bergot et al. which showed that ESR teaching of junior urologists allows them to acquire the skill to recognize the most frequently encountered stones quickly [START_REF] Bergot | The basis of endoscopic stones recognition, a prospective monocentric study[END_REF]. On the other hand, a computer-assisted approach delivers reproducible results and minimises operator dependency, as visual interpretation of stone images lacks a learning curve when the process is automated. It is therefore promising to train artificial intelligence algorithms with our confirmed endoscopic images. Black et al. have recently shown that artificial intelligence (deep learning) applied to in vitro surface and section images of stones represents a great asset for the automatic recognition of whewellite, weddellite, uric acid, brushite and struvite stones [START_REF] Black | Deep learning computer vision algorithm for detecting kidney stone composition[END_REF]. Such algorithms may be able to do tests with various predefined score and error levels. Combined with our confirmed endoscopic images, deep neural networks, such as convolutional neural networks (CNN), which are able to process efficiently a high number of specific images and videos, are good candidates for automatic ESR.

In our study, the SPIR examination detected 12 struvite cases that were not described using morphological analyses. The combination of morphological and SPIR examinations improved the diagnostic concordance of urolithiasis [START_REF] Daudon | Urinary Calculi: Review of classification methods and correlations with etiology[END_REF][START_REF] Daudon | Recurrence rates of urinary calculi according to stone composition and morphology[END_REF][START_REF] Cloutier | Kidney stone analysis: "Give me your stone, I will tell you who you are![END_REF][START_REF] Estrade | Why should urologist recognize urinary stone and how? The basis of endoscopic recognition[END_REF]. Our findings of the epidemiological distribution of the different crystalline types as well as the rate of pure/mixed stones are in line with already published results [START_REF] Daudon | Epidemiology of urolithiasis[END_REF]. Therefore, our cohort is representative of the stone distribution that any urologist might encounter. Our study confirmed the following six morphological types: Ia or Ib, IIa or IIb, IIIa or IIIb, IVa, IVb and IVd. These six morphological types cover 95% of the most common pure stones that urologists encounter in its practice daily [START_REF] Daudon | Comprehensive morpho-constitutional analysis of urinary stone improves etiological diagnosis and therapeutic strategy of nephrolithiasis[END_REF].

Our study showed that recognising the main morphological surface criteria is easier than the other stages of the analysis (i.e. section and nucleus). This result is crucial because the majority of the surface area of a urinary stone represents recent lithogenesis. Its destruction by a LASER de facto leads to a loss of aetiological information.

The enthusiasm of urologists for the LASER spraying of calculations in the "dusting" mode and the upcoming arrival of a new generation of very high frequency "super Thulium fibered" LASER [START_REF] Andreeva | Preclinical comparison of superpulse thulium fiber laser and a holmium:YAG laser for lithotripsy[END_REF][START_REF] Traxer | Thulium fiber laser: the new player for kidney stone treatment? A comparison with Holmium:YAG laser[END_REF] will positively affect the success of interventions to minimise the rate of residual fragments and improve results without fragments.

However, the urologist should not forget that the morpho-constitutional examination of a urinary stone, which is currently the only solution for observing the entire stone, is as important as the pathological examination of a surgical specimen for the therapeutic strategy of onco-urology cases.

Despite the variety of compositions and morphologies observed in urinary calculi, about 90% are composed of a limited number of crystalline species and morphological characteristics that are easily recognised through an endoscopic examination. The morpho-constitutional classification of urinary stones previously published is particularly suitable for this purpose.

Endoscopic examination of 399 stones revealed good concordance between endoscopic and microscopic typing of the stones. For example, concordance of the results was observed in 86.1% of type I stones, 85% of type II stones, 91% of type III stones and 79% of type IV stones made of calcium phosphate. Within stone types, it was possible to identify with good agreement specific subtypes related to a more accurate aetiological diagnosis, such as subtype Ia or Ib, which is mainly related to dietary hyperoxaluria due to oxalate-rich food intake or low diuresis, while subtype Ie is related to enteric hyperoxaluria. Subtype Id stones suggest stasis and anatomical confinement related to urological anomalies. The presence of a thin greyish layer on the surface of a type Ia stone is a marker for dietary hyperoxaluria in most cases. ESR may be more critical when mixtures of phosphates are present, as suggested by differences in concordance among the various subtypes. As shown in Table 2, concordance was high (about 80%) for subtype IVa1 and IVb stones. In contrast, concordance was lower (65%) for subtype IVd, which is a marker for brushite-containing stones. Concordance was significantly better when brushite was the main component of the stone, reaching 77.8% of cases. However, the concordance between the endoscopic and microscopic examinations was only 21% for IVa2 stones, which are mainly composed of carbapatite. Two reasons could explain this low concordance. First, subtype IVa2 is uncommon. Second, some morphological characteristics of these stones, such as the presence of tiny cracks within the structure, are reinforced by drying the stone before the microscopic examination, and thus could be less visible during the endoscopic analysis. Training based on an examination of a larger number of samples would improve the concordance for the stone types providing poor agreement between the endoscopic and microscopic examinations. Notably, the global concordance level was 81.6% in this first study. As Almeras et al described a higher incidence of IVa2 stones when intra papillary cristallization occured and others descriptions correlated to different lithogenesis mechanisms, the description of papillary abnormalities during flexible ureteroscopies would certainly add a diagnostic value to this endoscopic stone recognition [START_REF] Almeras | Classification of the renal papillary abnormalities by flexible ureteroscopy: evaluation of the 2016 version and update[END_REF].

It should be emphasised that it is essential to record the traceability of the morphology of the surface, the section and the nucleus of a urinary stone. The urologist must archive all images to allow retrospective expertise, if necessary, by a biologist, an expert urologist or machine learning.

CONCLUSION

This study provided a 81.6% concordance between endoscopic and microscopic urine stone characterisations. We created didactic boards of confirmed endoscopic images for the most frequently encountered urinary pure stones, including whewellite (Ia or Ib, Id, Ie), weddellite (IIa or IIb), uric acid (IIIa or IIIb), carbapatite (IVa or IVb), brushite (IVd) and cystine (Va or Vb).

These boards can be used by urologists to learn how to recognise stones in-situ using an endoscopic examination before they are destroyed. Thereby, urologists must be more involved in the aetiological diagnosis, as well as in the therapeutic nephrolithiasis strategy, to maintain an essential role in the decision tree of urinary stone management. This first series of endoscopic images will be supplemented by more stone images in future studies to increase the number of examples of rarer urinary stones. This task is already in progress at our institution and will rely on the epidemiological distribution of the occurrence of urinary stones to obtain a sufficiently large population for an opposable statistical approach.

Accurate recognition of the most frequently encountered kidney stone morphologies allows for the development of an endoscopic stone image data base designed for automatic computer-aided in-situ recognition. Unorganized porous section of same color as surface IIId Heterogeneous, embossed, rough and porous surface with a greyish to darkbrown color.

Alternative thick, brown and thin, beige layers with small porous zones.

IVa1

Whitish, rough homogenous surface.

Poorly organized structure, with loose concentric layers of same color as the surface.

IVa2 Very peculiar morphology characterized by a yellow-brown, smooth surface with a glazed appearance and small cracks Section showing irregularly arranged thin whitish and thick yellow-brown concentric layers.

IVb

Heterogeneous surface, both embossed and rough. Heterogeneous color from whitish to dark brown.

Alternate thick whitish and thin brown-yellow layers.

IVc

Aggregates of large crystals with blunt angles and edges. Color : Whitish color.

Diffuse, loose radial crystallization. Whitish color.

IVd

Large rod-shaped crystals thereafter evolving toward slightly rough surface or resembling cabbage. Color : whitish or beige.

Radial crystallization with locally concentric layers. Color : whitish to beige.

Va

Homogeneous, bumpy or rough, surface. Waxy aspect. Color: light brown yellow.

Diffuse radial organization or unorganized section. Color: light brown-yellow.

Vb

Homogeneous smooth or finely rough surface, Color: whitish to pale beige.

Heterogenous structure made of finely concentric microcrystalline beige organization in periphery around a compact, crystalline unorganized light brown-yellow core. VIa Soft matrices, light-brown in color in contrast with other types of protein-rich calculi.

Unorganized structure.

Color: light brown.

TABLE LEGENDS Table 1 :

 LEGENDS1 Morpho-constitutional classification of urinary stones.

	Type	Surface	Section
	Ia	Smooth or mammillary, dark-brown	Section showing compact concentric layers
		surface. Frequent papillary umbilication	with a radiating organization starting from a
		with a piece of Randall's plaque	nucleus (often a Randall's plaque)
	Ib	Mammillary, rough surface. Color :	Unorganized section. Absence of an
		brown to dark-brown,	umbilication
	Id	Very smooth surface, beige or pale	Compact section showing thin concentric lay-
		brown in color	ers without radiations
	Ic	Pale brown-yellowish budding surface	Loose, unorganized section. Light color
		(without any umbilication)
	Ie	Locally budding, mamillary or rough	Heterogeneous structure, with a mixture of
		surface. Heterogeneous color from	poorly organized brown-yellow pale areas and
		brown-yellow pale to dark brown	of locally concentric dark brown layers with
			radiating organization.
	IIa	Yellow or light-brown prickly,	Loose radial crystallization. Color : pale
		spiculated surface due to presence of	brown-yellow
		aggregated bipyramidal crystals with
		sharp angles and edges
	IIb	Yellowish or light-brown surface with	Compact poorly organized crystalline section.
		smooth, long bipyramidal crystals, thus	Color : pale brown-yellow
		resembling small desert roses
	IIc	Rough grey-beige	Diffuse concentric structure at the periphery
		to pale brown surface,
	IIIa	Homogeneous smooth surface. Color :	Concentric layers with a radiating organization
		typically orange	around a well-defined nucleus. Color : ochre to
			orange
	IIIb	Rough and porous surface with a heter-	Poorly organized, porous structure. Color :
		ogeneous, beige to orange-red color	Ochre to orange
	IIIc	Rough and locally porous beige or grey-
		ish surface.

Table 2. Match results between endoscopic and microscopic studies. p-values denoted by "/" corresponds to data that are insufficiently populated for statistical comparison. Nonsignificant differences between endoscopy and microscopy analysis are denoted by "NS". AUC: area under the ROC curve; PPV: positive predictive value; NPV: negative predictive value; FPR: false predictive rate; FNR: false negative rate. AUROC, Sensitivity, Specificity, PPV, NPV, FPR and FNR are reported only for non-significant differences between endoscopy and microscopy findings.

Section :

IIIa : Compact concentric layers with a radiating organization around a well-defined nucleus.

Color: ocher to orange IIIb : with a poorly organized, porous structure.

IIIa IIIa IIIb

IIIb IIIa IIIb IIIa IIIb

Common Etiology :

IIIa :Low urine pH, stasis prostate hypertrophy, metabolic syndrome, ammoniagenesis defect

IIIb :Insulin resistance, metabolic syndrome, type 2 diabetes, ammoniagenesis defect, low urine pH

Section :

Diffuse or unorganized radial, light brown-yellow

Common

Etiology : Cystinuria


Johannes R Björk
email: rbjork@nd.edu

Robert B O'hara

Marta Ribes
email: mribes@icm.csic.es

Rafel Coma
email: coma@ceab.csic.es

José M Montoya

The dynamic core microbiome: Structure, dynamics and stability

The long-term stability of microbiomes is crucial as the persistent occurrence of beneficial microbes and their associated functions ensure host health and well-being. Microbiomes are highly diverse and dynamic, making them challenging to understand. Because many natural systems work as temporal networks, we present an approach that allows identifying meaningful ecological patterns within complex microbiomes: the dynamic core microbiome. On the basis of six marine sponge species sampled monthly over three years, we study the structure, dynamics and stability of their microbiomes. What emerge for each microbiome is a negative relationship between temporal variability and mean abundance. The notion of the dynamic core microbiome allowed us to determine a relevant functional attribute of the microbiome: temporal stability is not determined by the diversity of a host's microbial assemblages, but rather by the density of those microbes that conform its core microbiome. The core microbial interaction network consisted of complementary members interacting weakly with dominance of comensal and amensal interactions that suggests self-regulation as a key determinant of the temporal stability of the microbiome. These interactions have likely coevolved to maintain host functionality and fitness over ecological, and even evolutionary time scales.

Introduction

Microbes form intricate relationships with most animals and plants, with symbiosis postulated as one of the driving forces behind diversifications across the tree of life ( [START_REF] Mcfall-Ngai | Animals in a bacterial world, a new imperative for the life sciences[END_REF]). Research on host-microbe symbioses are typically restricted to highly specialized reciprocal interactions with one or a few microbes interacting with a single host, resulting in mutual benefits for both parties ( [START_REF] Gil | Extreme genome reduction in Buchnera spp.: toward the minimal genome needed for symbiotic life[END_REF][START_REF] Spencer | The winnowing: Establishing the squid-vibrio symbiosis[END_REF]). However, more diverse and complex host-associated microbial communities (hereafter microbiomes) are increasingly found in different plant and animal species ( [START_REF] Mcfall-Ngai | Animals in a bacterial world, a new imperative for the life sciences[END_REF]). This poses a challenge because the pairwise specificity, coevolution and reciprocity of host-microbe interactions might not explain the structure, dynamics and functioning of microbiomes. The mere existence of multiple microbes interacting with a host suggests that interactions among microbes might also be an important driver regulating the overall abundance and composition of microbiomes and their associated ecosystem stability and functions.

The diversity, complexity and highly dynamic nature of microbiomes makes them challenging to understand. We thus require approaches that embrace the complexity and dynamics but still allow identifying meaningful ecological patterns within and across microbiomes. The quest for core microbiomes is a promising avenue. A core microbiome is typically defined cross-sectionally, rather than longitudinally ( [START_REF] Tschöp | Getting to the core of the gut microbiome[END_REF][START_REF] Peter | A core gut microbiome in obese and lean twins[END_REF][START_REF] Shade | Beyond the Venn diagram: the hunt for a core microbiome[END_REF][START_REF] Schmitt | Deep sequencing reveals diversity and community structure of complex microbiota in five Mediterranean sponges[END_REF][START_REF] Ainsworth | The coral core microbiome identifies rare bacterial taxa as ubiquitous endosymbionts[END_REF][START_REF] Astudillo-García | Evaluating the core microbiota in complex communities: A systematic investigation[END_REF][START_REF] Hernandez-Agreda | Defining the Core Microbiome in Corals' Microbial Soup[END_REF]), thus failing to capture temporal dynamics. However, a core microbiome characterized as a set of microbes consistently present over long periods of 2 time ( [START_REF] Björk | Specificity and temporal dynamics of complex bacteria-sponge symbiotic interactions[END_REF][START_REF] Faith | The long-term stability of the human gut microbiota[END_REF][START_REF] Flores | Temporal variability is a personalized feature of the human microbiome[END_REF]), is more likely to be important for the development, health and functioning of its host; for example, several aspects of human health, including autoimmune disorders ( [START_REF] Round | The gut microbiota shapes intestinal immune responses during health and disease[END_REF][START_REF] Scher | The microbiome and rheumatoid arthritis[END_REF]), diabetes ( [START_REF] Qin | A metagenome-wide association study of gut microbiota in type 2 diabetes[END_REF]) and obesity ( [START_REF] Ley | Microbial ecology: human gut microbes associated with obesity[END_REF][START_REF] Peter | A core gut microbiome in obese and lean twins[END_REF]) can be linked to severe shifts in the gut microbiome. Whether these disorders emerge as a consequence of perturbed core microbiomes (defined longitudinally) still remains to be seen. However, arguably, the long-term stability of the core microbiome is likely critical as the persistent occurrence of beneficial microbes and their associated ecosystem functions ensure host health and well-being ( [START_REF] Hartman | Human gut microbiome adopts an alternative state following small bowel transplantation[END_REF][START_REF] Lozupone | Diversity, stability and resilience of the human gut microbiota[END_REF][START_REF] David A Relman | The human microbiome: ecosystem resilence and health[END_REF][START_REF] Cho | The human microbiome: at the interface of health and disease[END_REF][START_REF] Missaghi | Perturbation of the Human Microbiome as a Contributor to Inflammatory Bowel Disease[END_REF]).

Despite the recent realizations that complex microbiomes pervade the tree of life, little is known about microbiome dynamics beyond humans. Here we study the structure, dynamics and stability of microbiomes from six coexisting marine sponges (Porifera) belonging to different orders that were sampled over 36 consecutive months.

Sponges are keystone species in marine coastal areas due to their filter-feeding activities: they regulate primary and secondary production by transferring energy between the pelagic and benthic zones ( [START_REF] Goeij | Surviving in a Marine Desert: The Sponge Loop Retains Resources Within Coral Reefs[END_REF][START_REF] Coppari | The role of Mediterranean sponges in benthic-pelagic coupling processes: Aplysina aerophoba and Axinella polypoides case studies[END_REF]). Despite their constant influx of water, they maintain highly diverse, yet specific microbiomes with very little intraspecific variation ( [START_REF] Thomas | Diversity, structure and convergent evolution of the global sponge microbiome[END_REF]). As Porifera is a sister-group to all other multicellular animals ( [START_REF] Simion | A Large and Consistent Phylogenomic Dataset Supports Sponges as the Sister Group to All Other Animals[END_REF]), their association with microbes is likely the oldest extant form of an host animal-microbe symbiosis ( [START_REF] Taylor | Sponge-associated microorganisms: evolution, evology, and biotechnological potential[END_REF][START_REF] Hentschel | Genomic insights into the marine sponge microbiome[END_REF][START_REF] Yin | Sponge grade body fossil with cellular resolution dating 60 Myr before the Cambrian[END_REF]).

We analyzed sponge species corresponding to two different groups that differ markedly in numerous traits illustrating their dependence upon their association with microbes. This classification is based on the diversity and abundance of microbes they harbor-High and Low Microbial Abundance (HMA and LMA) sponges. The classification pervades sponge host morphology and physiology; LMA hosts have an interior architecture fitted for pumping large volumes of water, whereas HMA hosts are morphologically adapted to harbor denser microbial assemblages within their tissue ( [START_REF] Weisz | Do associated microbial abundances impact marine demosponge pumping rates and tissue densities?[END_REF][START_REF] Lise | Evidence of nitrification and denitrification in high and low microbial abundance sponges[END_REF]). As a result, LMA hosts are more dependent on nutrient uptake from the water column ( [START_REF] Jiménez | Sponges as a source of dissolved inorganic nitrogen: Nitrification mediated by temperate sponges[END_REF][START_REF] Weisz | Do associated microbial abundances impact marine demosponge pumping rates and tissue densities?[END_REF][START_REF] Lise | Evidence of nitrification and denitrification in high and low microbial abundance sponges[END_REF][START_REF] Freeman | Complex interactions between marine sponges and their symbiotic microbial communities[END_REF][START_REF] Maldonado | Chapter three -Nutrient Fluxes Through Sponges: Biology, Budgets, and Ecological Implications[END_REF]) compared to HMA hosts that rely more heavily on nutrients produced by their microbial symbionts ( [START_REF] Freeman | Complex interactions between marine sponges and their symbiotic microbial communities[END_REF][START_REF] Maldonado | Chapter three -Nutrient Fluxes Through Sponges: Biology, Budgets, and Ecological Implications[END_REF][START_REF] Ribes | Functional convergence of microbes associated with temperate marine sponges[END_REF][START_REF] Freeman | Quality or quantity: is nutrient transfer driven more by symbiont identity and productivity than by symbiont abundance?[END_REF][START_REF] Freeman | Metabolic diversity and niche structure in sponges from the Miskito Cays, Honduras[END_REF][START_REF] Poppell | Sponge heterotrophic capacity and bacterial community structure in high-and lowmicrobial abundance sponges[END_REF]). These two sets of hosts provide an ideal system from which we can generalize whether the structure and temporal dynamics of complex microbiomes differ across hosts with different eco-evolutionary characteristics and lifestyles.

Our general aim is to understand the temporal dynamics of complex microbiomes. More specifically, we try to answer: (i) What is the diversity and community structure of each core microbiome and how does it differ from that of the many transient taxa passing through the host? (ii) What are the temporal dynamics and stability of each microbial taxon and their aggregated effect on the core microbiome? And (iii) what are the likely ecological processes that underpin the observed temporal dynamics and stability? We expect the answer to each of these questions to differ across hosts with different lifestyles (i.e., HMA vs. LMA), reflecting their different dependency on their respective micro-3 biomes. In particular, we expect the core microbiomes harbored by HMA hosts to be more diverse, compositionally similar and less variable over time.

Results

We analyzed six microbiomes-Three belonging to host species classified as HMA (Agelas oroides; Chondrosia reniformis; and Petrosia ficiformis), and three from hosts classified as LMA (Axinella damicornis; Dysidea avara; and Crambe crambe) ( [START_REF] Volker Gloeckner | The HMA-LMA Dichotomy Revisited: an Electron Microscopical Survey of 56 Sponge Species[END_REF][START_REF] Erwin | Stable symbionts across the HMA-LMA dichotomy: low seasonal and interannual variation in sponge-associated bacteria from taxonomically diverse hosts[END_REF])

Opportunistic taxa dominate the sponge microbiome

To better understand the commonness and rarity of the taxa occurring throughout the time series, we divided each microbiome into three different temporal assemblages based on the persistence of individual taxa over the 36 consecutive months. Core microbiomes were defined as those taxa that were present in more than (or equal to) 70% of each time series (i.e., persisting 26 months), whereas opportunistic assemblages comprised taxa present in less than (or equal to) 30% of each time series (i.e., persisting  11 months). Intermediately persistent taxa (i.e., those persisting between 12 and 25 months) formed transient assemblages. This resulted in three temporal assemblages each spanning 36 months, with individual taxa not restricted to consecutively occur over the course of the time series. While this categorization is rather arbitrary, it agrees with previous classifications based on temporal occurrences ( [START_REF] Grime | Benefits of plant diversity to ecosystems: immediate, filter and founder effects[END_REF][START_REF] Magurran | Explaining the excess of rare species in natural species abundance distributions[END_REF]). Grime (1998) for example, suggested a similar approach (the mass ratio hypothesis) to better understand the relationship between plant diversity and ecosystem properties, dividing species into three categories-dominants, subordinates and transients-reflecting their contribution to ecosystem biomass, stability and functioning.

As typically observed for macroecological systems (see e.g., [START_REF] Mcgill | Species abundance distributions: Moving beyond single prediction theories to integration within an ecological framework[END_REF]), we found that opportunistic assemblages made up the bulk of the total diversity (i.e., species richness) that resulted from the aggregation of all taxa observed throughout the time series with the core microbiomes only representing a small fraction (1.24%) of this diversity (Table 1).

Core taxa still represented an important fraction of the monthly diversity in some of the hosts; Compared to LMA hosts that on a monthly basis harbored 19 to 25 times more opportunistic taxa, core taxa were only 4 to 5 times less numerous in HMA hosts (Table 1). As common species are increasingly recognized as key contributors to ecosystem functions ([44, 45, 46]), the observed difference in the number of common taxa may mirror variation in the aforementioned symbiont dependency among HMA and LMA host species, respectively.

4

Core microbiomes are mainly driven by changes in abundance

Temporal turnover is an intrinsic property of our definition of core microbiomes, transient and opportunistic assemblages; for example, as opportunistic taxa persist less than 30% of the total time series, these assemblages are bound to heavily fluctuate in microbial composition. Nevertheless, in order to quantify the temporal turnover of each assemblage, we applied a measure that disentangles the two additive determinants of temporal turnover: change in species composition and change in total abundance ( [START_REF] Shimadzu | Measuring temporal turnover in ecological communities[END_REF]). We found that core microbiomes were mainly driven by changes in abundance, whereas transient and opportunistic assemblages were governed by changes in species composition (Fig- ures S1). Due to this large turnover, i.e., new taxa rapidly replace extinct taxa, the opportunistic assemblages displayed an overall higher average evenness (0.83±0.17), measured as Pielou's J over time, compared to the core microbiomes (0.64±0.25).

Microbiome composition and specificity differ across host species and lifestyles

Microbiome composition varied markedly across host species (Figure 1; Phylogenetic composition, Figure S2), with nearly no overlap amongst core microbiomes (Figure 2). The observed community structure was different from the expectation generated by a null-model (Figure S3; See Supplementary material for more information). In addition, HMA and LMA host species displayed very different taxonomic profiles (Figure 3). HMA hosts harbored three dominant phyla that accounted for roughly half of their diversity: Chloroflexi; Actinobacteria; and Acidobacteria.

Despite the dominance of these three phyla, each HMA host still kept a unique taxonomic fingerprint by harboring other phyla, such as Gemmatimonadetes, Nitrospira, and Bacteroidetes. The core microbiomes of LMA hosts, on the other hand, were largely dominated by taxa from a single phylum-Proteobacteria.

Sponges typically harbor certain taxa highly specific to the phylum Porifera. These "sponge-specific" 16S rRNA gene sequence clusters are monophyletic and span 14 known bacterial and archaeal phyla ( [START_REF] Hentschel | Molecular Evidence for a Uniform Microbial Community in Sponges from Different Oceans Molecular Evidence for a Uniform Microbial Community in Sponges from Different Oceans[END_REF][START_REF] Simister | Sponge-specific clusters revisited: A comprehensive phylogeny of sponge-associated microorganisms[END_REF]). Taxa that fall into these sponge-specific clusters are only detected at very low abundances outside of sponge hosts, e.g., in the seawater and sediment ( [START_REF] Taylor | Sponge-associated microorganisms: evolution, evology, and biotechnological potential[END_REF][START_REF] Taylor | Sponge-specific' bacteria are widespread (but rare) in diverse marine environments[END_REF][START_REF] Thomas | Diversity, structure and convergent evolution of the global sponge microbiome[END_REF]). Some of these taxa are transmitted vertically from parent to offspring, suggesting spongemicrobe coevolution and cospeciation ( [START_REF] Robert | Chapter two -Sponge-Microbe Symbioses: Recent Advances and New Directions[END_REF]). We found that the core microbiomes of HMA hosts harbored a larger proportion of taxa that corresponded to sponge-specific clusters than the cores of LMA hosts (Table 2; Figure S4). Furthermore, these taxa had a higher average monthly abundance in the core microbiomes compared to the transient and opportunistic assemblages (Table S1). Then, despite the fact that identifying what constitutes a core microbiome remains elusive ([6, 25, 9]), the definition of the core used here (i.e., taxa present in at least 70% of the time series) provided a clear distinction in species composition between the core microbiomes and the transient and opportunistic assemblages (Figure 2) that displayed prominent differences in the average monthly abundance of sponge-specific taxa (Table S1) and in their temporal turnover (Figures S1).

The innate immune defense of some sponge species can differentiate between pathogens, food bacteria and commensals in a manner similar to the adaptive immune system of vertebrates ([52, 53, 54, 55, 56, 57]). Thus, the above findings suggests that sponge-specific clusters, although present in the water column and sediment as part of the rare biosphere ( [START_REF] Webster | Deep sequencing reveals exceptional diversity and modes of transmission for bacterial sponge symbionts[END_REF][START_REF] Taylor | Sponge-specific' bacteria are widespread (but rare) in diverse marine environments[END_REF]), represent taxa important for host functioning.

Core taxa are more stable and abundant

The unusual temporal coverage of our dataset allowed us to address a number of fundamental questions related to the temporal variability within and across microbiomes. We computed the mean abundance and the coefficient of variation CV = µ for each taxon's abundance trajectory over time. What emerged for each microbiome was a negative relationship between the coefficient of variation and the log mean abundance: individual core and transient taxa were more stable and abundant than opportunistic taxa. However, in a few microbiomes, abundant conditionally rare taxa (CRTs) greatly outnumbered the abundance of core taxa. CRTs are extremely rare or below detection limit throughout most of the time series, but occasionally reach high abundance ( [START_REF] Shade | Temporal patterns of rarity provide a more complete view of microbial diversity[END_REF]). Our results suggest that the presence of these taxa can have a negative impact on the stability of the focal core microbiome (Figure 4 AB; Figure S5).

There are several non-mutually exclusive mechanisms that can underpin the observed relationship between temporal variability and mean abundance. First, in agreement with recent observations on a broad range of marine communities, neutrality alone is unlikely to explain the high abundance of the most common species ( [START_REF] Connolly | Commonness and rarity in the marine biosphere[END_REF]). Second, lower variability can imply the presence of self-regulation (i.e., density dependent processes) at higher population abundances ( [START_REF] Hanski | Density Dependence, Regulation and Variability in Animal Populations[END_REF][START_REF] Peter | Direct evidence that density-dependent regulation underpins the temporal stability of abundant species in a diverse animal community[END_REF]). Furthermore, on a log-log scale, our observed relationship between temporal variability (now variance 2 ) and mean abundance describes Taylor's power law ( [START_REF] Taylor | Aggregation, Variance and the Mean[END_REF]). The null expectation of the temporal version of Taylor's law states that the variance scales with the mean abundance following a power law with an exponent equal to 2. This implies that population or community variability is constant. Conversely, if the exponent is less than 2, then the variability in population abundance decreases with increasing mean population abundance, as seen in Figure 4 CD and Figure S6. Third, for large population sizes, as is characteristic of microbial communities like those studied here, environmental stochasticity is more likely than demographic stochasticity to reduce Taylor's power law exponent ( [START_REF] Anderson | Variability in the abundance of animal and plant species[END_REF][START_REF] Kilpatrick | Species interactions can explain Taylor's power law for ecological time series[END_REF]). Fourth, even weak interactions among species can reduce this exponent; for example, interspecific interactions, such as competition typically lead to smaller fluctuations of common versus rare species within communities ( [START_REF] Kilpatrick | Species interactions can explain Taylor's power law for ecological time series[END_REF])-a 6 result that could be extended to other types of interactions. Below, we explore the contribution of these mechanism to microbiome dynamics and stability in more detail.

Core microbiome density, not diversity, begets its stability

We aggregated all individual microbial population abundances within each temporal assemblage-core, transient, and opportunistic (Figure 5). This revealed two markedly different temporal dynamics across microbiomes: in the hosts A. oroides, C. reniformis and C. crambe, core microbiomes were very dense, i.e., they accounted for the majority of microbiome relative abundance. In contrast, the core microbiomes of hosts D. avara, A. damicornis and P. ficiformis were sparser, and instead transient and/or opportunistic assemblages dominated microbiome relative abundance. We found that dense cores were more stable over time than sparse cores, measured as community-level invariability (the inverse of variability, [START_REF] Haegeman | Resilience, invariability, and ecological stability across levels of organization[END_REF]; Figure 5). The association between core density and stability was also robust to a more stringent definition ( 85%) of the core microbiome (Figure S6). For some microbiomes, no taxa persisted on or above this threshold.

Community-level stability showed a weak (non-significant) relationship with diversity (species richness; Figure S8 A). Instead, we observed a significant positive relationship between community-level stability and the median relative abundance (Figure S8 B). This agrees with some studies that demonstrated a positive relationship between communitylevel stability and the relative abundance or biomass of common plant species ( [START_REF] Polley | Dominant species constrain effects of species diversity on temporal variability in biomass production of tallgrass prairie[END_REF][START_REF] Sasaki | Dominant species, rather than diversity, regulates temporal stability of plant communities[END_REF][START_REF] Ma | Climate warming reduces the temporal stability of plant community biomass production[END_REF]). Importantly, this suggests that the effect of diversity on stability can be constrained by the dynamics of only a few common and abundant species, as seen in the case of our high-density cores (Figure 5); for example, host species P. ficiformis harbored the second largest core diversity (S=40, Table 1) that still resulted in the least stable low-density core. On the other hand, sponge species C. crambe achieved a stable high-density core by harboring only a few core taxa (S=8, Table 1).

Furthermore, in agreement with the HMA-LMA dichotomy, the metabolic profiles of P. ficiformis and C. crambe yet match those of other archetypal HMA and LMA hosts, respectively ( [START_REF] Morganti | Trophic niche separation facilitates co-existence of high and low microbial abundance sponges is revealed by in situ study of carbon and nitrogen fluxes[END_REF]). While P. ficiformis harbored an unstable low-density core, its core was nonetheless comprised of the largest proportion of sponge-specific clusters (Table 2;

Figure S4), suggesting that microbiome stability may not always be a prerequisite for host functioning.

Intraspecific interactions drive core microbiome dynamics

In order to further disentangle the main drivers of core microbiome temporal dynamics, we used a model that decomposes temporal fluctuations in species abundances into three contributions-Interspecific interactions, intraspecific 7 interactions, and environmental variability. Environmental variability includes apart from ecological drift, the aggregated effects of the host on the microbial consortium, as well as the external environment acting on the host. While we included several environmental covariates that had been recorded over the 36 months, such as water temperature, salinity and nutrients, our model is able to capture unmeasured effects through the use of latent variables (see Methods for more information).

We found that intraspecific interactions explained the largest proportion of variation, suggesting that all core taxa on average, experienced strong self-regulation (Figure 6 A). There was a marked difference between HMA and LMA host species in terms of the variation explained by interspecific interactions: In the core microbiomes of the HMA hosts, interspecific interactions had a relatively large effect on the dynamics, whilst this effect was almost negligible in the LMA host species (Figure 6 A).

The main drivers we identified for core microbiome temporal dynamics somewhat differ from those commonly reported for other large species-rich systems; for example, studies assuming Lotka-Volterra dynamics in macroecological communities ( [START_REF] Crispin | A multispecies perspective on ecological impacts of climatic forcing[END_REF][START_REF] Almaraz | Size-mediated non-trophic interactions and stochastic predation drive assembly and dynamics in a seabird community[END_REF][START_REF] Martorell | Testing the roles of competition, facilitation and stochasticity on community structure in a species-rich assemblage[END_REF][START_REF] Crone | Contrasting effects of spatial heterogeneity and environmental stochasticity on population dynamics of a perennial wildflower[END_REF]) or other types of time series analyses applied on free-living marine microbial communities ( [START_REF] Gilbert | Defining seasonal marine microbial community dynamics[END_REF][START_REF] Fuhrman | Marine microbial community dynamics and their ecological interpretation[END_REF][START_REF] Bunse | Marine Bacterioplankton Seasonal Succession Dynamics[END_REF]) often report environmental variability as the major factor affecting population dynamics. While environmental variability was an important factor, intraspecific interactions was the single most important driver affecting the dynamics across all core microbiomes. Furthermore, the modeled environmental covariates only accounted for a small fraction of the variation explained by environmental variability (Figure 6 A), indicating that the sponge microbiome may experience a reduced influence of the external environment acting on the host.

A large body of theoretical and empirical literature suggests that self-regulation, which we equate here to intraspecific interactions, is a key determinant of temporal stability ( [START_REF] Hanski | Density Dependence, Regulation and Variability in Animal Populations[END_REF][START_REF] Peter | Direct evidence that density-dependent regulation underpins the temporal stability of abundant species in a diverse animal community[END_REF]). The simple premise of this regulatory mechanisms is that species abundances decrease per capita growth rates when population abundances are high and vice versa. Recently Barabás and colleagues (2017) showed that stability requires the majority of species in a community to experience negative self-effects. Although this work focuses on asymptotic stability (i.e., whether small perturbations of species' abundances away from an equilibrium point tend to be dampened, with the system returning to equilibrium) and not on temporal variability, their results in combination with ours suggest that strong self-regulation should also be the norm within microbiomes ( [START_REF] Barabás | Self-regulation and the stability of large ecological networks[END_REF]).

The core microbiomes are comprised of weak unilateral interactions

The nature and strength of interspecific interactions are pivotal to understand the dynamics of species interaction networks ( [START_REF] Coyte | The ecology of the microbiome: Networks, competition, and stability[END_REF][START_REF] Mougi | Food-web complexity, meta-community complexity and community stability[END_REF]). Our unique temporal series allows for determining how microbes interact with each other within 8 different host species. Because interspecific interactions were almost negligible in the LMA hosts, we analyzed the most credible network structure of each HMA core microbiome (Figure 7; Figures S9). For these three networks, only a small fraction among the possible interactions were likely to occur (Figures S10), resulting in low network connectance (5-7%). The networks had a skewed distribution of interaction strengths toward many weak and a few strong interactions (Figures S11). This pattern mimics that found in empirical food-webs [START_REF] Paine | Food-web analysis through field measurement of per capita interaction strength[END_REF][START_REF] Mark | Predator-prey body size, interaction strength and the stability of a real food web[END_REF]) and mutualistic networks ( [START_REF] Vazquez | The strength of plant-pollinator interactions[END_REF]). Theory shows that skewed distributions of interaction strengths beget stability, and arises during the assembly of persistent communities ( [START_REF] Mccann | Weak trophic interactions and the balance of nature[END_REF][START_REF] Allesina | Stability criteria for complex ecosystems[END_REF]). A. oroides and C. reniformis were largely dominated by commensal interactions, whereas the core microbiome of host P. ficiformis had a higher frequency of amensalism. It was also the only core network that had competitive interactions among its members, although at very low frequencies. Sponge-specific clusters, although prevalent in the core microbiomes of HMA hosts were not more connected than other nodes within these networks Certainly, opportunistic and transient taxa can influence the dynamics of the core microbiome. However, we focused on inferring interspecific interactions among core taxa for three main reasons: (1) Computation-Including non-core taxa would increase network size ten folds, and thus not be computationally feasible; (2) Information-There is very little information in taxa only occurring in one or a few timesteps, thus the inference for those taxa would be very poor; and (3) Ecology-Interspecific interactions require species to frequently co-occur. The many occasional taxa observed throughout the time series, especially those occurring at high densities (i.e., the CRTs) are likely acting together with ecological drift as sources of stochastic variability. In our model, the effects of unmodeled interspecific interactions are captured by the inclusion of the latent variables (see Methods for more information).

Microbial transmission mode can affect microbiome stability and functioning

High-density cores (i.e., those core microbiomes whose taxa accounted for the majority of abundance) were found in sponge species that transmit microbes vertically from adult to larvae ( [START_REF] Levi | Embryogenese de Chondrosia reniformis (Nardo), demosponge ovipare, et transmission des bacteries symbiotiques[END_REF][START_REF] María | Morphology and Ultrastructure of the Swimming Larvae of Crambe crambe (Demospongiae, Poecilosclerida)[END_REF][START_REF] Schmitt | Molecular microbial diversity survey of sponge reproductive stages and mechanistic insights into vertical transmission of microbial symbionts[END_REF]), whereas low-density cores represent sponge species with larvae largely deprived of microbes ( [START_REF] Lepore | The ultrastructure of the mature oocyte and the nurse cells of the ceractinomorpha Petrosia ficiformis[END_REF][START_REF] Riesgo | Differences in reproductive timing among sponges sharing habitat and thermal regime[END_REF][START_REF] Maldonado | Gametogenesis, embryogenesis, and larval features of the oviparous sponge Petrosia ficiformis (Haplosclerida, Demospongiae)[END_REF]).

9

Vertical transmission provides an evolutionary mechanism for preserving particular combinations of microbes, including their interaction structure and the ecosystem functions that emerge from it ( [START_REF] Thompson | The geographic mosaic of coevolution[END_REF]). Evidence from other ecological communities, including the human gut microbiome, suggest that priority effects-the order and timing of species arrivals-determine the interactions among species, and in turn also community assembly and stability ( [START_REF] Chase | Stochastic community assembly causes higher biodiversity in more productive environments[END_REF][START_REF] Tadashi | Historical Contingency in Community Assembly: Integrating Niches, Species Pools, and Priority Effects[END_REF][START_REF] Sprockett | Role of priority effects in the early-life assembly of the gut microbiota[END_REF]). The process of vertical inheritance of microbes likely has similar outcomes as priority effects. Furthermore, as pathogenic microbes use cooperative secretions that modify their environment to enhance their growth and expansion ( [START_REF] Nogueira | Horizontal Gene Transfer of the Secretome Drives the Evolution of Bacterial Cooperation and Virulence[END_REF][START_REF] Mcnally | Cooperative secretions facilitate host range expansion in bacteria[END_REF]), it is reasonable to assume that commensal microbes do too. We therefore hypothesize that the complementary set of microbes that are vertically transmitted from parent to offspring pre-empt the initial host niche by fast reaching carrying capacity, while simultaneously modifying it in their favor. This would inhibit the subsequent colonization of some microbes, while facilitating the establishment of certain others.

As previously mentioned, host species P. ficiformis harbored an unstable low-density core, yet it harbored the largest consortia of sponge-specific clusters within its core. This does not only suggest that microbiome stability may not always be a requirement for host functioning, it also indicates that host species P. ficiformis achieves HMA archetypal functional characteristics by means of horizontally selecting commensal microbes from the water column, many of which fall into sponge-specific clusters. These sponge-specific microbes are examples of conditionally rare taxa (CRTs) outside the sponge host, but through some recognition mechanism are allowed to flourish within the host.

Conclusion

Natural systems commonly work as temporal networks ( [START_REF] Li | The fundamental advantages of temporal networks[END_REF]). To increase our understanding of the processes governing microbiome assembly, stability and functioning, this study highlights the importance of defining the core microbiome longitudinally (i.e., analysing temporal dynamics) rather than cross-sectionally. Most studies to date have focused on the human microbiome (see e.g., [START_REF] Caporaso | Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample[END_REF][START_REF] Flores | Temporal variability is a personalized feature of the human microbiome[END_REF][START_REF] Backhed | Dynamics and stabilization of the human gut microbiome during the first year of life[END_REF][START_REF] Coyte | The ecology of the microbiome: Networks, competition, and stability[END_REF][START_REF] Sprockett | Role of priority effects in the early-life assembly of the gut microbiota[END_REF]), with only a few studies exploring temporal dynamics in other host systems ([11, 106, 107]). By focusing on the six most common sponge species of the temperate Mediterranean benthic community, we successfully characterized microbiome dynamics and stability under natural conditions.

The observed negative relationship between temporal variability and abundance that emerged from the analyses of each microbiome is conducive to the notion of the dynamic core microbiome. This notion allowed us to determine that irrespective of host's eco-evolutionary characteristics and lifestyles, it is the density of the core microbiome rather than its diversity that determines the stability of the sponge microbiome. We hypothesize that priority effects 10 mediated by vertical transmission underpins this pattern, which may further suggest that high-density cores confer hosts resistance against the establishment of occasional taxa that sponges are constantly exposed to through their filter-feeding activities. The core microbiome has been proposed as the common taxa shared among microbiomes in an habitat ( [START_REF] Peter | The human microbiome project: exploring the microbial part of ourselves in a changing world[END_REF][START_REF] Shade | Beyond the Venn diagram: the hunt for a core microbiome[END_REF]). Our finding reveals a relevant functional attribute that constitutes an step forward towards the identification and characterization of the so-called core microbiome.

We further found that intraspecific self-regulation is much more important for microbiome dynamics than environmental forcing. The most credible core interaction networks consisted of members interacting weakly with each other with a dominance of comensal and amensal interactions. Altogether, this suggests that host-associated microbiome dynamics and its emerging interaction structures differ from the temporal dynamics of free-living microbial communities. These interactions have likely coevolved to maintain host functionality and fitness over ecological, and even evolutionary time scales.

We have focused on compositional stability over time-a notion inclusive of variation in both species relative abundances and species composition. More generally, compositional and functional stability often show complex interdependencies ( [START_REF] Hillebrand | Decomposing multiple dimensions of stability in global change experiments[END_REF]), because, for example, compensatory dynamics between species may reduce functional loss within communities ( [START_REF] Connell | Resisting regime-shifts: the stabilising effect of compensatory processes[END_REF]). For most host species studied here, our results suggest that there may exist a positive relationship between compositional and functional stability. However, for host-associated microbiomes in general, it is difficult to disentangle the effects of compositional stability from the host's own ability to control the identity and abundance of its microbes on the overall functioning. For one pair of hosts, our findings of high-density and low-density cores did not match the notion of a positive relationship between compositional and functional stability. This suggests that, at least for some host species, functioning is achieved by other mechanisms that do not require compositional stability. The relationship between microbiome compositional stability and host functionality will benefit from further investigation that focuses on the temporal dynamics and functioning of the core microbiome.

Methods

Sponge collection

The sponge species Agelas oroides, Chondrosia reniformis, Petrosia ficiformis, Axinella damicornis, Dysidea avara

DNA extraction and sequencing

Following the manufacturers Animal Tissue protocol, 16S rRNA gene sequences were PCR-amplified from sponge samples (n=648) and seawater filters (n=108) using the DNeasy tissue kit (Qiagen, CA, USA), and subsequently submitted to the Research and Testing Laboratory (Lubbock, TX, USA) for gene amplicon pyrosequencing. Samples were amplified with primer 28F and amplicons were sequenced using 454 Titanium chemistry (Roche, CT, USA), producing 250 base pair read lengths in the 5'!3' direction.

Analysis of sequencing data

454 reads were processed in mothur v. 1.29.2 ([111]). Raw reads were pooled from replicates belonging to the same sponge species. Fasta, qual and flow files were extracted from binary sff files; sffinfo(...,flow=T). Flow files were then filtered based on barcodes to speed-up the proceeding de-noising process; trim.flow. Sequences were de-noised;

shhh.flows(..., lookup= LookUp_Titanium.pat). The LookUp-file is necessary and specific to the 454 technology used.

Next the barcode and primer sequences were removed together with sequences shorter than 200bp and/or contained homopolymers longer than 8bp; trim.seqs(..., pdiffs =2, bdiffs =1, maxhomop =8, minlength =200). In order to minimise computational effort, files were reduced to non identical sequences; unique.seqs. Non redundant sequences were aligned to SILVA 102 reference alignment with default kmer search and Needleman-Wunsch algorithm; align.seqs(..., flip =F). Non overlapping sequences were removed; screen.seqs(..., optimize= end, start= 1044, criteria = 95), in addition to empty columns that were introduced from the alignment process; filter.seqs(...,vertical =T, trump =.).

Aligned sequences were reduced to non redundant sequences; unique.seqs. To further reduce amplification errors, less abundant sequences were binned to more abundant sequences if they were within 2bp of a difference; pre.cluster(..., diffs =2). Chimeric sequences were identified; chimera.uchime(..., dereplicate =T) and removed; remove.seqs. Se-quences were classified using the RDP reference taxonomy; classify.seqs(..., template =trainset9_ 032012.pds.fasta, taxonomy =trainset9_ 032012.pds.tax, cutoff =80), and non bacterial lineages were removed; remove.lineage(..., taxon= Mitochondria-Chloroplast-Archaea-Eukaryota-unknown). We calculated pairwise distances between aligned sequences; dist.seqs(..., cutoff =0.050).

Due to an uneven sequence distribution across samples, we pooled sequences across the three host species replicates prior to OTU clustering. As we found a plate effect on the number of sequences, we sub-sampled 1500 sequences from the resulting monthly samples. This number corresponded to the average of the three lowest host-plate averages.

Sequences were thereafter clustered into OTUs defined at 97% similarity; classify.otu(..., label=0.030) and outputted to an OTU-table (.shared-file); make.shared(..., label=0.030).

Identification of sponge-specific clusters

A representative sequence from each OTU was taxonomically assigned using a BLAST 62 search against a curated ARB-SILVA database containing 178 previously identified sponge-specific clusters (SC) ( [START_REF] Simister | Sponge-specific clusters revisited: A comprehensive phylogeny of sponge-associated microorganisms[END_REF]). For each BLAST search, the 10 best hits were aligned to determine sequence similarities. The most similar OTU sequence to the respective reference sequence within the database was then assigned to an SC based on a 75% similarity threshold:

(i) a sequence was only assigned to any given SC if its similarity was higher to the members of the cluster than to sequences outside the cluster; and (ii) if its similarity to the most similar sequence within the cluster was above 75%.

A majority rule was applied in cases where the assignment of the most similar sequences was inconsistent, and the OTU sequence was only assigned to the SC if at least 60% of the reference sequences were affiliated with the cluster.

Microbiome stability and dynamics

Temporal turnover

We applied a measure of temporal turnover that describes the extent to which individual OTUs and consequently the microbiome changes over time ( [START_REF] Shimadzu | Measuring temporal turnover in ecological communities[END_REF]. Importantly, this measure decomposes abundance fluctuations into two additive contributions of change due to microbiome composition and total abundance.

13

Total turnover D between times t and u, (u > t) is defined as

D(t : u) = S X i=1 d i (t : u) = S X i=1 log i,u i,t ! p i,t (1) = S X i=1 log p i,t p i,u ! p i,t + u t ! (2) = D 1 (p t : p u ) + D 2 ( t : u ) (3)
where t = P S i=1 i,t represent the sum of the expected total abundance of each OTU in the microbiome. The expected abundance of OTU i in time t, i.e., i,t ,i = 1, 2,..., S is unknown and therefore needs to be estimated from an observed time series. p i,t represents the relative abundance of OTU i in time t, and is calculated as

p i,t = ( i,t t
). As such, the total turnover D can be decomposed into D 1 which is related to the amount of change in microbiome composition, and D 2 reflecting the amount of change in total abundance.

As noted above, the expected abundance needs to be estimated. We thus modeled each time series of time series of high-throughput DNA sequence N i,t assuming a Poisson log-linear model with a time-varying mean parameter i,t

N i,t ⇠ Pois( i,t ) (4)
log i,t = N c X j=1 X t,k k,j (5)
where X k,t is a time series of k = 1, 2,..., N c environmental covariates, and k,j is the corresponding regression coefficient that needs to be estimated. We included temperature, salinity, chlorophyll, bacterial cell density, nitrite (N O 2 ), ammonia (N H 4 ), and phosphate (PO 4 ) as the N c environmental covariates. All covariates where standardized to have zero-mean and unit variance.

Community-level stability

We applied a newly developed measures of temporal stability defined as invariability ( [START_REF] Haegeman | Resilience, invariability, and ecological stability across levels of organization[END_REF]) at the community-level, i.e., for each core microbiome, transient and opportunistic assemblages, respectively. Invariability at the communitylevel is defined as

I = 1 CV (N tot ) 2 = ( b N tot ) 2 V ar(N tot ) (6)
where b N tot and b N i denote the average of the total abundance and the average of the abundance of OTU i over the time series, respectively.

Core dynamics and ecological interactions

We adapted a multivariate first-order autoregressive model assuming Gompertz population dynamics ( [START_REF] Crispin | What drives community dynamics?[END_REF]) to model core dynamics and to infer ecological interactions among core taxa from time series of high-throughput DNA sequence counts. Compared to the commonly used Generalized Lotka-Volterra (GLV) model (see e.g., [START_REF] Stein | Ecological Modeling from Time-Series Inference: Insight into Dynamics and Stability of Intestinal Microbiota[END_REF]), the Gompertz model is linear on the log-scale and thus a good approximation of nonlinear dynamics. Furthermore, a severe problem with even the simplest GLV model is that the number of parameters that need to be estimated often is larger than the number of available data points. Mutshinda and colleagues (2009) addressed this problem by using Gibbs Variable Selection ( [START_REF] Robert | A review of bayesian variable selection methods: What, how and which[END_REF]) in order to induce sparseness: the N (N + 1) interspecific interactions are constrained so that most are shrunk to zero. We further adapted this framework to address two additional challenges of microbiome data.

First, to model covariances between a large number of taxa using a standard multivariate random effect is computationally challenging: the number of parameters that has to be estimated when assuming a completely unstructured covariance matrix increases quadratically as the number of taxa increases. We therefore parameterized the residual covariance matrix using a latent variable approach. Latent variables also capture unmeasured effects (e.g., host effects and unmodeled interspecific interactions) which if not accounted for may lead to erroneous inference ( [START_REF] Warton | So Many Variables: Joint Modeling in Community Ecology[END_REF]).

Second, high-throughput DNA sequencing produces compositional data, i.e., non-negative counts with an arbitrary sum imposed by the sequencing platform, which can produce spurious correlations if not properly accounted for (see e.g., [START_REF] Li | Microbiome,Metagenomics,and High-Dimensional Compositional Data Analysis[END_REF][START_REF] Matthew | Compositional data analysis of the microbiome: fundamentals, tools, and challenges[END_REF][START_REF] Gloor | Microbiome Datasets Are Compositional: And This Is Not Optional[END_REF]). We therefore used a log-linear model with a Poisson likelihood that includes a random effect accounting for sample size. This model produces a likelihood equivalent to that of the logistic normal multinomial model but is more convenient for computation and estimation ( [START_REF] Stuart | The Multinomial-Poisson Transformation[END_REF]).

Process model

To accommodate the time series, we refer to OTU i 2 {1,..., I} in time point t 2 {1,..., T }. Let N (µ, 2 ) denote a univariate normal distribution with mean µ and variance 2 , and analogously, let MV N (µ, ⌃) denote a multivariate normal distribution with mean vector and covariance matrix ⌃. If we denote n i,t ⇤ as the expectation of n i,t which is the natural logarithm of the observed time series N i,t , then on the natural logarithmic scale we have the expected number of 16S rRNA gene sequences from OTU i in time point t within a given core microbiome described by 15

n i,t ⇤ |n i,t 1 = n i,t 1 + r i  1 S X j=1 ↵ i,j n j,t 1 k i + N c X j=1 X t,k k,j + ✏ i,t (7) t = 2, 3,..., T ; k = 1, 2,..., N c
where we assume r i ⇠ N (0, 10) and K i ⇠ Exp(1). The coefficients measuring each taxon's response the k th environmental covariate are assumed k,j ⇠ N (0, 100). We assumed correlated residual responses to the environment by included temperature, salinity, chlorophyll, bacterial cell density, nitrite (N O 2 ), ammonia (N H 4 ) and phosphate (PO 4 ) as the N c environmental covariates potentially influencing the modeled core OTUs. All covariates where standardized to have zero-mean and unit variance. ✏ i,t ⇠ MV N (0, ⌃) represents the residual variance, where ⌃ is the residual covariance matrix which was parameterized using latent variables.

Latent variables

As previously stated, to improve statistical efficiency we parameterized the residual covariance matrix ⌃ using a latent variables approach, where the index q runs over the q = 1,..., 2 latent variables.

✏ i,t = 2 X q=1 ⌘ t,q q,i + i,j " i,t (8)
" i,t ⇠ N (0, 2 i )

⌃ = > + diag( 2 i ) (9)
were ⌘ denote the latent variables and their corresponding factor loadings. Both are assigned standard normal priors N (0, 1) with the assumption of zero mean and unit variance to fix the location (see chapter 5, [START_REF] Skrondal | Generalized Latent Variable Modeling: Multilevel, Longitudinal, and Structural Equation Models[END_REF]). i,j denote the Kronecker's delta such that i,i = 1 and i,j = 0 for i , j. Thus, the covariance matrix ⌃ can be computed from the factor loadings (Equation 9). The diagonal elements of the covariance matrix ⌃ quantify the amount of residual variation for OTU i not captured by the modeled environmental covariates.

Observation model

The

n i,t ⇠ MV N (n i,t ⇤, 2 i ) (12)
where log N t is an offset accounting for the observed total number of DNA sequences in time t, while µN t ⇠ N (0, 100) represent the random effect accounting for sample size in time t. Both thus represent the total abundance in time t, and account for the aforementioned compositional nature of high-throughput DNA sequencing data.

Gibbs Variable Selection

As mention above, we used Gibbs Variable Selection ( [START_REF] Robert | A review of bayesian variable selection methods: What, how and which[END_REF]) method in order to constrain the model to only use interspecific interaction coefficients ↵ i,j for which there were strong support in the data. This was achieved by introducing a binary indicator variable i,j for i , j, and assuming i,j ⇠ Bernoulli(p), such that i,j = 1 when OTU j is included in the dynamics of OTU i, and i,j = 0 otherwise. Where there was low support for ↵ i,j in the data, i,j = 0 and the corresponding interaction was excluded from the model. On the other hand, when i,j = 1, ↵ i,j was freely estimated from the data. The parameter p represents our prior belief about the proportion of realized interspecific interactions:

we set p=0.1, thus assuming that 90% of all interspecific interactions were zero.

We used Markov chain Monte Carlo (MCMC) simulation methods through JAGS ( [START_REF] Plummer | JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling[END_REF]) in R ( [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF]) using the runjags package ( [START_REF] Denwood | runjags: An R package providing interface utilities, model templates, parallel computing methods and additional distributions for MCMC models in JAGS[END_REF]) to sample from the joint posterior distribution of all the model parameters. We ran 10 independent chains with dispersed initial values for 5e6 iterations, discarding the first 2e6 samples of each chain as burn-in and thinned the remainder to every 50th sample. We evaluated convergence of model parameters by visually inspecting trace and density plots using the packages coda ( [START_REF] Plummer | CODA: Convergence Diagnosis and Output Analysis for MCMC[END_REF]) and mcmcplots ( [START_REF] Mckay | Create Plots from MCMC Output[END_REF]), as well as using the Geweke diagnostic ( [START_REF] Geweke | Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments[END_REF]). In addition, to ensure good mixing of the parameter ↵ i,j , we calculated the number of jumps the parameter i,j made between its two states, 0 and 1.

Variance partitioning

The total variance V i affecting the dynamics of core OTU i can be decomposed into additive sources reflecting interspecific interactions, intraspecific interactions, and environmental variability (i.e., measured environmental covariates and residual variation) as follows,

V i = Interspecific interactions z }| { " r i K i # 2 X j,i v j,j ↵ 2 i,j + Intraspecific interactions z }| { " r i K i # 2 v i,i + Environmental variability z }| { N c X i=1 2 i,q + 2 i ( 13
)
where v i,i represent the stationary variance for n i (Equation 11), 2 i,k the variance attributable to each k covariate, and 2 i correspond to the residual variance (the diagonal elements of ⌃, Equation 9). As a consequence of Equation 13, the proportion of variation attributed to e.g., interspecific interactions can be calculated as

2 i Inter = 8 > > > < > > > : " r i K i # 2 X j,i v j,j ↵ 2 i,i 9 > > > = > > > ; , V i ( 14
)
Finding the most credible network structure

To characterize each core microbiome network, we analyzed the interaction structure of the marginal posterior distribution for the interaction coefficient ↵ i,j . The parameter ↵ i,j is a probability distribution, thus containing the probability of OTU j having a per capita effect on the growth rate of OTU i. This means that we can decompose ↵ i,j into two marginal posterior distributions-one for interaction sign and another one for interaction strength. We constructed the core microbiome networks for each HMA host as means of visualizing the most credible network structures. This was done by mapping the marginal posterior distribution of the average number of links onto ↵ i,j , thus extracting the marginal posterior average number of links with the highest probability of non-zero interactions. As a way of further validating network structure, we compared the marginal posterior distribution of connectance to the empirical connectance of each constructed network. Notes: For each assemblage, the first row displays the total number of unique taxa (S), while the second row shows the monthly average (±SD) number of coexisting taxa. The last row shows the total species richness (S T ) of each microbiome. LMA HMA Figure 3: Taxonomic profiles of core microbiomes across host species and lifestyles (HMA vs. LMA). Taxonomic classification at the phylum level of core taxa and their relative contribution to diversity (species richness) for each core microbiome. The core microbiomes of HMA hosts harbored a larger taxonomic diversity than those of LMA hosts. Colors correspond to different bacterial phyla and dash and solid lines represent positive and negative interactions, respectively. Nodes marked with SC correspond to taxa that assigned to sponge-specific clusters. See Figure S10 for the corresponding link probabilities between OTU i and j, and Figure S9 for the networks belonging to the core microbiomes of HMA hosts C. reniformis and P. ficiformis. See Methods for how the networks were constructed.

Tables and Figures

and

  Crambe crambe were collected monthly from March 2009 until February 2012 close to the Islas Medas marine reserve in the NW Mediterranean Sea 42 3 0 0 00 N , 3 13 0 0 00 E by SCUBA at depths between 10-15 metres. Three replicates per sponge species (i.e., different individuals per sampling time) and ambient seawater samples were collected.11The collected sponge species belong to six different orders that represent common members of the Mediterranean benthic community. Each sponge species were identified based on distinct morphological features. Specimens were sublethally sampled and excised fragments were placed in separate plastic bottles and brought to the surface were they were frozen in liquid nitrogen until DNA extractions. Samples of the ambient water were taken at 5 m depth and poured into three separate 5 L jars. All samples were stored at -80 C until DNA extraction. Aliquots of seawater (300-500 mL each, 1 aliquot per sample jar) were concentrated on 0.2 µm polycarbonate filters, submerged in lysis buffer and stored at -80 C until DNA extraction.

33

 33

Figure 4 :Figure 5 :

 45 Figure4: The top panel shows the relationship between temporal variability (CV) and mean log abundance for each taxon's abundance trajectory over time, for host (A) A. oroides and (B) A. damicornis. Overlaying points have been separated jitter (random noise). Opportunistic, transient and core taxa are represented by stars, diamonds and circles, respectively. Individual core and transient taxa are more stable (Kruskal-Wallis test: H =2198, df=2, P<0.001 twotailed; Dunn's post-hoc test with bonferroni correction; P <0.001) and abundant (Kruskal-Wallis test: H =1694, df= 2, P<0.001; Dunn's post-hoc test with bonferroni correction; P< 0.001 two-tailed) than opportunistic taxa. The bottom panel shows the relationship between log variance and log mean abundance for each taxon's abundance trajectory over time (Taylor's power law), shown for host (C) A. oroides and (D) A. damicornis. Solid lines represent the null expectation of Taylor's power law, i.e., an exponent a=2, and the dashed lines correspond to the exponent of each host's microbiome: 1.52 for A. oroides and 1.77 for A. damicornis.

Figure 7 :

 7 Figure7: The most credible network structure for HMA host A. oroides's core microbiome. Nodes represent core taxa and links their inferred ecological interactions. Node size is scaled to their degree (i.e. in and out-going links). Colors correspond to different bacterial phyla and dash and solid lines represent positive and negative interactions, respectively. Nodes marked with SC correspond to taxa that assigned to sponge-specific clusters. See FigureS10for the corresponding link probabilities between OTU i and j, and FigureS9for the networks belonging to the core microbiomes of HMA hosts C. reniformis and P. ficiformis. See Methods for how the networks were constructed.

  Theory shows that reciprocal interactions such as exploitation {+/ }, cooperation {+/+} and competition { / } differ in their effects on community stability and ecosystem functions ([START_REF] Mougi | Diversity of Interaction Types and Ecological Community Stability[END_REF][START_REF] Mougi | Stability of competition-antagonism-mutualism hybrid community and the role of community network structure[END_REF][START_REF] Mougi | Adaptation in a hybrid world with multiple interaction types: a new mechanism for species coexistence[END_REF][START_REF] Lurgi | The effects of space and diversity of interaction types on the stability of complex ecological networks[END_REF][START_REF] Alix | Structure-stability relationships in networks combining mutualistic and antagonistic interactions[END_REF][START_REF] Coyte | The ecology of the microbiome: Networks, competition, and stability[END_REF][START_REF] Mougi | Food-web complexity, meta-community complexity and community stability[END_REF]). Specifically, communities consisting of a mixture of unilateral interactions are more stable than those with only reciprocal inter-

actions (

[START_REF] Mougi | Food-web complexity, meta-community complexity and community stability[END_REF]

). We found that amongst all possible interspecific interactions (Figures S12), unilateral interactions in the form of commensalism {+, 0} and amensalism { , 0} dominated (Figure

6 B

). Reciprocal interactions-cooperative {+, +}, competitive { , } and exploitative {+, } interactions were exceptionally rare. The core microbiomes of hosts

  time series of high-throughput DNA sequence counts y i,t were modeled as Poisson random variables with means

	y i,t ⇠ Pois( i,t )	(10)
	log t,i = n i,t + log N t + µN t	(11)
	t,i satisfying the following log-linear model
	16

Table 1 :

 1 Microbiome diversity by temporal assemblage ± 66.7 141.4 ± 30.2 126.8 ± 78.1 126.3 ± 42.8 160.1 ± 60.4

			HMA			LMA
		A. oroides C. reniformis P. ficiformis A. damicornis	D. avara	C. crambe
	Core	45 38.4 ± 11.7	33 27.4 ± 11.1	40 32.1 ± 12.3	8 6.7 ± 2.1	6 6.4 ± 2	8 6.4 ± 2.3
	Transient	90 41.3 ± 14.4	54 26.4 ± 12.5	140 67.9 ± 30.8	31 12.7 ± 5.9	44 21.4 ± 7.4	41 17.7 ± 6.2
	2658 140 ± 77.7 129.9 Total biome 2436 Opportunistic 2793 2523	2580 2760	2443 2482	2763 2813	3465 3514

Table 2 :

 2 Percentage of taxa within each assemblage and host assigning to sponge-specific clusters. Dark gray denotes transient and light gray opportunistic assemblages, respectively.

			HMA			LMA
		A. oroides C. reniformis P. ficiformis A. damicornis D. avara C. crambe
	Core	42.2	45.6	60	25	0	12.5
	Transient	43.3	40.7	47.1	25.8	9.1	9.8
	Opportunistic	22.8	33.9	35.6	22.2	9.5	11.8

Figure

1

: Microbiome compositional similarity across host species. Non-metric multidimensional scaling (NMDS) calculated on Jaccard distances among all 36 monthly samples for each host species (n=36x6=216). Colors and shapes denote all monthly samples from a given host species surrounded by an ellipse showing the 95% confidence interval. Red circles A. oroides; blue triangles C. reniformis; pink diamonds P. ficiformis; green stars A. damicornis; yellow crosses (+) D. avara; and orange crosses (⇥) C. crambe. ANOSIM: R=0.Figure

2

: Microbiome compositional similarity across assemblages and host species. Non-metric multidimensional scaling (NMDS) calculated on Jaccard distances among all 36 monthly samples for each host and assemblage (n=36x6x3=648). Colors and shapes denote all monthly samples from a given host and assemblage, respectively. Different shapes for each host species: circles A. oroides; triangles C. reniformis; diamonds P. ficiformis; stars A. damicornis; crosses (+) D. avara; crosses (⇥) C. crambe. Different colors for each core microbiome: red A. oroides; blue C. reniformis; pink P. ficiformis; green A. damicornis; yellow D. avara; and orange: C. crambe.

Acknowledgements

JRB. was supported by an FPI Fellowship from the Spanish Government (BES-2011-049043), JMM. was supported by the French Laboratory of Excellence Project 'TULIP' (ANR-10-LABX-41; ANR-11-IDEX-002-02) and by a Region

Code and Data Availability

All data and code will eventually be available at the Open Science Framework.

Author contributions

JRB and JMM conceived the study. JRB performed all the analyses and wrote the manuscript. JRB and RBO adapted the MAR(1) model. JRB and JMM refined the manuscript. RC and MR collected the data. All authors commented on later versions of the manuscript. Figure 6: Processes explaining core microbiome dynamics and the frequency of of all possible interaction types within HMA core microbiomes. Panel A shows the relative contribution of interspecific (yellow) and intraspecific (blue) interactions and environmental variability (greens) to temporal variation in microbial population abundances across core microbiomes. Environmental variability is decomposed into residual variation (light green, Env) and variation attributed to the modeled environmentalcovariates (dark green, Cov). In all hosts, core microbiome dynamics were mainly driven by intraspecific interactions. While the modeled environmental covariates explained relatively little variation across core microibomes, an important driver of the dynamics in HMA cores was interspecific interactions which were almost negligible in LMA cores. Panel B shows the relative frequency of all inferred possible interaction types within HMA cores. This was calculated from the marginal sign posterior distribution of alpha i,j (see Methods).

Commensalism {+, 0} and amensalism { , 0} were the most frequent interaction types across HMA cores. Competitive { , } and exploitative {+, } interactions were exceptionally rare. Noteworthy, cooperative interactions {+, +} was never inferred.