Thibault Bonnemain

Thierry Gobron

Denis Ullmo

Universal behavior in non stationary Mean Field Games

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Mean Field Games provide a powerful framework to analyze the dynamics of a large number of controlled objects in interaction. Though these models are much simpler than the underlying differential games they describe in some limit, their behavior is still far from being fully understood. When the system is confined, a notion of "ergodic state" has been introduced that characterizes most of the dynamics for long optimization times. Here we consider a class of models without such an ergodic state, and show the existence of a scaling solution that plays a similar role. Its universality and scaling behavior can be inferred from a mapping to an electrostatic problem.

Mean Field Games are a powerful framework introduced about a decade ago by Lasry and Lions [1] as an alternative approach to differential game theory when the number of agents becomes large (see [2] for a suitable introduction to the field for physicists). Their applications are numerous, ranging from finance [3,4] and economy [START_REF] Guéant | Mean field games and applications[END_REF][START_REF] Achdou | [END_REF] to engineering sciences [START_REF] Mériauxi | Conference Record of the Forty Sixth Asilomar Conference on Signals, Systems and Computers[END_REF][START_REF] Kizilkale | Control of Complex Systems[END_REF], and wherever one has to deal with optimization issues for many coupled subsystems. On a quite general basis, a Mean Field Game involves a set of N players (or agents) which are characterized by a continuous "state variable" X i ∈ R d , i = 1 . . . N , which, depending on the context, may represent a physical position, the amount of resources owned by a company, the house temperature in a network of controlled heaters, etc.. These state variables evolve on a time interval [0, T ] according to some controlled dynamics, which we assume here to be described by a linear, d-dimensional, Langevin equation, dX i t = a i t dt + σdW i t , where each component of W i is an independent white noise of variance one, σ is a constant, and the "control parameter" is the velocity a i t . This control is adjusted in time by the agent i in order to minimize a cost functional over the remaining time to play, which in the simplest case can be assumed of the form

c[a i ](x i t , t) = E T t µ 2 (a i t ) 2 -V [m τ ](X i τ ) dτ + c T (X i T ) X i t = x i t . (1)
In ( 1), c T (x) is the terminal cost, an additional cost player i would pay if his state at the end of the game, X i (T ), is equal to x, µ is positive constant, akin to a mass in physics, and V [m t ](x) is a functional of the empirical density at time t m t (x) = 1 N j δ(x -X j (t)), through which the agents' optimization problems are coupled. We shall assume V [m t ](x) takes the simple form

V [m t ](x) = U 0 (x) + g m t (x) . (2)
For a very large number of players, like in a mean field theory, the fluctuations of the empirical density are neglected and m t (x) becomes a deterministic quantity governed by a Fokker Planck equation. Furthermore, the optimization problems decouple and the optimal value of the cost (1) (the "value function") for the agent i becomes a function of the variable x i t , solution of an Hamilton-Jacobi Bellman equation [START_REF] Bertsekas | Dynamic Programming and Optimal Control[END_REF]. The resulting model is called a Mean Field Game and can be defined as a pair of coupled equations, describing the (forward in time) evolution of the density of players m(x, t) and the (time backward) evolution of the value function u(x, t). In the simple case we consider here the system of deterministic equations takes the form [2, 10]

       ∂ t m - 1 µ ∇ [m∇u] - σ 2 2 ∆m = 0 (3)
∂ t u - 1 2µ [∇u] 2 + σ 2 2 ∆u = V [m] . (4)
The coupling between the two PDE's comes from two parts: in the Fokker-Planck equation (3), the optimal velocity appears in the drift term and here is proportional to the gradient of the value function, a = -1 µ ∇u; in the Hamilton-Jacobi Bellman equation (4) the term V [m] reflects the dependence of the cost functional (1) on the density. This structure also induces rather atypical boundary conditions: the (forward) Fokker-Planck equation is associated with an initial condition m(0, x) = m 0 (x) specifying the initial distribution of agents, while the terminal cost C T in Eq. (1) imposes a final condition for the value function, u(T, x) = c T (x). This forward-backward structure together with mixed initialfinal boundary conditions leads to new challenges when trying to characterize, either analytically or numerically, solutions to this system of equations.

For a large class of settings, which includes in particular the case of repulsive interactions when the system is confined either by a stabilizing external potential or because of a geometry with fixed spatial extension, such a system exhibits an ergodic (stationary) state, independent on the boundary conditions, which can be rigorously defined in the limit T → ∞ as an hyperbolic fixed point. The importance of this result, as proven by Cardialaguet et al. [11] is twofold: for finite but long enough optimization time T , the game will stay very close to this ergodic (time independent) state except possibly in its initial and final parts. Furthermore, the transient dynamics near t = 0 and t = T completely decouple one from the other, and the mixed boundary problem simplifies accordingly: they describe how the initial and final boundary conditions match with the ergodic state and are characterized by two possibly different time scales, τ i and τ f , independent on T .

The problem we want to address in this letter is the behavior of such a system when an ergodic state cannot exist, for instance when the state space is unbounded and the interactions repulsive, so that any initially localized configuration will expand forever. The natural question in such conditions is whether some kind of limiting (non-stationary) regime still exists for such systems that would, at some level, play a similar role as the ergodic state in confined systems.

Considering more specifically a one-dimensional problem without "one-body" term (i.e. U 0 (x) ≡ 0), a first indication that this is indeed the case is brought out by numerical solutions for the system (3)-( 4) for repulsive interactions (that is V [m t ](x) = gm(t, x), g < 0), starting from a localized initial density. Numerical results (see Figure 1) show clearly that except for transient times near t = 0 and near t = T , the density can be accurately fitted on an inverted parabolic shape

m(x, t) =    1 N (t) (z(t) 2 -x 2 ) if |x| ≤ z(t) 0 otherwise , (5)
where the prefactor is found to be N (t) = 4 3 z(t) 3 , as total mass conservation requires. Assuming this shape as an ansatz, the evolution reduces to that of the scaling factor z(t), which is for most times found to grow as a power law (see left insert in Figure 1)

z(t) ∼ t 2/3 . ( 6
)
These features appear for sufficiently long optimization time T , and are essentially independent of the initial and final boundary conditions provided the initial distribution has a bounded extension and the final cost is close to zero everywhere. By analogy with the physics of Bose condensate, let us introduce a characteristic length-scale, the healing length ν = µσ 2 /|g|, which represent the typical distance on which interactions balances quantum pressure (or diffusion in our case) [12]. A necessary condition for the behavior observed in Fig. 1 to appear is the smallness of ν with respect to other length-scales. This condition will however eventually be fulfilled at one point as the extension of density distribution grows in time, possibly inducing a time shift in [START_REF] Achdou | [END_REF].

What those numerical results tell us is that in this setting, the notion of ergodic state has been replaced by the next best thing, namely a universal scaling solution. The goal of this paper is to understand this puzzling result, in particular the 2 3 scaling exponent, and give some arguments in favor of its universality.

For this purpose, we now introduce a few formal transformations which allow to show that this result is rather natural and intuitive, once the problem is cast in its proper language. In turn, we will also gain a better understanding of this regime

I. HYDRODYNAMIC REPRESENTATION

The main idea behind the approach below comes from the deep link between Mean Field Games and the nonlinear Schrödinger equation, discussed at lengths in [2], which allows for the use of several techniques developed to study Bose-Einstein condensates. One of these techniques, the Madelung substitution, is particularly well suited to deal with the small ν regime [START_REF] Pethick | Bose-Einstein Condensation in Dilute Gases[END_REF]. It consists in defining a velocity field v(t, x)

v(x, t) = - ∇u µ - σ 2 ∇m 2m , (7)
which maps the evolution (3) for the density m to a simple transport/continuity equation

∂ t m + ∇(mv) = 0 . (8)
The evolution for the velocity field v derives from the HJB equation ( 4) and reads

∂ t v + ∇ σ 4 2 √ m ∆ √ m + v 2 2 + g µ m = 0 , (9)
and involves a O(σ 4 ) term. As in the context of cold atoms, this term can be neglected as long as the charac-teristic length of the system is large in front of the healing length ν = µσ 2 /|g| [START_REF] Pethick | Bose-Einstein Condensation in Dilute Gases[END_REF], leading to what is referred to as the Thomas Fermi approximation. This weak noise limit is also one of the requirements for the appearance of the parabolic behavior described in Figure 1. In this approximation, the O(σ 4 ) term in ( 9)is dropped and the equations read

   ∂ t m + ∇(mv) = 0 ∂ t v + ∇ v 2 2 + g µ m = 0 . ( 10
)
These equations are formally very close to those studied for instance in the field of cold atoms [12], the main differences being the (negative) sign of g, which makes the system elliptic rather than hyperbolic, and the nature of the boundary conditions. Within this approximation, it can easily be verified that the ansatz Eq. ( 5) with

z(t) = 3 |g| 4µ 1/3 t 2/3 , (11)
is a particular solution of the Thomas-Fermi-like equations Eq. ( 10). The real mystery is therefore not that such a particular solution does exist, but to understand why this t 2/3 behavior shows up at all intermediate times for a large class of boundary conditions, and in other terms, why it is universal in the large optimization time limit.

II. RIEMANN INVARIANT AND HODOGRAPH TRANSFORM

To answer this question, we shall turn to an approach developed in the context of non-linear waves [START_REF] Kamchatnov | Nonlinear Periodic Waves and Their Modulations: An Introductory Course[END_REF], which relies on the notions of Riemann invariants and hodograph transform. Riemann's method can be considered an extension of the method of characteristics. It amounts to finding curves (characteristics) on which some quantities (Riemann invariants) are conserved. Here, one can show that there exists a pair of Riemann invariants, (λ + (x, t), λ -(x, t)), namely λ ± = v ± 2i |g|m/µ, so that (10) reads

       ∂ t λ + + 3 4 λ + + 1 4 λ -∂ x λ + = 0 ∂ t λ -+ 1 4 λ + + 3 4 λ -∂ x λ -= 0 . ( 12
)
Though characteristics do not exist in this context (they are curves in the complex plane C 2 ), this change of variables still allows us to linearize these equations using an hodograph transformation [START_REF] Kamchatnov | Nonlinear Periodic Waves and Their Modulations: An Introductory Course[END_REF]. Taking the pair (λ + , λ -) as independent variables, we express x and t as functions of them, so that the system (12) transforms into a linear one

∂ -x -β + ∂ -t = 0 ∂ + x -β -∂ + t = 0 , ( 13
)
where β ± = 3 4 λ ± + 1 4 λ ∓ , and where we have introduced the notation ∂ ± ≡ ∂ ∂λ± . This system can be readily integrated once as

x -β + t = ω + x -β -t = ω - , (14)
with ω ± solution of

