Hugo Saulnier
email: hsaulnier@lirmm.fr

Olivier Strauss
email: strauss@lirmm.fr

Ines Couso
email: couso@uniovi.es

Comparing Interval-Valued Estimations with Point-Valued Estimations

Keywords: Interval-valued data, Imprecise probabilities, Engineering

In the last decade, numerous proposals have been made to deal with imprecision in estimation problems. Those approaches, many of which involve dealing with interval-valued outputs, deal with the subtle difference between uncertainty and imprecision. One of the crucial pointswhich to our knowledge has never been addressedis "how to compare an interval-valued method with a precise valued method?"

The usual way to compare two estimation methods is to use benchmark data with ground truths and to compute a distance between the estimates of each method and the ground truth. However, most of the mathematical available extensions of distances are either biased in favor of a precise approach or in favor of an imprecise approach.

This paper proposes a new tool, the weighted variation of the midpoint distance (WVD), that is more suitable to achieve this kind of comparison, dealing with imprecision with a particular semantic. After reviewing existing distances, we introduce the WVD, first from an intuitive perspective, then from a more mathematical point of view. Its very satisfactory properties are highlighted through an experiment.

Introduction

Scientists willing to consider imprecise data and methods in their analysis (such as [START_REF] Beer | Imprecise probabilities in engineering analyses[END_REF][START_REF] Aughenbaugh | The value of imprecise probabilities in engineering design[END_REF] or [START_REF] Oberguggenberger | Classical and imprecise probability methods for sensitivity analysis in engineering: a case study[END_REF]) face the problem of comparing methods with interval-valued output with methods with precisely valued outputs. In this paper, we consider imprecise valued regression methods producing interval valued estimates of a precise reference. One of the most usual ways for assessing a preference between one method and another is to use a set of benchmark data to compute a distance between the output of each method and the known ground truth, the current fashion in engineering problems being to rather use the L 1 distance (see [START_REF] Afonso | Image reconstruction under multiplicative speckle noise using total variation[END_REF] for example). However, if comparing two interval-valued methods or two precise valued methods is straightforward, a comparison between an intervalvalued method and a precise valued method is more intricate since any existing extension of the L 1 distance is either biased towards or against imprecision: the supremum (or Hausdorff1 ) distance tends to disfavor imprecise valued estimates while the infimum distance promotes immoderately imprecision.

In order to establish a preference between methods we need to compute the distance between the reference (precise) and the estimate (precise or imprecise). This preference must express as wisely as possible the intricate pros and cons of using an imprecise estimation versus using a precise estimation. Figure 1 illustrates the kind of situation we could fall into. A reference (plain line) is estimated (dotted line) in four different settings.

(1)

(2) Probably (2) would be preferred to [START_REF] Beer | Imprecise probabilities in engineering analyses[END_REF]. The estimation of ( 2) is imprecise but informative. It describes the possible variations of the reference. On the other hand the precise estimation of (1) gives an inaccurate description of the data. Again (3) is better than [START_REF] Aughenbaugh | The value of imprecise probabilities in engineering design[END_REF], because it has narrower bounds while still being informative. But (4) should be preferred among all, because it is both precise and accurate.

(3) (4)
In this paper, we propose a nice candidate to achieve this kind of ordering. The rest of the paper will be organized in three main sections: after reviewing existing metrics, we introduce the WVD, first from an intuitive perspective, then from a more mathematical point of view. Finally some experiments highlight its very nice properties.

A Review of Existing Distances Between Points and Intervals

Notations

Let R be the real line and IR be the set of closed intervals of R. X will denote a closed interval of IR, x ∈ R its lower element, x ∈ R its upper element, x = x+x 2 its center, and r = x+x 2 its radius. Here we are interested in defining distances between a finite sequence of points (scalar vector) and a finite sequence of intervals (interval vector). In the following we will consider distances from a reference vector of scalars s = (s 1 , . . . , s N ) to an indexed collection (vector) of intervals X = (X 1 , . . . , X N ) where ∀i ∈ {1, . . . , N}, X i = [x i , x i ]. We also denote by x the vector of the center values of X. The literature includes different alternatives to compute the distance between pairs of intervals of the real line. Based on each of those proposals, one can define the distance from a point s ∈ R to an interval X ∈ IR as the distance between the singleton {s} (which is in turn an element in IR) and X. We can therefore define the distance from s = (s 1 , . . . , s N ) to X as the arithmetic mean of the distances between their components. We next review different proposals of distances between pairs of real intervals from the literature, and construct their associated distances from vectors of points to vectors of intervals.

Hausdorff Distance

Let (U, d) be a complete metric space. Let K(U ) denote the family of non-empty compact subsets of U . The Hausdorff distance between two non-empty compact subsets A, B ∈ K(U ) is

d H (A, B) = max{sup a∈A inf b∈B d(a, b), sup b∈B inf a∈A d(a, b)}.
In particular, we can define the Hausdorff distance from s ∈ R to X ∈ IR as the Hausdorff distance2 between {s} and X, i.e.:

d H (s, I) = max{sup x∈I d(s, x), inf x∈I d(s, x)} = sup x∈I d(s, x), (1)
The Hausdorff distance is one the most widely used distances from points to sets [START_REF] De Souza | Clustering of interval data based on city-block distances[END_REF]. If d is the L 1 distance, Eq. ( 1) can be simplified in:

d H (s, I) = |s -x| + r, (2)
where x and r respectively denote the mid-point and the radius of the closed interval X = [x-r, x+r]. Computation of the distance from s ∈ R N to X ∈ IR N is obtained by averaging Expression (2) on their components:

d H (s, X) = 1 N N i=1 (|s i -xi | + r i ). (3)

Mid-Spread Distance

Expression (2) considers the interval X as a pair composed of its centre x and its radius r like in the mid-spread distance approach [START_REF] Trutschnig | A new family of metrics for compact, convex (fuzzy) sets based on a generalized concept of mid and spread[END_REF], where the distance is computed as a linear combination of the distance between the two centers and the distance between the two radii (the radius of s being equal to 0). The midspread distance between two intervals

X = [x -r, x + r] and Y = [ỹ -r , ỹ + r ] is defined as: d γ ms (X, Y ) = |x -ỹ| + γ|r -r |.
According to it, we can define the mid-spread distance from s to X = [x+r, x+r] as:

d γ ms (s, X) = |s -x| + γr.
Based on it we can calculate the mid-spread distance from s to X as:

d γ ms (s, X) = 1 N N i=1 (|s i -xi | + γr i ), (4)
where γ ∈ R + is a weight given to the radius w.r.t. the mid. This approach can be seen as a generalization of the Hausdorff distance since d H = d 1 ms .

Mid-Point Distance

The mid-point distance leads to a less technical distance. It is just the L 1 distance between s and x. Within this approach the intervals are reduced to points and the imprecision is simply ignored. This approach is easy both conceptually and computationally. For this reason it is used a lot in practice [START_REF] Augustin | Statistical inference[END_REF]. It is defined as:

d m (s, X) = 1 N N i=1 |s i -xi |. ( 5
)
It can also be considered either as a special case of the mid-spread distance, where γ = 0 or as a two step Hausdorff distance. In this last interpretation, the method consists of choosing, as a representative single point of the interval X, the precise value that minimizes the Hausdorff distance of the L 1 distance, i.e. its center x, and then compute the usual distance between s and this representative.

Infimum Distance

The classical way to define the distance from a point s to a set X in topology is to consider the shortest distance from s to any point in X, i.e.:

d inf (s, X) = inf x∈X d(s, x). (6)
When, in particular, d denotes the L 1 distance from a point s to a closed interval X = [xr, x + r], Expression (7) reduces to:

d inf (s, X) = 0 i f x ∈ X |x -s| -r if x / ∈ X (7)
Based on this definition, we can compute the distance from s to X by averaging expression [START_REF] Augustin | Statistical inference[END_REF] on the N components:

d inf (s, X) = 1 N N i=1 d inf (s i , X i ). ( 8
)
This definition captures the notion of imprecision, but does not have enough separating power. As a long as a point is included in the interval, the distance will be zero. It means that an estimation consisting of intervals of the form (-∞, +∞), i.e. a completely vacuous estimation, would always have a distance of zero to any precise estimation, because it contains them all. This distance has exactly the opposite flaw than the Hausdorff distance: the wider the interval, the lower the distance. Thus an interval-valued estimation will always be considered as less distant from a reference value than an equivalent precise-valued estimation: ∀y ∈ X, d(s, y) ≥ d inf (s, X).

3 Weighted Variation of the Mid-Point Distance

Definition

A perfect candidate would be an extension which would not penalize imprecision when it conveys information, for example when it reflects coherently the variability of the quantities under consideration. This is in line with the guaranteed approach of [START_REF] Jaulin | Applied Interval Analysis[END_REF] and the confidence interval interpretation of imprecision [START_REF] Mauris | Possibility distributions: a unified representation of usual directprobability-based parameter estimation methods[END_REF]. The weighted variation of the mid-point distance is simply defined as the mean of deviations from the center of the intervals weighted by the inverse of the radius of the interval:

d w (s, X) = 1 N N i=1 r r i |s i -xi |, ( 9
) with r = 1 N N i=1 r i .
It can be seen as a mid-point distance computed on a space where the imprecision is uniform (usually unitary). It is then back-projected to the original space by multiplying the result by the mean radius of the intervals. The weights of the weighted extension are proportional to the inverses of the radii of the intervals. Thus for a certain index, the wider the radius, the lower will be the impact of an important deviation from the center of the interval. It is a very straightforward way of translating the idea that if an interval estimation suitably describes the variability of the reference, then it should be less penalized. Still, having wide radii will be penalized by increasing the mean radius r.

Obviously the WVD does not satisfy metric properties. The main problem is that if the vector X contains both intervals of strictly positive radii and degenerated intervals (i.e. singletons) then the WVD we propose is not defined. It thus does not formally stand as an extension of the L 1 distance. However we shall insist that the aim of this tool is to propose a practical solution to the unsolved problem of computing distances between points and intervals. When we restrict its use to cases where X contains only non-degenerated intervalswhich happens most of the time in practical casesthen the WVD fulfills its role in a way that is both very simple computationally and that makes sense from a formal point of view.

Formal Interpretation

Let us consider a regression problem where X X : Ω → R n denotes the vector of attributes and Y : Ω → R represents the response variable. Let us consider a sample of size N ((x 1 , y 1 ), . . . , (x N , y N )) . Let us consider a regression model f : R N → R and let ŷi = f (x i ) be the punctual estimation of y i based on it. The average distance d = 1

N |y iŷi | can be used to estimate the degree of goodness of our model, in terms of L 1 loss function. Let now X i,α = [ŷ ic α r i , ŷi + c α r i ] represent a prediction interval, with (exact) confidence level 1α, based on x i for the value of the response variable of another individual, randomly picked in the subpopulation of those whose vector of attributes coincides with the vector of attributes of the i th individual of our initial sample, x i . Due to the variability of this subpopulation, such a response value does not necessarily coincide with the observed y i . The length of prediction intervals in linear regression problems, under the usual normality assumptions, takes this form of c α r i . The average radius of the N prediction intervals calculated from the sample is c α r. Let us now consider another test sample of size N of pairs of the form (x i , y i ). An unbiased estimation of 1α based on this test sample is:

1 -α = 1 N #{i| |ŷ i -y i | ≤ c α r i } = 1 N #{i| |ŷ i -y i | r i ≤ c α }.
Let the random variable D denote the average L 1 distance between the true values y i and their estimations ŷi based on the x i 's. Let us prove that the WVD from the vector of responses y i to the sequence of confidence intervals X i,α is an unbiased estimation of the expectation of D. According to above assumptions, for every j = 1, . . . , N, the probability that D j is less than c α •r j is 1-α. Furthermore, as we have shown, the relative cumulative frequency 1 N #{i|

|y i -ŷi| ri ≤ c α }
is an unbiased estimation of 1α. In other words, F (c α ) = 1 N #{i|

Ê(D) = 1 N N i=1 |y i -ŷi | r i r.

Experiment

Experimental Procedure

This section aims at comparing the behavior of the proposed WVD of the L 1 distance to the four other extensions we have presented in Sect. 2. The experiment is based on a data set composed of K = 3900 8×8 subsampled patches extracted from a 2117×3006 high resolution image of the painting La Joconde by Leonardo Da Vinci. The reduction factor is set to 5. The subsampling procedures aims at simulating acquisitions of images of this painting using different imagers having the same numerical resolution but different point spread functions.

To do so, we used the so-called imprecise filtering sub-sampling method [START_REF] Strauss | Non-additive interval-valued f-transform[END_REF]. It consists of replacing the smoothing anti-aliasing kernel used to transform a high resolution image into a low resolution image by a capacity that represents a convex set of bell-shaped smoothing kernels. Filtering a patch with this method leads to an interval-valued subsampled patch that represent the convex set of all the patches that would have been obtained by subsampling the original patch with all the kernels belonging to the so-defined convex set (here a set of unimodal centered smooth kernels whose support is lower than 8). Practically, each interval-valued patch is composed of an upper patch P k and a lower patch P k (for the k th patch).

How Consistent Is the Weighted Variation of the Mid-Point

Distance?

This experiment aims at illustrating the fact that the extension we propose behaves consistently. Ideally an imprecise patch contains the information of several precise patches. Therefore a distance between imprecise and precise patches which wisely uses the information of imprecision would behave the same way as the L 1 distance would when used to compare precise patches drawn randomly inside the imprecise patches. The comparison is achieved by computing a ratio of imprecise-to-precise L 1 distances on various configuration of patches. If the distribution of this ratio is centered around 1 for a large number of samples, it means that the considered extension is consistent, because it reflects the information conveyed by the imprecise patches in a way that is consistent with how the L 1 distance would behave on a set of precise patches. Moreover, the lower the variance, the more consistent is the extension. The experiment is carried as follows. We divide the set of 3900 patches into two subsets of 1950 patches. The motivation for separating our data into two subsets is to compare patches having non-null distances (even if some patches bear some similarities).

We take the first subset (k = 1, . . . , 1950) to define references by considering the mid patches Pk = 1 2 (P k + P k ). We take the second subset (j = 1951, . . . , 3900) to define the imprecise estimations defined by their upper (P j ) and lower (P j ) patches. For each imprecise patch of the second subset, we draw 300 precise patches P n j (n = 1, . . . , 300) included in the imprecise patch [P j , P j ] (300 offers a good tradeoff between a statistically meaningful sampling of the interval and a reasonable computation time). We then compute, for each (k, j, n) ∈ [START_REF] Beer | Imprecise probabilities in engineering analyses[END_REF]1950] × [1951, 3900] × [1, 300], a ratio-distance which is defined as:

r e =
d e ( Pk , [P j , P j ])

d ( P n j , Pk ) , ( 10
)
where d is the (precise) L 1 distance between two patches and d e is one of the extensions of the L 1 distance (e ∈ {Hausdorff, infinimum, mid, weighted }. The distribution of the ratios for the different extensions are presented as a box-plot in Fig. 2. Table 1 shows the mean, median and standard deviation of the ratios distributions for each extension.

As might have been expected, the Hausdorff distance always leads to an over-evaluated distance, while, on contrary, the infimum distance leads to an under-evaluated distance. WVD and mid-point distance seem to provide distances that are consistent in that the means of their ratios are close to 1. For these extensions, comparing the imprecise patch with a precise one is statistically equivalent to comparing this precise patch with any patch contained in the imprecise patch. Although the ratios computed for the WVD have a slightly higher standard-deviation (which can be explained by the divergent behavior it has when some intervals tend to points) a Wilcoxon test applied on all the extensions showed the distribution of ratios of the WVD to be the closest to 1.

Conclusion

In this article we have introduced a new mathematical tool, the weighted variation of the mid-point distance, that allows to compare the performances of an interval-valued method with those of a conventional precise-valued method. We have presented its behavior through an experiment, where we compared it with other possible alternatives, namely the Hausdorff, mid-point and infimum distances. The WVD has some problems, such as its divergent behavior when some of the intervals radii tend to zero. However when considering a case of imprecise estimations where the imprecision has to be informative, meaning that the radius of the intervals should reflect the quality of the information provided, the new tool we proposed proved itself to have the best tradeoff between informativeness and consistency with the L 1 distance.

Fig. 1 .

 1 Fig. 1. Precise estimations (1,4) and imprecise estimation (2,3). Plain line is the reference, dotted lines are estimations.

  |y i -ŷi| ri ≤ c α } is an unbiased estimation of the cumulative probability F D j r j (c α ). Therefore the relative frequency p(x) = 1 N N i=1 1 |y i -ŷi | r i (x) is an unbiased estimation of the probability P Dj rj = x = P (D j = r j • x). Thus, the expectation of D j is unbiasedly estimated by Ê(D j ) = 1 N N i=1 |y i -ŷi| ri r j . Therefore, the expectation of D, E(D) = 1 N N j=1 E(D j ) can be unbiasedly estimated by:

Fig. 2 .

 2 Fig. 2. Box-plot results of the first experiment computed on the four extensions. Weighted extension is WVD.

Table 1 .

 1 Moments of the ratios.

		Hausdorff Infinimum Mid-points Weighted
	mean	1.54	0.503	0.95	0.99
	median	1.71	0.32	0.99	1.01
	standard deviation 0.36	0.27	0.14	0.16

The Hausdorff distance is also called the Pompeiu-Hausdorff distance. Here we simply refer to it as the Hausdorff distance.

We will use the same notation dH in order to denote the mapping defined on R × IR derived from dH : K(R) × K(R) → R. Obviously this new mapping does not satisfy metric properties. It is not even applied to objects of the same kind.

Acknowledgements. This work is partially supported by TIN2014-56967-R (Spanish Ministry of Science and Innovation), FC-15-GRUPIN14-073 (Regional Ministry of the Principality of Asturias) and Merimee exchange program.


Xiaomei Xu

Raphaël Rachedi

Maryline Foglino

Dr Emmanuel Talla
email: talla@imm.cnrs.fr

Amel Pr

Latifi
email: latifi@imm.cnrs.fr

