Konstantin Khalturin

Isabelle M L Billas

Yassmine Chebaro

Adam M Reitzel

Ann M Tarrant

Vincent Laudet

Gabriel V Markov
email: gabriel.markov@sb-roscoff.fr

NR3E receptors in cnidarians : a new family of steroid receptor relatives extends the possible mechanisms for ligand binding

Keywords: A-ring aromatic steroid, aromatization, steroid receptor, cnidarian

Cnidarians do not have vertebrate-type estrogen signaling.

• A cytoplasmic factor triggers nuclear translocation of hERα-GFP in Hydra epithelium.

• Medusozoan cnidarians have a specific type of steroid-related receptor, called NR3E.

• Anthozoan cnidarians have lost the NR3E receptor but can produce aromatic steroids.

• NR3E has a cnidarian-specific anchor that can bind to aromatic steroids in a novel way.

" /" 1 29

Introduction

Nuclear receptors are important regulators of life-history transitions in various phyla of bilaterian animals. They are well documented ligand-activated transcription factors in vertebrates, arthropods and nematodes. However, their role is largely unknown in other metazoan groups such as molluscs or annelids [START_REF] Handberg-Thorsager | The ancestral retinoic acid receptor was a low-affinity sensor triggering neuronal differentiation[END_REF]. In cnidarians, following initial identification through PCR screens [START_REF] Escriva | Ligand binding was acquired during evolution of nuclear receptors[END_REF][START_REF] Grasso | The evolution of nuclear receptors: evidence from the coral Acropora[END_REF], some nuclear receptors have been implicated as regulators of nervous system development [START_REF] Kostrouch | Retinoic acid X receptor in the diploblast, Tripedalia cystophora[END_REF][START_REF] Gauchat | The orphan COUP-TF nuclear receptors are markers for neurogenesis from cnidarians to vertebrates[END_REF] and one " /" of them was recently shown to be involved in the molecular cascade triggering the polyp-tojellyfish transition in the scyphozoan medusa Aurelia aurita [START_REF] Fuchs | Regulation of polyp-to-jellyfish transition in Aurelia aurita[END_REF]. However, the endogenous ligands, if any, for these cnidarian nuclear receptors are still unclear.

Strictly speaking, cnidarians do not have a bilaterian-like endocrine system, with signaling molecules transported by means of an internally circulating body fluid [START_REF] Tarrant | Endocrine-like Signaling in Cnidarians: Current Understanding and Implications for Ecophysiology[END_REF]. However, they do regulate their physiology in response to nutritional and environmental inputs, which also makes them sensitive to environmental pollution, and diffusible signaling molecules are involved in such processes [START_REF] Tarrant | Hormonal signaling in cnidarians: do we understand the pathways well enough to know whether they are being disrupted?[END_REF][START_REF] Bosch | How do environmental factors influence life cycles and development? An experimental framework for early-diverging metazoans[END_REF].

Historically, numerous efforts have been made to search for vertebrate-like signaling features in cnidarians, especially components of the estrogen signaling pathway, which, in vertebrates, is primarily mediated through the estrogen nuclear receptor belonging to the NR3 family. In particular, aromatization, the final step in the synthesis of estrogens from testosterone (Figure 1A), has been studied in multiple species. Early attempts to detect aromatase activity in cnidarian tissues (scleractinian corals and octocorals), using vertebrate precursors, were unsuccessful [START_REF] Gassman | Cytochrome P-450 content and xenobiotic metabolizing enzyme activities in the scleractinian coral, Favia fragum (Esper)[END_REF][START_REF] Slattery | Steroid metabolism in Antarctic soft corals[END_REF][START_REF] Tarrant | Metabolism of estrogens and androgens by scleractinian corals[END_REF], but subsequently, aromatization activity has been reported in anemones and scleractinian corals [START_REF] Twan | Sex steroids in scleractinian coral, Euphyllia ancora: implication in mass spawning[END_REF][START_REF] Rougee | Molecular reproductive characteristics of the reef coral Pocillopora damicornis[END_REF][START_REF] Goodall | Characterization of reproductive hormonal pathways in cnidarians[END_REF].

Aromatization of exogenous androgen precursors by tissue homogenates from cnidarians or molluscs has been interpreted as evidence for the ability to endogenously synthesize vertebrate-type estrogens [START_REF] Armoza-Zvuloni | Trioecy, a unique breeding strategy in the sea anemone Aiptasia diaphana and its association with sex steroids[END_REF]. However, analyses of genomic data indicate that cnidarians cannot synthesize estrogens through the same pathways as vertebrates because essential components are lacking in cnidarians. Specifically the CYP19 aromatase is chordate-specific, and the CYP11A1 enzyme responsible for the early cholesterol side-chain cleavage is vertebrate-specific [START_REF] Markov | Independent elaboration of steroid hormone signaling pathways in metazoans[END_REF]. Moreover, there is evidence that corals can take up vertebrate steroids from the water column, especially around human polluted sites, so estrogen-like compounds detected in cnidarian tissues may have exogenous origins [START_REF] Tarrant | Uptake of estrone from the water column by a coral community[END_REF][START_REF] Armoza-Zvuloni | Reproductive characteristics and steroid levels in the scleractinian coral Oculina patagonica inhabiting contaminated sites along the Israeli Mediterranean coast[END_REF]. These observations are also in line with the reinterpretation of the presence of vertebrate-type steroids in molluscs [START_REF] Scott | Do mollusks use vertebrate sex steroids as reproductive hormones? Part I: Critical appraisal of the evidence for the presence, biosynthesis and uptake of steroids[END_REF][START_REF] Schwarz | Mussels (Mytilus spp.) display an ability for rapid and high capacity uptake of the vertebrate steroid, estradiol-17beta from water[END_REF]. Arthropod-like ecdysteroids and nematode-like dafachronic acids have also been isolated from corals [START_REF] Lafont | Diversity of Ecdysteroids in Animal Species in Ecdysone: Structures and Functions[END_REF][START_REF] Saini | [END_REF], but have not been viewed as evidence for ecdysozoan-type hormonal signaling in these cnidarians. In case of ecdysteroids, a defensive role and a dietary origin are the main interpretations that have been used to explain the presence and high abundance of those compounds [START_REF] Lafont | Diversity of Ecdysteroids in Animal Species in Ecdysone: Structures and Functions[END_REF]. Similarly, the presence and activity of " /" gonadotropin-like neuropeptides in cnidarians have been used to draw parallels between vertebrate and cnidarian reproductive signaling cascades [START_REF] Twan | Corals have already evolved the vertebrate-type hormone system in the sexual reproduction[END_REF]. However, here again, phylogenetic analyses indicate that the cnidarian peptides are equally similar to arthropod neuropeptides that do not act in reproductive signaling [START_REF] Tsai | Gonadotropin-releasing hormone by any other name would smell as sweet[END_REF].

To date, ten distinct A-ring aromatized steroids have been identified from four different octocoral species belonging to different genera (Alcyonium gracillimum, Capnella sp., Dendronephtya studeri, Scleronephthya pallida). Among these steroids are aromatized pregnanes and aromatized C9-C10 secosteroids, which are compounds that have core structures similar to progesterone and vitamin D, respectively [START_REF] Sarma | Marine metabolites: the sterols of soft coral[END_REF]. A third group consists of A-ring aromatized steroids bearing a variety of side chains. These have been found in the Taiwanese octocoral Dendronephthya studeri [START_REF] Yan | Steroids with aromatic A-rings from the Hainan soft coral Dendronephthya studeri Ridley[END_REF], see also Figure 1. Based on the co-occurrence of A-ring aromatic steroids and corresponding dienones (see compounds ( 1) and (3) on Figure 1), that bear an additional unsaturation (in yellow on Figure 1) compared to testosterone, we hypothesize that the aromatization reaction in those corals is distinct from that in vertebrates, as proposed for the two other groups [START_REF] Sarma | Marine metabolites: the sterols of soft coral[END_REF]. norcholesta-1,3,5(10)-trien-3-ol) has been proposed as an ancestral steroid [START_REF] Markov | Origin of an ancient hormone/receptor couple revealed by resurrection of an ancestral estrogen[END_REF]. We have named « dienone aromatization » the proposed aromatization reaction to stress the difference from the vertebrate case where there is no delta1-2 double bond (highlighted in yellow) on the A ring.

" /"
Among the steroids present in Dendronephthya studeri, one compound has an identical structure to what we hypothesize to be the ancestral NR3 ligand in chordates. This compound indeed binds an ancestral steroid receptor with low affinity [START_REF] Markov | Origin of an ancient hormone/receptor couple revealed by resurrection of an ancestral estrogen[END_REF]. Due to the combination of data discussed above, we have previously interpreted paraestrol A (19-norcholesta-1,3,5(10)-trien-3-ol) and other aromatized steroids from corals as defensive compounds [START_REF] Markov | Origin of an ancient hormone/receptor couple revealed by resurrection of an ancestral estrogen[END_REF]. The diversity of steroids in cnidarians led us to reevaluate the distribution and potential signaling functions of NR3 in the phylum Cnidaria. Our data suggest that NR3 orthologs are present in three of the four major cnidarian lineages and that an endogenous ligand may be present based on heterologous expression of vertebrate ERα in Hydra. Together, our data suggest that some ring-A aromatic steroids might function as endogenous signaling molecules in cnidarians, particularly paraestrol A that has the potential to act as a ligand for cnidarian NR3 receptors.

" /"

Material and Methods

Sequence analysis

Cnidarian NR3 sequences were identified by tBLASTn and reciprocal BLASTp queries at NCBI GenBank using human ERα and Trichoplax adhaerens NR3 (previously identified [START_REF] Baker | Trichoplax, the simplest known animal, contains an estrogen-related receptor but no estrogen receptor: Implications for estrogen receptor evolution[END_REF]). In addition to NCBI data, putative NHRs were identified by tBLASTn searches of the de novo assembled transcriptomes of Aurelia, Morbakka, Chironex, Tripedalia and Copula (Khalturin, unpublished data). We required that identified sequences had to be complete or nearly complete predicted nuclear receptor proteins (i.e., both DNA and ligand binding domains) to facilitate phylogenetic analyses. All sequences were downloaded and, if necessary, translated into the correct reading frame. Collected sequences were aligned using Clustal Omega [START_REF] Sievers | Clustal Omega, accurate alignment of very large numbers of sequences[END_REF] and alignments were checked by eye and edited with Seaview [START_REF] Gouy | SeaView version 4: A multiplatform graphical user interface for sequence alignment and phylogenetic tree building[END_REF]. A phylogenetic tree (Figure 2) was made using PHYML [START_REF] Guindon | A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood[END_REF], a fast and accurate maximum likelihood heuristic method, using the LG model [START_REF] Le | Accounting for solvent accessibility and secondary structure in protein phylogenetics is clearly beneficial[END_REF] with a gamma law. Reliability of nodes was assessed by likelihood-ratio test [START_REF] Anisimova | Approximate likelihood-ratio test for branches: A fast, accurate, and powerful alternative[END_REF].