∂ ± ω ∓ = -(∂ ± β ∓ )t = - 1 4 t . ( 15
)
Thus the functions ω ± can be expressed as derivatives of a potential ω ± = ∂ ± χ, where χ(λ + , λ -) is solution of an Euler-Poisson-Darboux equation:

∂ 2 +-χ - 1 2(λ + -λ -) (∂ + χ -∂ -χ) = 0 . ( 16
)
The main difference with the traditional treatment of NLS that we have closely followed until now is that here the Riemann invariants (λ + , λ -) are complex conjugates. In terms of its real and imaginary parts, ξ =

1 2 (λ + + λ -) = v and η = 1 2i (λ + -λ -) = 2 |g|m/µ, equation (16) becomes ∂ ξξ χ + ∂ ηη χ + 1 η ∂ η χ = 0 , (17)
which is the Laplace equation in cylindrical coordinates (with no angular dependence) with η and ξ as radial and axial coordinates, respectively. Equations ( 14) now read

ηt = -E η 2(x -ξt) = -E ξ , (18)
with E η and E ξ the radial and axial components of the electric field, E = -∇χ. Note that even if (17) is originally a two-dimensional problem, a clear connection with electrostatics emerges when considering it a three dimensional one with axial symmetry.

III. POTENTIAL REPRESENTATION

Through the hodograph transform we have shown that for any potential χ, solution of the Laplace equation ( 17), there is a solution to the Thomas Fermi equation [START_REF] Swiecicki | [END_REF] provided that the relations (18) between x, t and the electric field E hold. The linear Laplace equation (and the related electrostatic problem) is clearly significantly simpler than the original non-linear hydrodynamic equations. The price to pay for that simplification is that taking into account the boundary conditions becomes highly non trivial since the locus of the curves t(ξ, η) = 0 or t(ξ, η) = T on which these conditions are expressed actually depend on the particular potential χ(ξ, η) considered.

However, since the dynamics we are interested in is associated with the spreading of the density of agents, the curves t(ξ, η) = const. are contracting as t increases, smaller time curves including every larger time ones. If we consider χ as generated by a distribution of charge ρ(ξ, η), Eq. ( 17) implies that ρ(ξ, η) = 0 between the curves t(ξ, η) = 0 and t(ξ, η) = T but can be non-zero either near the origin (for times larger than T ) or at large distance (corresponding to negative times). If the optimization time is long enough we can assume that there exists a range of times [ tmin , tmax ], 0 tmin , tmax T , so that for any curve t(ξ, η) ∈ [ tmin , tmax ] the distributions of charges both at the origin and "at infinity" are far enough so that the effects of their moments are essentially negligible. Since the total charge at infinity contributes only by an irrelevant constant, we can give a good approximation of χ in this range of time as the potential created by a point charge Q 0 located at the origin

χ(η, ξ) ≈ Q 0 η 2 + ξ 2 , ( 19
)
with a relation between Q 0 and the boundary conditions of the problem yet to be determined.

The main result of this paper is that the approximation of the charge distribution by a monopole centered at the origin is precisely the observed universal behavior expressed by Eqs. ( 5)-( 6). In addition the conditions under which this approximation is valid provide the regime of validity of the scaling form Eqs. ( 5)- [START_REF] Achdou | [END_REF]. Indeed, inserting Eq. ( 19) into Eq. ( 18) and inverting the relation between (ξ = v, η = |g|m/µ) and (x, t) readily gives

         m(t, x) = 3 (µQ 0 /2g) 2/3 z(t) 2 -x 2 4z(t) 3 v(t, x) = - z (t) z(t) x . ( 20
)
for |x| < z(t), with z(t) given by Eq. (11). It suffices now to show that independently of the boundary condition we have Q 0 = 2g/µ (which is in any case required by the normalization of m(x, t)), to prove that Eq. ( 20) is in fact the numerically observed scaling form Eqs. ( 5)-( 6).

To show this, we can simply apply Gauss's law

St ( E • n)dS = 4πQ 0 (21)
on a surface S t such that t(η, ξ) = t = const., for any 0 ≤ t ≤ T . Parametrizing this surface by (x, θ), with dS = ηj(x, t)dθdx (with j(x, t) ≡ (∂ x ξ) 2 + (∂ x η) 2 ), and the normal n to the surface S t given by n = (n ξ , n η , n θ ) = j(x, t) -1 (∂ x ξ, -∂ x η, 0), we get

Q 0 = 1 4π 2π 0 R η 2(x -ξ t)∂ x η -η t∂ x ξ dθdx = 1 2 R -t∂ x (η 2 ξ) + 2xη∂ x η dx . (22)
If we assume η to decrease sufficiently fast with x, the first, time dependent, term integrates to zero. This was to be expected because no matter the time t, as long as 0 < t < T , the total charge included in S t is the same: Q 0 is by construction a constant of the motion. Integrating by part and recalling that η = 2 |g|m/µ, Eq.( 22) yields Q 0 = 2g µ R mdx, which, because of the normalization condition on m(x, t), is the required result.

IV. CONCLUSION

Introducing a potential representation for the [1d potential-free] Mean Field Game problem Eqs. ( 3)-( 4), we have shown that the remarkable, and a priori quite puzzling universal scaling form we have observed numerically for large optimization times can be derived in a very natural way. In this representation, it appears as the simple fact that the related potential χ, sufficiently far away from the charge distribution which creates it, can be well approximated by the potential generated by a single monopole. The condition of being sufficiently far from both the charges near the origin and the one near infinity is reflected in the MFG problem by considering times far from both t = 0 and t = T , which is of course only possible in the limit of very long optimization time. In some sense, we have thus been allowed to replace the notion of ergodic state by what could be thus called "an ergodic scaling form" in a situation where a genuine ergodic state cannot exist.

At a more general level, the potential representation of this Mean Field Game underlies the integrability of the hydrodynamical equations [START_REF] Swiecicki | [END_REF]. Beyond the simple monopolar approximation Eq. ( 19), we can construct a complete multipolar expansion for the potential χ, and it can be seen easily that each "charge" of this expansion correspond to a conserved quantity of the dynamics. The mapping between the boundary conditions and these charges being thus equivalent to the mapping between boundary conditions and constants of the motion. This consideration emphasizes the fact that the deep reason behind the scaling law characterizing the 1d potential free MFG that we consider in this paper is actually their integrable character. A more in depth discussion of this question will appear in a separate publication [START_REF] Bonnemain | [END_REF].

Figure 1 .

 1 Figure 1. Numerical simulation of the time evolution described by the system (3)-(4) in one space dimension with T = 200, g = -2, σ = 0.5 and µ = 1. The initial distribution is a Gaussian of variance 0.1 centered at the origin, and the terminal cost is cT (x) = 0. The left inset shows a log-log plot of the time evolution of the density at the origin m(0, t) (full line), showing a t 2/3 scaling behavior (dashed line) at almost all times. On the right inset the rescaled density z(t)m(x/z(t), t) is shown to have an invariant parabolic shape (with z(t) = 4/(3m(0, t))).


Dr Alexandre Sadate
email: alexandre.sadate@chu-nimes.fr

Bob V Occean

Jean-Paul Beregi

Aymeric Hamard

Takieddine Addala

Hélène De Forges

Pascale Fabbro-Peray

Julien Frandon

Systematic review and meta-analysis on the impact of lung cancer screening by low-dose computed tomography

Keywords: Lung cancer, screening, low-dose computed tomography, meta-analysis List of abbreviations CI: Confidence Interval CXR: Chest X-ray LC: Lung Cancer LDCT: Low-Dose Computed Tomography NLST: National Lung Screening PLCO: Prostate, Lung, Colorectal and Ovarian Cancer Screening

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Lung cancer shows the highest cancer mortality worldwide and a poor prognosis with a 5-year survival rate of 16% in the United States [START_REF] Siegel | Cancer statistics[END_REF]. It is the fourth most common cancer in France; around 30 000 lung cancer specific deaths in 2012 for around 40 000 diagnosed diseases (2). The incidence rates have been stable among men whilst they are rising rapidly in women in France and Europe. One cause may be that the main risk factor of lung cancer is tobacco consumption, that increases in women yet decreases in men [START_REF] Francim | Service de biostatistique des Hospices civils de Lyon, Institut de veille sanitaire[END_REF][START_REF] Bray | Lung cancer mortality trends in 36 European countries: secular trends and birth cohort patterns by sex and region 1970-2007[END_REF]. This confirms the promotion of smoking cessation is the most effective strategy to reduce LC mortality. One major issue is early detection, as nearly half of patients are diagnosed at an advanced stage, with various prognosis depending on cancer stage at diagnosis [START_REF] Walters | Lung cancer survival and stage at diagnosis in Australia, Canada, Denmark, Norway, Sweden and the UK: a population-based study, 2004-2007[END_REF].

Two main trials have evaluated the benefits of chest radiography with or without sputum cytology, the PLCO Cancer Screening Trial and the Mayo Lung Project, with no efficacy in reduction of cancer-specific mortality [START_REF] Gohagan | Colorectal and Ovarian Cancer Screening Trial Project Team. The Prostate, Lung, Colorectal and Ovarian (PLCO) Cancer Screening Trial of the National Cancer Institute: history, organization, and status[END_REF][START_REF] Fontana | Lung cancer screening: the Mayo program[END_REF]. Several years later, a Japanese trial showed that low-dose computed tomography (LDCT) was superior to chest radiograph for the detection of peripheral lung cancer in a population highly exposed to tobacco [START_REF] Kaneko | Peripheral lung cancer: screening and detection with low-dose spiral CT versus radiography[END_REF]. The Early Lung Cancer Action Project (ELCAP) program then showed that LDCT not only detected more lung cancers than chest radiography but also at an earlier and more resectable stage, which could potentially reduce cancer-specific mortality [START_REF] Henschke | Early Lung Cancer Action Project: overall design and findings from baseline screening[END_REF].

Radiation risk is one of the main obstacles preventing a more wide-spread lung cancer screening in populations highly exposed to tobacco. However, LDCT allows dose reduction without altering image quality [START_REF] Macri | Ultra-low-dose chest CT with iterative reconstruction does not alter anatomical image quality[END_REF], thanks to tube voltage and current reduction and iterative reconstructions.

Based on these results, several randomized controlled trials (RCT) were launched to evaluate the efficacy and cost effectiveness of a potential lung cancer screening in high-risk populations.

Following the results of a large randomized American trial, the National Lung Screening Trial (NLST) [START_REF] Lung | Reduced lung-cancer mortality with low-dose computed tomographic screening[END_REF], which demonstrated a significant relative reduction in all-cause and lung cancer-specific mortality, the US Preventive Services Task Force recommended screening with LDCT all adults between 55 and 80 years with a 30 pack-year smoking history [START_REF] Moyer | Screening for lung cancer: U.S. Preventive Services Task Force recommendation statement[END_REF]. Other international guidelines have also published recommendations supporting lung cancer screening in various eligible populations [START_REF] Shen | Low-dose CT for lung cancer screening: opportunities and challenges[END_REF].