Interaction network among factors involved 1 in heterocyst-patterning in cyanobacteria 2 3

Keywords: Cell differentiation, Cyanobacteria, Genomics, Pattern establishment, Protein-19 Protein interaction

The genetically regulated pattern of heterocyst formation in multicellular cyanobacteria represents the 31 simplest model to address how patterns emerge and are established, the signals that control them, and 32 the regulatory pathways that act downstream. Although numerous factors involved in this process have 33 been identified, the mechanisms of action of many of them remain largely unknown. The aim of this 34 study was to identify specific relationships between 14 factors required for cell differentiation and 35 pattern formation by exploring their putative physical interactions in the cyanobacterium model Nostoc 36 sp. PCC 7120 and by probing their evolutionary conservation and distribution across the cyanobacterial 37 phylum. A bacterial two-hybrid assay indicated that 10 of the 14 factors studied here are engaged in 38 more than one protein-protein interaction. The transcriptional regulator PatB was central in this network 39 as it showed the highest number of binary interactions. A phylum-wide genomic survey of the 40 distribution of these factors in cyanobacteria showed that they are all highly conserved in the genomes 41 of heterocyst-forming strains, with the PatN protein being almost restricted to this clade. Interestingly, 42 eight of the factors that were shown to be capable of protein interactions were identified as key elements 43 in the evolutionary genomics analysis. These data suggest that a network of 12 proteins may play a 44 crucial role in heterocyst development and patterning. Unravelling the physical and functional 45 interactions between these factors during heterocyst development will certainly shed light on the 46 mechanisms underlying pattern establishment in cyanobacteria.

Introduction 54

Spatial patterns of cellular differentiation occur in many developmental organisms across the tree of life, 55 and in most cases, key aspects of the molecular interactions that are involved in patterning are poorly 56 understood. In this work, we questioned the possible involvement of protein-protein interactions in 57 patterning. Our study model is the development of heterocysts in multicellular and diazotrophic 58 cyanobacteria. Indeed, in prokaryotes, heterocyst formation in multicellular cyanobacteria represents 59 the simplest model to address one-dimensional patterns based on repeated peaks of diffusible 60 morphogens. These peaks are generated by the activity of a transcriptional regulator that activates the 61 expression of the gene encoding its own inhibitor (the diffusible morphogen); thus, this loop has been 62

called the reaction-diffusion model (Fig. 1a) [START_REF] Turing | The chemical basis of morphogenesis[END_REF]. 63

Heterocysts are nearly anoxic cells that form semi-regularly along the filaments and function to host the 64 oxygen-sensitive nitrogenase. Therefore, the developmental program in these heterocyst-forming 65 bacteria allows the occurrence of two antagonistic metabolic activities (i.e., oxygenic photosynthesis 66 and nitrogen fixation) in the same organism [START_REF] Meeks | Regulation of cellular differentiation in filamentous cyanobacteria in free-659 living and plant-associated symbiotic growth states[END_REF][START_REF] Kumar | Cyanobacterial heterocysts[END_REF]Herrero et al. 67 2016). Much of the knowledge about cell differentiation in cyanobacteria stems from studies of the 68 model bacterium Anabaena/Nostoc sp. PCC 7120 (termed Nostoc hereafter). Combined nitrogen 69 starvation triggers the activation of the expression of the hetR gene, which encodes the master and 70 specific transcriptional activator of cell differentiation [START_REF] Buikema | Characterization of a gene controlling heterocyst differentiation in 559 the cyanobacterium Anabaena 7120[END_REF]. HetR is essential 71 for cell differentiation, as its absence abolishes differentiation and its overexpression induces the 72 differentiation of multiple contiguous heterocysts (MCH) under combined nitrogen starvation and 73 allows differentiation even under non-permissive conditions [START_REF] Buikema | Characterization of a gene controlling heterocyst differentiation in 559 the cyanobacterium Anabaena 7120[END_REF]. The 74 regulation of hetR transcription is complex and involves the global regulator NtcA and a positive 75 autoregulation loop mediated by HetR itself [START_REF] Black | Spatial expression and autoregulation of hetR, a gene involved in the 552 control of heterocyst development in Anabaena[END_REF][START_REF] Frias | Requirement of the regulatory protein NtcA for the expression of 583 nitrogen assimilation and heterocyst development genes in the cyanobacterium Anabaena 584 sp. PCC 7120[END_REF]Rajagopalan and 76 Callahan 2010). The inhibition of HetR activity in the neighboring cells of the heterocysts initiates and 77 maintains the pattern of differentiation throughout the growth of the filaments. Three genes (patS, patX, 78 and hetN), the expression of which is induced in developing cells, negatively control this process (Black 79 and Wolk 1994;[START_REF] Yoon | Heterocyst pattern formation controlled by a diffusible peptide[END_REF][START_REF] Yoon | PatS and products of nitrogen fixation control heterocyst pattern[END_REF][START_REF] Higa | The RGSGR 597 amino acid motif of the intercellular signalling protein, HetN, is required for patterning of 598 heterocysts in Anabaena sp. strain PCC 7120[END_REF]Corrales-Guerrero 80 et al. 2014b;[START_REF] Elhai | Ancient association of cyanobacterial multicellularity with the regulator 573 HetR and an RGSGR pentapeptide-containing protein (PatX)[END_REF][START_REF] Khudyakov | Inactivation of Three RG(S/T)GR Pentapeptide-Containing 618 Negative Regulators of HetR Results in Lethal Differentiation of Anabaena PCC 7120[END_REF]. Deletion of the patS gene also leads 81 to the formation of multiple contiguous heterocysts (Fig. 1b). A common feature of the PatS, PatX, and 82

HetN proteins is that they contain a conserved hexapeptide (ERGSGR for PatS and HetN,and HRGTGR 83 for PatX). The RGSGR peptide inactivates the DNA-binding activity of HetR [START_REF] Huang | HetR homodimer is a DNA-binding protein required for heterocyst 603 differentiation, and the DNA-binding activity is inhibited by PatS[END_REF]84 Feldmann et al. 2011). These peptides are considered as diffusible morphogens that ensure patterning 85 by lateral inhibition [START_REF] Corrales-Guerrero | Functional dissection and evidence for 564 intercellular transfer of the heterocyst-differentiation PatS morphogen[END_REF]. Recently, HetR was shown to interact with a 86 pentapeptide-repeat-containing protein called HetL at the same interface used by the ERGSGR peptide 87 [START_REF] Xu | HetL, HetR and PatS form a reaction-711 diffusion system to control pattern formation in the cyanobacterium nostoc PCC 7120[END_REF]. Moreover, the transcription of hetL, which is a heterocyst-specific gene, is activated 88 by HetR [START_REF] Xu | HetL, HetR and PatS form a reaction-711 diffusion system to control pattern formation in the cyanobacterium nostoc PCC 7120[END_REF]. hetL overexpression suppresses the inhibitory effect of PatS (Liu and Golden 89 2002) and the interaction of HetL with HetR provides protection against PatS [START_REF] Xu | HetL, HetR and PatS form a reaction-711 diffusion system to control pattern formation in the cyanobacterium nostoc PCC 7120[END_REF], thus 90 explaining how the developing cell becomes immune to PatS/HetN/PatX inhibition in the heterocyst 91 (Fig. 1a). 92