Heterologous hERα-GFP expression in Hydra

Nucleotide sequence of human estrogen receptor (hERα) was adjusted to fit the codon usage of Hydra vulgaris with the average GC content reduced to 33%. A synthetic gene (shERα) produced by GeneArt (Regensburg, Germany) was inserted into Hydra expression vector (ligAC-6) between the actin promoter and EFGP coding region using the SbfI and PacI restriction sites (see Figure 3a).

The resulting plasmid ligAL-5 was propagated in E. coli and sequenced. DNA for microinjection was prepared using the Plasmid Midi Kit (QIAGEN). Construct DNA was microinjected into fertilized eggs of Hydra vulgaris strain AEP as described previously [START_REF] Wittlieb | Transgenic Hydra allow in vivo tracking of individual stem cells during morphogenesis[END_REF]. In about 70% of injected embryos EGFP-positive nuclei were detectable 48 hours after microinjection. Polyps started to hatch 2-3 weeks after microinjection and contained small patches of ectodermal epithelial cells with clear nuclear localisation of the shERα-EFGP protein. Polyps were cultivated under intensive feeding conditions and transgenic patches were enriched by selection of EGFP-positive buds and removal of non-transgenic areas. Within 3 months of selection a polyp culture where ~80% of ectodermal cells were transgenic was obtained (see Figure 3). " /"

Homology modeling of NR3E and paraestrol A docking

Homology modeling was performed using the program Modeller [START_REF] Šali | Comparative protein modelling by satisfaction of spatial restraints[END_REF] (version 9.18). The selected structural templates for the homology modeling were the crystal structure of ERα bound to 17βestradiol (PDB id: 1ERE) [START_REF] Brzozowski | Molecular basis of agonism and antagonism in the oestrogen receptor[END_REF], ERRγ in the constitutively active conformation (PDB id: 1KV6) [START_REF] Greschik | Structural and functional evidence for ligand-independent transcriptional activation by the estrogenrelated receptor 3[END_REF] and apo ERRα in its active conformation (PDB id: 3D24) [START_REF] Greschik | Communication between the ERRalpha homodimer interface and the PGC-1alpha binding surface via the helix 8-9 loop[END_REF]. Missing loops in the ERα structure, namely between residues 331-336 and 462-464, were built prior to homology modeling [START_REF] Fiser | Modeling of loops in protein structures[END_REF]. Sequence alignment of target and template was performed using the align2d procedure in Modeller, taking into account structural information from the template when constructing the alignment, resulting in 26.2, 29.5 and 25.8% of sequence identity for the modeling using ERα, ERRγ and ERRα respectively. A total of 50 models were generated and the structure with the best DOPE score [START_REF] Shen | Statistical potential for assessment and prediction of protein structures[END_REF], a statistical potential optimized for model assessment in Modeller, was selected for ligand docking. Examination of the best homology models based on ERα, ERRγ and ERRα revealed that the C-alpha traces nicely follow their respective crystal structure templates with an RMSD of 0.172, 0.289 and 0.293Å respectively after sequence-structure superposition using Pymol software and with 96.48, 95.87 and 96.36% of residues estimated in the Ramachandran favored space using Molprobity server [START_REF] Chen | MolProbity: all-atom structure validation for macromolecular crystallography[END_REF]. We noticed that homology modelling of NR3E LBD based on the crystal structure of ERRγ, led to a slight local deformation in a helix turn belonging to helix H10. This is to maintain the structure/sequence conservation in the rest of the model, but this does not affect the ligand-binding pocket and thus does not impact the docking results. Docking was performed with Autodock4 [START_REF] Morris | Autodock4 and AutoDockTools4: automated docking with selective receptor flexiblity[END_REF] using the Lamarckian genetic algorithm and consisted of 10 runs per search, with a maximum of 2.5 million energy evaluation per run and a population size of 150. Affinity maps for the receptor were computed using Autogrid4 prior to docking, using the default grid map spacing of 0.375Å. Paraestrol A was superimposed on the 17β-estradiol ligand in the 1ERE structure, and the homology models on the ERα protein in the same structure. Based on this, residues in the models displaying steric clashes were selected as flexible residues in the docking procedure. In order to compare the estimated free energies of binding, a docking run was performed using the 17β-estradiol ligand in the ERα crystal structure with the same docking " /" parameters. Interaction diagram was drawn using LIGPLOT v.4.5.3 [START_REF] Wallace | LIGPLOT: a program to generate schematic diagrams of protein-ligand interactions[END_REF]. All structure figures were prepared using PyMOL (The PyMOL Molecular Graphics System, Version 1.7.4.0 , Schrödinger LLC).

Results

NR3E is a novel nuclear receptor type present in three of the four main cnidarian lineages

We identified new nuclear receptor NR3 sequences from cnidarians that group together as a separate family at the base of the bilaterian NR3s, suggesting they are not particularly related to one of the bilaterian NR3 families (i.e., ER, ERR, or SR). This is also true for the Trichoplax sequence initially described as an ERR [START_REF] Baker | Trichoplax, the simplest known animal, contains an estrogen-related receptor but no estrogen receptor: Implications for estrogen receptor evolution[END_REF]. We have previously suggested, based on the uncertainty in the phylogenetic relationships between vertebrate ER and sequences from molluscs and annelids described as ERs to rename the protostome clade as NR3D, to avoid misleading inferences that estrogen would be their natural physiological ligand [START_REF] Markov | Origin and evolution of the ligand-binding ability of nuclear receptors[END_REF]. The most recent analyses of bilaterian NR3s, based on extensive sampling, are in agreement with this view, putting the molluscan and annelid NR3Ds at the base of a clade grouping the vertebrate NR3A and NR3C [START_REF] Bridgham | Vestigialization of an allosteric switch: genetic and structural mechanisms for the evolution of constitutive activity in a steroid hormone receptor[END_REF][START_REF] Katsu | A second estrogen receptor from Japanese lamprey (Lethenteron japonicum) does not have activities for estrogen binding and transcription[END_REF]. Here, NR3A and NR3D group together in a weakly supported node, so we do not think this grouping contradicts these previous studies. Most importantly, the cnidarian sequences follow the known subdivisions between the three clades in which they are present : cubozoans, scyphozoans and hydrozoans (Figure 3). Therefore we propose to name this clade as NR3E. For the same reason, we propose that the Trichoplax ERR should be renamed NR3F. The tree was calculated using maximum-likelihood. Likelihood ratio test branch support values upper than 0.70 are indicated. Bilaterian sequences are shown in blue and cnidarian sequences in brown. The Hydra NR3E is highlighted in red. Mollusc and annelid NR3Ds correspond to proteins that were previously labelled as ERs [START_REF] Bridgham | Vestigialization of an allosteric switch: genetic and structural mechanisms for the evolution of constitutive activity in a steroid hormone receptor[END_REF]. Similar analyses using the LBD alone gave largely similar results and did not affect the placement of the cnidarian nuclear receptors.

" /"

Heterologous Human ERα-GFP construct translocates into the Hydra cytoplasm

In the context of functional studies based on transgenesis in Hydra [START_REF] Wittlieb | Transgenic Hydra allow in vivo tracking of individual stem cells during morphogenesis[END_REF], we made an inducible reporter system for Hydra based on human ERα. This was based on the assumption that the vertebrate estrogen signaling system does not naturally function within Hydra. The approach was to make transgenic animals where the translocation of transcription factors into the nucleus is " /" controlled by addition of estrogen (in the construct hERα would be fused to the gene of interest).

Initially, a construct was made containing hERα fused to green fluorescent protein (hERα-GFP), which theoretically should stay in the cytoplasm in the absence of exogenous estrogen. Surprisingly, hERα-GFP was translocated to the nucleus in Hydra even without addition of estrogen (Figure 3).

At the same time pure GFP protein never concentrates in the nuclei of transgenic Hydra. There might be several explanations for the translocation of the hERα-GFP construct. For example, hsp90 and hsp70 of Hydra may not bind hERα-GFP properly and thus ERα would go into the nucleus by default. However, the high level of conservation between human and cnidarian Hsp70/90 speaks against "default" translocation. Another plausible scenario is the presence of a cytoplasmic factor in Hydra which can specifically bind the ligand-binding pocket of the human estrogen receptor and thus trigger nuclear translocation.

" /"

Docking of paraestrol A inside the ligand-binding pocket of ERα-, ERRγ-and ERRα-based homology models of NR3E

Docking experiments were performed to assess the potential binding of paraestrol A to Hydra NR3E. This was done using homology models for Hydra NR3E that were developed using the crystal structures of ERα, ERRγ and ERRα as templates. For comparison, we employed the same parameters (see Materials and Methods) for docking the 17β-estradiol in the ERα crystal structure, estimated the corresponding binding free energy and compared it to that of paraestrol A in the different homology models. In the case of 17β-estradiol in ERα, the conformation obtained for 17βestradiol was nearly identical to the one observed in the crystal structure (rmsd of 0.7Å) and the calculated binding free energy was estimated to be -10.32 kcal/mol.

Considering the Hydra NR3E ERα-based homology model, superimposition of paraestrol A onto the steroid core of 17β-estradiol in the ERα-based homology model shows that paraestrol A fits inside the ligand-binding pocket without any steric clash with the surrounding residues. As a consequence, the initial conformations of the latter residues were maintained throughout the whole docking procedure. The best docking pose obtained had a binding free energy of -11.28 kcal/mol after 10 runs of docking attempts. These observations together with a comparison of the free energy with that of 17β-estradiol in ERα suggest that paraestrol A would nicely fit inside the ligand-binding pocket of the ERα-based homology model with interactions as described below.

For the ERRγ-based model, ligand superimposition of paraestrol A onto the steroid core of 17βestradiol indicates that residues W80 and F35 are in steric clash with the ligand. We therefore considered them as flexible in the docking procedure, while the conformation of the other residues inside the binding cavity was kept as modeled initially. The binding free energy of the best docking pose was estimated to be -9.38 kcal/mol. As in the case of the ERα-based homology model, the results suggest that paraestrol A could be well-accommodated inside the ligand-binding pocket of the ERRγ-based homology model of Hydra NR3E LBD with an interaction pattern rather similar to the one observed for the docking in the ERα-based homology model We further applied the same procedure to the ERRα-based homology model and observed that several residues, especially aromatic ones, were in steric clash with the paraestrol A initial conformation. Despite the fact that the residues in steric clash were allowed to be flexible in the " /" docking procedure, no reasonable docking pose was obtained, suggesting that successful docking of paraestrol A strongly depends on the nature of the homology model used and further supports the validity of the docking predictions obtained in the ERα-and ERRγ-based homology models.

Interactions of paraestrol A with residues inside the ligand-binding pocket