The objective of our systematic review and meta-analysis was to evaluate the efficacy of lung cancer screening by LDCT in populations highly exposed to tobacco on cancer specific-and overall mortality.

Materials and methods

The systematic review protocol was registered with PROSPERO on March 27 th , 2018 (CRD42018091720). Results are reported here following the PRISMA Equator Network guidelines [START_REF] Moher | Preferred reporting items for systematic reviews and meta-analyses: the PRISMA statement[END_REF].

Search strategy

The search was made using the Medline and Cochrane Library databases. The research terms were covering lung cancer, screening by LDCT compared with any other intervention, randomized controlled trials (RCT) and their synonyms, from inception until the February 19 th , 2018, as follows:

(((("Lung Neoplasms"[Mesh] OR lung neoplasm*[TIAB] OR pulmonary neoplasm*[TIAB] OR bronchopulmonary neoplasm*[TIAB] OR broncho-pulmonary neoplasm*[TIAB] OR bronchial neoplasm*[TIAB] OR lung cancer*[TIAB] OR pulmonary cancer*[TIAB] OR bronchopulmonary cancer*[TIAB] OR broncho-pulmonary cancer*[TIAB] OR bronchial cancer*[TIAB] OR lung carcinoma*[TIAB] OR pulmonary carcinoma*[TIAB] OR bronchopulmonary carcinoma*[TIAB] OR broncho-pulmonary carcinoma*[TIAB] OR bronchial carcinoma*[TIAB] OR bronchogenic carcinoma*[TIAB] OR lung blastoma*[TIAB] OR pulmonary blastoma*[TIAB] OR bronchopulmonary blastoma*[TIAB] OR broncho-pulmonary blastoma*[TIAB] OR bronchial blastoma*[TIAB] OR lung tumor*[TIAB] OR pulmonary tumor*[TIAB] OR bronchopulmonary tumor*[TIAB] OR bronchopulmonary tumor*[TIAB] OR bronchial tumor*[TIAB])) AND ("Mass Screening"[Mesh] OR "Early Detection of Cancer"[Mesh] OR screen*[TIAB] OR test[TI] OR testing[TI] OR detection*[TI])) AND (random*[TI] OR randomly[TIAB] OR randomized[TIAB] OR placebo[TIAB] OR Random Allocation[MH] OR Double-Blind Method[MH] OR Single-Blind Method[MH] OR Cross-Over Studies[MH] OR randomized controlled trial[PT] OR controlled clinical trial[PT]))))

It was completed by a Medline search with the names of the trials identified, the reference list of the eligible articles, a clinicaltrials.gov database search, a free-text Internet search and a regular update Medline search. The search was limited to articles published in English and French. The initial search was updated on February 12 th , 2020.

Study selection criteria

Inclusion criteria for the systematic review and meta-analysis were topics about lung cancer screening, RCT study design, LDCT compared to any other intervention, population who reported an average smoking history over 15 pack-years (corresponding to the lowest criteria of the European RCTs on lung cancer screening) and the report of data on all-cause mortality or lung cancer-specific mortality.

Two double-blind reviewers (AS, QD) selected the publications by screening the titles and abstracts first, and then on the full-text articles. Discrepancies were resolved by consensus between the two readers. If a consensus was not made, an epidemiologic expert (PFP), it was consulted.

Quality assessment and data extraction

The critical appraisal of each eligible RCT was made by two reviewers (AS, JF) using a CONSORT checklist including random assignment, complete diagnostic work-up planned, respect of inclusion criteria, valid measurement of mortality and blinded outcomes assessment. The contents of the checklist and the inclusion or exclusion of each study for the meta-analysis were then discussed with two epidemiologic experts (PFP, BO). Data were collected and managed using REDCap electronic data capture tools.

The extraction of the data for the meta-analysis was made independently by two double-blind reviewers (AS, JF), and the differences were resolved by an epidemiologic expert (PFP). When several articles were available for a single study, data from the most recent publication were used, unless unique data was found in a previous publication. The following data were extracted: title of the article, sample size, participants, characteristics (age, sex, smoking history), type of intervention in the control arm, median follow-up, number of screening rounds, time between rounds, number of any-cause and lung cancerassociated deaths.

Statistical analysis

Results on all-causes and lung cancer-specific mortality of the studies included were combined for meta-analysis. In the collection of studies considered the effect size and the true effect size were estimated. A general model was then specified by both standard and random effects models to incorporate the between-study variance [START_REF] Higgins | Quantifying heterogeneity in a meta-analysis[END_REF]. The inverse variance weighting method was used for pooling, on the hypothesis of common "effect size" normality and weights assigned to individual studies. This method allows maximizing the accuracy of the estimate of the common effect, minimizing the variance of the pooled result. The binary outcome data as mortality rate was extracted in a 2X2 table for all selected studies. A relative risk and 95% confidence interval (95% CI) were used as summary statistics. The heterogeneity between trials was taken into account using chi-square tests. I 2 statistic was used to estimate the percentage of total variation across studies arising from heterogeneity rather than chance. I 2 values of 25%, 50%, and 75% represented a low, moderate or substantial heterogeneity [START_REF] Fleiss | The statistical basis of meta-analysis[END_REF].

The reasons for substantial heterogeneity were explored and subgroup analyses were performed to consider possible variations across studies. Supplementary analyses excluding studies with heterogeneity and low contribution to the overall effect, then including only those with the highest weighting, were performed. A subgroup analysis was performed to take into account variability across studies during follow-up. Possible publication biases were explored through visual analysis of funnel plots. Statistical analysis was performed using General Package for Meta-Analysis "meta" v4.9-1 with R software v3.5.1.

Results

Study selection

Of the 891 publications identified with the search strategy in the Medline and Cochrane Library databases, 799 were excluded on titles and abstracts because of topic or study design: 449 were not related to lung cancer screening, 142 reported results of lung cancer screening without using LDCT, 40

were not RCTs and 168 were systematic review or general considerations. The additional search identified 22 other publications which were added to the 92 studies. Reading the full-text of these 114 publications, 107 studies were excluded: 11 were pilot study, 31 did not report mortality data in the control group, 59 did not assess mortality as primary endpoint, and 6 presented intermediate results. A total of 7 RCTs were finally included, including a total of 84,558 patients (Figure 1) [START_REF] Infante | Long-Term Followup Results of the DANTE Trial, a Randomized Study of Lung Cancer Screening with Spiral Computed Tomography[END_REF][START_REF] Wille | Results of the Randomized Danish Lung Cancer Screening Trial with Focus on High-Risk Profiling[END_REF][START_REF] Paci | Mortality, survival and incidence rates in the ITALUNG randomised lung cancer screening trial[END_REF][START_REF] Pastorino | Prolonged Lung Cancer Screening Reduced 10-year Mortality in the MILD Trial[END_REF][START_REF] Pinsky | The National Lung Screening Trial: results stratified by demographics, smoking history, and lung cancer histology[END_REF][START_REF] Becker | Lung cancer mortality reduction by LDCT screening -results from the randomised German LUSI trial[END_REF][START_REF] De Koning | Reduced Lung-Cancer Mortality with Volume CT Screening in a Randomized Trial[END_REF].

Study characteristics and quality assessment

The characteristics of the 7 RCTs included are detailed in Table 1. For all studies included in the meta-analysis, randomization was performed at inclusion. The intervention in the control arm was chest X-ray (CXR) at inclusion followed by annual clinical examination in the DANTE study, annual CXR in the NLST study, primary prevention for smoking cessation in four RCTs (NELSON, LUSI, ITALUNG, MILD) and annual clinical examination for the Danish Randomized Lung Cancer CT Screening Trial (DLCST). Five studies included both men and women; the DANTE and NELSON studies only included men.

All 7 trials presented a complete diagnostic work-up planned in the protocol with well-described and respected inclusion criteria. Validity of mortality measurement was checked in all studies. Quality assessment results of all 7 RCTs are shown in Table 2.

Overall and lung cancer-specific mortality

Two trials reported statistically significant results on mortality data. The NLST trial showed a relative reduction of lung cancer-specific mortality of 20% (95% CI: 6.8-26.7; p=0.004) and a relative reduction of overall mortality of 6.7% (95% CI: 1.2-13.6) in the LDCT group compared to the control annual CXR group. The NELSON study reported a relative reduction of lung-cancer-specific mortality of 24% (95% CI: 6-39; p=0.01) compared with standard primary prevention; there was no difference of overall mortality in the two groups.

The final results of the ITALILUNG, DANTE, MILD and LUSI trials did not find any difference for both overall and lung cancer-specific mortality in the LDCT groups compared to the control groups.

Publication biases

The risk of publication bias was evaluated using a funnel plot (Figure 2). The symmetric distribution of the relative risk across the global effect paired with the standard deviation of the screening effect confirmed the studies included did not present major biases.

Results of the meta-analysis

A total of 84,558 participants were included in the meta-analysis. There was no heterogeneity in the data (I 2 =0%, tau 2 =0, p=0.67). A relative reduction of overall mortality of 4% was observed in the experimental screening group versus control group (RR=0.96, 95% CI: 0.92-1.00) (Figure 3).

Concerning lung cancer-specific mortality, a significant relative reduction of 17% was observed in the experimental screening group (RR=0.83, 95% CI: 0.76-0.91) (Figure 4). To prevent one lung-cancer related death, 294 patients needed to be screened.

Discussion

Our meta-analysis on the impact of lung cancer screening using LDCT is, to our knowledge, the first meta-analysis including the final mortality results of the recently published NELSON, DLCST, ITALUNG, MILD and LUSI studies [START_REF] Wille | Results of the Randomized Danish Lung Cancer Screening Trial with Focus on High-Risk Profiling[END_REF][START_REF] Paci | Mortality, survival and incidence rates in the ITALUNG randomised lung cancer screening trial[END_REF][START_REF] Pastorino | Prolonged Lung Cancer Screening Reduced 10-year Mortality in the MILD Trial[END_REF][START_REF] Becker | Lung cancer mortality reduction by LDCT screening -results from the randomised German LUSI trial[END_REF][START_REF] De Koning | Reduced Lung-Cancer Mortality with Volume CT Screening in a Randomized Trial[END_REF]. Our results show a significant relative reduction of lung cancer-specific mortality of 17% in the screening group compared to the control group, but no significant difference in the mortality of all causes.