Correct patterning requires the function of several genes that have been termed pat genes (Fig. 1b). The 93 phenotype of mutants of these genes provides an important indication of their role and establishes a 94 starting point for the investigation of how the pattern is established. These genes include patA and patL, 95 the mutation of which results in the formation of rare heterocysts only at the ends of the filaments (Liang 96 et al. 1992;[START_REF] Liu | Mutations in genes patA and patL of Anabaena sp. strain PCC 7120 result in 653 similar phenotypes, and the proteins encoded by those genes may interact[END_REF]. PatA is a two-domain protein with a C-terminal region similar to the 97 CheY-receiver domain and a conserved N-terminal domain called PATAN [START_REF] Makarova | Cyanobacterial response regulator PatA 656 contains a conserved N-terminal domain (PATAN) with an alpha-helical insertion[END_REF]. PatL 98 is a pentapeptide-repeat-containing protein with two transmembrane domains [START_REF] Liu | Mutations in genes patA and patL of Anabaena sp. strain PCC 7120 result in 653 similar phenotypes, and the proteins encoded by those genes may interact[END_REF]. In 99 addition to sharing similar mutant phenotypes, the PatA and PatL proteins interact, as shown by yeast-100 two-hybrid assays [START_REF] Liu | Mutations in genes patA and patL of Anabaena sp. strain PCC 7120 result in 653 similar phenotypes, and the proteins encoded by those genes may interact[END_REF]. A cluster of three genes (hetZ, patU3, and patU5) that are 101 expressed in the heterocysts and encode proteins of unknown function modulates heterocyst formation 102 and pattern establishment [START_REF] Zhang | A gene cluster that regulates 721 both heterocyst differentiation and pattern formation in Anabaena sp. strain PCC 7120[END_REF]. Inactivation of patU3 leads to an MCH phenotype and 103 restores heterocyst formation in a patA mutant background [START_REF] Zhang | A gene cluster that regulates 721 both heterocyst differentiation and pattern formation in Anabaena sp. strain PCC 7120[END_REF]. A hetZ mutant of Nostoc 104 is unable to differentiate heterocysts [START_REF] Zhang | A gene cluster that regulates 721 both heterocyst differentiation and pattern formation in Anabaena sp. strain PCC 7120[END_REF][START_REF] Videau | The hetZ gene indirectly regulates heterocyst development at the level of pattern 699 formation in Anabaena sp. strain PCC 7120[END_REF]. HetZ action has been 105

proposed to be involved in its interaction with PatU3 and, to a lesser extent, with HetR (Videau et al. 106 2018;[START_REF] Zhang | Functional Overlap of hetP and hetZ in Regulation of 719 Heterocyst Differentiation in Anabaena sp. Strain PCC 7120[END_REF][START_REF] Du | Expression from DIF1-motif 570 promoters of hetR and patS is dependent on HetZ and modulated by PatU3 during 571 heterocyst differentiation[END_REF]). In N. punctiforme, the deletion of the patN gene, which encodes 107 a putative membrane protein, increases heterocyst frequency and decreases the number of vegetative 108 cells between two heterocysts, thus leading to a multiple singular heterocyst (MSH) phenotype (Risser 109 et al. 2012) (Fig. 1b). A similar phenotype has been observed in a patN mutant of Nostoc (Abed et al. 110 2009). PatN was proposed to modulate the competency of a vegetative cell to differentiate into a 111 heterocyst [START_REF] Risser | Biased inheritance of the protein PatN frees vegetative cells to 676 initiate patterned heterocyst differentiation[END_REF]; therefore, it also affects patterning. Interestingly, the MSH phenotype 112 has been observed in a mutant of the patB (cnfR) gene, the product of which is predicted to be a 113 transcriptional regulator with a DNA-binding helix-turn-helix at its C-terminal extremity and two 114 putative ferredoxin domains at its N terminus [START_REF] Jones | Heterocyst-specific expression of patB, a gene required 609 for nitrogen fixation in Anabaena sp. strain PCC 7120[END_REF]. In Anabaena variabilis, the homolog 115 of PatB (CnfR1) activates the transcription of the heterocyst specific nitrogenase-encoding genes (Pratte 116 and Thiel 2016). Furthermore, in the non-heterocyst-forming cyanobacterium Leptolyngbya boryana, 117

CnfR is essential for the transcription of the nif genes and for the diazotrophic growth under micro-oxic 118 conditions [START_REF] Tsujimoto | Transcriptional regulators ChlR and CnfR are essential for 686 diazotrophic growth in nonheterocystous cyanobacteria[END_REF]. In Nostoc, the transcription of the nif genes has been shown to be 119 abolished in a patB (cnfR) mutant [START_REF] Kurio | The CRP-family transcriptional regulator 636 DevH regulates expression of heterocyst-specific genes at the later stage of differentiation in 637 the cyanobacterium Anabaena sp. strain PCC 7120[END_REF]) and the DNA-binding sites of this regulator has 120 been predicted in the promoter of the nif operon [START_REF] Brenes-Alvarez | 556 Elements of the heterocyst-specific transcriptome unravelled by co-expression analysis in 557 Nostoc sp. PCC 7120[END_REF]. The function of CnfR/PatB 121 as activator of the nif genes seems therefore conserved in cyanobacteria, but the manner in which PatB 122 exerts this negative effect on heterocyst spacing in Nostoc is unknown. Two other pat genes, patC and 123 patD, both encoding proteins of unknown functions, seem to negatively regulate patterning, as their 124 deletion increases heterocyst frequency (Corrales-Guerrero et al. 2014a;[START_REF] Wang | patD, a Gene 704 Regulated by NtcA, Is Involved in the Optimization of Heterocyst Frequency in the 705 Cyanobacterium Anabaena sp. Strain PCC 7120[END_REF]. 125

In addition to the Pat proteins, other factors involved in the early steps of the differentiation process can 126 affect pattern formation. The hetC gene, which encodes a putative bacterial ABC protein transporter, is 127 involved in the regulation of the early steps of heterocyst differentiation [START_REF] Khudyakov | hetC, a gene coding for a protein similar to bacterial ABC protein 621 exporters, is involved in early regulation of heterocyst differentiation in Anabaena sp. strain 622 PCC 7120[END_REF]128 Corrales-Guerrero et al. 2014a;[START_REF] Videau | ABC Transporter Required for Intercellular Transfer 693 of Developmental Signals in a Heterocystous Cyanobacterium[END_REF]. hetC inactivation mutants exhibit a pattern of spaced 129 cells that are still able to divide [START_REF] Khudyakov | hetC, a gene coding for a protein similar to bacterial ABC protein 621 exporters, is involved in early regulation of heterocyst differentiation in Anabaena sp. strain 622 PCC 7120[END_REF][START_REF] Xu | Role for hetC in the transition to a nondividing state during heterocyst 713 differentiation in Anabaena sp[END_REF]. These cells are likely 130 to represent an early stage of heterocyst differentiation; thus, HetC has been suggested to be required 131 for the transition of heterocysts to their terminal non-dividing state [START_REF] Xu | Role for hetC in the transition to a nondividing state during heterocyst 713 differentiation in Anabaena sp[END_REF]. hetF, which 132 encodes a putative protease, is another gene that is essential for heterocyst formation, as a hetF mutant 133 does not form heterocysts and has altered morphology, whereas hetF overexpression results in an MCH 134 phenotype [START_REF] Wong | The hetF gene product is essential to heterocyst differentiation and 707 affects HetR function in the cyanobacterium Nostoc punctiforme[END_REF][START_REF] Risser | HetF and PatA control levels of HetR in Anabaena sp. strain PCC 7120[END_REF]. 135 Shortly after the initiation step and pattern establishment, the decision to differentiate is irreversibly 136 fixed, which is known as commitment. The HetP protein (encoded by the alr2818 gene) and its 137 homologs (encoded by the asl1930, alr3234, and alr2902 genes) have been proposed as important 138 players in the regulation of this step [START_REF] Videau | The 695 heterocyst regulatory protein HetP and its homologs modulate heterocyst commitment in 696 Anabaena sp. strain PCC 7120[END_REF]. The inactivation of the hetP gene affects the 139 differentiation process downstream of initiation and patterning, as the mutant develops heterocysts, 140 albeit at a lower frequency and at a later time than the wild type strain [START_REF] Fernandez-Pinas | A third genetic locus required for the formation of 578 heterocysts in Anabaena sp. strain PCC 7120[END_REF]141 Higa and Callahan 2010). Unlike hetP, inactivation of the genes encoding the 3 HetP homologs does 142 not impair cell differentiation, indicating that they are not required for heterocyst formation or patterning 143 (40). Overexpression of hetP bypasses the absence of HetZ and vice versa, which indicates a functional 144 overlap between these factors [START_REF] Zhang | Functional Overlap of hetP and hetZ in Regulation of 719 Heterocyst Differentiation in Anabaena sp. Strain PCC 7120[END_REF]. This observation highlights the connection between 145 patterning and commitment and places HetZ at the intersection between the two steps. 146

In summary, pattern establishment has clearly been shown to be genetically controlled, and the 147 involvement of several factors in this process has been demonstrated (Fig. 1b). However, their function 148 and their potential molecular interactions remain poorly known. Because protein-protein interaction 149 approaches have been shown to be useful for identifying protein functions and understanding several 150 cellular pathways, we used a bacterial two-hybrid assay (BACTH) to address the possible binary 151 interactions among the 14 factors involved in the early, commitment, and patterning stages in Nostoc. 152

In addition, we performed a phylum-wide analysis of the conservation of these factors in 212 153 cyanobacterial genomes. A comparative genomics analysis of the genes encoding these factors across 154 cyanobacteria supported the contention that they likely were present in the common ancestor of 155 cyanobacteria and that they evolved through gene loss in clades of cyanobacteria that do not form 156 heterocysts. Taken together, our data suggests that a set of 12 proteins (DevH,HetC,HetF',HetL,HetP,157 HetR,HetZ,PatA,PatB,PatL,PatN,and PatU3), which encoding genes are conserved in heterocyst-158 forming cyanobacterial genomes, are associated through multiple cross-talk interactions. Our findings 159 provide a framework for a better understanding of heterocyst pattern establishment. 160

Materials and methods

161

Bacterial strains and growth media 162

The E. coli K-12 strain DH5 alpha (Invitrogen) was used in all the cloning steps. BACTH 163 complementation assays were carried out with the E. coli cya strain BTH101 araD139,164 galE15,galK16,rpsL1 (Strr ),hsdR2,mcrA1,mcrB1) [START_REF] Karimova | A bacterial two-hybrid system based on a 614 reconstituted signal transduction pathway[END_REF]. Bacteria were grown at 165 30°C in Luria-Bertani (LB) broth supplemented with ampicillin at 100 µg/ml and kanamycin at 50 µg/ml 166 when needed. Screening for the ability to complement the cya mutation was performed on LB agar plates 167 supplemented with 5-bromo-4-chloro-3-indolyl--D-galactoside (X-Gal; 40 µg/ml), isopropyl--D-168 galactopyranoside (IPTG; 0.5 mM), and the appropriate antibiotics. The E. coli K-12 Lemo21(DE3) 169 strain (Biolabs) was used for protein production and purification. Spectinomycin (50 µg/ml) was used 170 for the selection of pCDF Duet plasmid (Novagen) and derivatives. Kanamycin (100 µg/ml) for the 171 selection of pLIC07 plasmid [START_REF] Michel-Souzy | Direct 662 interactions between the secreted effector and the T2SS components GspL and GspM reveal 663 a new effector-sensing step during type 2 secretion[END_REF]) and derivatives. 172