In the ERα-based homology model, paraestrol A is globally positioned in a similar manner as is 17β-estradiol in its cognate receptor (Figure 4a). In fact, the aliphatic side chain of paraestrol A is oriented towards the helices H7 and H11, while the steroid core is positioned between helices H3 and H5 and rather close to the first strand of the β-sheet (Figure 4a). Several hydrophobic residues, including aromatic ones form stabilizing interactions with the ligand. In addition, a remarkable π-π interaction is observed between Trp80 in helix H5 and the steroid core (Supplemental Figure 1); this residue and the corresponding stabilizing interaction are specific to NR3E (Figure 5). We further observe that the 3-OH group of paraestrol A is hydrogen-bonded to Asn79 (H5) with a donor-acceptor distance of 2.9Å, and examination of the surroundings of this interaction shows a potential additional hydrogen bond with Gln42 (H3) (Supplemental Figure 1). Both residues could form a clamp that would strongly position the ligand inside the pocket. To test this hypothesis, we started with the paraestrol A docked conformation in the ERα-based model and performed additional docking runs where both Gln42 and Asn79 were allowed to be flexible. Interestingly, the best docked conformation in this case showed structural adaptation of these two residues, which are then both implicated in hydrogen bond formation with the 3-OH group of paraestrol A (Figure 4b), with estimated binding free energy of -11.89 kcal/mol, and a similar π-π interaction with Trp80 in H5 as in the previous docking result. The helices from the ligand-binding domain are mapped according to the structure of human ERα [START_REF] Tanenbaum | Crystallographic comparison of the estrogen and progesterone receptor's ligand binding domains[END_REF]. Calculations of identity percentages for both domains are also provided in Supplemental Table 1.

" /"

" /"

" /"

The relevance of these observed interactions is emphasized by the remarkable conservation of residues W80, Q42 and N79 in cnidarian NR3 sequences (Figure 5). In fact, W80 is strictly conserved and specific to the cnidarian NR3E, while being replaced by M/L/V/I/C/A in bilaterian steroid receptors and ERRs as well as in Trichoplax NR3F. Furthermore, the anchoring residues Q42 and N79 are strictly conserved in all cnidarian NR3E sequences. Q42 is also present in human oxosteroid receptors, but replaced by E in human ERs and ERRs, as well as in Trichoplax NR3F.

N79 is also lineage-specific, being strictly conserved in cnidarian NR3E, and replaced by G/D/A/Q/ N in bilaterian steroid receptors and in Trichoplax NR3F. In human ERα bound to 17β-estradiol, the 3-OH group is stabilized by an electrostatic anchor formed by two conserved residues E353 (H3) and R394 (H5) (Figure 4c). These residues are crucial for the positioning and the stabilization of the natural 17β-estradiol ligand inside its cognate ligand-binding pocket. In the case of cnidarian NR3E, the anchor is made of Q42 and N79. In this case, the 3-OH group of paraestrol A acts both as an acceptor and a donor for these residues. An inspection of the multiple sequence alignment together with the superimposition of the homology model with the crystal structure indicate that Q42 in H3 is the cnidarian equivalent of human E353 that forms the first half of the electrostatic anchor. On the other hand, the second half of the anchor differs in the two cases. Residue N79 of Hydra NR3E is located one helix turn upstream of the position occupied by R394 and replaced by G in human ERα (Figure 4b and4c). Altogether, we observed that the 3-OH functional group of paraestrol A is in an adequate environment to form complementary interactions with the conserved electrostatic Q(H3)/N(H5) anchor and further stabilized by π-π interaction with the cnidarianspecific W80 residue located in H5.

It is noteworthy that in our second set of docking experiments, the two best docked structures display hydrogen-bonding with Q42 and N79, however in the third one, the distance between 3-OH group of paraestrol A and the amide group of Q42 increases to 4.8Å, leading to a severe displacement of the ligand in the pocket, with the aliphatic side chain replacing the position of the aromatized A-ring. This structure supports the crucial role of the anchor in positioning of the ligand in the pocket, as strongly suggested by the sequence conservation of these residues.

The positioning of paraestrol A inside the ligand-binding pocket of the ERRγ-based homology model is similar as is its general interaction pattern with surrounding residues in the best docked " /" conformation with an estimated binding free energy of -9.38 kcal/mol. Even though the exact position of paraestrol A is slightly more towards the outside of the pocket and slightly shifted between helices H3 and H5, W80 still can form a π-π interaction and a hydrogen bond with Q(H3) (Supplemental Figures 3a and3b). We then performed a similar docking experiment where this latter structure was used and both residues Q42 and N79 are allowed to be flexible in the search.

Similarly to the previous results, adjustment of the side chains of these residues permits in the best model hydrogen bonding with the 3-OH group of the paraestrol, leading to a more favorable binding free energy of -11.50 kcal/mol (Supplemental Figure 3b). The three conserved residues that were shown to be crucial for the interaction of the ligand with the ERα-based homology model hence fulfill the same role in the ERRγ-based homology model. Altogether, our docking studies of paraestrol A in NR3E homology models suggest that paraestrol A can be readily accommodated in the ERα-based model and depicts a binding free energy that is of the same order of magnitude as the one calculated for 17β-estradiol in ERα. Paraestrol A can also fit inside the ligand-binding pocket of ERRγ, with comparable energy. In both cases, we have highlighted cnidarian-specific residues that interact with the ligand, among which two polar residues, Q(H3)/N(H5), that define a novel anchor that stabilizes paraestrol A inside the pocket.

Discussion

Novel cnidarian-specific candidate players for a steroid signaling pathway

All evidence to date has shown that the vertebrate-type estrogen-ER signaling pathway is not present in cnidarians [START_REF] Markov | Origin of an ancient hormone/receptor couple revealed by resurrection of an ancestral estrogen[END_REF]. However, recent findings, in the literature and the data presented in this study indicate the presence both of aromatic steroids and nuclear receptors within a novel subclade of the NR3 family (NR3E). Previous surveys of nuclear receptor diversity in cnidarian have largely studied anthozoans due to availability of sequence data as well as interest in possible effects of environmental pollution on reef-building corals. Our analysis here used the availability of sequence " /" data from throughout the phylum and revealed NR3 members in the other three cnidarian classes that form the medusozoan clade. The presence of these genes helps to better resolve the history of the NR3 family in animals because NR3 was present in the placozoan Trichoplax, absent from surveyed ctenophores and sponges, and until this study, inferred to absent from cnidarians. Our hypothesis is thus that NR3 evolved in the common ancestor of the placozoan and cnidarianbilaterian with lineage-specific loss in the anthozoans. Despite the presence of NR3 genes in cnidarians, it has not yet been demonstrated whether these aromatic steroids or any other compounds can serve as ligands for NR3E receptors.

Interestingly, aromatic steroids have only been documented within a few species of anthozoans but NR3E genes are only present within the medusozoa (Hydrozoa, Cubozoa, Scyphozoa). Similarly, NR2B (RXR) homologs have been lost from the Anthozoa and are only present within Medusozoa.

The loss of these two receptors from the anthozoan lineage is of considerable interest and may be related to the loss of the medusa (jellyfish) stage. This step implied the loss of metamorphosis between polyp and medusa, and/or the reduction of sensory organs (eyes). Regarding aromatic steroids in anthozoans, data about their quantitative variation under physiological conditions will be necessary to determine if they play a defensive function, or if they may be involved in endogenous intercellular signaling. If they do play any role in endogenous signaling, it could be mediated through binding to different nuclear receptors from other families than NR3, or to membrane receptors [START_REF] Schwartz | Rapid steroid hormone actions via membrane receptors[END_REF]. The translocation data on the shERα-GFP heterologous construct expressed in Hydra epithelial cells suggests that there may be a cytoplasmic factor that binds to shERα-GFP and triggers its translocation into the nucleus. Because Hydra, an hydrozoan, as other medusozoans, has an endogenous NR3E, which is the closest relative to ER and other vertebrate steroid receptors, we thought it would be important to explore its ability to bind an aromatic steroid which is present in another cnidarian, even from a distantly related species.

Cnidarian-specific molecular mechanisms for ligand binding in the NR3E pocket

Docking experiments with paraestrol A inside the ligand-binding pocket of Hydra NR3E indicate that this receptor might indeed bind an aromatic steroid. In this case, the binding would involve π -π interactions with a tryptophan residue, which is conserved in all cnidarian sequences, and a glutamine residue which, to date, was mainly correlated to the binding of oxosteroids, like " /" progesterone, androgens and corticoids. Thus, even the most estrogen-like aromatic steroid from cnidarians shows binding properties that are distinct from the mechanism of estradiol binding to the vertebrate estrogen receptor. This strongly suggests that identification of endogenous ligands for the cnidarian NR3E receptors and their concomitant functional characterization will enlarge our knowledge of the diversity ligand-binding interactions in play between steroid hormones and their cognate nuclear receptors. Over the last few years, a great deal of information has accumulated about physiologically relevant alternative ligands for the human estrogen receptors [START_REF] Baker | The promiscuous estrogen receptor : evolution of physiological estrogens and response to phytochemicals and endocrine disruptors[END_REF]. In particular, steroids without an aromatic A-ring were also reported to be able to bind ER. So far, all those described alternative steroid ligands have a hydroxy group on carbon 3, just like the hydroxy group on carbon 3 of estradiol. Interestingly, the distinction between binding of steroids with a 3hydroxy group, like estrogens, and binding of steroid with a 3-oxo group, like androgens, progesterone and corticoids, has long been correlated to the distinction between the E and Q ligandbinding residues in the H3 part of the ligand-binding pocket, and interpreted as an evidence for more ancient origin of estrogen as a ligand relative to oxosteroids [START_REF] Thornton | Evolution of vertebrate steroid receptors from an ancestral estrogen receptor by ligand exploitation and serial genome expansions[END_REF]. Docking experiments with NR3E and paraestrol A indicate that, structurally speaking, nothing prevents a 3-hydroxysteroid from binding the pocket through a Q residue, which is classically related to 3-oxosteroid binding.

Therefore, the distinction between the two classes of molecules may not be so clear-cut, and even in vertebrates, we could find unexpected new selective estrogen receptor modulators by looking at endogenous 3-oxosteroids. Many of them are indeed present in the human body as precursors in bile acid synthesis [START_REF] Russell | The Enzymes, Regulation, and Genetics of Bile Acid Synthesis[END_REF].

Long distance chemical communication in the cnidarian body

Because they lack an internal circulating body fluid, cnidarians cannot have a vertebrate-like hormonal system. This does not necessarily preclude steroid-mediated intercellular signaling in those animals. Indeed, even in vertebrates, some steroids do not act through the canonical hormonal pathway. Classically, hormones are defined as internally circulating molecules, which excludes vertebrate bile acids based on this definition [START_REF] Norris | Vertebrate endocrinology[END_REF]. However, this view is changing with the acknowledgment of their functional role not only as facilitators of lipid digestion, but also as signaling molecules, that allow the coupling of the nutritional state with various parameters, from digestive physiology to behavioural traits such as the regulation of appetite [START_REF] Thomas | Bile acids and the membrane bile acid receptor TGR5connecting nutrition and metabolism[END_REF]. This example " /"