A previous meta-analysis [START_REF] Fu | A meta-analysis: is low-dose computed tomography a superior method for risky lung cancers screening population?[END_REF] based on the pooled analysis of four RCTs (DANTE, MILD, NLST and the intermediate results of DLCST) showed a significant decrease of lung cancer-specific mortality in the control group (OR=0.84 [0.78, 0.96]); no difference was reported for overall mortality between the two groups (OR=1.04 [0.72, 1.51]). However, this pooled analysis included RCTs with both short (intermediate results of MILD with a 4.4-year follow-up) or long follow-up (DANTE: 8.3year follow-up), one reason of the heterogeneity of the study (Q test p-value = 0.13). In our metaanalysis, we have included more studies as recent RCTs have been published since, including the final results of the MILD study, i.e. we included only final results for the RCTs included, and thus mostly long follow-up studies.

One major issue regarding lung cancer risk and screening is defining target populations at risk. Among studies included in the meta-analysis, only the NLST study included older participants with higher smoking history. Subgroups analyses divided the population into quintiles according to the risk of lung-cancer associated death [START_REF] Kovalchik | Targeting of[END_REF]. Lung-cancer specific mortality increased significantly with each quintile and the proportion of false positive results decreased significantly. Overall, 88% of the prevented lung cancer-associated deaths in the LDCT group corresponded to the 60% of participants within the highest risk subgroups, whereas 1% of the prevented lung cancer deaths accounted for 20% of participants in the lowest risk subgroups. Furthermore, the number of screened patients required to prevent one death from lung cancer decreased from 5276 in the lowest risk subgroup to 161 in the highest risk subgroup of patients. Another study with post-hoc analysis of the targeted population was published on the DLCST cohort [START_REF] Fu | A meta-analysis: is low-dose computed tomography a superior method for risky lung cancers screening population?[END_REF]. In the highest risk subgroup (patients with chronic obstructive pulmonary disease and a smoking history of more than 35 pack-years), there was a non-significant a 20% decrease in the hazard ratio for LC mortality in the LDCT group was reported although it did not reach significance (26). The recent NELSON study included patients with relatively small smoking history (less than 36 packs.years) as compared with other studies and especially the NLST study [START_REF] De Koning | Reduced Lung-Cancer Mortality with Volume CT Screening in a Randomized Trial[END_REF]. Overall, these results suggest the screening program may be restricted to patients with the highest risk, i.e. patients with a significant smoking history.

Our meta-analysis included the eagerly awaited recent NELSON study [START_REF] De Koning | Reduced Lung-Cancer Mortality with Volume CT Screening in a Randomized Trial[END_REF]. This well-designed study, the biggest European study, reported results after a long follow-up (10 years). The results were significant regarding the specific overall mortality, but did not show any difference between the two groups for the all-causes mortality. This was explained by the authors by the study design, as the study included around 13,000 patients, a number, although consequent, which was not sufficient to show a difference. Also, a strength of the NELSON study was the very low number of false-positive (1.2% only). This may be due to the study design and the fact that the positive screening tests were assessed by nodule volume and not size as many other studies.

One strength of our meta-analysis is the design and quality of our study; we conducted a systematic review following the PRISMA guidelines with two double-blind independent reviewers to avoid missing publications. Quality assessment of each eligible trial and data extraction were also made by two double-blind independent readers to prevent data errors. Other strengths included inclusion in our meta-analysis of only controlled randomized studies with precise inclusion and exclusion criteria. Also, the studies included had long-term follow-up and we only included final analyses results (no intermediate or preliminary results). The symmetric distribution of the relative risk in the funnel plot confirms the minimal publication bias among the studies included. Last, our meta-analysis included a large final total number of patients (n=84,558).

Limitations of our study include the partial heterogeneity of the protocols studied, in particular the interventions in the control arm, either prevention or clinical examination in all the studies except the NLST (annual CXR), and DANTE (CXR at inclusion) studies. Also, heterogeneity among studies concerned the smoking history of patients included, much higher in the NLST study than in other RCTs included. low-dose CT screening according to the risk of lung-cancer death. N Engl J Med. [START_REF] Wille | Results of the Randomized Danish Lung Cancer Screening Trial with Focus on High-Risk Profiling[END_REF]

Figure 1 :Figure 2 :Figure 3 :Figure 4 :

 1234 Figure 1: PRISMA flowdiagram of the meta-analysis Figure 2: Funnel plot reporting relative risk for overall mortality Figure 3: Forrest plot of relative risk for overall mortality Figure 4: Forrest plot of relative risk for lung cancer specific mortality

  X-ray / 95% CI: 95% confidence interval / NR: not reported i 5 for the annual arm and 3 for the biennial arm; p-yr: person years; yrs: years

Table 1 :

 1 Details of the 7 Randomized Clinical Trials included

	juill

Table 2 :

 2 Quality assessment results

Table 1

 1 Randomized Controlled Trials included in the meta-analysis

	Characteristics	DANTE	DLSCT	ITALILUNG	MILD	LUSI	NLST	NELSON
		2015	2015	2017	2019	2019	2013	2020
		(1)	(24, 29)	(4)	(5)	(6,7)	(8-10)	(11)
	Country	Italy	Denmark	Italy	Italy	Germany	USA	Belgium /
								Netherland
	Control arm	CXR at inclusion	Annual clinical	Primary	Primary	Primary	Annual CXR	Primary
			examination	prevention for	prevention for	prevention		prevention for
				smoking	smoking	for smoking		smoking cessation
				cessation	cessation	cessation
	Total sample size	2450	4104	3206	4099	4052	53452	13195
	LDCT arm	1264	2052	1613	2376	2029	26722	6583
	Control arm	1186	2052	1593	1723	2023	26730	6612
	Age, yrs	60-74	50-70	55-69	≥49	50-69	55-74	50-74
	Men, %	100	55	65	66	65	59	100
	Median follow-up, yrs	8.35	Intermediate: 4.8	9.3	4.4	3	6.5	10
			Final: 9.8
	Screening, p/yr	10875	19439	14658	6449.5	/	144103	/
	Control, p/yr	10104	19547	14247	5556.7	/	143368	/
	Smoking history	Mean [95% CI]:	Mean [95% CI]:	Mean:	Median:	NR	Mean:	Mean:
	LDCT, pack-yrs	47.3 [45.7 -49]	36.4 [23 -49.8]	40	39		56.04	38.0 [29.7-49.5]
	Controls, pack-yrs	47.2 [45.5 -49]	35.9 [22.5 -49.3]	38	38		55.93	38.0 [29.7-49.5]
	Current smoker, n (%)	1395 (57%)	3124 (76%)	2077 (65%)	3176 (77%)	2506 (62%)	25779 (48%)	7254 (55%)
	Former Smoker, n (%)	1055 (43%)	980 (24%)	1129 (35%)	923 (23%)	1545 (38%)	27677 (52%)	5941 (45%)
	Screening rounds, n	5	5	4	5 or 3 i

Table 2 :

 2 Quality assessment (CONSORT Check List) of the seven randomized controlled trials included in the meta-analysis.

	Quality criteria	DANTE DLSCT ITALILUNG MILD LUSI	NLST	NELSON
	Recruitment	Yes	Yes	Yes	Yes	Yes	Yes	Yes
	Strategies well
	defined
	Complete	Yes	Yes	Yes	Yes	Yes	Yes	Yes
	diagnostic work-
	up planned
	Random	Yes	Yes	Yes	Yes	Yes	Yes	Yes
	assignment
	Inclusion criteria	Yes	Yes	Yes	Yes	Yes	Yes	Yes
	described and
	respected
	Measure of all-	Yes	Yes	Yes	Yes	Yes	Yes	Yes
	cause mortality
	Measure of	Yes	Yes	Yes	Yes	No	Yes	Yes
	specific mortality
	Valid	Yes	Yes	Yes	Yes	Yes	Yes	Yes
	measurement of
	mortality
	Blinded of	Yes	NR	Yes	Yes	NR	Yes	Yes
	outcome
	Assessment
	Long enough	Yes	Yes	Yes	No	Yes	Yes	Yes
	follow-up
	NR : not reported

Acknowledgements

We would like to thank Sarah Kabani for valuable editorial assistance and Quentin Durand for the study selection. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Conclusions

Our meta-analysis is the first systematic review to include all recent RCTs including the recent NELSON study. Our results confirm that of the NELSON study, showing an impact of lung cancer screening on lung cancer-specific mortality, reduced in the LDCT group. No impact of such screening on all-cause mortality was reported.

Conflict of Interest statement

None declared


Didier Lomet

Xavier Druart

David Hazlerigg

Massimiliano Beltramo

Hugues Dardente
email: hugues.dardente@inrae.fr

Circuit-level analysis identifies target genes of sex steroids in ewe seasonal breeding

Keywords: Biological rhythms, breeding, circannual clock, estradiol, GnRH, kisspeptin, Kiss1, melatonin, pars tuberalis, photoperiod, pituitary, progesterone, RFRP3, seasonality, sheep, tanycytes, thyrotropin estradiol; IdF: Ile-de-France breed;

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The pars tuberalis (PT) of the pituitary is the key target tissue of melatonin for the control of seasonal breeding [START_REF] Dardente | Thyroid hormone and seasonal rhythmicity[END_REF]Dardente et al 2019a;[START_REF] Wood | The pars tuberalis: The site of the circannual clock in mammals?[END_REF]. In the PT, a short duration of nightly melatonin production, typical of longer spring/summer days, triggers the expression of Tshb through TEF/EYA3/SIX1, a triad of transcription factors and co-activators [START_REF] Dardente | A molecular switch for photoperiod responsiveness in mammals[END_REF]. This specific PT-derived TSH acts on neighboring tanycytes, which line the walls of the third ventricle, to induce the expression of deiodinase II (Dio2) and repress expression of Dio3, leading to increased local T3 production [START_REF] Nakao | Thyrotrophin in the pars tuberalis triggers photoperiodic response[END_REF][START_REF] Hanon | Ancestral TSH mechanism signals summer in a photoperiodic mammal[END_REF][START_REF] Ono | Involvement of thyrotropin in photoperiodic signal transduction in mice[END_REF]. This opposite photoperiodic control over Dio2 and Dio3 is responsible for the tight seasonal control of T3 levels within the medio-basal hypothalamus (MBH), which governs physiological switches in the annual breeding cycle of multiple species, including sheep [START_REF] Karsch | Involvement of thyroid hormones in seasonal reproduction[END_REF][START_REF] Dardente | Melatonin-dependent timing of seasonal reproduction by the pars tuberalis: pivotal roles for long daylengths and thyroid hormones[END_REF]).