Plasmid constructions 173

For BACTH analysis: the genes involved in heterocyst differentiation and patterning were amplified 174 from genomic DNA of Nostoc using the primers listed in Table S1. The recombinant plasmids were 175 constructed by inserting gene sequences of interest in pKT25 (T25 at N-terminus), pKNT25 (T25 at the 176 C-terminus), pUT18 (T18 at the C-terminus) and PUT18C (T18 at the N-terminus) (Karimova et al. 177 1998). All the open reading frames (ORFs) of the analyzed genes, except patN, were inserted between 178

Pst1 and EcoRI restriction sites of the BACTH plasmids. patN ORF was inserted between PstI and Kpn1 179 restriction sites. 180

For pull-down assays: the hetR, hetL or hetF' sequences fused at their 3' extremities to the sequence 181 coding the 89-106 amino acids of the human influenza hemagglutinin protein (HA-tag) were synthesized 182 by Genewiz (https://www.genewiz.com/en-GB/). The obtained genes were inserted between NdeI and 183 XhoI restriction sites of the pCDF Duet plasmid. The patB coding sequence was fused to a 6xHis tag at 184 its 5'extremity. For this, the patB ORF was amplified by PCR using pLIC07-patB-Fwd and pLIC07-185 patB-Rev primers (Table S1) and inserted between the BsaI restriction sites of the pLIC07 plasmid using 186 the SLIC method [START_REF] Li | SLIC: a method for sequence-and ligation-independent cloning[END_REF]. 187

All the constructs were analyzed by sequencing. The plasmids used in this study are listed in Table S2. 188

SDS-PAGE and immunodetection 189

Proteins were separated on 4-20% polyacrylamide gradient gels (NuSep) and transferred onto 190 nitrocellulose membranes using a semidry blotting apparatus. Membranes were blocked with 3% bovine 191 serum albumin (BSA) in TBST (10 mM Tris, 150 mM NaCl, 0.05% Tween 20). Membranes were 192 incubated with monoclonal anti-HA (Enzo) or anti-His (Qiagen) antibodies, diluted at 1:5000 and 193 1:2000, respectively in TBST, 3% BSA. After three-step washes in TBST, the membranes were 194 incubated with peroxidase-coupled anti-mouse antibodies (1:5000, Promega). The detection was 195 performed using the enhanced chemiluminescence reagents (Cytiva). The membranes were scanned 196 using ImageQuant TL analysis software (GE Healthcare). 197

Coproduction of proteins in E. coli and Nickel-affinity chromatography 198

Lemo21 (DE3) competent cells were co-transformed with pCDF Duet, pLIC07, and their derivatives. 199

The recombinant clones were grown until reaching an optical density at 600 nm of 0.5 at 37 °C in LB 200 medium with appropriate antibiotics (50 µg/ml spectinomycin, 100 µg/ml kanamycin). Induction was 201 performed with 0.5 mM IPTG for 3 hours at 30 °C. Bacteria were collected by centrifugation and broken 202 using a French press (10,000 p.s.i.) in cold buffer (100 mM Tris-HCl, pH 7.3, 400 mM NaCl, 0.5 mg/ml 203 lysozyme, 20 µg/ml DNase, 1X cOmplete™ EDTA-free protease inhibitor cocktail (Roche). The lysate 204 was cleared by ultracentrifugation. The cleared lysate was loaded onto a 1-ml HiFliQ Ni-NTA FPLC 205 column (Generon) using an ÄKTA Prime apparatus. After a washing step (10 ml of buffer: 100 mM 206 Tris-HCl, pH 7.3, 400 mM NaCl), the immobilized proteins were eluted in buffer: 100 mM Tris-HCl, 207 pH 7.3, 400 mM NaCl, 500 mM imidazol. The loaded flow-through, and elution fractions were analyzed 208 by SDS-PAGE and immunodetection. 209 210

Bacterial two hybrid assays 211

Bacterial two-hybrid assays were performed following the procedure described by Karimova et al 212 (1998) [START_REF] Karimova | A bacterial two-hybrid system based on a 614 reconstituted signal transduction pathway[END_REF]. Briefly, after co-transforming the BTH101 strain with the two plasmids 213 expressing the T18-and T25-fusions, LB plates containing ampicillin and kanamycin were incubated 214 at 30° C for 2 days. For each assay, 10 independent colonies were inoculated in 3 ml of LB medium 215 supplemented with ampicillin, kanamycin, and 0.5 mM IPTG and incubated at 30 °C overnight. 216 ß-galactosidase activity was then determined as previously described [START_REF] Zubay | Detection and isolation of the repressor protein for 724 the tryptophan operon of Escherichia coli[END_REF].

217

Bioinformatics analysis 218

Datasets 219

The genome data of 212 cyanobacterial strains available in March 2019 were downloaded from NCBI 220 ftp site (ftp:/ftp.ncbi.nih.gov/genomes/) and constituted the primary data source (Table S3, Sheet 1). 221

This data source includes all complete genomes in the highest levels of assembly, their refseq category 222 (reference or representative genome, or none) as well as their annotation features and taxonomy lineages. 223

The ability of each strain to form heterocysts was provided by the literature [START_REF] Rippka | Generic assignments, strain stories and 672 properties of pure cultures of cyanobacteria[END_REF]Henson 224 B.J 2004;[START_REF] Shih | Improving the coverage of the 683 cyanobacterial phylum using diversity-driven genome sequencing[END_REF][START_REF] Komárek | Taxonomic classification of cyanoprokaryotes 627 (cyanobacterial genera) 2014, using a polyphasic approach[END_REF][START_REF] Hammerschmidt | The Order of Trait 588 Emergence in the Evolution of Cyanobacterial Multicellularity[END_REF] S4. 230

Identification of homologs 231

HMMER package [START_REF] Mistry | Challenges in homology search: HMMER3 and 665 convergent evolution of coiled-coil regions[END_REF] and HMM domain profiles were first used to establish the list of 232 seed functional domains associated to reference seed proteins (Table S4). Alignments with a score 233 higher than the Pfam trusted thresholds were considered as significant seed domains. Then, HMMER3 234 [START_REF] Mistry | Challenges in homology search: HMMER3 and 665 convergent evolution of coiled-coil regions[END_REF]) and self-written Perl scripts were used to search for protein homologs (with reference 235 seed domains) in the complete genomes. For each reference protein, the presence of one or more seed 236 domains was a requisite to identify putative homologs (Table S4). As described above, alignments with 237 a score higher than the Pfam thrusted thresholds were also considered significant. Note that HetP protein 238 and its 3 homologs in Nostoc were considered as a unique seed protein HetP since they share an identical 239 functional domain (Table S4). In contrary, the two seed proteins HetF and HetF' were analyzed 240 independently since their domain organizations are different (Table S4). For each reference protein, 241 putative homologs were subsequently analyzed with the same software in order to determine the 242 potential presence of additional functional domains. In-house Perl scripts were used to define the domain 243 organization of each homolog. Putative homologs were then selected for the strictly presence of seed 244 domains (i.e., in the right order and without additional domain), and therefore led to the set of Hmm-245 homologs. To select specific homologs, additional criteria were added for the selection of HetC, HetL 246 and PatL. For this, we considered that to be a HetC homolog, the gene must be located in the same 247 genomic context as hetP (based on the co-localization of these two genes in the genome of Nostoc, see 248

Table S4) with a maximal distance of 5 genes upstream or downstream from hetP. HetL and PatL 249 homologs were limited to protein homologs sharing similar features with their seed sequences [a protein 250 length between 0.8 and 1.25 times when compared to the reference seed proteins (Table S4) and a 251 number of repeats fixed to 3-5, and 4-6, respectively]. Since PatN and PatU3 do not exhibit a known 252 functional domain, a Blast [START_REF] Altschul | Basic local alignment search tool[END_REF]) analysis search was performed using all reference seed 253 proteins as queries and an E-value of 10 -5 . Blast alignments with a MinLrap [Alignment size / Minimun 254 size (Query, Subject)] >= 0.8 and MaxLrap [Alignment size / Minimun size (Query, Subject)] >= 0.8 255 were considered as significant and constituted the set of Blast-homologs. For PatU3 and PatN, for which 256 no seed functional domain was found, MinLrap and MaxLrap were set to 0.7. Hmm-homolog and Blast-257 homolog sets were finally combined in order to yield the final set of homologs for each reference protein. 258The location of transmembrane domains in the proteins was performed using the TMHMM software 259 [START_REF] Krogh | Predicting transmembrane protein 631 topology with a hidden Markov model: application to complete genomes[END_REF]. 260