shows that there is no objective reason to strictly limit long-distance intercellular communication to the internal milieu. We thus hypothesize that in cnidarians, the digestive surfaces may be a major vehicle for signaling molecules, and that the steroids identified in these animals, if they really have a role in distant intercellular communication, could be transported in this way. The gut lumen is sometimes viewed as the continuation of the external milieu, because both are communicating through the oral openings. Thus, secretions of molecules in the gut are considered as exocrine secretions. The real chemical composition of some parts of the digestive tract is distinct from the external environment, including secretions and selective absorption processes from the digestive surfaces of the animal. The gut goes throughout the body for some cnidarians and is restricted to particular regions in others (e.g. mesenteries) and could be a carrier of signaling molecules between distant body parts or different ramets in colonial species (Figure 6). Anyways, caution is needed before extrapolating anatomical knowledge from bilaterians to cnidarians, because the homology between the germ layers is debated [START_REF] Steinmetz | Gut-like ectodermal tissue in a sea anemone challenges germ layer homology[END_REF].

" /" The pink tissue represents the mesoglea, which consists of a gelatinous matrix that contains collagen fibers and usually some cells. The mesoglea forms a hydrostatic skeleton, but does not contain a circulating body fluid or play a known role in circulation.

Regulation of reproduction and metamorphosis have already been proposed as potential roles for cnidarian nuclear hormone receptors. Regulation of energetic metabolism and growth is another prominent function of many nuclear receptors that has not yet been investigated in cnidarians. Indeed, within the NR3 family, ERR plays important and diverse roles in energetic homeostasis in both protostomes and deuterostomes [START_REF] Giguere | Transcriptional control of energy homeostasis by the estrogen-related receptors[END_REF][START_REF] Tennessen | The Drosophila estrogenrelated receptor directs a metabolic switch that supports developmental growth[END_REF]. In cnidarians, regulation of energetic metabolism is of considerable interest, but the mechanisms are still poorly understood. Many species are tightly regulating their body size based on nutrition intake. For example, many anthozoan anemone species can survive several months without food, during which they can exhibit negative growth, consume stored lipids and reduce their metabolic rates [START_REF] Hand | Asexual reproduction by transverse fission and some anomalies in the sea anemone Nematostella vectensis[END_REF][START_REF] Chia | Development and juvenile growth of the sea anemone Tealia crassicornis[END_REF][START_REF] Sebens | The regulation of asexual reproduction and indeterminate body size in the sea anemone Anthopleura[END_REF][START_REF] Fitt | Effects of starvation, and light and dark on the energy metabolism of symbiotic and aposymbiotic sea anemones, Anthopleura elegantissima[END_REF]. Studies with the hydrozoan Hydra suggest an important role of autophagy during starvation. Starving hydras do not die, but simply shrink proportionally [START_REF] Chera | Autophagy in Hydra: A response to starvation and stress in early animal evolution[END_REF]. For both anemones and Hydra, growth rapidly resumes following feeding " /" (references cited in two preceding sentences and personal observations of authors). Hence, there must be some regulators and nuclear receptors are among the best candidates.

Conclusions

While molecular studies of steroid signaling in cnidarians are still in early stages, the fragmentary data are already sufficient to state that, despite some similarities to bilaterians, cnidarians have lineage-specific receptors and potential ligands. A detailed look at the mechanisms of steroid signaling and receptor function is certain to give new insights on the many different possible ways to achieve interactions between a receptor and a steroid ligand. Obviously, this also holds true for other chemical classes of possible ligands. In comparison to anthozoans, very few chemical analyses have targeted the steroid composition of medusozoan cnidarians [START_REF] Cimino | G, Cholest-4-en-4,16β,18,22R-tetrol-3-one 16,18-diacetate, a novel polyhydroxylated steroid from the hydroid Eudendrium sp[END_REF][START_REF] Milkova | Sterols from black sea invertebrates-I. Sterols from Scyphozoa and Anthozoa (Coelenterata)[END_REF][START_REF] Rosa | Sterol composition of the Black Sea Hydrozoan, Obelia longissima (Pallas 1766)[END_REF]. There is growing awareness among the global metabolomics community that « time is ripe to focus on model organism metabolomes » [START_REF] Edison | The Time Is Right to Focus on Model Organism Metabolomes[END_REF]. Along these lines, metabolomic approaches have already proven extremely useful in understanding phenotypic plasticity in a mediterranean zoanthid coral [START_REF] Cachet | Metabolomic profiling reveals deep chemical divergence between two morphotypes of the zoanthid Parazoanthus axinellae[END_REF]. We hope that our results will stimulate interest of natural product chemists in identifying potential endogenous ligands for cnidarian nuclear receptors and particularly in determining whether A ring aromatic sterols are present in Hydra and other medusozoan cnidarians, as well as other molecules that could function as endogenous ligands to cnidarian nuclear receptors. This will not only illuminate the field of cnidarian zoology but, more widely, the field of nuclear receptor pharmacology to uncover new mechanisms of action and novel signaling molecules.

Supplemental

Figure 1 .

 1 Figure 1. A model for A-ring aromatic steroid synthesis pathway in the octocoral cnidarian Dendronephthya studeri, based on the molecules described in Yan et al., 2011, and the knowledge of enzymatic pathways known from other organisms. (1): Cholesta-1,4,22-trien-3one; (2): (22E)-19-Norcholesta-1,3,5(10),22-tetraen-3-ol; (3): 24-Methylenecholesta-1,4,22trien-3-one; (4): 24-Methylene-19-norcholesta-1,3,5(10),22-tetraen-3-ol. Paraestrol A (19-

Figure 2 .

 2 Figure 2. A new NR3E family at the base of the bilaterian ER/ERR/SR clade.

Figure 3 .

 3 Figure 3. An endogenous cytoplasmic factor induces hERα-GFP translocation from cytoplasm to nucleus in Hydra epithelial cells. a) Schematic representation of transgenic constructs with GFP only and with human ERα-GFP driven by 1.4 kb of Hydra ß-actin promoter. b) One week old Hydra embryo injected with the control construct. GFP positive blastomeres are visible. c) Hydra embryo injected with actin::ERα-GFP construct. Several GFP positive nuclei are visible and the cytoplasm of blastomeres is not GFP positive. d) Head of transgenic Hydra polyp. e) Body column of a transgenic polyp. Green spots are the nuclei of ectodermal epithelia cells which are filled with hERα-GFP without addition of any ligand.

Figure 4 .

 4 Figure 4. Paraestrol A docking into the ligand binding pocket of an estrogen receptor α homology-based model of Hydra NR3E. (a) Superposition of the NR3E structural model using ERα structure (in violet) where the paraestrol A is represented in pink, the ERα structure used as template (PDB id: 1ERE) is represented in cyan and the 17β-estradiol in blue. (b) Detailed view of the binding pocket of the NR3E ERα-based model, highlighting interactions (in yellow dashed lines) between the paraestrol and residues Q42 and N79. W80 involved in the cnidarian-specific π-π interaction is also shown. (c) Detailed view of the binding pocket of ERα, highlighting interactions (in yellow dashed lines) between the 17β-estradiol and residues E353, R394 and H524 (corresponding to V201 in NR3E model).

Figure 5 .

 5 Figure 5. Alignment of H3, H4-H5 and H10-H11 helices in the NR3 family, compared to NR3E from Hydra. Residues involved in estrogen binding in ER are highlighted in green. Homologous residues that are shared by the cnidarian NR3E and the vertebrate oxosteroid receptors (SRs) are highlighted in pink. The cnidarian-specific N79-W80 anchor is highlighted in brown. An alignment for the entire DNA-binding domain and ligand-binding domain is shown in Supplemental Figure 2.

Figure 6 .

 6 Figure 6. A possible pathway for steroid enterocrine signaling in cnidarians. Steroids, in green, could be produced from dietary sterols in the gut (in yellow) and go to the gonad (in orange) along with other nutrients. Gonads may be of ectodermal or endodermal origin, depending on the species.

Supplemental Figure 3 .

 3 Paraestrol A docking into the ligand binding pocket of estrogen related receptor γ homology-based model of Hydra NR3E. (a) Superposition of the NR3E structural model using ERRγ structure (in dark green) where the paraestrol A is represented in light green, the ERRγ structure used as template (PDB id: 3D24) is represented in red and for comparison the paraestrol A position found in the NR3E ER model in pink. (b) Detailed view of the ligand binding pocket of NR3E ERRγ-based model, highlighting interactions (in yellow dashed lines) between the paraestrol and residues Q42 and N79. Results of the first round of docking experiments are represented in grey, showing one hydrogen bond (in red dashed line) between paraestrol and residue Q42. (c) Detailed view of the binding pocket of NR3E ERα-based model, highlighting interactions (in yellow dashed lines) between the paraestrol and residues Q42 and N79. Results of the first round of docking experiments are represented in grey, showing one hydrogen bond (in red dashed line) between paraestrol and residue N79.

Table 1 . Identity percentages of DBDs and LBDs when compared to NR3E from Hydra.

 1

	10/11

Supplemental Figure 2 (continued)

Acknowledgments

The authors wish to express their sincere thanks to Jörg Wittlieb (CAU, Kiel) for microinjection of transgenic construct into Hydra AEP embryos. KK is supported by grant from Japan Society for the Promotion of Science (JSPS 17K07420). I.M.L.B and Y.C. acknowledge the support and the use of resources of the French Infrastructure for Integrated Structural Biology FRISBI ANR-10-INBS-05 and of Instruct-ERIC. AMR was supported by NIH Award R15GM114740. AMT was supported by an Internal Research and Development Award from the Woods Hole Oceanographic Institution. The work from the VL laboratory is supported by CNRS, Sorbonne Université, ANSES and EMBRC-" /"

Supplemental information for

anchor is highlighted in brown.


Alexandros Alexakis
email: alexakis@lps.ens.fr

3D instabilities and negative eddy viscosity in thin-layer flows

Keywords: numbers: 47.20.-k, 47.11.St, 47.11.Kb, 47.15.Fe

The stability of flows in layers of finite thickness H is examined against small scale three dimensional (3D) perturbations and large scale two-dimensional (2D) perturbations. The former provide an indication of a forward transfer of energy while the latter indicate an inverse transfer and the possibility of an inverse cascade. The analysis is performed using a Floquet-Bloch code that allows to examine the stability of modes with arbitrary large scale separation. For thin layers the 3D perturbations become unstable when the layer thickness H becomes larger than

, where U is the rms velocity of the flown, U is the correlation length scale of the flow, ν the viscosity and Re = U U/ν is the Reynolds number. At the same time large scale 2D perturbations also become unstable by an eddy viscosity mechanism when Re > c2, where c1, c2 are order one non-dimensional numbers. These relations define different regions in parameter space where 2D and 3D instabilities can (co-)exist and this allows to construct a stability diagram. Implications of these results for fully turbulent flows that display a change of direction of cascade as H is varied are discussed.