Our recent RNAseq analysis of the MBH in ovariectomized estradiol-implanted (OVX+E2) ewes revealed that ~10% of protein-coding genes display differential expression levels between May and November, which illustrates the potent transcriptome-wide impact of photoperiod (Lomet et al 2018;[START_REF] Dardente | Photoperiod and thyroid hormone regulate expression of l-dopachrome tautomerase (Dct), a melanocyte stem-cell marker, in tanycytes of the ovine hypothalamus[END_REF]. We validated seasonal changes in the expression level for the top ~30 differentially expressed genes (DEG). Most DEG were found either in the PT (e.g. Tshb,Fam150b,Vmo1,Ezh2,Suv39H2,Vcan, or in tanycytes of the beta subtypes, located at the bottom of the third ventricle (e.g. Dio2, Tmem252, Shh, SlcO1c1, Dct and NpSR1).

Kiss1 and Npvf (encoding the RF-amide peptides KISS1 and RFRP3), which are expressed by discrete clusters of neurons within the MBH parenchyma [START_REF] Smith | Sex steroid control of hypothalamic Kiss1 expression in sheep and rodents: comparative aspects[END_REF][START_REF] Lomet | The impact of thyroid hormone in seasonal breeding has a restricted transcriptional signature[END_REF][START_REF] Angelopoulou | Functional Implications of RFRP-3 in the Central Control of Daily and Seasonal Rhythms in Reproduction[END_REF], were also identified as DEG. The MBH is key to the central control of breeding and appears to host the long sought-after GnRH pulse generator. There is now strong evidence that KNDy (KISS1/NeurokininB/Dynorphin) neurons located in the arcuate nucleus are the GnRH pulse generator [START_REF] Lehman | Minireview: kisspeptin/neurokinin B/dynorphin (KNDy) cells of the arcuate nucleus: a central node in the control of gonadotropin-releasing hormone secretion[END_REF][START_REF] Han | Selective optogenetic activation of arcuate kisspeptin neurons generates pulsatile luteinizing hormone secretion[END_REF][START_REF] Moore | KNDy Cells Revisited[END_REF][START_REF] Herbison | The Gonadotropin-Releasing Hormone Pulse Generator[END_REF].

We found that Kiss1 displays a marked (~10-fold) seasonal rhythm in expression in OVX+E2 ewes, leading to heightened expression during short winter days [START_REF] Lomet | The impact of thyroid hormone in seasonal breeding has a restricted transcriptional signature[END_REF], consistent with prior findings [START_REF] Smith | Variation in kisspeptin and RFamide-related peptide (RFRP) expression and terminal connections to gonadotropin-releasing hormone neurons in the brain: a novel medium for seasonal breeding in the sheep[END_REF][START_REF] Wagner | Redefining the limits of day length responsiveness in a seasonal mammal[END_REF][START_REF] Urias-Castro | Hypothalamic concentrations of kisspeptin and gonadotropin-releasing hormone during the breeding season and non-breeding season in ewes[END_REF]. In intact ewes, this increase in KISS1 presumably leads to enhanced GnRH secretion, which triggers breeding [START_REF] Messager | Kisspeptin directly stimulates gonadotropin-releasing hormone release via G protein-coupled receptor 54[END_REF][START_REF] Caraty | Kisspeptin synchronizes preovulatory surges in cyclical ewes and causes ovulation in seasonally acyclic ewes[END_REF][START_REF] Caraty | Gonadotrophin-releasing hormone release into the hypophyseal portal blood of the ewe mirrors both pulsatile and continuous intravenous infusion of kisspeptin: an insight into kisspeptin's mechanism of action[END_REF]. Furthermore, in sheep as in other mammals a vast majority of KNDy neurons express receptors for E2 and P4 (ERa and PGR, respectively ; [START_REF] Goubillon | Identification of neurokinin B-expressing neurons as an highly estrogen-receptive, sexually dimorphic cell group in the ovine arcuate nucleus[END_REF][START_REF] Dufourny | Progesterone receptor, estrogen receptor alpha, and the type II glucocorticoid receptor are coexpressed in the same neurons of the ovine preoptic area and arcuate nucleus: a triple immunolabeling study[END_REF][START_REF] Franceschini | Kisspeptin immunoreactive cells of the ovine preoptic area and arcuate nucleus co-express estrogen receptor alpha[END_REF][START_REF] Smith | KiSS-1 messenger ribonucleic acid expression in the hypothalamus of the ewe is regulated by sex steroids and season[END_REF][START_REF] Foradori | Colocalization of progesterone receptors in parvicellular dynorphin neurons of the ovine preoptic area and hypothalamus[END_REF][START_REF] Campbell | A molecular census of arcuate hypothalamus and median eminence cell types[END_REF][START_REF] Chen | Single-Cell RNA-Seq Reveals Hypothalamic Cell Diversity[END_REF], which confers them direct responsiveness to sex steroid feedback. This is highly relevant since the current model for the photoperiodic control of seasonal breeding, established in the early 80's [START_REF] Legan | The endocrin control of seasonal reproductive function in the ewe: a marked change in response to the negative feedback action of estradiol on luteinizing hormone secretion[END_REF][START_REF] Karsch | Neuroendocrine basis of seasonal reproduction[END_REF]Goodman & Inskeep, 2015), emphasizes the role of "a seasonal switch in the ability of E2 to inhibit the GnRH pulse generator" as the key event responsible for timely onset and offset of breeding. Indeed, OVX+E2 ewes display overt seasonal differences in the frequency and amplitude of GnRH/LH pulses, which are lost in OVX ewes [START_REF] Karsch | Seasonal changes in gonadotropin-releasing hormone secretion in the ewe: alteration in response to the negative feedback action of estradiol[END_REF]. The neural mechanism(s) driving the seasonal switch in the ability of E2 to exert negative feedback remain unclear even though the MBH appears to be the key neural site [START_REF] Blache | Ventromedial hypothalamus as a target for oestradiol action on proceptivity, receptivity and luteinizing hormone surge of the ewe[END_REF][START_REF] Caraty | Evidence that the mediobasal hypothalamus is the primary site of action of estradiol in inducing the preovulatory gonadotropin releasing hormone surge in the ewe[END_REF]. Since OVX+E2 ewes undergo properly timed seasonal switches in LH/FSH, P4 appears dispensable for these [START_REF] Legan | The endocrin control of seasonal reproductive function in the ewe: a marked change in response to the negative feedback action of estradiol on luteinizing hormone secretion[END_REF]. However, P4 exerts negative feed-back onto the MBH during the luteal phase, which is mandatory for normal estrus cycles and fertility during the breeding season [START_REF] Dardente | Melatonin-dependent timing of seasonal reproduction by the pars tuberalis: pivotal roles for long daylengths and thyroid hormones[END_REF]Goodman & Inskeep, 2015).

Even though KNDy neurons are strong contenders as a common conduit for photoperiod and sex steroids towards GnRH seasonal control, we lack a comprehensive circuit-level analysis of the cellular and molecular targets of E2 (and P4) responsible for the "seasonal switch". Indeed several non-mutually exclusive hypotheses can be put forth: E2 may impair the photoperiodic read-out mechanism in the PT, and/or modulate transduction of the signal towards tanycytes and/or modulate signaling directly within the hypothalamic parenchyma. The latter possibility would likely involve Kiss1 neurons, but might also include an impact on Npvf neurons. Here, we used our validated set of seasonal DEG to test these possibilities in ewes sampled at different times of the year and under various sex steroid states (intact, OVX+E2 and OVX ewes).

Materials & Methods

Ethics statement