Statistical analysis. 261

-galactosidase activities: the values presented are means of 3 independent assays containing each 10 262 independent colonies. The results are expressed as the means ± standard variation. 263

Genomics analysis: one-sample t-test as well as boxplot, heatmap and distribution graphs were done 264 using the R (R Core Team, 2018). For boxplot analysis, taxonomic orders with low number of genomes 265 (i.e., Chroococcidiopsidales, Gloeobacterales, Gloeoemargaritales, Pleurocapsales, Spirulinales, and 266 unclassified Cyanobacteria) were fused into a single group (named Others) before computation (Table 267 S3, Sheet 4). 268 269

Results

270

To characterize the physical interactions between proteins that are involved in heterocyst differentiation 271 and pattern establishment in Nostoc, 13 full-length proteins (i.e., DevH, HetC, HetF, HetL, HetP, HetR, 272 HetZ, PatA, PatB, PatD, PatL, PatN, and PatU3) were tested systematically for pairwise interactions 273 using BACTH. We also included the All1730 ORF, which potentially encodes a HetF homolog (named 274 HetF') and the expression of which (i.e., all1730) has been shown to be induced in the same manner as 275

HetF after combined nitrogen starvation, as demonstrated by RNA-deep-sequencing experiments 276 [START_REF] Flaherty | Directional RNA deep sequencing sheds 580 new light on the transcriptional response of Anabaena sp. strain PCC 7120 to combined-581 nitrogen deprivation[END_REF]. 277

As for the BACTH assay the T18 and T25 domains must be in the cytosol, we first determined the 278 functional and structural domain composition of the studied proteins, as described in the Materials and 279

Methods. This allowed the definition of a subset of membrane proteins (HetC, HetF, PatL, and PatN) 280 and cytosolic factors (DevH, HetF', HetL, HetP, HetR, HetZ, PatA PatB, PatD, and PatU3) (Fig. 2, 281 Table S4). Each protein was fused to a part of the catalytic domain of the chimeric adenylate cyclase 282 (T25 or T18) of Bordetella pertussis, at the C and N termini, when possible. According to the predicted 283 topology of the membrane proteins, insertions were made at one terminus to have the T18 and T25 284 domains in the cytosol. Therefore, membrane proteins (HetC, HetF, PatL, and PatN) were fused to T18 285 and T25 at the extremity of their domain, which was predicted to be in the cytoplasm (Fig. 2). Direct 286 interactions were screened by transforming all possible pairwise combinations into the E. coli BTH101 287 strain (500 in total). Protein interactions were probed as explained in the Materials and Methods section. 288

In addition to the negative T25-T18 negative-control combination, all constructs were tested using a 289 relevant control vector (as an example: T25-HetR+ NT18 and NT25 + T18-HetR served as negative 290 controls for HetR) and only results above the control values were considered as being positive in the 291 remainder of the study (Fig. 3 and4). This control is important for the elimination of self-activation 292 values. In addition to the BACTH T18-Zip and T25-Zip positive control, the well-established 293 interactions between HetR subunits and HetL subunits, and between HetR and HetL, served as positive 294 internal controls. The pairwise combinations that yielded positive values are displayed in histograms 295 (Fig. 3 and4), whereas negative results are grouped in Table 1. 296

Network analysis of early factors 297

The major regulator of heterocyst differentiation, HetR, displays a specific fold that has been proposed 298 to allow its interaction with other factors [START_REF] Kim | Structure of transcription factor HetR required for heterocyst 625 differentiation in cyanobacteria[END_REF][START_REF] Hu | Structural insights into HetR-601 PatS interaction involved in cyanobacterial pattern formation[END_REF]. In fact, HetR has been shown 299 to interact with homologs of HetP [START_REF] Videau | The 695 heterocyst regulatory protein HetP and its homologs modulate heterocyst commitment in 696 Anabaena sp. strain PCC 7120[END_REF], with the serine/threonine kinase Pkn22 (Roumezi 300 et al. 2019), and with HetL [START_REF] Xu | HetL, HetR and PatS form a reaction-711 diffusion system to control pattern formation in the cyanobacterium nostoc PCC 7120[END_REF]. The potential interaction between HetR and other factors 301 acting early after the initiation step and during pattern establishment was analyzed. Among the factors 302 that are induced early, HetR, HetL, and HetF' were found to interact with the PatB transcription factor, 303 as the -galactosidase activities obtained were 10-to 40-fold higher than those of the negative controls 304 (Fig. 3a,b). In addition, interactions between HetL (Fig. 3b) and HetF' and multimerization of HetF' 305 (Fig. 3a) were observed. A weaker interaction (80 -galactosidase units) between HetF' and DevH was 306 also identified (Fig. 3a). Because HetF' has not been studied to date, these results suggest the possible 307 implication of this protein in the differentiation process, which is worth exploring. Among this group of 308 factors, the assays involving HetC and HetF all yielded negative results (Table 1) 309

Network analysis of commitment and patterning factors 310

Among the three factors (HetP, HetZ, and PatU3) that formed the "commitment group", a moderate but 311 significant interaction was observed between HetP and PatB (Fig. 4a). The highest -galactosidase 312 values (400 units) were obtained for the HetP-T25/PatB-T18 combination. In addition, BACTH data 313 confirmed the previously published interaction between HetZ and PatU3 [START_REF] Du | Expression from DIF1-motif 570 promoters of hetR and patS is dependent on HetZ and modulated by PatU3 during 571 heterocyst differentiation[END_REF], as two 314 combinations (T25-HetZ/PatU3-T18 and PatU3-T25/T18-HetZ) showed average -galactosidase 315 activities of 300 units (Fig. 4a). 316

The patterning factors consist of PatA, PatD, PatB, PatL, and PatN. Our analysis also included DevH in 317 this group, because it is required for patB transcriptional expression [START_REF] Kurio | The CRP-family transcriptional regulator 636 DevH regulates expression of heterocyst-specific genes at the later stage of differentiation in 637 the cyanobacterium Anabaena sp. strain PCC 7120[END_REF]). In addition to 318 the interactions with HetR, HetL, HetF', and HetP mentioned above, PatB formed multimers and 319 interacted with PatL in the PatL-T25/PatB-T18 combination (Fig. 4b). Moreover, PatA interacted with 320 PatL, which is in agreement with published data obtained using the yeast-two-hybrid approach (Liu and 321 Wolk 2011). PatD and PatN did not show any interaction in our assays (Table 1). 322

The patterning protein network of Nostoc as defined by BACTH 323 Among the 14 factors analyzed here, 10 were engaged in protein-protein interactions. In total, 16 324 interactions were detected, including self-associations (i.e. multimerization) (Fig. 5). Five proteins 325 (HetR, HetL, HetF', PatB, and PatL) were able to associate with multiple partners (at least two). The 326 interactions identified form a network that connects the elements involved in the initiation stage of 327 differentiation, commitment, and patterning. As expected, early and commitment factors mostly interact 328 among themselves, whereas the patterning factors, in particular PatB (with five interactions), appear to 329 be the central network and, therefore, serve as the master link between the three developmental stages 330 (Fig. 5). In addition, HetL and HetF' displayed three interactions each, suggesting that these proteins 331 are also important key players in the interaction network. 332

Validation of the interaction between PatB and the early factors 333

To validate our BACTH procedure, we sought to confirm some of the protein-protein interactions 334 described above by another approach. Since PatB appears to be central in the network, we decided to 335 analyze its ability to interact with the early factors by pull-down assays. PatB was produced as a His6-336 tagged fusion and its capacity to bind nickel column was verified by affinity chromatography 337 purification (Fig. S1a). HetR, HetL and HetF' fused to the HA-tag were produced as shown by the 338 immunoblots of Fig. S1b. To proceed with the pull-down assay, His6-PatB and each of the HA-tagged 339 proteins were coproduced in E. coli and protein purifications were performed using Nickel-affinity 340 chromatography. The analysis of the elution fractions by SDS-PAGE followed by immunoblotting with 341 anti-His6 and anti HA antibodies showed a specific co-elution of HetR, HetL and HetF' with PatB (Fig. 342 6, left panels). To control for non-specific binding of the HA-tagged proteins, HetR-HA, HetL-HA and 343

HetF'-HA produced in the absence of His6-PatB were purified following the same procedure. The data 344 presented in Fig. 6, right panels show that none of them was detected in the elution fractions, indicating 345 that they were unable to bind to the column. Altogether these results reveal direct interactions between 346

PatB and the three early factors and thus confirm the data obtained by BACTH. 347

Evolutionary conservation of genes underlying heterocyst formation and patterning 348

To address the general relevance in cyanobacteria of the interaction network described above, a large-349 scale (212 complete genomes in total) in silico identification and analysis of the set of the 14 protein 350 factors studied was performed (see Materials and methods and Table S4). Note that PatS, PatX and PatC 351 proteins were not included in our in silico study. The reasons are that these proteins are too short, do not 352 contain any known functional domain and in particular, patC was only present in Nostoc PCC7120 and 353

Trichormus variabilis NIES-23 genomes. These proteins are thus not compatible with significant 354 sequence comparisons based on functional domains search. 355

Interestingly, this genome survey detected 63 (for PatN) to 779 (for HetL) homologs within the genomes 356 (Table S3, Sheets 2 and 3). This large variation in the total number of homologs for each reference 357 protein reflects the wide variation in the number of copies of each gene per genome. In fact, the 358 distribution of the number of copies within a genome (Table S3, Sheet 3; Fig. S2) showed that: (i) PatN, 359 HetP, PatB, HetF', HetR, PatU3, HetZ, and HetF did not have genes encoding the corresponding 360 homologs in a high number of genomes (80); (ii) when present, they were on average present in ~1 361 copy (for PatN, PatB, HetF', HetR, and PatU3) or ~2 copies (for HetZ and HetF); and (iii) HetL group 362 of homologs displayed up to 4 copies per genome on average, with 13 copies encoding orthologs for 363 HetL in Microcoleus sp. PCC 7113. Moreover, HetL, which is encoded by a gene that is present in 364 multiple copies within most of the cyanobacterial genomes, harbored a low copy number of homologs 365 in Prochloraceae, with no copies detected in some Prochlorococcus marinus strains (e.g., P. marinus 366 str. NATL2A, P. marinus str. MIT9515, and P. marinus str. MIT 9211). Interestingly, genes encoding 367 HetR homologs were present in 2 or 3 copies in some genomes (e.g., Nodosilinea nodulosa PCC 7104 368

or Prochlorothrix hollandica PCC 9006, respectively). More surprisingly, most of the Nostococaceae 369 and Aphanizomenonaceae strains displayed homolog(s) for each group of reference proteins, except for 370 few organisms (e.g., Cylindrospermopsis raciborskii CS-505 and Dolichospermum circinale 371 AWQC310F) for which no homolog of PatU3, PatB, and HetF was found. 372

Next, we asked whether a relationship exists between the number of homologs in an organism and its 373 taxonomic lineage and, therefore, its ability to form heterocysts. For this purpose, a heatmap graphical 374 representation of the number of homologs within organisms, according to the taxonomic orders, was 375 performed (Fig. 7). Our analysis showed that except for HetC, nearly all the Nostocales contain at least 376 one gene encoding a homolog of each reference group, with a variable number of homologs per 377 organism, as described above. Interestingly, PatN seemed to be specifically present in Nostocales 378 exclusively. To a lesser extent, HetP and PatB also appeared to be poorly represented outside the 379 Nostocales group (Fig. 7). The distributions of the genes encoding HetC, HetF ',HetP,HetZ,PatU3,380 PatB and PatL homologs over the cyanobacterial lineage were similar to the one of hetR, suggesting that 381 which suggests that they may be involved in biological processes that are common to all these 387 prokaryotes. Unveiling the mechanisms beyond their speciation for heterocyst formation in the 388

Nostocales will be a significant progress in this field. 389

Another important question that needed to be addressed was whether the distribution of the studied 390 genes in Nostocales was significantly distinct compared with other taxonomic orders. For this, we 391 computed for each reference group of proteins the average number of genes encoding homologs per 392 genome for each cyanobacterial order (Materials and methods). These data are shown in the boxplots of 393 Fig. 8a (see also Table S3, Sheet 4). With the exception of patN, the other analyzed genes exhibited a 394 normal distribution among cyanobacterial orders, according to the Shapiro-Wilk test. As expected, HetL 395 and DevH exhibited a high number of homologs on average, regardless of the cyanobacterial taxonomic 396 orders. Interestingly, the presence of HetC, HetP, PatA, PatB, PatL and of PatN, was significantly higher 397 (P-value < 0.01) in Nostocales compared with other taxonomic orders (Fig. 8a). The high prevalence of 398 these six studied genes in Nostocales indicates their functional importance in heterocyst-forming 399 cyanobacteria, and are also considered as key factors genes. In addition, the data presented in this Fig. 400 highlight again the specific occurrence of PatN in genomes of Nostocales strains. 401

Altogether, when combined the two comparative genomics analyses (Fig. 7 In this study, we combined protein-protein interaction and genomics analyses to determine the physical 408 and evolutionary relationships that exist between the factors involved in heterocyst development and 409 patterning. BACTH has proven to be a powerful approach for probing interactions among the proteins 410 involved in the same process. For instance, it allowed the deciphering of the interactions that occur 411 between the membrane proteins involved in cell division in E. coli [START_REF] Karimova | Interaction network among Escherichia coli membrane 611 proteins involved in cell division as revealed by bacterial two-hybrid analysis[END_REF], in cell wall 412 synthesis in Corynebacterium glutamicum [START_REF] Jankute | Elucidation of a protein-protein interaction 606 network involved in Corynebacterium glutamicum cell wall biosynthesis as determined by 607 bacterial two-hybrid analysis[END_REF], and in spore coat formation in Bacillus 413 subtilis [START_REF] Krajcikova | Exploring the interaction network of the Bacillus 629 subtilis outer coat and crust proteins[END_REF]. In eukaryotes, a yeast-two-hybrid system was used to determine the 414 physical interactions between proteins involved in vulva development in Caenorhabditis elegans 415 [START_REF] Walhout | 701 Protein interaction mapping in C. elegans using proteins involved in vulval development[END_REF], and this study was the starting point for the elucidation of the global interactome 416 of this organism [START_REF] Li | A map of the 646 interactome network of the metazoan C. elegans[END_REF]). In a similar way, our study may help reveal the global interactome in 417

Nostoc. The fact that the interactions observed by BACTH for two of the analyzed proteins (HetL and 418 HetR) were confirmed by a BioLayer interferometry [START_REF] Xu | HetL, HetR and PatS form a reaction-711 diffusion system to control pattern formation in the cyanobacterium nostoc PCC 7120[END_REF]) and that pull-down assay confirmed 419 the interaction of PatB with 3 factors (Fig. 6) validates the BACTH approach. One of the interesting 420 perspectives of this work will be the study of all these possible interactions in Nostoc throughout the 421 differentiation process. The development of techniques such as bimolecular fluorescence 422 complementation (BifC) [START_REF] Kerppola | Design and implementation of bimolecular fluorescence complementation (BiFC) 616 assays for the visualization of protein interactions in living cells[END_REF]) and mass spectrometry [START_REF] Gingras | Analysis of protein complexes using mass 586 spectrometry[END_REF]) should help 423 achieve this objective. This study represents an important step-forward future works aiming to reveal 424 the interactions occurring in the heterocyst. We propose to refer to the protein network involved in 425 pattern formation as the "patternome". 426