I. INTRODUCTION

In many systems in nature there is a transfer of energy to both large and small scales [START_REF] Celani | Turbulence in More than Two and Less than Three Dimensions[END_REF][START_REF] Marino | Inverse cascades in rotating stratified turbulence: Fast growth of large scales[END_REF][START_REF] Pouquet | Geophysical Turbulence and the Duality of the Energy Flow Across Scales[END_REF][START_REF] Seshasayanan | On the edge of an inverse cascade[END_REF][START_REF] Deusebio | Dimensional transition in rotating turbulence[END_REF][START_REF] Seshasayanan | Critical behavior in the inverse to forward energy transition in two-dimensional magnetohydrodynamic flow[END_REF][START_REF] Sozza | Dimensional transition of energy cascades in stably stratified forced thin fluid layers[END_REF][START_REF] Marino | Resolving the Paradox of Oceanic Large-Scale Balance and Small-Scale Mixing[END_REF][START_REF] Santiago | Critical transitions in thin layer turbulence[END_REF][START_REF] Musacchio | Split energy cascade in turbulent thin fluid layers[END_REF]. This transfer at large Reynolds numbers leads to a split cascade of energy such that part of the energy cascades to the small scales and gets dissipated while at the same time part of the energy cascades to the large scales typically forming large scale coherent structures [START_REF] Byrne | Robust inverse energy cascade and turbulence structure in three-dimensional layers of fluid[END_REF][START_REF] Seshasayanan | Condensates in rotating turbulent flows[END_REF][START_REF] Alexakis | Rotating taylor-green flow[END_REF][START_REF] Antonio M Rubio | Upscale energy transfer in three-dimensional rapidly rotating turbulent convection[END_REF][START_REF] Shats | Turbulence Decay Rate as a Measure of Flow Dimensionality[END_REF][START_REF] Xia | Upscale energy transfer in thick turbulent fluid layers[END_REF][START_REF] Francois | Inverse energy cascade and emergence of large coherent vortices in turbulence driven by faraday waves[END_REF][START_REF] Frishman | Turbulence statistics in a 2d vortex condensate[END_REF]. For a recent review and classification of such split cascades see [START_REF] Biferale | Cascades and transitions in turbulent flows[END_REF]. A typical example of such a behavior is observed in the Earth's atmosphere where, due to the combined effects of confinement, rotation and stratification, there is this split cascade leading to both small scale turbulence and fast mixing as well as large-scale structures like zonal flows and hurricanes. This split cascade has been quantified by in situ aircraft measurements in the hurricane boundary layer [START_REF] Byrne | Height-dependent transition from 3-D to 2-D turbulence in the hurricane boundary layer[END_REF]. In the presence of a split cascade the amplitude of the inverse cascade depends on a control parameter, whose variation, makes the system transition from a state exclusively cascading energy to the small scales to a state that the energy cascade is split or strictly inverse. The simplest example perhaps is turbulence in layers of finite thickness H such that at scales L much larger than the thickness, L H, the flow looks like two dimensional (2D) while at scales much smaller than the thickness, H

, the flow behaves like a three dimensional (3D) flow. Depending on the relative value of H with respect to the forcing scale U the system can behave like 3D cascading energy forward or like 2D cascading energy inversely. The amplitude of the inverse cascade appears to decrease as H/ U is increased until a critical value is met H c / U such that the inverse cascade becomes exactly zero.

The mechanisms however under which such a coexistence of counter directed cascades exist are not yet understood. In this work we try to unravel some of the mechanisms involved by examining the stability properties of simple laminar flows to both large and small scale perturbations. Although such an approach can not be used to study the full nonlinear dynamics of a split cascade it provides a much simpler setup that can shed light on the mechanisms involved in the transfer of energy to small or larger scales. On the one hand it is known that 3D instabilities of laminar flows generate in general smaller scales transferring energy to them. On the other hand 2D instabilities can couple the forced modes to large scale 2D modes transferring energy to larger scales. This non-local interaction of scales can be quantified with the use of an eddy-viscosity that can change sign depending on the flow parameters. A possible way to understand the mechanisms involved in the transition from a forward to an inverse cascade is to examine when 3D instabilities dominate and when there is a change of sign for the eddy viscosity of the flow.

The notion of eddy viscosity has been introduced in turbulence very early by Taylor [START_REF] Taylor | Eddy motion in the atmosphere[END_REF]. First attempts for its calculation were made by Kraichnan [START_REF] Robert H Kraichnan | Eddy viscosity in two and three dimensions[END_REF] in an attempt to quantify the loss or gain of energy to the small scales. In this framework the evolution of a weak large scale flows v due to small scale turbulent fluctuations U can be described as

∂ t v i = - j,k,m νδ i,j δ l,m + ν i,j,l,m eddy ∇ l ∇ m v j ( 1
)
where ν is the regular viscosity and δ i,j stands for the Kronecker delta. The tensor ν i,j,l,m eddy is the eddy viscosity that models the effect of the small scale fluctuations U on v. Since the work of Kraichnan there have been many attempts to calculate an eddy viscosity for turbulent flows [START_REF] Forster | Large-distance and long-time properties of a randomly stirred fluid[END_REF][START_REF] Fournier | Remarks on the renormalization group in statistical fluid dynamics[END_REF][START_REF] Yakhot | Renormalization group analysis of turbulence. i. basic theory[END_REF][START_REF] Zhou | Renormalization-group theory for the eddy viscosity in subgrid modeling[END_REF][START_REF] Yakhot | Renormalization group formulation of large-eddy simulations[END_REF][START_REF] Zhou | Renormalized eddy viscosity and kolmogorovs constant in forced navierstokes turbulence[END_REF][START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF]. However the cascade processes of turbulent flows that excites a continuous spectrum of scales prevents from having a closed expression for the eddy viscosity and even the notion of eddy viscosity for turbulent flows can be questioned. One case where the eddy viscosity tensor ν i,j,l,m eddy can be rigorously defined is in the presence of large scale separation between the forced small scale field U that evolves at scales U and the flow velocity u that evolves at scales L U . It can then be calculated in the low Reynolds number limit using homogenization theory [START_REF] Dubrulle | Eddy viscosity of parity-invariant flow[END_REF][START_REF] Gama | Negative eddy viscosity in isotropically forced two-dimensional flow: linear and nonlinear dynamics[END_REF]. Here the Reynolds number is defined as Re ≡ U U /ν where U is the rms of the small scale velocity field and U is the typical length-scale of the flow. In the small Reynolds number limit and for isotropic flows it takes the form

ν i,j,l,m eddy = c 1 ν U 2 2 U δ i,j δ l,m + O (Re) for Re 1, (2)
where c 1 is a non-dimensional number that depends on the detailed structure of the small scale field U. For large Re a viscosity independent value of ν eddy is expected to be reached

ν i,j,l,m eddy = c 2 U U δ i,j δ l,m + O 1 Re for Re 1. (3)
If -ν i,j,l,m eddy has positive eigenvalues that are larger than the viscosity then the system can develop large scale instabilities [START_REF] Cameron | Large-scale instabilities of helical flows[END_REF].

An other possibility for the effect of small scales on the large is the anisotropic kinetic alpha effect [START_REF] Frisch | Large-scale flow driven by the anisotropic kinetic alpha effect[END_REF]. This effect is similar to the alpha effect in dynamo theory, and leads to a faster growth-rate compared to a negative eddy viscosity.

However, unlike the dynamo problem that requires only the presence of helicity, in hydrodynamics in order for the anisotropic kinetic alpha effect to be present the flow at small scales needs to satisfy much stricter conditions that are not met in the examined flows in this work. Fot this reason this work is only going to consider the effect of negative eddy viscosity as a large scale instability mechanism.

In practice it is not always feasible to calculate ν eddy analytically, and calculations are limited to either very simple flows [START_REF] Ld Meshalkin | Investigation of the stability of a stationary solution of a system of equations for the plane movement of an incompressible viscous liquid[END_REF][START_REF] Aa Nepomnyashchy | On the stability of the secondary flow of a viscous fluid in an infinite domain[END_REF][START_REF] Sivashinsky | Negative viscosity effect in large-scale flows[END_REF][START_REF] Gregory | On negative eddy viscosity under conditions of isotropy[END_REF] or in the low Re limit where the eddy-viscosity is sub-dominant to the regular viscosity. In this work a different path is followed by calculating the growth-rate of small and large scale instabilities using Floquet-Bloch theory from which the value of the eddy viscosity can be extracted. By doing so we can determine for the examined flows when the regions in the parameter space that the eddy viscosity is negative and large scale instabilities are present, and where the flow is dominated by three dimensional small scale instabilities instead.

II. FORMULATION

We begin by considering the flow of a unit density fluid in a triple periodic box of size (2πL, 2πL, 2πH) with L H so that the layer has a small height compared to its length in the other directions. The triple periodic box although unphysical it simplifies significantly the following calculations. No significant differences are expected from a flow in layers with free slip (stress-free) boundary conditions. For no-slip boundary conditions however the drag force from the boundary can have a significant effect on the dynamics. The flow satisfies the Navier-Stokes equation that is given by:

∂ t u + u • ∇u = -∇P + νF + ν∆u. (4)
Here, u is the fluid velocity, ν is the viscosity and F is an external body force that maintains the flow. We assume the F is such that it supports a laminar time independent solution for the fluid U.

F = -ν∆U + U • ∇U -∇P.
For simplicity a simple family of flows is going to be considered here, given by

U = U cos rπ 2   sin(k U y) sin(k U x) 0   + √ 2U sin rπ 2   sin(k U y) sin(k U x) 0   sin(k U z), (5)
where k U = 1/ U is the wavenumber of the flow and U = U is the L 2 norm of the flow. Note that U is independent of the parameter r. For the remainder of this investigation we will set k U = 1 and U = 1. With this normalization the Reynolds number of the flow is always given by

Re = U U ν = 1 ν . (6)
For r = 0, the flow is a two dimensional (2D) cell flow and is a solution of the Euler equations. For r = 0, U is three dimensional (3D) as it varies in all three dimensions. For r = 1 the flow reduces to a Taylor-Green flow [START_REF] Taylor | Mechanism of the production of small eddies from large ones[END_REF] a prototypical flow that was used to study the small scale generation. Note that if r is not zero U is not a solution of the Euler equations and F has to have a more complex behavior than U to sustain it against viscosity and the non-linearities. Furthermore we note that for r = 0, H can take any value, while for r = 0, H is restricted to be an integer multiple of the flow period H = n U with n an integer.

The main objective of this work is to examine the stability of this flow and examine the evolution of an infinitesimal perturbation v to the laminar flow so that u = U + v. The linear evolution equation for v reads