All experimental procedures were performed in accordance with international (directive 2010/63/UE) and national legislation (décret n° 2013-118) governing the ethical use of animals in research (authorization n° E37-175-2 and n°A38 801). All procedures used in this work were evaluated by a local ethics committee (Comité d'Ethique en Expérimentation Animale Val de Loire) and approved by the Ministry of Higher Education and Research (project n°00710.02). All surgeries were performed after sodium thiopental anesthesia (Nesdonal®, 1g/80kg) under constant isoflurane administration (Vetflurane®) and all efforts were made to minimize suffering. Following surgery, animals received an injection of antibiotics (oxytetracycline, Terramycine LA®, 1ml/10kg) and an injection of a non-steroidal anti-inflammatory drug (Finadyne®, flumixin megumine, 2ml/50kgs). Animals were taken care of daily throughout the experiment and killed by decapitation under deep barbiturate anesthesia (Nesdonal®, 5mL). To minimize potential issues linked to time-of-day fluctuation in gene expression (i.e. circadian rhythms) all animals were killed in the early day (ZT1-6, with ZT0 being the time of lights on, or sunrise).

Animals & Experimental procedures

Experiments were conducted on adult Ile-de-France (IdF) ewes (3-5 years old; weight 60-80 kg) and adult Romanov (Rom) ewes (3-5 years old; weight 60- S1.

In situ hybridization

ISH was performed as described previously [START_REF] Lomet | The impact of thyroid hormone in seasonal breeding has a restricted transcriptional signature[END_REF]Dardente et al 2019b).

Hypothalamic blocks for were cut into 20µm sections using a cryostat (CryoStar NX70, ThermoScientific) and thaw-mounted onto SuperFrost Plus slides (ThermoScientific). Specific probes were generated and validated (cloning & sequencing) as described above, using a mix of MBH cDNA as template. Details of the inserts used to generate radioactive probes are provided in Table S1. Radioactive cRNA riboprobes were prepared by plasmid linearisation and in vitro transcription (Riboprobe System, Promega) including 35 S-UTP (Perkin-Elmer). Probes were purified using Illustra Probe Quant G50 micro-columns (Fisher) and counted with a liquid scintillation counter (Tri-Carb 2900TR, Packard). Slides were post-fixed at 4°C for 20 min in 4% PFA, 0.1 M PB, rinsed with 0.1 M PB (2 X 5min), acetylated with 3.75% v/v of acetic anhydride in 0.1 TEA, 0.05 N NaOH and finally rinsed with 0.1 M PB (2 X 5min). Slides were then dehydrated through graded ethanol solutions (50%, 70%, 95% and 100%; 3min each) and dried under vacuum for 60 min. Sections were hybridized overnight at 58ºC with 10 6 cpm of probe per slide in hybridization buffer (50% deionized formamide, 10% dextran sulfate, 1 X Denhardt's solution, 300 mM NaCl, 10 mM Tris, 10 mM DTT, 1 mM EDTA, 500 µg/mL tRNA). Sections were then rinsed in 4 X SSC (3 X 5 min) and subjected to RNase-A digestion (20 µg/mL) in a buffer containing 500 mM NaCl, 1 mM Tris, 1 mM EDTA for 30 min at 37°C. Stringency washes in SSC (with 1mM DTT) were performed to remove non-specific probe hybridisation: 2 X SSC (2 X 5 min), 1 X SSC (10 min), 0.5 X SSC (10 min), 0.1 X SSC (30 min at 60°C), 0.1 X SSC (5 min).

Slides were then dehydrated through graded ethanol solutions, dried under vacuum for 60 min and exposed to an autoradiographic film (BioMax MR, Kodak). Exposure duration was optimized for each gene by repeated film exposures, depending on labeling intensity (from 3 days to 4 weeks).

Films were scanned on a transmittance image scanner (Amersham, UK) along with a calibrated optical density (OD) transmission step wedge (Stouffer, USA). Calibrated Integrated OD (IOD) measurements of gene expression within the MBH were performed using ImageJ software.

Hormonal profile in OVX ewes

Plasma levels of LH in OVX ewes (experiment 3) were assayed by RIA [START_REF] Pelletier | Seasonal variation in LH and testosterone release in rams of two breeds[END_REF]. All samples were included in a single assay and every sample was measured in duplicate. The assay standard was 1051-CY-LH (equivalent to 0.31 NIH-LH-S1). Intra-and inter-assay coefficients of variation averaged 9% and 15%, respectively with an assay sensitivity of 0.1 ng/mL.

Data analysis

Data were analysed using GraphPad Prism 6. ISH data were analysed by t-test (Mann-Whitney).

Data for qRT-PCR were analysed by one-way ANOVA (followed by post-hoc Tukey test when applicable; Expts 2 and 4) or two-way ANOVA (Expt 3) using time (May and November) and treatment (intact, OVX+E2 and OVX) as variables. The Sidak's post-hoc test was used for multiple comparisons.

Results

We first used ISH to compare the seasonal expression pattern of candidate genes in the MBH of intact ewes to the pattern we previously reported for OVX+E2 ewes (Lomet et Intact and OVX+E2 ewes display roughly similar basal levels of E2 throughout the year [START_REF] Yuthasastrakosol | Luteinizing hormone, oestrogen and progesterone levels in peripheral serum of anoestrous and cyclic ewes as determined by radioimmunoassay[END_REF] but only intact ewes are further exposed to P4 during the breeding season (Goodman & Inskeep 2015). We therefore went on to consider the importance of P4 for the seasonal pattern of expression of Kiss1 and Npvf, reasoning that the impact of cycling P4 might somehow superimpose to the tonic E2-dependent photoperiodic impact on their expression level.

We attempted to get circuit-level insights on the potential impact of P4 by probing expression of tanycytes themselves harbor circannual clocks can not be ruled out [START_REF] Milesi | Downregulation of Deiodinase 3 is the earliest event in photoperiodic and photorefractory activation of the gonadotropic axis in seasonal hamsters[END_REF]Lomet et al 2018). The circannual clock is predicted to be sheltered from cues other than melatonin in order to consistently provide reliable timing information to downstream circuits.

We and others previously showed that gene expression in the PT is unresponsive to TRH and TH (Bockman et al 1997 ;[START_REF] Dardente | Melatonin-dependent timing of seasonal reproduction by the pars tuberalis: pivotal roles for long daylengths and thyroid hormones[END_REF][START_REF] Lomet | The impact of thyroid hormone in seasonal breeding has a restricted transcriptional signature[END_REF]. We now show that seasonal gene expression in the PT is also unaffected by sex steroids. Therefore, neither TH nor E2, the two key hormones in the progression of the seasonal breeding cycle, have any noticeable impact on photoperiod decoding in the PT. We believe this is a strong argument in favor of the view that melatonin-responsive PT-specific thyrotropes (Klosen et Tmem252,Shh,SlcO1c1,Dct) or repressed (NpSR1) by lengthening photoperiods. We speculate that functional relevance of this P4 feedback at the level of tanycytes, if any, may be more obvious when the strength of the photoperiodic input is comparatively weak, which would correspond to physiological transitions to anestrus or breeding. In the latter case, since P4 can not increase unless ovulation occurs, P4 is unlikely to play any role in timing breeding onset. However, P4 may modulate the timing of anestrus. However, we acknowledge that only Tmem252 ans Shh in IdF ewes do show a statistically significant decrease upon P4 exposure. Once more, since OVX+E2 ewes undergo properly timed seasonal switches in LH/FSH, P4 is clearly dispensable for these [START_REF] Legan | The endocrin control of seasonal reproductive function in the ewe: a marked change in response to the negative feedback action of estradiol on luteinizing hormone secretion[END_REF]. Therefore, great caution is required in the interpretation of the general trends observed here, and the notion that P4 may exert a modulatory role will have to be tested further.

We notice that expression of Shh is affected by E2 (see and PGR are expressed in cells lining the third ventricle [START_REF] Lehman | Distribution of estrogen receptor-immunoreactive cells in the sheep brain[END_REF][START_REF] Scott | The distribution of progesterone receptor immunoreactivity and mRNA in the preoptic area and hypothalamus of the ewe: upregulation of progesterone receptor mRNA in the mediobasal hypothalamus by oestrogen[END_REF].

Whether these cells are tanycytes is unknown, and co-localization with tanycyte-specifc markers would be required. In mouse, tanycytes appear to express ERa and PGR and display some responses to P4 [START_REF] Parkash | Semaphorin7A regulates neuroglial plasticity in the adult hypothalamic median eminence[END_REF]. However, a recent scRNAseq study failed to detect ERa or PGR in murine tanycytes [START_REF] Campbell | A molecular census of arcuate hypothalamus and median eminence cell types[END_REF]. expression of these two genes is governed by photoperiod and sex steroids. There is strong evidence that "final processes governing the onset of puberty in the lamb and the onset of the annual season in the ewe share a common mechanism", based on a decrease in the ability of E2 to inhibit the GnRH pulse generator [START_REF] Legan | The endocrin control of seasonal reproductive function in the ewe: a marked change in response to the negative feedback action of estradiol on luteinizing hormone secretion[END_REF][START_REF] Foster | Endocrine mechanisms governing transition into adulthood: a marked decrease in inhibitory feedback action of estradiol on tonic secretion of luteinizing hormone in the lamb during puberty[END_REF][START_REF] Karsch | Neuroendocrine basis of seasonal reproduction[END_REF]Goodman & Inskeep, 2015). In this respect, Kiss1-expressing neurons of the arcuate nuclei are prime candidate [START_REF] Beltramo | Cellular mechanisms and integrative timing of neuroendocrine control of GnRH secretion by kisspeptin[END_REF][START_REF] Pinilla | Kisspeptins and reproduction: physiological roles and regulatory mechanisms[END_REF]Simonneaux 2018).

We show that Kiss1 is expressed at basal low levels throughout the year in intact ewes, which contrasts with the large May/Nov difference seen in OVX+E2 ewes. While E2 is present at roughly similar levels throughout the year in intact ewes, not overlooking the functionally critical but modest 2-3 fold rise during the late follicular phase, P4 is produced exclusively during the breeding season [START_REF] Yuthasastrakosol | Luteinizing hormone, oestrogen and progesterone levels in peripheral serum of anoestrous and cyclic ewes as determined by radioimmunoassay[END_REF]Goodman & Inskeep 2015). Therefore, sex steroid status during the anestrus season -presence of E2 and lack of P4 -can be considered similar in intact and OVX+E2 ewes. This suggests that it is the presence of P4 during the breeding season, which actively down-regulates the photoperiodic drive upon Kiss1 expression. On the other hand, OVX ewes display intermediate-to-high levels of Kiss1 in both May and November, consistent with prior findings [START_REF] Pompolo | Colocalization of kisspeptin and gonadotropin-releasing hormone in the ovine brain[END_REF][START_REF] Smith | KiSS-1 messenger ribonucleic acid expression in the hypothalamus of the ewe is regulated by sex steroids and season[END_REF][START_REF] Smith | Variation in kisspeptin and RFamide-related peptide (RFRP) expression and terminal connections to gonadotropin-releasing hormone neurons in the brain: a novel medium for seasonal breeding in the sheep[END_REF][START_REF] Smith | Kisspeptin neurons in the ovine arcuate nucleus and preoptic area are involved in the preovulatory luteinizing hormone surge[END_REF][START_REF] Merkley | KNDy (kisspeptin/neurokinin B/dynorphin) neurons are activated during both pulsatile and surge secretion of LH in the ewe[END_REF]Goodman & Inskeep, 2015). As assesssed by c-Fos labeling, OVX is accompanied by activation of KNDy neurons -but not of Kiss1 neurons of the preoptic areawhich suggests that only Kiss1 neurons in the arcuate nucleus are important for the E2 negative feedback [START_REF] Merkley | KNDy (kisspeptin/neurokinin B/dynorphin) neurons are activated during both pulsatile and surge secretion of LH in the ewe[END_REF]. Collectively, seasonal expression of Kiss1 reflects E2-mediated inhibition during anestrus, while E2-independent increase during the breeding season would be masked by P4 inhibition during the luteal phase, and unmasked during the follicular phase. There is prior evidence that both E2 and, to a lesser extent, P4 exert negative feedback on Kiss1 expression [START_REF] Smith | KiSS-1 messenger ribonucleic acid expression in the hypothalamus of the ewe is regulated by sex steroids and season[END_REF][START_REF] Smith | Variation in kisspeptin and RFamide-related peptide (RFRP) expression and terminal connections to gonadotropin-releasing hormone neurons in the brain: a novel medium for seasonal breeding in the sheep[END_REF][START_REF] Goodman | Evidence that the arcuate nucleus is an important site of progesterone negative feedback in the ewe[END_REF]Goodman & Inskeep, 2015).

We acknowledge that our experimental paradigm is not suitable to reveal a potential negative impact of P4 upon Kiss1 expression (Fig 4), which is decreasing at the transition towards anestrus due to lengthening photoperiod.

Following a line of reasoning similar to the one developed above for Kiss1, we conclude that seasonal expression of Npvf reflects E2-mediated activation during anestrus, while neither E2 nor P4 impinge on Npvf expression during the breeding season. Consistent with our findings in intact IdF ewes, we previously reported that photoperiod had little impact on Npvf expression in intact Soay ewes [START_REF] Dardente | RFamide-related peptide and its cognate receptor in the sheep: cDNA cloning, mRNA distribution in the hypothalamus and the effect of photoperiod[END_REF]. Furthermore, Npvf expression is not impacted by the transition to long days in either intact Soay ewes [START_REF] Dardente | RFamide-related peptide and its cognate receptor in the sheep: cDNA cloning, mRNA distribution in the hypothalamus and the effect of photoperiod[END_REF] or OVX+E2 IdF ewes [START_REF] Lomet | The impact of thyroid hormone in seasonal breeding has a restricted transcriptional signature[END_REF]. Since such acute photoperiodic transition almost immediately impacts the gonadotropic axis, we conclude that changes in Npvf transcription are not necessary for this physiological response, a conclusion which holds true also for Kiss1 [START_REF] Lomet | The impact of thyroid hormone in seasonal breeding has a restricted transcriptional signature[END_REF]. Interestingly, the number of RFRP3-ir cells is 40% higher during the non-breeding season in OVX+E2 ewes [START_REF] Smith | Variation in kisspeptin and RFamide-related peptide (RFRP) expression and terminal connections to gonadotropin-releasing hormone neurons in the brain: a novel medium for seasonal breeding in the sheep[END_REF]. However, the expression level of Npvf was similar across breeding and non-breeding seasons, in both OVX and OVX+E2 ewes [START_REF] Smith | Variation in kisspeptin and RFamide-related peptide (RFRP) expression and terminal connections to gonadotropin-releasing hormone neurons in the brain: a novel medium for seasonal breeding in the sheep[END_REF], which is not congruent with our current findings. The reasons for this are unclear but may relate to the methods used (nonradioactive ISH vs qRT-PCR).