The possibility that proteins with transmembrane segments may not be correctly inserted in the 427 membrane of E. coli constitutes a limitation of the BACTH assays. We did not observe any interactions 428 for the HetC, HetF, and PatN proteins, which are predicted to be membrane proteins. In the case of 429 HetC, in addition to this prediction, its localization at the inner membrane of the heterocyst has been 430

proposed, as assessed using fusion to the green fluorescent protein (Corrales-Guerrero et al. 2014a). 431

Similarly, membrane localization has been indicated for PatN in Nostoc punctiforme [START_REF] Risser | Biased inheritance of the protein PatN frees vegetative cells to 676 initiate patterned heterocyst differentiation[END_REF]. 432

The possibility that the lack of interactions of these proteins could be related to an incorrect topology of 433 the transmembrane domains in E. coli cannot be ruled out. HetF has recently been shown to be part of 434 the divisome of the vegetative cells under specific light conditions [START_REF] Xing | HetF Protein Is a New Divisome 709 Component in a Filamentous and Developmental Cyanobacterium[END_REF]. The fact that no 435 interaction was found with the proteins involved in heterocyst formation and patterning might therefore 436 reflect that it does actually not belong to the protein network formed by these factors. 437

Fourteen of the 16 protein interactions observed in our BACTH analysis were detected for the first time. 438

These new data provide additional clues about the molecular mechanisms underlying patterning, as they 439

show that the three steps of the process are likely to be connected through physical interactions between 440 key factors (Fig. 5). Strikingly, our results indicate that the transcriptional regulator PatB interacts with 441 multiple partners, which renders it a central factor in the network (Fig. 3, 4, and6). Among the cytosolic 442 proteins, PatD alone was shown to not interact with any other factor (Fig. 5). This result agrees with the 443 hypothesis that the effect of PatD on patterning might occur indirectly through general regulators, such 444 as ppGpp [START_REF] Wang | patD, a Gene 704 Regulated by NtcA, Is Involved in the Optimization of Heterocyst Frequency in the 705 Cyanobacterium Anabaena sp. Strain PCC 7120[END_REF]. The BACTH analysis also revealed that PatA interacts with two proteins 445 that are required for cell division, ZipN and SepF [START_REF] Valladares | Interactions of PatA with the Divisome during 691 Heterocyst Differentiation in Anabaena[END_REF]. As PatA also interacts with 446