∂ t v + U • ∇v = -v • ∇U -∇P + ν∇ 2 v. (7)
The linearity and homogeneity in time of the problem implies that asymptotically at large times v will have an exponential behavior with time and the goal is to determine the growth rate γ of the perturbation in terms of the parameters of the system. Because U is periodic in space and we are interested in the behavior of scales much larger than the scale of U it is convenient to use Floquet (or Bloch) theory. Floquet theory states that if U is periodic the perturbation field v can be decomposed as

v(x, t) = e iq•x ṽ(x, t) + c.c. (8)
where ṽ(x, t) is a complex vector field that has the same spatial periodicity as the velocity field U, and q is an arbitrary wave number. The linear evolution equation for the field ṽ(x, t) then becomes

∂ t ṽ + U • ∇ṽ + iU • qṽ = -ṽ • ∇U -∇ P + ν(∇ + iq) 2 ṽ, (9)
where the pressure P is such so that it guarantees the incompressibility condition ∇ṽ + iq • ṽ = 0. The advantage of studying eq. ( 9) numerically as opposed to eq. ( 7) is that one can consider arbitrary large scale separations (determined by the vector parameter q) with no additional computational cost. Furthermore, the Floquet formulation gives a clear distinction between small scale and large scale instabilities. For q = |q| < k U , the volume average ṽ over the spatial period (2π U ) 3 gives the amplitude of v at large scales L ∝ 2π/q (ie v = ṽ e iq•x + c.c.). Fields with q = 0 and/or ṽ = 0 correspond to purely small scale fields, and if these modes are unstable it amounts to a pure small scale instability ie : an instability that involves only wavenumbers that are integer multiples of the basic period k U .

Eq. ( 9) is solved numerically using a pseudospectral code. Details of the code can be found in [START_REF] Cameron | Large-scale instabilities of helical flows[END_REF][START_REF] Cameron | Fate of alpha dynamos at large r m[END_REF]. The resolutions examined varied from 32 3 grid points for most of the results here although resolutions of 64 3 grid points and 128 3 grid points were also used for large Re and to test convergence.

III. TWO DIMENSIONAL FLOWS r = 0

We begin by examining the case for which r = 0 so that the laminar flow is two-dimensional. In this case the laminar flow can be written in terms of a stream function

Ψ = cos(k U x) -cos(k U y) as U = [∂ y Ψ, -∂ x Ψ, 0], where the stream function is connected to the vertical vorticity W as W = -∇ 2 Ψ = k 2 U Ψ.
The last equality holds because U contains modes of single wavenumber k U . The stability of this 2D flow will be examine in the next subsections, first against 2D perturbations and second against 3D perturbations. The stability for 2D perturbations has been examined before in [START_REF] Thess | Instabilities in two-dimensional spatially periodic flows. part ii: Square eddy lattice[END_REF] however in the presence of large scale linear friction that does not allow to directly compare with these results.

A. Two dimensional perturbations

For velocity perturbations that are also 2D (ie ∂ z v = 0) the vertical component of the perturbation velocity decouples and follows a passive advection diffusion equation

∂ t v z + U • ∇v z = ν∇ 2 v z . ( 10
)
As a consequence its L 2 norm v z 2 follows

∂ t v z 2 = -ν ∇v z 2 ≤ -ν 2π L 2 v z 2
and therefore the norm v z 2 decays monotonically with time. The remaining two components v 2D = [v x , v y , 0] can be written in terms of a stream function ψ as v x = ∂ y ψ, v y = -∂ x ψ. The evolution of these components can be written in terms of their vorticity w = -∇ 2 ψ as

∂ t w + U • ∇w = -v 2D • ∇W -ν∇ 2 w. (11)
Multiplying eq. ( 11) by w and space averaging we obtain after some manipulation (and using the laminar flow property

W = -∇ 2 Ψ = k 2 U Ψ) the enstrophy evolution equation ∂ t w 2 = k 2 U wU • ∇ψ -ν ∇w 2 . ( 12
)
Multiplying eq. ( 11) by ψ and space averaging we obtain after similar manipulations the 2D energy evolution equation

∂ t v 2D 2 = -ψU • ∇w -ν w 2 . ( 13
)
Multiplying the energy evolution equation by k 2 U and subtracting from the enstrophy equation we obtain

∂ t ( w 2 -k 2 U v 2D 2 ) = -ν( ∇w 2 -k 2 U w 2 ) ( 14
)
In general, the terms in parenthesis on the left and right hand side of the equation above 14 can have either sign and the equation above can not exclude the growth of w. If however k U = 2π/L so that the domain size L coincides with the velocity scale 2πk -1 U (ie no scales larger than the forcing scale are allowed) both (

w 2 -k 2 U v 2D 2
) and ( ∇w 2 -k 2 U w 2 ) are non-negative (by Poincaré inequality). The right hand side of eq. ( 14) is then negative and the positive quantity (

w 2 -k 2 U v 2D
2 ) has to decrease monotonically. This implies that single scale 2D flows such that

∇ 2 Ψ = -k 2 U Ψ
where k U is the largest scale of the domain are linearly stable to all 2D perturbations. In fact eq. 12 and 13 also hold for the nonlinear evolution of w and v 2D and the flow is thus non-linearly stable to all 2D perturbations as shown for the first time in [START_REF] Marchioro | An example of absence of turbulence for any reynolds number[END_REF]. In terms of the Floquet decomposition given in eq. ( 8) this conclusion translates to: the flow is stable for q = 0 ie there are no small scale instabilities, and any unstable mode has to excite simultaneously scales both smaller and larger than the forcing.

Given this restriction we investigate numerically eq. ( 9) using a pseudo spectral code (described in [START_REF] Cameron | Large-scale instabilities of helical flows[END_REF][START_REF] Cameron | Fate of alpha dynamos at large r m[END_REF]) and calculate the growth-rate γ of the unstable modes as a function of q and the viscosity ν. The procedure to calculate this growth-rate is as follows. We chose the wave-vector q = [q x , 0, 0]. In this case since we are interested in 2D perturbations q is restricted in the x,y plane. For simplicity, results are presented only for q = [q x , 0, 0] that was found (but not proven) to be the most unstable from the general cases examined q = [q x , q y , 0]. For this chosen wave-vector q = [q x , 0, 0] the complex velocity field ṽ is initialized using random initial conditions. Then the complex field ṽ is evolved based on eq. 9 until a clear exponential growth of the energy E = |ṽ| 2 ∝ e 2γt is observed. The growth-rate γ is then the calculated by fitting. This procedure is then repeated for different values of q and different values of viscosity. This allows to calculate the growth-rate of the most unstable mode as a function of q and ν. It is noted that this procedure reveals only the fastest growing unstable mode and not all unstable modes in the system. It is worth keeping this in mind in the following sections where general 3D perturbations will be considered. The resulting growth rate γ as a function of q for different values of the viscosity is shown in the left panel of figure 1 in linear scale. The right panel shows the absolute value of the same data in logarithmic scale. As discussed in the beginning of this section for q = 0 there can be no instabilities and the growth rate becomes zero as q → 0. The growth rate γ is approaching zero as q → 0 following a quadratic power-law. This is demonstrated more clearly in logarithmic scale shown in the right panel of figure 1 where it is fitted as

γ = a q 2 (15)
shown by the dotted lines. The proportionality coefficient a is negative for large values of the viscosity, while bellow a critical value it becomes positive. For larger q higher even powers (due to symmetry) of q are needed to describe γ = a q 2 +βq 4 +. . . and as ν is decreased the maximum growth-rate moves to larger q. At ν sufficiently small the most unstable mode is found at q = 1/2, indicating that the wavenumbers close to the forcing one are the most unstable. Nonetheless, the onset of the large scale instability is determined by the sign of the pre-factor a of the quadratic power-law in eq. 15. This particular scaling of the growth-rate γ with the wavenumber q given in eq. 15 implies that the instability is of the form of an eddy-viscosity. The eddy viscosity is a tensor and the proportionality coefficient a in eq. ( 15) indicates the maximum eigenvalue of -νδ i,j + ν i,j,x,x eddy . For simplicity a + ν will be denoted as -ν eddy and referred as eddy viscosity although it must not be forgotten that the eddy viscosity is truly a tensor. The relation of ν eddy with a is a = -(ν + ν eddy ). ( 16) The eddy viscosity depends on the small-scale flow examined and the value of viscosity. Figure 2 shows the value of ν + ν eddy as a function of ν for the examined flow. The dotted lines indicates the prediction (2) for c 1 = 0.5. The fit to the power-law is surprisingly good even up to ν 0.05 which is far from the asymptotic regime Re 1. Negative values of ν + ν eddy imply that the flow amplifies large scale flows. For large values of ν the value of the eddy-viscosity is such that ν + ν eddy is positive although still smaller than ν. So although ν eddy is negative and the flow reduces the decay rate caused by regular viscosity it is not sufficiently strong to drive any large scale scale instabilities. Bellow some critical value of the viscosity ν = ν * the value of ν + ν eddy becomes negative. For values of ν smaller than ν * the large scales of the flow are unstable. The critical value of the viscosity ν * can be estimated to be

ν * 0.7 U U ( 17
)
where we have kept U = 1 and U = 1 in the expression to recover the dimensions.

Figure 3. Energy spectrum of the unstable modes e iq•x ṽ(x, t) for different values of q and ν = 0.2, and k U = 1 is the wavenumber of the laminar flow. The value of q is indicated by the furthest point to the left.

Before concluding this section it is worth examining the spectral shape of the unstable modes. Figure 3 shows the energy spectra of the unstable modes as a function of the wavenumber for different values of q. The spectrum is composed by the large scale component located at the wave number q plus the part localized at small scales k ≥ k U = 1. The amplitude of the energy at large scales is independent on the value of q. This is seen in the figure 4 where the ratio of the energy contained in the largest scale E 0 = 1 2 | ṽ | 2 to the total energy E = 1 2 |ṽ| 2 is shown. Note that E 0 /E is shown also for modes with negative growth rate. The amplitude of this ratio is initially decreased as q is decreased but asymptotes to a q independent value for sufficiently small values of q. The fact that the q = 1/2 modes have a larger value of E 0 /E than small q also indicates that the q = 1/2 are more efficient at extracting energy from the forcing scale. This asymptotic value at small q decreases as ν decreases, which can be shown in the asymptotic expansion in [START_REF] Dubrulle | Eddy viscosity of parity-invariant flow[END_REF].

B. Three dimensional perturbations

In the previous section we examined the stability of the 2D flow against large scale 2D perturbations. In this section we investigate the stability of the flow against 3D perturbations. This question is important because it determines when 3D variations enter the system that can transfer energy to the small scales.

For r = 0 the flow U has no dependence on the z direction. We can therefore consider a layer with arbitrary thickness compared to the flow length-scale k -1

U . The homogeneity in the z-direction also implies that each q z mode can be considered independently for the linear problem with the perturbation field ṽ being a function of x, y, t only and not of z. This section is restricted to 3D perturbations with wavenumbers q = [0, 0, q z ]. For this choice the field ṽ satisfies the incompressibility condition