Compared to KNDy neurons, the evidence for a direct control of Npvf-expressing neurons by sex steroids is scarce (Goodman & Inskeep, 2015). In hamsters, about 40% of Npvf neurons express ERa and E2 injection in OVX animals induces c-Fos in a subset of Npvf cells [START_REF] Kriegsfeld | Identification and characterization of a gonadotropin-inhibitory system in the brains of mammals[END_REF]. In mouse, however, only 20-25% of Npvf neurons co-express ERa as assessed by ISH [START_REF] Molnar | Estradiol down-regulates RF-amide-related peptide (RFRP) expression in the mouse hypothalamus[END_REF][START_REF] Poling | Development, sex steroid regulation, and phenotypic characterization of RFamide-related peptide (Rfrp) gene expression and RFamide receptors in the mouse hypothalamus[END_REF], while scRNAseq provides no evidence for co-expression [START_REF] Chen | Single-Cell RNA-Seq Reveals Hypothalamic Cell Diversity[END_REF]. In sheep, we also found no evidence for co-expression (our unpublished observations) while others reported that ~20% of Npvf neurons may express ERa (unpublished observations, cited in [START_REF] Clarke | New concepts of the central control of reproduction, integrating influence of stress, metabolic state, and season[END_REF]. These data are consistent with the overall distinct distribution patterns of Npvf- [START_REF] Dardente | RFamide-related peptide and its cognate receptor in the sheep: cDNA cloning, mRNA distribution in the hypothalamus and the effect of photoperiod[END_REF] and ERa-expressing cells in the ovine MBH [START_REF] Lehman | Distribution of estrogen receptor-immunoreactive cells in the sheep brain[END_REF][START_REF] Blache | Oestrogen receptors in the preoptico-hypothalamic continuum: immunohistochemical study of the distribution and cell density during induced oestrous cycle in ovariectomized ewe[END_REF]. Npvf neurons are most numerous at the level of the VMH/DMH, where they tend to cluster medially, close to the third ventricle, while ERa-ir cells are few at this level and rather scattered in lateral regions, farther from the third ventricle. In sheep, the distribution of PGR does not provide support for co-expression with Npvf, since labeled nuclei are not found in the DMH/VMH [START_REF] Goubillon | Identification of neurokinin B-expressing neurons as an highly estrogen-receptive, sexually dimorphic cell group in the ovine arcuate nucleus[END_REF][START_REF] Dufourny | Progesterone receptor, estrogen receptor alpha, and the type II glucocorticoid receptor are coexpressed in the same neurons of the ovine preoptic area and arcuate nucleus: a triple immunolabeling study[END_REF][START_REF] Smith | KiSS-1 messenger ribonucleic acid expression in the hypothalamus of the ewe is regulated by sex steroids and season[END_REF][START_REF] Foradori | Colocalization of progesterone receptors in parvicellular dynorphin neurons of the ovine preoptic area and hypothalamus[END_REF]. These neuroanatomical considerations suggest that Npvf-expressing neurons are not primary targets of sex steroids, even though the possibility that ERa/PGR are expressed at low levels (i.e. below the detection threshold) can not be excluded.

Overall, the opposite expression patterns for Kiss1 and Npvf in intact, OVX+E2 and OVX ewes is striking. We suggest two non-mutually exclusive scenarios to explain this : sex steroids have independent and opposite actions on both genes, or the two neuronal populations are somehow connected such that the impact of sex steroids on one population leads to inhibition of expression in the other. In both scenarios, the interplay between T3 and E2 is required to govern proper seasonal timing [START_REF] Lincoln | Seasonal breeding: nature's contraceptive[END_REF][START_REF] Karsch | Neuroendocrine basis of seasonal reproduction[END_REF][START_REF] Karsch | Involvement of thyroid hormones in seasonal reproduction[END_REF][START_REF] Dardente | Melatonin-dependent timing of seasonal reproduction by the pars tuberalis: pivotal roles for long daylengths and thyroid hormones[END_REF][START_REF] Dardente | Thyroid hormone and seasonal rhythmicity[END_REF]. The idea that Npvf-expressing neurons might be the primary target for photoperiodic control has already been put forth in hamsters [START_REF] Henningsen | RFRP Neurons -The Doorway to Understanding Seasonal Reproduction in Mammals[END_REF][START_REF] Angelopoulou | Functional Implications of RFRP-3 in the Central Control of Daily and Seasonal Rhythms in Reproduction[END_REF]. However, we deem this unlikely for several reasons. As discussed above, KNDy cells do directly respond to sex steroids (presence of ERa and PGR)

and also contain the thyroid hormone receptor TRa [START_REF] Dufourny | The Content of Thyroid Hormone Receptor alpha in Ewe Kisspeptin Neurones is not Season-Dependent[END_REF] ; they are therefore well-placed to mediate the T3-and E2-dependent switches in seasonal breeding. Anatomically, KNDy cells are also likely to receive contacts from beta-tanycytes processes, which traverse the arcuate nucleus and reach the median eminence [START_REF] Rodriguez | Hypothalamic tanycytes: a key component of brain-endocrine interaction[END_REF]Prévot et al 2018 ;[START_REF] Rodriguez | Tanycytes: A rich morphological history to underpin future molecular and physiological investigations[END_REF]. Indeed, KNDy cells might also receive signals, including T3, directly from PT cells through the peculiar micro-vasculature of the arcuate nucleus region [START_REF] Clarke | Hypothalamus as an endocrine organ[END_REF]. Finally, the evidence that KISS1 is the key activator of reproduction in mammals is overwhelming [START_REF] Beltramo | Cellular mechanisms and integrative timing of neuroendocrine control of GnRH secretion by kisspeptin[END_REF][START_REF] Pinilla | Kisspeptins and reproduction: physiological roles and regulatory mechanisms[END_REF]Simonneaux 2018), while evidence implicating RFRP3 is weak [START_REF] Leon | Dissecting the Roles of Gonadotropin-Inhibitory Hormone in Mammals: Studies Using Pharmacological Tools and Genetically Modified Mouse Models[END_REF][START_REF] Angelopoulou | Functional Implications of RFRP-3 in the Central Control of Daily and Seasonal Rhythms in Reproduction[END_REF].

In sheep, KISS1 and analogues consistently trigger GnRH/LH release and ovulation, even during the anestrus season [START_REF] Messager | Kisspeptin directly stimulates gonadotropin-releasing hormone release via G protein-coupled receptor 54[END_REF][START_REF] Caraty | Kisspeptin synchronizes preovulatory surges in cyclical ewes and causes ovulation in seasonally acyclic ewes[END_REF][START_REF] Beltramo | Rational design of triazololipopeptides analogs of kisspeptin inducing a long-lasting increase of gonadotropins[END_REF]Decourt et al 2016b;[START_REF] Beltramo | Towards new strategies to manage livestock reproduction using kisspeptin analogs[END_REF]. In comparison, we could not find any evidence for a role of RFRP3 in GnRH/LH secretion (Decourt et al 2016a), which contrasts with earlier findings [START_REF] Clarke | Potent action of RFamide-related peptide-3 on pituitary gonadotropes indicative of a hypophysiotropic role in the negative regulation of gonadotropin secretion[END_REF]. Genetic ablation of the RFRP3 receptor (NPFFR1) in mouse has negligible impact upon reproduction [START_REF] Leon | Physiological roles of gonadotropin-inhibitory hormone signaling in the control of mammalian reproductive axis: studies in the NPFF1 receptor null mouse[END_REF], while absence of KISS1R (a.k.a.

GPR54) leads to infertility [START_REF] Messager | Kisspeptin directly stimulates gonadotropin-releasing hormone release via G protein-coupled receptor 54[END_REF][START_REF] Pinilla | Kisspeptins and reproduction: physiological roles and regulatory mechanisms[END_REF][START_REF] Beltramo | Cellular mechanisms and integrative timing of neuroendocrine control of GnRH secretion by kisspeptin[END_REF].

Finally, there are only minimal projections from RFRP3 cells towards the arcuate nucleus, which makes it unlikely that RFRP3 neurons make synaptic input onto KISS1 neurons in sheep [START_REF] Qi | Projections of RFamide-related peptide-3 neurones in the ovine hypothalamus, with special reference to regions regulating energy balance and reproduction[END_REF]. However, some evidence for such a unidirectional control has been found in mouse [START_REF] Poling | Kisspeptin neurones do not directly signal to RFRP-3 neurones but RFRP-3 may directly modulate a subset of hypothalamic kisspeptin cells in mice[END_REF]. Therefore, we can not exclude a primary role for Npvf-expressing neurons in the seasonal response, but current evidence in sheep rather favor a model in which these neurons are located downstream of KNDy neurons.

In conclusion our findings reveal the circuit-level logic of the impact of sex steroid of seasonal breeding (Fig 5). The PT behaves as would be expected for a true circannual oscillator, sheltered from the impact of sex steroids and TH, but acutely sensitive to photoperiod through melatonin.

Tanycytes display gene-specific regulation by both TH, E2 and P4, which we believe might help fine-tune the timing for the onset and offset of the breeding season. Finally, neurons expressing Kiss1 emerge as a key node for the integration of environmental (photoperiod, through T3 responsiveness) and internal (feedback to sex steroid) cues, which together drive seasonal breeding. Our analysis strongly suggests that Kiss1 neurons of the arcuate nuclei are the cellular substrate for the "seasonal switch in the ability of E2 to inhibit the GnRH pulse generator". Finally the opposite transcriptional control of Kiss1 and Npvf expression over distinct states of sex steroid feedback suggest some level of interconnection between these neuronal populations, the nature of which remains to be identified.

Supplemental Materials

Figures S1-S5

  photoperiod-driven DEG in intact ewes at different phases of the estrous cycle (follicular vs luteal)

  Fig 4A), which implies that E2 and P4 impact upon tanycytes in a gene-specific manner. The

  Further experiments are required to clarify expression of ERa and PGR in tanycytes and define whether they play any physiologically relevant role. Considering the impact of tanycyte-derived T3 on the expression of Kiss1 and Npvf in sheep (Lomet et al 2018) and hamsters (Henson et al 2013 ; Klosen et al 2013), a potential role for sex steroids in regulating DEG within tanycytes is of interest and warrants further investigation. Amongst the DEG we characterized, Kiss1 and Npvf are the only ones expressed by hypothalamic neurons. Kiss1 and Npvf are also the only DEG for which seasonal expression is abolished in both intact and OVX ewes, compared to OVX+E2 ewes (Fig 1, Fig 2 & Fig 3). Therefore, seasonal