PatB and PatL, taken together, these data show that the physical protein interactions described here may 447 impact the evolution of the heterocyst toward its terminal differentiation state (i.e., the arrest of cell 448 division). Contrary to what has been published [START_REF] Du | Expression from DIF1-motif 570 promoters of hetR and patS is dependent on HetZ and modulated by PatU3 during 571 heterocyst differentiation[END_REF], we did not observe any interaction 449 between HetR and HetZ (Table 1). The facts that we systematically tested combinations of the fused 450 constructs with the BACTH control vectors as negative controls (Fig. 4a) and that we considered as 451 being positive only values that were higher than those obtained for these controls may explain the 452 discrepancy between our data and the published ones. The genomics analysis indicated that the 14 genes 453 analyzed here are highly conserved in the genomes of heterocyst-forming cyanobacteria (Fig. 7); 454 therefore, an evolutionary history of these genes can be drawn. In fact, the location of almost all genes 455 over the taxonomic lineages suggests that they were present within the genome of the cyanobacteria 456 common ancestor. Throughout the evolution process, some of these genes were lost and/or functional 457 specialization might have occurred for those present in multiple copies (e.g., HetL). The process of gene 458 loss appears to be faster in some lineages than others, with Chroococcaceae, Gloeomargaritaceae, 459

Gloeobacteraceae, and Synechococcaceae being the clades in which the gene losses are the most 460 pronounced. Thus, gene loss occurs in clades encompassing strains that do not exhibit cell 461 differentiation, whereas Nostococcales genomes have conserved this set of genes (Fig. 7). Regarding 462 functional specialization, the most emblematic example is that of HetF and HetF', as the present copy 463 of HetF would probably be closer to the ones of the common ancestor. If future investigation unveils a 464 role for HetF' in the differentiation process, then one might suggest that its functional specialization 465 during evolution would have consisted in a loss of the transmembrane region, followed by the 466 acquisition of the Forkhead-associated domain (FHA). As the FHA domain has been maintained during 467 evolution, its importance in the differentiation process and in the interaction of HetF' with PatB, DevH, 468 and HetL may be questioned. 469

The heatmap generated here revealed that various levels can be proposed for the evolutive and functional 470 importance of the analyzed proteins in each step of the differentiation process. Among the early factors, 471 HetF', HetR, and HetC are the major evolutionarily conserved proteins; hetF' and hetR for their single 472 copy in the genome, and hetC for a very significant presence bias in favor of the Nostocales compared 473 with the other orders (Fig. 7). The essential role of HetR in Nostoc differentiation is fully understood. 474

As mentioned above, both the protein interaction analysis and the genomics survey suggest that the 475 HetF' protein is an interesting factor that deserves investigation. Genetic analyses revealed the 476 involvement of HetC and its N-terminal domain in heterocyst differentiation (Khudyakov and Wolk 477 1997;[START_REF] Xu | Role for hetC in the transition to a nondividing state during heterocyst 713 differentiation in Anabaena sp[END_REF]Corrales-Guerrero et al. 2014a). However, whether HetC actually acts as an 478

ABC transporter remains to be demonstrated. A comparative and evolutionary analysis of HetP, HetZ, 479 and PatU3 showed that all of them are key players in the commitment step, even if there is a significant 480 presence bias in favor of HetP in the Nostococales. It has been proposed that cell-fate commitment in 481

Nostoc is modulated by a hierarchy of four homologous HetP proteins, with HetP promoting cell 482 differentiation, whereas the HetP-like proteins seem to exert an opposite control [START_REF] Videau | The 695 heterocyst regulatory protein HetP and its homologs modulate heterocyst commitment in 696 Anabaena sp. strain PCC 7120[END_REF]. 483

These observations combined with the presence of multiple copies of HetP-encoding genes in 484

Nostocales clearly demonstrate that, during the evolutionary time, some HetP paralogs underwent a 485 interacts with factors belonging to the three defined groups and is mostly present in heterocyst-forming 494 cyanobacteria. Little is known about the mechanism of this regulator and about the genes that are 495 controlled by it in Nostoc, although it appears from our study that it deserves additional attention. Among 496 all the key players defined here, 8 out of 10 exhibited interactions with each other and with other 497 partners. This shows that protein-protein interactions may play pivotal roles in heterocyst 498 differentiation. Therefore, unveiling their dynamics through this process is a determining direction for 499 the understanding of development in these prokaryotes. 500

Given that three factors in the interaction network are transcriptional regulators (HetR, DevH, and PatB), 501 interactions occurring inside the network likely impact the transcriptional program of the developing 502 cell. To obtain insights into this interconnection, we first referred to the literature to build the 503 Each rectangle of the matrix represents a tested two hybrid interactions (from horizontal and vertical 732 axes). The β-galactosidase activities correspond to the functional complementation of the T18/T25 733 adenylate cyclase domains following interaction of the two hybrid proteins. β-galactosidase activity 734 lower than 20 Miller units, or lower than the self-activation controls (ie: the fused-construct tested with 735 the control vector) was considered as negative. Those whose interactions were detected as negative by 736 all possible combinations are marked as dark grey. Those whose interactions were detected at least by 737 one combination are marked as blank (see Fig. 3 and4). The repeat copy of each protein pair in this 738 table is marked as light grey. The size of the gray cycles is relative to the total number of biological links, including protein-protein 820 interactions and transcription regulatory actions. Key player proteins (as defined by in silico analysis) 821 are underlined. The transcriptional regulations include direct and indirect actions (see Table S4 for 822 details). As described within the main text, HetN, PatS and PatS were not included in our analyses. 823 824

  they might share a similar evolutionary history from the ancestor of cyanobacteria. Considered this 382 observation and the central and essential role of HetR within the heterocyst differentiation process, it 383 can be suggested that HetC, HetF'; HetP, HetZ, PatU3; PatB, PatL, and PatN might also act as key 384 factors for the early, commitment, and patterning stages, respectively, in Nostocales. Fig. 7 also shows 385 that HetL, DevH, PatA, and PatD are conserved over the entire cyanobacterial taxonomic lineages, 386

  and 8a), the overall key 402 factors consist of a network encompassing 10 proteins (HetC, HetR, HetF', HetP, HetZ, PatU3, PatA, 403 PatB, PatL, and PatN) that subsequently refers to the patterning factor network of heterocyst-forming 404 cyanobacteria from an evolutionary point of view. 405 406 Discussion 407

  neofunctionalization process, leading to new functions. The interactions of HetP with PatB and PatU3 486 with HetZ indicate that the protein complexes formed might explain the role of these proteins in 487 connecting the initiation step with subsequent steps. From an evolutionary point of view, PatA, PatB, 488 PatL and PatN appear as the key players among the patterning factors. The presence of PatN exclusively 489 in Nostocales seems to indicate that its action is important for the formation of heterocysts in all 490 organisms of this clade. Therefore, the elucidation of the molecular mechanism of PatN is crucial for 491 the understanding of this process. The significant enrichment of PatA in the Nostocales also makes it a 492 key factor that deserves further attention. In the protein network, PatB occupies a central position as it 493

Figure 7 807

 7 Figure 5

Figure

  Figure

  /ftp.ncbi.nlm.nih.gov/hmm/4.0/ ftp sites, respectively. Proteins involved in cell differentiation in 228 Nostoc used as reference seed proteins (retrieved from Uniprot database, www.uniprot.org), are 229 indicated in Table

	225					. Hidden Markov Models
	226	(HMMs) of protein family profiles (Pfam (version 33.0, February 2020) and Pgap (version 4, November
	227	2020)	were	downloaded	from	http://ftp.ebi.ac.uk/pub/databases/Pfam/releases/	and
		ftp:/

Table 1 : Matrix of BACTH data highlighting negative interactions between proteins 728 involved in heterocyst differentiation within Nostoc.

 1

	729
	730

Acknowledgments: 532

The authors thank D r Badreddine Douzi and D r Romé Voulhoux for helpful discussions regarding pull-533 down assays. They also thank Anais Scholivet for technical assistance. The domains predicted by Pfam are indicated by blue boxes. For HetR, the domains correspond to those 752 identified by the resolution of the structure [START_REF] Hu | Structural insights into HetR-601 PatS interaction involved in cyanobacterial pattern formation[END_REF]. The topology of the membrane protein was 753 predicted by the TMHMM software [START_REF] Krogh | Predicting transmembrane protein 631 topology with a hidden Markov model: application to complete genomes[END_REF].

754

Protein sizes are in parenthesis (aa: amino acids).

Please read the important information on page 4 before you begin

This form should be used by authors to request any change in authorship including changes in corresponding authors. Please fully complete all sections. Use black ink and block capitals and provide each author's full name with the given name first followed by the family name.

Please note: In author collaborations where there is formal agreement for representing the collaboration, it is sufficient for the representative or legal guarantor (usually the corresponding author) to complete and sign the Authorship Change Form on behalf of all authors.

Section 1: Please provide the current title of manuscript (For journals: Please provide the manuscript ID, title and/or DOI if available.) (For books: Please provide the title, ISBN and/or DOI if available.)

First name(s)

Family name (this name will appear in full on the final publication and will be searchable in various abstract and indexing databases) The change in the author-list consists in including Raphaël Rachedi as co-first author The reason for this addition is thar R. Rachedi performed all the co-expression and copurification and immunoblotting experiments that were performed to answer the reviewers and editor requests (Figure 6 and Figure S1)