∂ x ṽx + ∂ y ṽy + iq z ṽz = 0. ( 18
)
Note that unlike the previous case of 2D perturbations (with q z = 0) the incompressibility condition implies that the z-component does not act independently. Figure 5 displays the growth rate γ as a function of q z for different values of the viscosity ν. Bellow a critical value of ν < ν 0 0.1 modes with positive growth rates appear with q z ∝ k U = 1. As ν is decreased further the range of unstable wave-numbers extends to larger and large values of q z . For ν ν 0 the unstable wavenumbers are found in the range q z < q * z , here q * z is the largest unstable wavenumber. This critical wavenumber q * z can be estimated from the graph in the left panel of figure 5 as the wavenumber for which the growth rate intersects the zero growth rate γ = 0 line. It depends on ν and it is plotted in the right panel of figure 5 in a logarithmic scale. It is shown to follow the scaling q * z ∝ ν -1/2 for small ν. This scaling is obtained from a balance of the viscous term νq 2 ṽ in eq. ( 9) with the stretching term v • ∇U. A best fit results in

q * z 0.8 U ν U 1/2 (19)
where we have kept U = 1,and U = 1 in the expression to recover the dimensions and the result is valid for small values of ν (large q * z ). For layers of finite height H the presence of these unstable modes depends on H. This is because the periodicity of the domain in the z direction imposes that the wavenumber of a 3D perturbation satisfies q z = n/H where n ≥ 1 is an integer. Therefore even for ν < ν 0 the flow will be stable to 3D perturbations if the layer thickness H is thin enough so that 1/H is larger than the largest unstable wavenumber q * z . It is already possible to draw a first conclusion for 2D flows. From the analysis in section III A it was shown that the large 2D modes are unstable if the viscosity is small enough

ν < ν * 0.7U U (20)
so that the eddy viscosity for 2D perturbations satisfies ν eddy + ν < 0 (see eq. ( 17)). At the same time in section III B it is shown that 3D modes are stable if the layer is thin enough so that 1/H is larger than maximum unstable wave number q * z (see eq. ( 19)). This implies stability of the 3D modes if

H ≤ 1 q * z ν 1/2 0.8(U k U ) 1/2 . ( 21
)
Therefore there is a region is the parameter space (H, ν) determined by the conditions 20,21 where 3D modes are stable while 2D modes are unstable to large scale perturbations by an eddy viscosity mechanism. This is summarized in the diagram in figure 6 where condition 20 is given by the solid line while the condition 21 is marked by a dashed line. The region of 3D stable and 2D unstable modes is marked by gray. In this region a transfer of energy to the large scales is expected and possibly an inverse cascade can develop. Of course whether an inverse cascade builds up or not depends on the nonlinear evolution of the unstable modes and can not be determined by the stability of the laminar flow alone. Furthermore in the region that both 3D and 2D instabilities are present we can not certify nor exclude a transfer of energy to the large scales. This is because we can not a priori conclude if the nonlinear evolution of the unstable 3D modes will suppress or not the inverse transfer by the negative eddy viscosity mechanism.

IV. THREE DIMENSIONAL FLOWS

This next section considers the stability of 3D laminar flows given in eq. ( 5) with in general r = 0. For these flows the minimum layer thickness H one can consider is given by the periodicity of U in the z direction H = U . Since the interest in this work is on thin layer for the remaining of this section we will fix H to this minimum value. This limits the values of q z in the Floquet expansion to zero q z = 0 since all perturbations can only have vertical scales smaller than H. The investigation therefore will be restricted to large scale 2D perturbations with q of the form q = [q x , 0, 0]. It is noted that setting q z = 0 does not imply that v = e iqx ṽ is independent of the coordinate z because the field ṽ depends in principle in all coordinates x, y, z having the same periodicity as U. Therefore 3D perturbations can still be unstable and as it is shown later they play an important role. Figure 7 shows the measured growth rate as a function of the wavenumber q x for a flow with r = 0.4 for different values of ν. For large values of ν the growth rate shows a similar behavior as for the r = 0 case shown in figure 1. For ν below a critical value the growth rate is negative and with a quadratic dependence in q x as in eq. [START_REF] Shats | Turbulence Decay Rate as a Measure of Flow Dimensionality[END_REF]. Above this critical value the growth rate retains its quadratic behavior (eq.15) but with a positive value of the proportionality coefficient a. This behavior can then again be interpreted as an effect of negative eddy viscosity such that γ = -(ν + ν eddy )q 2

x . However, if viscosity is further decreased the quadratic dependence on q x is lost and the growth rate attains a finite value at q x = 0.

A finite growth rate at q = 0 can be a bit surprising since the growth of a non zero mean v = ṽ = 0 field would violate momentum conservation. The contradiction is resolved by examining the relevant projection of the unstable mode to the large scales by looking at the ratio of E 0 /E = | ṽ | 2 / |ṽ| 2 that is shown in figure 8. This figure demonstrates that while for modes for which γ → 0 as q → 0 the ratio E 0 /E remains finite, for modes with γ = 0 as q → 0 the projection to the large scales decreases as E 0 /E ∝ q 4 . At q = 0 therefore these last modes have a zero projection to the large scale modes, and correspond to purely small scales instabilities [START_REF] Cameron | Large-scale instabilities of helical flows[END_REF][START_REF] Cameron | Fate of alpha dynamos at large r m[END_REF].

Therefore the linear evolution of the large scale modes can be expressed in terms of an negative eddy viscosity up to this second critical value of ν = ν 3D bellow which 3D small scale instabilities begin. Large scale instabilities therefore are limited to values of the viscosity ν in the range ν * > ν > ν 3D , where ν * is the value above which large scales are stable, and for ν smaller than ν 3D small scale instabilities are present. In principle for values of ν < ν 3D large scale modes (0 < q 1) with positive growth rates γ = -a(ν + ν eddy )q 2 > 0 still exist but their growth rate can be smaller than the growth rate of small scale instabilities of the q = 0 modes. As a result for ν < ν 3D we can not predict if the flow would be able to transfer energy to the large scales with the linear evolution model given in 7 and a nonlinear theory for the upscale transfer of energy is required. In the range ν > ν 3D , the value of the eddy viscosity can be evaluated by fitting to a parabola as was done in section III A. The values of the eddy viscosity for different values of ν and different values of r were calculated and are plotted in the left panel of figure 9. This allows us to estimate the value of ν * at which large scale 2D modes become unstable. Similarly we can calculate the value ν 3D by looking at the growth rate of small scale instabilities q = 0 and determining by extrapolation at which value of ν the growth rate γ(q = 0) of small scale modes becomes positive. This is shown in the right panel of figure 9 for different values of r.

The end result is shown in figure 10 where the stability diagram of the parameter space is shown in the parameter space (ν, r). The space is split in three regions. One (marked by white) for which the growth rate of all instabilities is negative, one (marked by light gray) for which large scale modes are unstable by a negative eddy viscosity instability but 3D instabilities have negative growthrate, and finally one region (marked by dark gray) for which both small-scale and large-scale instabilities are present. For 2D flows r = 0 the region for which the effective viscosity is negative, is larger. As r is increased, the value of ν 3D (dashed line) is increased and the value of ν * (solid line) is decreased and as a result the range of values for which large scale instabilities can be observed is shrunk. Nonetheless even for r = 1 for which the laminar flow does not contain any 2D part there is still a small range of values of ν for which only the large scale modes are unstable. This implies that transfer of energy to large scales can be observed in 3D flows as well.

V. CONCLUSIONS

This work examined the linear stability of some basic linear flows against 2D large scale instabilities 3D small scale instabilities. The results showed that for the examined flows, 2D instabilities appeared first as the viscosity was decreased transferring energy to larger scales. The 3D instabilities appeared afterwards at smaller values of the viscosity and transfer energy to smaller scales. At sufficiently small viscosity therefore both types of instabilities are present and thus transfer to both large and small scales is possible.

In more detail, 2D laminar flows in thin layers become unstable to the large scale 2D modes due to an negative eddy viscosity type instability above a critical Re independent of the layer height H. The onset of this instability is determined solely by the value of the eddy viscosity. The onset Reynolds number for 3D instabilities on the other hand increased as the layer thickness decreased. This result implies that for sufficiently thin layers the 2D part of the flow can be unstable and transfer energy to the large scales while 3D instabilities are suppressed. This result already indicates that depending on the layer thickness the flow can have different behaviors. For small H it can behave like a 2D flow cascading energy inversely while if H is large the flow can be dominated by 3D instabilities cascading energy forward. The stability diagram for the examined 2D flows is given in figure 6.

For 3D flows it was shown that the negative eddy viscosity type of instabilities persist. Therefore, even for 3D flows there are unstable modes that move energy to larger scales. The large scale instabilities compete with the 3D small scale instabilities if Re is above their onset. For the examined flows is was shown that 2D instabilities have in general a smaller onset Re than the 3D instabilities, and thus a transfer to the large scales is still possible. However, as the flow becomes more 3D the difference between the two critical values decreases leaving little room to observe a large scale instability. These results are summarized in figure 10.

Finally we need to give a word of caution. The present results deal only with the stability of the laminar flows. Therefore we can only draw conclusions for the early stages of evolution of these instabilities. First, if the flow is 3D unstable and 3D turbulence develops then the transfer of energy to the large scales is not necessarily through an instability. A new nonlinear definition of an eddy-viscosity would be required if the notion of eddy-viscosity is at all valid in this case. Second, if the laminar flow is only 2D unstable then at the non-linear stage the flow can start to develop an inverse cascade [START_REF] Celani | Turbulence in More than Two and Less than Three Dimensions[END_REF][START_REF] Santiago | Critical transitions in thin layer turbulence[END_REF][START_REF] Musacchio | Split energy cascade in turbulent thin fluid layers[END_REF] and form condensates [START_REF] Van Kan | Condensates in thin-layer turbulence[END_REF]. The stability properties these new states however will change and can be unstable to 3D perturbations. Interestingly, the shape of the condensate formed in a 2D inverse cascade in a square box is similar to the 2D flows examined here [START_REF] Gallet | A two-dimensional vortex condensate at high reynolds number[END_REF][START_REF] Frishman | Jets or vortices -what flows are generated by an inverse turbulent cascade?[END_REF] and the stability against 3D perturbations can be applied but they need to be examined for vortexes the size of the domain. In general however the co-existence of an inverse and forward cascade can not be neither confirmed or excluded just from the linear stability analysis. This would require a study of the nonlinear state of the system that can be studied in some cases with bounds [START_REF] Gallet | Exact two-dimensionalization of low-magnetic-Reynolds-number flows subject to a strong magnetic field[END_REF][START_REF] Gallet | Exact two-dimensionalization of rapidly rotating large-reynolds-number flows[END_REF] or with the use of numerical simulations, and where the notion of eddy viscosity is not well established.