Figure legends Figure 1 :

 legends1 Figure legends

Figure 2 :

 2 Figure 2 : qRT-PCR for a set of seasonal genes expressed in the MBH of intact ewes in November,

Figure 3 :

 3 Figure 3 : Direct comparison of expression levels (as assessed by qRT-PCR) for a set of seasonal

Figure 4 :

 4 Figure 4 : qRT-PCR for a set of seasonal genes expressed in the MBH of intact ewes of the Ile-

Figure 5 :

 5 Figure 5 : A revised circuit-level model for the photoperiodic control of seasonal breeding (see

  Fig.1Intact ewes -ISH

RNA extraction and qRT-PCR analysis

  80 kg) maintained under normal husbandry at the research station of the Institut National de la Recherche Agronomique in Nouzilly

	killed (n=10 per breed and time point) on February 13-16 th , March 5-6 th and March 15-16 th , a
	narrow seasonal time window during which anestrus develops in both breeds. The 1 st two groups
	(i.e. IdF and Rom breeds) were killed during a "natural" follicular phase (display of estrus behavior
	when presented with a teaser ram). The remaining 40 ewes received a progestogen treatment
	(PAO, (vaginal sponge) on February 16-17 th , which was removed on March 3-5 th and followed by INRA, 2018. Animal Physiology Experimental Facility, DOI:
	10.15454/1.5573896321728955E12). We used either intact, ovariectomized (OVX) or pregnant mare serum gonadotropin (PMSG) injection. Ewes were killed 2 or 12 days later, when
	ovariectomized, estradiol-implanted (OVX+E2; 1cm silastic implant) ewes. This OVX+E2 model displaying an "induced" follicular or luteal phase, respectively.
	normalizes the level of circulating E2, which uncovers the well-documented central seasonal shift
	RNA extraction was performed on each MBH block using TriReagent (Sigma). Concentration and
	purity of individual samples were determined with Nanodrop 2000 (ThermoScientific) and
	integrity was checked by standard agarose gel electrophoresis. All procedures, including qPCR
	all other experiments RNA was extracted from a MBH block, which assays for all target genes, have been thoroughly validated (Lomet et al 2018; Dardente and Lomet
	comprises the PT, the median eminence, the arcuate nucleus and the dorsomedial and ventromedial 2018). cDNA was synthesised from each individual RNA sample using Omniscript RT kit
	hypothalamic nuclei. The RNA samples were used to perform qRT-PCR; some of these samples (Qiagen) and Oligo-dT primers (synthesized by Eurofins, Germany). The optimal cDNA dilution
	have been used in prior studies (see below and Lomet al 2018, Dardente & Lomet 2018). and calibration curves for each PCR primer pair (Eurofins) were established using cDNA
	synthesized from an equimolar mix of all individual RNA from each experiment. As a negative
	In experiment 1, intact IdF ewes were kept in an open barn and killed in May (non-breeding season, control, the same mix with water instead of RT was prepared. Quantitative PCR (qPCR) was
	n=6) and November (breeding season, n=6). In experiment 2, intact IdF ewes were kept in an open performed using CFX-96 Real-Time PCR Detection System (Bio-Rad) and Sso Adv Universal
	barn and killed in November (n=6), February (when breeding stops, n=5), May (n=6) and August SYBR Green Supermix (Bio-Rad). All primer pairs had efficiencies ranging from 85% to 110%.
	(when breeding resumes, n=6). For experiment 3 (IdF breed), a direct comparison of expression All samples (unknown, standard curves) were assessed in triplicate and Rplp0 (ribosomal subunit
	levels between May and November was performed across three distinct conditions with respect to P0, a.k.a 36B4) was used as a housekeeping gene. The quantification of mRNA level was obtained
	sex steroid status: intact ewes (groups from Expt 2), OVX+E2 ewes (samples from Lomet et al by the 2 -∆CT method. Data are normalized and presented as fold-increase compared to the condition
	2018) and OVX ewes (n=5 for May and n=4 for November; see Fig S3A for details). In with the lowest expression level. For Expt 4, all cDNA were re-synthesized in a single run to allow
	experiment 4, we used intact ewes of IdF (n=30) and Romanov (Rom; n=30) breeds. Ewes were direct comparison across conditions. The list of qPCR primers is provided in Table

in the negative feedback action of E2 on gonadotropin secretion

[START_REF] Legan | The endocrin control of seasonal reproductive function in the ewe: a marked change in response to the negative feedback action of estradiol on luteinizing hormone secretion[END_REF][START_REF] Karsch | Neuroendocrine basis of seasonal reproduction[END_REF] Goodman & Inskeep, 2015)

.

For experiment 1, in situ hybridization (ISH) was performed on brain sections (see next section

and Lomet et al 2018)

. For

  and OVX+E2 ewes[START_REF] Lomet | The impact of thyroid hormone in seasonal breeding has a restricted transcriptional signature[END_REF][START_REF] Dardente | Photoperiod and thyroid hormone regulate expression of l-dopachrome tautomerase (Dct), a melanocyte stem-cell marker, in tanycytes of the ovine hypothalamus[END_REF], notwithstanding possible differences in amplitude, which are difficult to assess from comparing independent experiments. This situation in PT and tanycytes differed markedly for the hypothalamic markers Kiss1 and Npvf, which did not display any difference in expression levels in intact ewes killed in May and November (Fig1C). This result is in stark contrast with what we found in OVX+E2 ewes, which display Kiss1 mRNA levels that are ~10-fold higher in November than in May, while Npvf displays Following removal of the E2 implant in March, LH levels increased rapidly, reaching a plateau after 3-4 weeks. For the November group, LH levels had already increased to their near-maximal level when the E2 implant was removed in August. Kiss1 and Npvf levels were clamped to low and intermediate-high levels, respectively. Conversely, in OVX ewes Kiss1 and Npvf levels were clamped to intermediate-high and low levels, respectively. Therefore, across the distinctive experimental conditions, Kiss1 and Npvf display inverse patterns of expression.

	E2 was found to have little, if any, impact on the expression of the 9 PT markers investigated (Fig
	3A and Fig S3B): all DEG showed high-amplitude rhythms in intact, OVX+E2 and OVX ewes
	(2-way ANOVA, Time effect p<0.0001 for all DEG; see Table S2).
	To clarify the impact of E2 on molecular seasonal markers expressed in the three compartments
	of the photoperiod transduction cascade -PT, tanycytes and hypothalamus -we employed OVX

al 2018; Dardente and Lomet 2018). The PT markers Tshb, Fam150b, Vmo1, Ezh2, Klkb1, Suv39H2, Vcan and Elovl6 were expressed at higher levels in May than in November, while Chga and Aa-Nat showed the opposite seasonal pattern, with higher expression in November than in May (Fig 1A & Fig S1, see Table

S2

for stats). The markers of tanycytes Dio2, Tmem252, Shh and SlcO1c1 showed increased seasonal expression in May, while NpSR1 was more expressed in November (Fig

1B

). Therefore we conclude that seasonal molecular markers of PT and tanycytes showed similar expression patterns in intact the opposite pattern with mRNA levels ~10-fold higher in May than in November

[START_REF] Lomet | The impact of thyroid hormone in seasonal breeding has a restricted transcriptional signature[END_REF]

.

We then used qRT-PCR to track seasonal changes of DEG in intact ewes with a more quantitative approach and a finer temporal resolution (every 3 months). Seasonal high-amplitude rhythms of expression were revealed for all PT markers (>50-fold for Tshb; Fig 2A & Fig S2, see Table S2

for stats) and for all tanycytic markers (>45-fold for Tmem252; Fig 2A, see Table S2 for stats), whether peak expression level occurs in May, or in November as in the case of Chga and NpSR1. Dio2 and Dio3 displayed characteristic opposite patterns of seasonal expression, which predict increasing and decreasing levels of T3 levels at the transition to anestrus and breeding, respectively (deduced from the ratio Dio2/Dio3). Here again, expression levels for Kiss1 and Npvf did not display statistically significant seasonal rhythmicity (Fig 2C and Table S2). ewes, and used qRT-PCR to allow a direct side-by-side comparison of relative expression levels in May and November between intact, OVX+E2 and OVX ewes. Efficacy of OVX was validated by RIA for LH (Fig S3A). A slightly different conclusion emerged for the 7 tanycytic markers (Fig 3B and Fig S3C). All

genes displayed statistically significant differences with time, even though the amplitude of the rhythms appeared differentially affected across genes (2-way ANOVA, Time effect p<0.005 for all DEG; see Table

S2

).. All May/Nov pairwise comparisons were statistically significant with two exceptions: Dio3and Shh both displayed an increase in expression in OVX+E2 ewes, but neither in intact nor in OVX ewes. Nevertheless, the Dio2/Dio3 ratio remained higher in May for all 3 groups of animals: 4.6-fold for OVX+E2, 2.8-and 3.2-fold for intact and OVX, respectively. Consistent with our earlier findings using these samples from OVX+E2 ewes

[START_REF] Lomet | The impact of thyroid hormone in seasonal breeding has a restricted transcriptional signature[END_REF]

, we found a ~10-fold increase in expression for Kiss1 between May and November. Conversely, Npvf expression decreased by ~10-fold between May and November. Most importantly, these large opposite changes in Kiss1 and Npvf expression were lost not only in intact ewes, but also in OVX ewes see Table

S2

for stats). In intact ewes,

  Yamzaki et al 2000 ;[START_REF] Damiola | Restricted feeding uncouples circadian oscillators in peripheral tissues from the central pacemaker in the suprachiasmatic nucleus[END_REF], which are all able to phase-shift peripheral clocks[START_REF] Dibner | The mammalian circadian timing system: organization and coordination of central and peripheral clocks[END_REF][START_REF] Mohawk | Central and peripheral circadian clocks in mammals[END_REF].

	We show that P4 has a modest and gene-specific impact in tanycytes. Considering that tanycytes
	receive photoperiodic information through PT-derived products (TSH and others ; Lomet et al
	2018 ; Dardente et al 2019a), it leaves open the possibility that P4 somehow modulates the impact
	of these entraining signal(s). Similar trends in gene expression patterns were seen in two different
	ovine breeds (IdF and Rom) and for multiple molecular markers, which are either induced (Dio2,
	al 2002 ; Dardente et al 2003) are indeed
	circannual clocks. In this regard, a parallel between the PT and the SCN (master circadian clock)

may be drawn. The SCN is entrained by the light-dark cycle, while being insensitive to multiple resetting cues including glucocorticoids

[START_REF] Balsalobre | Resetting of circadian time in peripheral tissues by glucocorticoid signaling[END_REF]

, temperature

[START_REF] Buhr | Temperature as a universal resetting cue for mammalian circadian oscillators[END_REF]

or feeding schedule (

Table S1 :

 S1 Word document

Table S2 :

 S2 Excel file