Figure 1 .

 1 Figure1. The growth rate γ as a function of the wavenumber q for five different values of the the viscosity ν. The dotted lines show fits to a quadratic power-law γ = aq 2 . Left panel: linear scale. Right panel: The absolute value of γ in log log scale.

Figure 2 .

 2 Figure 2. The value of ν eddy + ν as a function of ν. When this quantity takes negative values there is a large scale instability. The dashed line indicates where ν eddy = 0. The doted line shows the fitting to the asymptotic result ν eddy = -0.5/ν. Large scales become unstable when ν < ν * 0.7.

Figure 4 .

 4 Figure 4. The ratio of the energy in the large scales mode E0 = 1 2 | ṽ | 2 to the total energy E = 1 2 |ṽ| 2 as a function of q for different values of ν. Left: on a linear scale, right in log-log scale.

Figure 5 .Figure 6 .

 56 Figure 5. Left panel: The energy growth rate as a function of qz. Right panel the critical wave number q * z bellow which qz modes are unstable as a function of the viscosity.

Figure 7 .

 7 Figure 7. The growth rate γ as a function of the wavenumber qx in linear scale (left panel) and logarithmic scale (right panel) for different values of ν and for r = 0.4.

Figure 8 .

 8 Figure 8. The ratio E0/E of the large energy E0 = 1 2 | ṽ | 2 to the total energy of the perturbation E = 1 2 |ṽ| 2 as a function of q, for different values of ν in linear scale (left panel) and logarithmic scale (right panel). The color-shade of the lines is the same for both panels.

Figure 9 .

 9 Figure 9. Left panel: Effective viscosity as a function of ν for the different values of r Right panel: Small scale instability growth rate γ(q = 0) as a function of ν for different values of r.

Figure 10 .

 10 Figure 10. Stability diagram of the parameter space.

ACKNOWLEDGMENTS

This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche and the HPC resources of GENCI-TGCC-CURIE & GENCI-CINES-JADE (Project No. A0010506421 & A0030506421) where the present numerical simulations have been performed. This work has also been supported by the Agence nationale de la recherche (ANR DYSTURB project No. ANR-17-CE30-0004). The author would also like to thank A. van Kan and S. Benavides for their critical reading of this manuscript.


Francois Renard
email: francois.renard@ujf-grenoble.fr

Dominique Bernard
email: bernard@icmcb-bordeaux.cnrs.fr

Jacques Desrues
email: jacques.desrues@hmg.inpg.fr

Erwan Plougonven

Audrey Ougier-Simonin

François Renard

-Erwan Plougonven

Characterisation of hydraulic fractures in limestones using X-ray microtomography

Keywords: X-ray computed microtomography, hydraulic fracture, permeability tensor MOTS-CLÉS: microtomographie, fracturation hydraulique, tenseur de perméabilité

Introduction

The growth of a fracture in rocks under the driving effect of a propagating fluid arises in various geological systems. Natural fractures involve the propagation of intrusive magmatic dykes and sills in the Earth's lithosphere, and sometimes the release of magma at the Earth's surface in volcanoes. Hydraulic fracturing is used in the oil industry to stimulate oil and gas production in low permeability underground reservoirs. Other industrial applications involve the extraction of heat in geothermal reservoirs or induced caving in the mining industry.

Usually, models of hydraulic fractures consider a 2D planar fracture plane inside a homogeneous elastic medium [START_REF] Bunger | Asymptotic solution for a penny-shaped near-surface hydraulic fracture[END_REF][START_REF] Zhang | Propagation of a hydraulic fracture parallel to a free surface[END_REF][START_REF] Garagash | Plane-strain propagation in a fluid-driven fracture during injection and shutin: Asymptotics of large toughness[END_REF]. However, rocks are heterogeneous at all scales. It can be shown that among other properties, the cohesion heterogeneity has a crucial effect on the fracturing pattern. Discrete simulations using elastic spring networks with disorder demonstrate that hydraulic fracture can even acquire fractal geometries, with strong spatial and temporal correlations as the fracture propagates into successive bursts clustered in time (Tzshichholz and Herrmann, 1995).

In field hydraulic experiments, two types of measurements can be performed. On one hand the pressure fluctuations can be monitored at the injection pump, and on the other hand acoustic emission signals can be recorded at the surface. Depending on the far-field stress, the injection rate and the rock permeability, either shear fractures or tension fractures are the dominant mechanism for hydraulic fracture propagation [START_REF] Lockner | Hydrofracture in Weber sandstone at high confining pressure and differential stress[END_REF]. As hydraulic fractures migrate parallel to the main compressive stress, they represent a measure of the local main stress direction in the Earth's crust [START_REF] Talebi | Analysis of the microseismicity induced by a fluid injection in a granitic rock mass[END_REF].

What is less known is how such fractures propagate in a porous medium and how the fracture tip interacts with pores and grains at small scale. Moreover, heterogeneities in the fluid flow related to local variations of the permeability can modify the stress field at the fracture tip during its propagation. Here we use a small hydraulic fracture cell to generate fractures into porous limestones. We image the core samples, using high resolution X-ray tomography, before and after fracturation. Using these 3D data, we can 1) calculate the permeability tensor of the core samples; 2) isolate the fracture path and characterise its roughness property, based on the initial heterogeneous distribution of pores in the rock. In the following sections we describe the experimental procedure and present some results on the calculation of the permeability and the roughness of hydraulic fractures.

Experimental hydrofractures in limestones

Sample preparation and experimental set-up

A hydraulic cell was especially designed to allow the measurement of fracture onset in limestone core samples (Figure 1). Four samples (Table 1) were drilled from blocks of porous limestones extracted from two different quarries: 1) The Estaillades limestone contains debris of shells cemented by calcite. This is a 95% calcite rock, with 30% porosity, and permeability between 0.2 and 0.4 Darcy. Its mechanical resistance to compression is 13 to 18 MPa. 2) The Anstrude rock is an oolitic limestone with 20% porosity, a low permeability of 2 to 4 mD, and a uniaxial compression strength between 30 and 40 MPa.

In each core samples a hole was carefully drilled along the cylinder axis. The resulting sample is a hollow cylinder, with an external diameter of 9 mm fixed by the range of the CCD camera used for microtomography, and a central 2 mm diameter hole. The advantage of this geometry is that, as water is injected through the central hole, the pressure gradient had a cylindrical symmetry in the sample.

Each sample was placed between two brass holders and a rubber joint ensured a flat contact between the sample and the holder. The sample was free to deform radially. This setup allowed a homogeneous contact of the sample with the holder, to avoid fluid leakage out of the central hole (Figure 1).

Synchrotron X-ray tomography data

The 3D microtomographic data were acquired from beamline ID19 at the European Synchrotron Radiation Facility under energy of 25 keV and at a pixel size of 5.1 micrometers. An X-ray beam was sent onto the sample, which was rotated over a range of 180° in 1000 steps around a vertical axis, and a series of radiographs was taken using a low-noise 2048x2048 FRELON CCD camera (Fast REadout LOw Noise). Based on the sets of projections, a reconstruction algorithm was used to obtain a tomogram representing the 3D microstructure of the sample. Wellcontrasted images were obtained, resulting from the difference in X-ray absorption between calcite grains and empty pores, before and after fracturation (Figure 2b,c). Table 1. Four core samples were broken in the hydraulic cell. For all of them, a microtomography scan was acquired before and after fracturation.

Data analysis

The permeability was estimated directly from the pressure drop and flux measured during the experiments. An alternative is to calculate the permeability tensor by solving the closure problem associated with the volume averaging of Stokes equations, as described by [START_REF] Whitaker | The Forchheimer equation: A theoretical development[END_REF]. Although computationally demanding, the solution of this closure problem can be calculated using established methods [START_REF] Bernard | Permeability assessment by 3D interdentritic flow simulations on microtomography mappings of Al-Cu alloys[END_REF].

Due to the large difference between the X-ray attenuation in the pores and in the grains, the grey-level distribution of the 3D images was almost bimodal. A simple threshold operation was applied to generate binary 3D images. This 3D image can be directly used as input geometry to perform permeability calculations. Only the pores forming a connected network between the two sides of the sample were kept for computing porosity. This technique was employed on three subsamples (see Table 1), where the full permeability tensor could be calculated. For each sample, the trace of the tensor gives the permeability, which is shown to be of the same order of magnitude than the measured permeabilities. Using a Lattice Boltzmann technique to solve Stokes equations, [START_REF] Fredrich | Predicting macroscopic transport properties using microscopic image data[END_REF] have also shown that the permeability can be calculated within a factor 5 of error in a 3D porous domain. In our approach, the complete permeability tensor can be obtained. In order to highlight the fracture path in the rock, the 3D images of the two halves of the fractured sample were registered in the initial non fractured sample coordinate system (figure 2). This registration process first requires disconnecting the two halves of the broken sample (figure 2c) and processing them separately.

In the solid, some denser grains attenuate X-rays more than calcite grains, resulting in high intensity regions in the images. These areas can easily be segmented, and their centres of gravity are precisely located. As a result, a set of reference points was obtained for the unbroken sample, and for each half of the broken sample. Using these points and an iterative closest points algorithm, it was possible to calculate the translation and rotations that superimposes the point set (and thus the microtomography) of the initial and broken samples (Figure 2d). Using this technique it was possible to localise with precision the path the fracture had followed. Preliminary results on sample Anstrude4M show that the fracture had connected the pores, and therefore broken the weak bonds between the grains. Note also that some grains have been broken.

Conclusions

A 3D simulation of permeability in porous limestones based on 3D real microstructure obtained by X-ray tomography has been performed. There is a reasonable agreement between the trace of the calculated permeability tensors and the experimentally measured permeability. The deviations can be attributed to the tendency for preferential flow channels in the sample compared to the subsamples used to calculate the permeability.

The microtomography data also show that the hydraulic fractures have propagated by connecting individual pores. Numerical simulations based on the deformation of a heterogeneous elastic solid under the effect of a pressure gradient should represent a future task for fracture propagation prediction.

Figure 1 .

 1 Figure 1. a) Limestone core sample 2S after fracturation. b) Miniature hydraulic cell. c) Experimental setup, with the fluid circuit and the position of the digital measurement devices. The hydraulic pump flux was increased gradually by the operator until fracturation occurred.

Figure 2 .

 2 Figure 2. Microtomography data. a) Connected porosity in a subsample of Anstrude4M, used to calculate the permeability. The solid is made transparent. b) Sections in sample Anstrude1S before (top) and after (bottom) fracturation. c) 3D view of a subsample of Estaillade3M after fracturation. d) Broken part of sample Anstrude4M relocated in the same coordinate system as the sample before fracturation. The fracture has propagated by connecting the pores.

2006 GeoX Conference.

Acknowledgements:

The authors would like to thank Elodie Boller at ESRF, Robert Guiguet for designing the hydraulic cell, and Christophe Rousseau for his help with the hydraulic load control system and the data acquisition software.