Q Poterek

P.-A Herrault

G Forestier
email: germain.forestier@uha.fr

D Schwartz

REVEALING LONG-TERM PHYSIOLOGICAL TRAJECTORIES OF GRASSLANDS FROM LEGACY B&W AERIAL PHOTOGRAPHS

Keywords: Archive Aerial Photographs, Radiometric Correction, Colorization, Digital Greenness Metrics, Phenology, Time Series

Landscape reconstruction is crucial to measure the effects of climate change or past land use on current biodiversity. In particular, retracing past phenological changes can serve as a basis for explaining current patterns of plant communities and predict the future extinction of species. Old spatial data are currently used to reconstruct vegetation changes, both morphologically (with landscape metrics) and semantically (grasslands to crops for instance). However, poor radiometric properties (single panchromatic channel, illumination variation, etc.) do not offer the possibility to compute environmental variables (e.g. NDVI and color indices), which strongly limits long-term phenological reconstruction. In this study, we propose a workflow for reconstructing phenological trajectories of grasslands from 1958 to 2011, in the French central Vosges, from old aerial black and white (B&W) photographs. Noise and vignetting corruptions were first corrected in B&W photographs with non-local filtering algorithms. Panchromatic scans were then colorized with a Generative Adversarial Network (GAN). Based on the predicted channels, we finally computed digital greenness metrics (Green Chromatic Coordinate, Excess Greenness) to measure vegetation activity in grasslands. Our results demonstrated the feasibility of reconstructing long-term phenological trajectories from legacy photographs with insights at different levels : (1) the proposed correction methods provided radiometric improvements in old aerial missions; (2) the colorization process led to promising and plausible colorized historical products; (3) digital greenness metrics were useful for describing past vegetation activity.

INTRODUCTION

It is widely recognized that present day biodiversity may reflect past land use or past climate [START_REF] Jansson | Global variation in diversification rates of flowering plants: energy vs. climate change[END_REF] because of a possible delay in the response of certain species to habitat perturbations [START_REF] Kuussaari | Extinction debt: a challenge for biodiversity conservation[END_REF]. This interval depends on changes themselves (e.g. intensity and cyclicity) as well as species traits, namely mobility dispersion abilities. In order to assess time-lag between changes in a landscape, and changes in the associated populations, studies rely most of the time on old spatial data that help highlighting spatio-temporal trajectories over large extents and long time periods [START_REF] Proenc ¸a | Global biodiversity monitoring: from data sources to essential biodiversity variables[END_REF]. Among these spatial data, legacy aerial photographs are largely under-exploited while they offer unique opportunities to monitor landscapes at a very high spatial resolution (≤ 1m) from up to the early 1930s [START_REF] Morgan | Aerial Photography: A Rapidly Evolving Tool for Ecological Management[END_REF][START_REF] Morgan | Historical Aerial Photography for Landscape Analysis[END_REF]. Despite these advantages, their use remains problematic for many reasons. First, they feature heterogeneous specifications and quality (e.g. noise and vignetting), mainly due to the properties of the acquisition system and the scanning procedure. Second, they lack inherent exploitable attributes, especially in the case of panchromatic pictures, making the development of automatic processing chains a complex task [START_REF] Paine | Aerial Photography and Image Interpretation[END_REF][START_REF] Aber | Small-Format Aerial Photography: Principles, Techniques and Geoscience Applications[END_REF].

Another major problem refers to the type of information that can be retrieved from these data. Landscape reconstruction * Corresponding author is generally carried out one of two different ways : (1) morphologically by monitoring landscape metrics (e.g. area and connectivity) over the considered time period [START_REF] Franco | Using historical records, aerial photography and dendroecological methods to determine vegetation changes in a grassy woodland since European settlement[END_REF][START_REF] Herrault | Combined effects of area, connectivity, history and structural heterogeneity of woodlands on the species richness of hoverflies (Diptera: Syrphidae)[END_REF]; (2) semantically by describing changes in the nature of spatial objects (grasslands to crops for instance) [START_REF] Treitz | Remote sensing for mapping and monitoring land-cover and land-use change-an introduction[END_REF]. Recently, analysis of radiometric properties in aerial photographs provided useful information to monitor vegetation health or forest phenology [START_REF] Franco | Using historical records, aerial photography and dendroecological methods to determine vegetation changes in a grassy woodland since European settlement[END_REF][START_REF] Reid | Using excess greenness and green chromatic coordinate colour indices from aerial images to assess lodgepole pine vigour, mortality and disease occurrence[END_REF]. Indeed, metrics such as the Green Chromatic Coordinate (GCC) or Excess Greeness (EG) are invariant to illumination conditions and outperform conventional indices such as the NDVI when processing near-sensing images [START_REF] Nijland | Monitoring plant condition and phenology using infrared sensitive consumer grade digital cameras[END_REF]. Unfortunately, inherent characteristics of the oldest aerial missions (i.e. B&W photographs) do not offer the possibility to retrieve these indices [START_REF] Morgan | Aerial Photography: A Rapidly Evolving Tool for Ecological Management[END_REF]. Thus, long term phenological reconstructions are strongly restricted temporally.

In this study, we develop a generic framework to reconstruct phenological trajectories of grasslands from 1958 to 2011, in the central Vosges from old aerial B&W photographs. We propose a workflow divided in three major steps :

1. Correction of noise and vignetting effects in old aerial B&W photographs;

2. Spectralization of old aerial B&W photographs;

3. Reconstruction of changes in grasslands with greenness indices.

STUDY SITE

The study site is located in the central Vosges, 30 km from Strasbourg in the East of France (Figure 1. It covers approximately 15,000 ha and includes 10 municipalities. This is a hilly region (100 to 400m) dissected by north-south valleys. The climate is semi-continental with oceanic influences. It comprises a majority of grasslands and crops (wheat, corn) with remnant traditional orchards delineated by hedges. A positive gradual landscape openness is directly observable from North to South, indicating changes in management practices, from traditional to intensive agriculture. Thus, heterogeneous landscape trajectories were provoked by this management diversity and might explain part of the spatial distribution of grasslands plant species over the study site. In the scope of this paper, 34 grasslands were selected to be monitored according to two criteria: (1) grasslands must persist without discontinuity over time, and (2) they vary in terms of openness throughout the years.

METHODOLOGY

Preparation of spatial data

Six air missions acquired between 1958 and 2011 were used to reconstruct long term phenological trajectories of grasslands (Table 1). For testing out the proposed methodology, only 1 single photograph was included in our database for each mission. All photographs that were selected overlapped. Each picture was acquired in spring (from May to June). This means that the proposed photographs captured the signal of grasslands just before the first peak of annual greenness or after the first mowing event of the year. Therefore, for specific years, some grasslands exhibited an almost bare soil.

Panchromatic scans were manually georeferenced and rectified using around 10 GCPs and by applying a second-degree polynomial transformation between image and ground coordinates (RM SE = 1.21 ± 0.33 m). The most recent RGB orthophotograph acquired in 2011 systematically served as a reference for the registration of B&W aerial photographs. 1. Specifications for the selected aerial photographs. * For scanned films, scale is followed by a pixel size estimate that was computed using calibration certificates.

In order to demonstrate the limitations imposed by historical photos, as well as the techniques available to address them, a processing chain has been developed to reconstruct long-term phenological trajectories (Figure 2). The workflow can be divided in three steps : (1) correction and harmonization of old aerial photographs, (2) spectralization of panchromatic photographs and (3) phenological time series analysis. [START_REF] Jalobeanu | Estimation of blur and noise parameters in remote sensing[END_REF]. The consequence is a random variation in pixels values, independent of the original data, making interpretation and processing complicated tasks [START_REF] Jalobeanu | Estimation of blur and noise parameters in remote sensing[END_REF][START_REF] Corner | Noise estimation in remote sensing imagery using data masking[END_REF]. In the context of old aerial photographs, Gaussian corruption is the most prominent, mainly due to illumination, sensor temperature and film scanning. Since aerial stills were not acquired following a standard and non-evolving scheme, the degree of corruption may vary across the entire time series [START_REF] Morgan | Aerial Photography: A Rapidly Evolving Tool for Ecological Management[END_REF]. Thus, noise analysis and correction are mandatory for harmonizing the available photographs.

In this paper, we assess noise with a pyramid-based method.

The standard deviation of the Gaussian noise distribution was first obtained by computing the median absolute deviation of the wavelet detail coefficients [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] at 6 different scales. Each photograph in the series was resampled into lower resolution overviews, by factors of 1, 2, 3, 4, 6 and 12. The outputs were then averaged for each scale and year so as to retrieve an estimate of the standard deviation for each photograph.

After estimation, noise in photographs was corrected with the non-local means (NLM) algorithm [START_REF] Buades | Non-local means denoising[END_REF]. This technique has been used in a significant number of remote sensing applications, where it was able to smooth surfaces while retaining fine details, such as roads and edges, as in [START_REF] Huang | A parallel nonlocal means algorithm for remote sensing image denoising on an Intel Xeon Phi platform[END_REF] for example. Unlike other filtering techniques, NLM performs denoising at patch-level. The algorithm recursively takes one pixel to process, along with its neighborhood. It then retrieves all the patches available in the image and measures their similarity to the processed area. An average of the patch values is then calculated, each contribution being weighted by the degree of similarity to the processed pixel. The standard deviation estimates that we previously computed were passed to the algorithm so as to guide the weighting process. Finally, the pixel value is updated with the weighted and averaged contributions.

Vignetting correction

In aerial photographs, spatially organized variations in brightness can be observed, most of the time following a radial gradient [START_REF] Yu | Vignetting distortion correction method for high quality digital imaging[END_REF]. This results in a bright principal point and dark borders around the fiducial marks. This phenomenon is called vignetting and can be troublesome when the photographs in question are used in photogrammetry or remote sensing applications. Indeed, the gradual darkening of the image results in a non-homogeneous signal for similar surfaces. This can notably disturb the matching of homologous points for the relative orientation of several photographs [START_REF] Kim | Robust Radiometric Calibration and Vignetting Correction[END_REF], reduce image classification performance, or compromise the analysis of time series [START_REF] Kelcey | Sensor correction of a 6-band multispectral imaging sensor for UAV remote sensing[END_REF].

We propose a simple but efficient technique for removing not only vignetting, but also other local variations in brightness.

After masking fiducial marks, the estimation of a correction map was carried out following a pyramid scheme similar to the one previously described for measuring the standard deviation of Gaussian noise. Each photograph was resampled into lower resolutions. All scales were sized back to the original photograph dimensions using bicubic interpolation, and later averaged, in order to compute a multi-scale mean illumination map.

It was then subtracted from the initial still, thus eliminating any illumination effect, whether due to exposure or vignetting for example. Finally, pixel values from the output were rescaled so as to match the distribution of the original.

Spectralization of B&W photographs

After correction and in order to ensure spatial coherence, all of the available photos were resampled to a 1m resolution. It was then necessary to enrich the spectral content of the available series before calculating phenological metrics and characterizing vegetation activity at each date. Indeed, panchromatic and single-infrared channels provided by old aerial missions are insufficient for the calculation of environmental variables (NDVI or color indices). Thus, we recommend the development of a spectralization technique that would help predict new channels in the visible domain, based on initial panchromatic information.

In essence, spectralization corresponds to colorization. Several techniques are available for colorizing black and white photographs. Recently, various techniques based on deep learning have been proposed, mostly non-convolutional networks [START_REF] Cheng | Deep colorization[END_REF], CNNs [START_REF] Zhang | Colorful image colorization[END_REF] and GANs [START_REF] Isola | Image-toimage translation with conditional adversarial networks[END_REF]. Unlike other methods that rely on human intervention for placing color scribbles [START_REF] Yatziv | Fast image and video colorization using chrominance blending[END_REF], guiding transfer algorithms [START_REF] Welsh | Transferring Color to Greyscale Images[END_REF] or computing appropriate attributes for standard machine learning (Desh-pande et al., 2015), the advent of deep learning has helped alleviate the colorization process, as models automatically learn feature maps to solve the grayscale-to-color mapping problem.

In the scope of this work, we propose a deep colorization model based on a generative adversarial network (GAN). It is built around two neural networks pitted against one another. The first corresponds to a generator G that learns how to produce samples that could belong to a reference distribution. The second corresponds to a discriminator D, whose role is to distinguish true samples from generated ones [START_REF] Goodfellow | Generative adversarial nets[END_REF]. The end objective is to train a generator capable of deceiving the discriminator. In this paper, we conditioned G by giving it grayscale samples as inputs, so as to learn color channels [START_REF] Mirza | Conditional Generative Adversarial Nets[END_REF]. It should be pointed out that, as of the submission of this paper, the colorization of aerial photographs time series has not been tackled in any other work. The only other papers that come close were either focused on processing singe-date products, with (1) the colorization of singlepolarization SAR images to obtain full-polarization samples [START_REF] Song | Radar Image Colorization: Converting Single-Polarization to Fully Polarimetric Using Deep Neural Networks[END_REF], or (2) the colorization of modern high resolution multispectral satellite imagery [START_REF] Liu | Single satellite imagery simultaneous super-resolution and colorization using multi-task deep neural networks[END_REF]. Thus, the proposed methodology is meant to address different shortcomings, with the colorization of archive spatial data that were compiled into long-term time series.

Since there is no proper aerial photographs database available at the moment for deep learning development, we built our own training and validation data sets. To minimize training time, and since the aim was only to colorize old photographs for a small extent, we randomly extracted a total of 3, 000 128 × 128 samples from the 05/05/2011 color reference, spatially independent from one another. Samples were split into training and validation sets given a 10:1 ratio.

In order to learn a panchromatic-to-color mapping, all samples were converted from the RGB color space to CIE Lab. This resulted in images with 3 channels : L, a and b. L contains luminance information, and approximately corresponds to a black and white photograph. Its values range from 0 to 100. The last two channels, a and b, are uncorrelated and contain color information. Their values range from -128 to +127. This step allowed separating grayscale and chromatic information. Thus, the proposed model only had to learn how to predict a and b, while preserving spatial information provided by L after concatenation of the input and the two predicted channels. This step also made it possible to learn only two channels (a and b) instead of three (R, G and B), making the training step less costly in terms of time and resources. At last, pixel values from the Lab samples were scaled to a [-1; +1] range as a prerequisite for traditional GAN training [START_REF] Goodfellow | Generative adversarial nets[END_REF].

The generator was based on a fully-convolutional UNet architecture (Ronneberger et al., 2015). It was fed with a batch of black and white (L) samples that were first downsampled, and later upsampled after passing through a bottleneck. The output corresponded to a and b channels for the input samples. Features learned during downsampling were concatenated in a symmetric fashion to their corresponding upsampled counterparts. This technique helped recognize spatial semantics at various scales, while maintaining spatial information. During training, samples passed to the generator were augmented by applying random horizontal and vertical flips, rotation, blur and variation in lightness. Data augmentation was necessary for simulating the look and feel of old aerial photographs taken under different conditions with specific sensors. In the end, the generator was trained to minimize the probability of the discriminator labeling generated data as not being part of the reference color distribution [START_REF] Goodfellow | Generative adversarial nets[END_REF]. The cost function for the generator was defined as:

LG

= E [log(D(G(p), p))] (1)
where

G = generator D = discriminator p = panchromatic sample G(p) = colorized sample
The discriminator was based on the PatchGAN architecture proposed by [START_REF] Isola | Image-toimage translation with conditional adversarial networks[END_REF]. It was fed alternatively with real color and colorized samples and was tasked with the evaluation of generated chroma, in regards to the reference distribution.

The discriminator was trained to maximize its probability of distinguishing between both real and generated samples [START_REF] Goodfellow | Generative adversarial nets[END_REF]. The cost function for the discriminator was defined as:

LD = E [log(D(c, p))] + E [log(1 -D(G(p), p))] (2)
where

D = discriminator G = generator p = panchromatic sample c = color sample G(p) = colorized sample
The model was trained for 1, 000 epochs, using a free Tesla K80 GPU provided by the Google Colab platform.

Colorized photographs from the validation set were evaluated with image quality metrics such as peak signal to-noise ratio (PSNR), and structural similarity index (SSIM). PSNR provides information regarding the quality of reconstruction after colorization. Its values range from 0 db (different) to +∞ (identical).

The SSIM was computed in order to provide information regarding the similarity in lightness, contrast and color fidelity. Unlike PSNR, it is said to correlate with human vision assessment [START_REF] Wang | Image Quality Assessment: From Error Visibility to Structural Similarity[END_REF]. Its values range from 0 (different) to 1 (identical).

Finally, after training and evaluation, the old aerial photographs (L) were colorized using the generator. The input L channel, along with the predicted a and b channels, were concatenated and converted to RGB, in order to compute phenological indices.

Time series analysis

Based on the corrected and harmonized photographs, we were then able to analyze changes for grasslands since 1958. It is important to note that the reconstructed trajectories cannot be used to evaluate long-term productivity trends because of irregular acquisition dates between air missions, and the use of one single date per year. Our objective was to assess the possibility of obtaining plausible trajectories after correction and spectralization of old aerial photographs.

First, one 10 × 10 m quadrat was randomly placed in each of the 34 selected grasslands. All pixels entirely contained in each quadrat were extracted in order to compute phenological indices.

Then, two digital phenological metrics were calculated. Recent works showed the usefulness of computing color indices for measuring vegetation activity in near-sensing images [START_REF] Reid | Using excess greenness and green chromatic coordinate colour indices from aerial images to assess lodgepole pine vigour, mortality and disease occurrence[END_REF]. Differences in scene illumination is a major issue when processing old aerial photographs, mainly due to vignetting and variations in exposition [START_REF] Aber | Small-Format Aerial Photography: Principles, Techniques and Geoscience Applications[END_REF]. Taking this issue into account, two metrics were computed in the scope of this study : (1) Green Chromatic Coordinate (GCC, see Equation 3), which is well suited for suppressing variability in scene illumination, and (2) Excess Greenness (EG, see Equation 4), which provides better differentiation between plant material and soil background compared to other color indices [START_REF] Nijland | Monitoring plant condition and phenology using infrared sensitive consumer grade digital cameras[END_REF][START_REF] Reid | Using excess greenness and green chromatic coordinate colour indices from aerial images to assess lodgepole pine vigour, mortality and disease occurrence[END_REF].

GCCi = Gi (Ri + Gi + Bi) (3) EGi = 2 × Gi -(Ri + Bi) (4)
where Gi, Ri and Bi are the predicted green, red and blue pixel values respectively.

In order to reconstruct greenness trajectories from 1958 to 2011, the mean and standard deviation of GCC and EG were computed for each quadrat and year. We also reported these metrics both for the original and the corrected photograph at each date to assess the effects of the proposed correction workflow on the greenness metrics.

Temporal profiles for greenness were finally clustered so as to highlight different groups of grasslands trajectories over time.

An agglomerative clustering technique was used. The number of clusters n was set empirically after trying out multiple values. Samples were then merged recursively using cosine distance affinity. More robust time series clustering techniques, such as the ones based on dynamic time warping [START_REF] Berndt | Using dynamic time warping to find patterns in time series[END_REF], were not explored in this study due to non-monotonic intervals between air missions [START_REF] Petitjean | Dynamic time warping averaging of time series allows faster and more accurate classification[END_REF].

RESULTS AND DISCUSSION

In this section, we first demonstrate the effectiveness of the proposed methodology for correcting and harmonizing old aerial photographs. Results from the pipeline are presented for 1958 in Figure 3.

For all available samples in the time series, we were able to estimate standard deviation of the Gaussian noise corruption. The estimated values based on a 1 byte pixel depth were the following for each year: 1.18 (1958), 3.01 (1976), 1.79 (1986), 2.06 (1992) and 2.44 (2000). Standard deviations were then used to denoise all of the available photographs, based on the NLM algorithm. A visual inspection of the corrected sample shows this technique has managed to retain details, while smoothing the provided photograph. This is further highlighted by the corresponding error map, which is highly textured and contains the removed noisy component (Figure 3). Even though there is no literature regarding the correction of additive Gaussian noise in archive photographs, other papers have come to the same findings. Indeed, the NLM algorithm was successfully applied in Error map are also provided, and correspond to the difference between the original and the corrected picture. Finally the colorization result is presented for the studied area.

both spatial [START_REF] Zhao | Satellite image deconvolution based on nonlocal means[END_REF] and wavelet [START_REF] Iqbal | Satellite image resolution enhancement using dual-tree complex wavelet transform and nonlocal means[END_REF][START_REF] Kang | Push-broomtype very high-resolution satellite sensor data correction using combined wavelet-fourier and multiscale non-local means filtering[END_REF] domains, and was able to remove noise while preserving edges and texture, at the cost of some blurring.

The proposed methodology for correcting local and non-local variations in brightness also proved to be efficient. In Figure 3, vignetting is indicated in the original photograph by a strong intensity fall-off around the edges. Meanwhile, the brightest area does not correspond to the center of the picture, nor to its principal point, as would normally be the case. The estimation of a multi-scale mean illumination map made it possible to account for these specific features, as can be seen in the corresponding error map. In an ideal case, the proposed method should only provide an illumination field. However, the error map shows that it also takes into account large spatial structures, such as forests and hedges. One technique to consider for estimating only vignetting, for example, would be to estimate it directly using a time series. The sample to be corrected would be compared to a reference photograph, taken at an earlier or later date. However, the difference between the two pictures should be small enough to minimize possible differences, such as changes in land use, which could affect the estimation of vignetting. Even though the proposed technique was not the same as ours, (da Silva, Candeias, 2012) have obtained similar results, with an increase in mean brightness value and a decrease in standard deviation after vignetting correction.

The distribution of pixel values for each year is shown in Figure 4. Mean pixel values along with standard deviations are given in Table 2 for both the original and corrected samples. Despite visually noticeable changes before and after correction on the photographs (Figure 3), the distribution of pixel values varied only slightly (Figure 4 and Table 2). This is encouraging, as the objective is to preserve the initial signal, while cleaning it of any imperfection that might disturb subsequent analyses. From 1976 to 2000, the mean of pixel values increased slightly, by an average of 0.0525 points. The standard deviation decreased by 0.4425 points on average for the same dates. Similar trends are observed for the 1958 photograph, but with a greater amplitude, as the mean of pixel values increased by 7.67 points, and the standard deviation decreased by 0.9 points. This picture was the most deteriorated, especially due to a strong vignetting, which explains the significance of the corrections that were made.

After correction, photographs were colorized using a generative adversarial network that was conditioned and pre-trained with a 2011 color reference. After 1000 epochs, results from the validation data set indicate robust performances for the colorization task, with P SN R = 39.17 ± 4.79 and SSIM = 0.93 ± 0.05.

Our model performs as well as or better than approaches proposed in other deep learning articles. In comparison, [START_REF] Varga | Twin deep convolutional neural network for example-based image colorization[END_REF] and [START_REF] Deshpande | Learning diverse image colorization[END_REF] obtained mean SSIM values of 0.89 and 0.93 respectively. Regarding PSNR, [START_REF] Larsson | Learning representations for automatic colorization[END_REF] and [START_REF] Liu | Single satellite imagery simultaneous super-resolution and colorization using multi-task deep neural networks[END_REF] obtained scores of 24.45 db and 25.05 db respectively. After training, chroma was predicted for all of the available black and white photographs. A colorization example is shown in Figure 3 for the 1958 sample. No color reference is available for old photographs. However, visual analysis of the results testify to the model's ability to process old photographs, although it was trained with recent pictures. Despite local errors, colors were properly predicted for all types of surfaces, including urban areas, crops, grasslands and bare soils (Figure 5). Most errors are located around image corners, or come in the form of large brown, blue or cyan patches. This behavior could imply that the model is sensitive to edge effect, and does not perform well over homogeneous surfaces larger than training images. Increasing the height and width of training samples could thus help overcome this issue. Given the decent performances of the model and the adequacy of the colorizations, the latter served as a medium for the calculation of phenological metrics.

We first observe similar relative greenness trends calculated before and after correction (Figure 6). It was an expected result since brightness correction was not meant to drastically change the radiometric distribution of aerial photographs, but rather correct local and global anomalies. Then, mean and standard deviations at each date showed comparable absolute results before and after correction (Table 3). Differences between GCC and EG means are lower than 0.05 between the two steps, indicating minor effects from brightness correction on the calculation of phenological metrics. These results might be explained by several reasons. First, considering observations made from Figure 4, noise correction had no major impact on the gray levels distribution except decreasing the amplitude of the distribution values. This effect is particularly noticeable in 1986, 1992 and 2000. Nonetheless, patterns in the texture of grasslands remained heterogeneous even if they were attenuated.

The NLM algorithm successfully preserved the radiometric heterogeneity in grasslands since it was not considered Gaussian noise. Second, the proposed metrics for measuring phenology are supposed to be insensitive to variations in illumination. It is particularly true for GCC for which the results after correction were similar to that of the original photographs. Outcomes are more nuanced for EG since a significant decrease in the standard deviation was observed after correction on the 1958 mission (from 0.072 to 0.066) for which photographs were the most corrupted by vignetting and dark outliers (Figure 4). Corrections allowed enhancing low brightness values, which reduced the EG variability among grasslands. Our results confirm findings from [START_REF] Nijland | Monitoring plant condition and phenology using infrared sensitive consumer grade digital cameras[END_REF], who recommend the use of GCC instead of EG in images where strong brightness variations occur.

Secondly, GCC and EG temporal profiles showed similar trends over the studied time period (Figure 6). Average profiles only indicated a reverse trend between 1986 and 1992, where GCC 1 9 6 0 1 9 7 0 1 9 8 0 1 9 9 0 2 0 0 0 2 0 1 was increasing (+0.3%) while EG was decreasing (-17%). One explanation might be that EG is more sensitive to extreme values than GCC. Therefore, an excessive number of bare soils in 1992 (19 out of 34 grasslands) contributed to the decrease of mean EG, leading to a trend reversal between those two dates. These results reinforce the necessity of using EG to maximize the differentiation between plant materials and background soils [START_REF] Nijland | Monitoring plant condition and phenology using infrared sensitive consumer grade digital cameras[END_REF][START_REF] Reid | Using excess greenness and green chromatic coordinate colour indices from aerial images to assess lodgepole pine vigour, mortality and disease occurrence[END_REF]. However, problems might occur in case of non-corrected illuminations variations, typical of the oldest B&W air missions.

High standard deviations were observed for each year, especially for 1958, 1992 and 2011, demonstrating a high greenness variability amongst the grasslands monitored in this study (see Figures 6 and7). This variability might be explained by several factors. First, species assemblages in grasslands vary strongly, leading to important phenological time-lags between grasslands. Furthermore, recent works showed that the intraannual phenological cycle is a powerful proxy to predict species diversity in grasslands from satellite image time series [START_REF] Rapinel | Mapping grassland plant communities using a fuzzy approach to address floristic and spectral uncertainty[END_REF][START_REF] Fauvel | Prediction of plant diversity in grasslands using Sentinel-1 and -2 satellite image time series[END_REF]. Secondly, the diversity of management practices provoked variations in the timing of mowing events. It results in a large contrast of greenness between grasslands at each date because of the presence of bare soils and grasslands advanced in their phenological cycle.

Last, we tested the possibility of clustering phenological grasslands trajectories (Figure 7). We remind that the retrieved trajectories cannot be used to analyze long-term grasslands productivity because of the non-monotonic acquisition dates, and the use of one single air mission per year. Nevertheless, clustering results showed the feasibility of identifying distinct groups of grasslands (6 clusters in total) and formulating assumptions.

In particular, clusters 1 (n = 7) and 3 (n = 6) exhibit similar trends over time, but grasslands from the third cluster show higher differences between dates, namely between 1958 and 1976. These grasslands might correspond to ancient hay meadows, with a function that has evolved over time. Individuals from the first cluster showed a higher relative stability, without the presence of bare soils. They might have not been used as a forage resource for a long time. Similarly, clusters 2 (n = 12) and 5 (n = 6) displayed similar trends, but grasslands from the latter group showed higher greenness values in the most recent photograph. Unlike cluster 2, the species composition of these grasslands might need more time to reach its peak of greenness, leading farmers to mow later in the season.

Figure 7. Clusters of trajectories obtained from the GCC on corrected products. In green, the average profile for the cluster.

In gray, the individual trajectories for the grasslands. Each square corresponds to one of the 10 × 10 m quadrats that were used for sampling. Background images correspond to the colorized aerial photographs.

These preliminary results are promising but should be pursued in order to retrace long-term productivity of grasslands. Several options must be considered. Fist, long-term grasslands trajectories must be retrieved from entire annual cycles to be reconstructed. However, only a handful of air missions are available each year in the best case scenario, with a maximum of 2 to 3 acquisitions for the studied area. One solution would be to densify these series by learning new acquisition dates based on present day satellite time series. Second, these reconstructed phenological series should be analyzed with caution and require further validation due to high uncertainty. Ex situ phenological metrics, such as satellite NDVI at lower resolutions, are available from up to early 1980s. In situ long-term data, such as herbariums, might also be interesting to analyze species position and vegetation signal for the available dates.

CONCLUSIONS

In this work, we presented a novel methodology for retrieving long-term trends in grasslands phenology. In order to go as far back in time as possible, and to promote the use of national photo libraries, a time series was created with six air missions from 1958 to 2011. Five of which were carried out with panchromatic sensors. Due to the heterogeneity of archive aerial products, the photographs were corrected from noise and vignetting, and then colorized with a conditional GAN to provide consistent multispectral information on the remotely sensed surfaces. The proposed techniques for correction and harmonization showed good performances for the entire series. The colorized photos then served as a basis for the computation of greenness indices, namely GCC and EG, which can be used as a proxy for assessing biodiversity. The output time series showed specific trends for different grasslands, allowing to group them by means of an agglomerative clustering technique. We assume that these clusters reflect different trends in grassland management or species composition in particular. However, the results must be considered with caution, due to the irregular frequency of photo acquisition and the lack of in situ and ex situ data for validation. Nonetheless, this work demonstrates the usefulness of archive aerial photographs for retrieving long-term and largescale data for environmental monitoring and evaluating changes in grasslands phenology.

Figure 1 .

 1 Figure 1. Study site location and surroundings. Infrared false color (bottom left panel) and natural color (right panel) images are provided for visualization.

Figure 2 .

 2 Figure 2. Proposed workflow for reconstructing phenological trajectories of grasslands from old aerial photographs

Figure 3 .

 3 Figure 3. Step-by-step results for the proposed correction and harmonization pipeline. Results are given for the 1958 photograph, as it was the most deteriorated. The original photograph and correction results are presented so as to show the contribution of each step.Error map are also provided, and correspond to the difference between the original and the corrected picture. Finally the colorization result is presented for the studied area.

Figure 4 .

 4 Figure 4. Distribution of pixel values, before and after correction, for each of the studied years.

Figure 5 .

 5 Figure 5. Colorization results for the corrected archive aerial photographs.

Figure 6 .

 6 Figure 6. Average temporal profiles of GCC and EG for the sampled grasslands.

Table

	Date	Scale/Pixel size(m)	Bands	Type	Provider
	01/05/1958	1:24,613 / 0.94*	Pan	Scanned film	IGN
	29/05/1976	1:19,707 / 0.40*	Pan	Scanned film	IGN
	29/06/1986	1:16,918 / 0.34*	Pan	Scanned film	IGN
	29/06/1992	1:28,108 / 0.62*	Pan	Scanned film	IGN
	08/06/2000	1:30,995 / 0.60*	Pan	Scanned film	IGN
	05/05/2011	/ 0.16	RGB	Digital	CIGAL

Table 2 .

 2 Average and standard deviation of pixel values, before and after correction, for each of the studied years.

	Date	Before correction	After correction
	1958	137.54 ± 26.83	145.21 ± 25.93
	1976	137.25 ± 12.42	137.26 ± 11.98
	1986	200.72 ± 16.65	200.76 ± 16.49
	1992	163.16 ± 18.22	163.25 ± 18.06
	2000	157.16 ± 11.39	157.26 ± 10.38

Table 3 .

 3 Mean and standard deviation for GCC and EG, before and after correction, for each of the studied years.

		Before correction	After correction
	Date	GCC	EG	GCC	EG
	1958	0.36 ± 0.02	0.08 ± 0.08	0.35 ± 0.02	0.07 ± 0.06
	1976	0.36 ± 0.01	0.13 ± 0.05	0.36 ± 0.01	0.12 ± 0.05
	1986	0.35 ± 0.01	0.12 ± 0.06	0.35 ± 0.01	0.12 ± 0.06
	1992	0.36 ± 0.02	0.10 ± 0.08	0.35 ± 0.02	0.09 ± 0.08
	2000	0.35 ± 0.01	0.07 ± 0.06	0.34 ± 0.01	0.05 ± 0.06
	2011	0.37 ± 0.02	0.11 ± 0.06	-	-

This contribution has been peer-reviewed. The double-blind peer-review was conducted on the basis of the full paper.https://doi.org/10.5194/isprs-annals-V-3-2020-549-2020 | © Authors 2020. CC BY 4.0 License.

ACKNOWLEDGEMENTS

We thank the French National Institute of Geographic and Forest Information (IGN) for providing us with the archive aerial photographs, along with the corresponding calibration certificates.


Keywords: Melatonin, pineal gland, circadian system, biological rhythms, suprachiasmatic nuclei, receptors, quantitative PCR

Objectives :

The pineal gland transduces photoperiodic informations to the neuroendocrine axis through the nocturnally melatonin secretion. This hormonal message plays a major role in the biological rhythm regulation. By autoradiography, more than 130 melatonin putative targets have been reported in the central nervous system (CNS) and in peripheral tissues. However, cross-species consensus concern only a few of them like the suprachiasmatic nuclei (SCN), the master circadian clock, and the pars tuberalis of the pituitary. Recently, MT1 melatonin receptor cDNA have been cloned in several mammals providing us with new tools to investigate its tissular location at the gene level. In the present study, we report a screening for MT1 mRNA by RT-PCR amplification of numerous tissue mRNA. Methods : mRNA were extracted from a large variety of rat tissues. To semi-quantify the melatonin receptor mRNA expression level, each cDNA was amplified concomitantly with both ß-actin and MT1 specific primers. Results : In central and peripheral tissues previously reported to bind melatonin, strong PCR signals were logically observed. More surprisingly, a vast majority of studied tissues express MT1 mRNA and then might be responsive to melatonin. Conclusion : Numerous biological functions express diurnal rhythmicity and internalsynchronization. As, most of them apparently do not receive any out-coming neuronal message from the SCN, endocrine communication was proposed to support biological rhythm synchronization. Our present data strengthen the idea that the nocturnally restricted melatonin secretion could be one internal zeitgeber that putatively distributes the endogenous circadian rhythmicity to all tissues expressing melatonin receptors.

INTRODUCTION :

Rhythmicity is a fundamental endogenous property of living organisms as it provides the living beings with the possibility to anticipate on incoming periodical events. Photoperiodic daily and seasonal variations are the major cues implicated in the synchronization of endogenous rhythmicity with environmental changes. In the pineal gland, the ambient lighting condition changes are transduced into an endocrine message : the secretion of melatonin [START_REF] Reiter | The melatonin rhythm: both a clock and a calendar[END_REF]. Melatonin synthesis is nocturnally restricted, and its duration is proportional to the length of the night, reflecting therefore annual photoperiodic changes [START_REF] Reiter | The melatonin rhythm: both a clock and a calendar[END_REF][START_REF] Reiter | Melatonin synthesis: multiplicity of regulation[END_REF]. Therefore, the melatonin secretion pattern provides both a daily and seasonal endocrine message, and it is now well established that these messages are directly involved in the regulation of both mammalian circadian and seasonal biological rhythms [START_REF] Reiter | The melatonin rhythm: both a clock and a calendar[END_REF][START_REF] Reiter | Melatonin synthesis: multiplicity of regulation[END_REF][START_REF] Cassone | Effects of melatonin on vertebrate circadian systems[END_REF].

As an hormone, melatonin is supposed to deliver its endocrine message to any tissue expressing its own receptors. Over the past 15 years, melatonin receptors, and especially their tissular detection, have been extensively studied. At first, autoradiographic methods to detect melatonin receptors were developed by using 2-125 iodo-melatonin as a radioligand [START_REF] Vanecek | Hypothalamic melatonin receptor sites revealed by autoradiography[END_REF]. However, this melatonin agonist is up to date the only radioligand available, which considerably limited achievement of exhaustive mappings of the melatonin receptor tissular distribution. Indeed, even if more than 130 melatonin binding site exhibiting areas have been detected in the central nervous system (CNS), the most striking evidence raised by comparing the published mappings is the great variability among species [START_REF] Weaver | Melatonin receptors and signal transduction in photorefractory Siberian hamsters (Phodopus sungorus)[END_REF][START_REF] Stankov | Characterization and mapping of melatonin receptors in the brain of three mammalian species: rabbit, horse and sheep. A comparative in vitro binding study[END_REF][START_REF] Masson-Pevet | An attempt to correlate brain areas containing melatonin-binding sites with rhythmic functions: a study in five hibernator species[END_REF][START_REF] Morgan | Melatonin receptors: localization, molecular pharmacology and physiological significance[END_REF][START_REF] Pang | Neuroendocrinology of melatonin in reproduction : recent developments[END_REF]. Two structures only have been reported to be almost consistently labeled in all mammalian species: the pars tuberalis (PT) of the pituitary, which exhibits the highest density of melatonin receptors [START_REF] Masson-Pevet | An attempt to correlate brain areas containing melatonin-binding sites with rhythmic functions: a study in five hibernator species[END_REF][START_REF] Morgan | Melatonin receptors: localization, molecular pharmacology and physiological significance[END_REF][START_REF] Morgan | The pars tuberalis: the missing link in the photoperiodic regulation of prolactin secretion?[END_REF] and is implicated in mediating melatonin photoperiodic effects on the neuroendocrine system [START_REF] Morgan | The pars tuberalis: the missing link in the photoperiodic regulation of prolactin secretion?[END_REF][START_REF] Lincoln | Evidence that melatonin acts in the pituitary gland through a dopamine-independent mechanism to mediate effects of daylength on the secretion of prolactin in the ram[END_REF] and the suprachiasmatic nuclei (SCN), which is the site of the circadian biological clock.

In 1994 and 1995, two mammalian receptor subtypes with high-affinity for melatonin have been cloned and characterized [START_REF] Reppert | Cloning and characterization of a mammalian melatonin receptor that mediates reproductive and circadian responses[END_REF][START_REF] Reppert | Molecular characterization of a second melatonin receptor expressed in human retina and brain: the Mel1b melatonin receptor[END_REF] and are now referenced as MT1 and MT2 receptors. When expressed in transfected cells, these subtypes exhibit the pharmacological and functional characteristics of the 2-iodo-melatonin binding sites described so far [START_REF] Reppert | Melatonin receptors step into the light: cloning and classification of subtypes[END_REF], and MT1 mRNA was clearly reported to be expressed in the PT and the SCN [START_REF] Reppert | Cloning and characterization of a mammalian melatonin receptor that mediates reproductive and circadian responses[END_REF][START_REF] Schuster | Photic regulation of mt1 melatonin receptors in the Siberian hamster pars tuberalis and suprachiasmatic nuclei: involvement of the circadian clock and intergeniculate leaflet[END_REF]. On the other hand, targeted molecular disruption of the MT1 receptor in mice resulted in the total disappearance of 2-iodo-melatonin binding in brain tissues, in particular in the SCN [START_REF] Liu | Molecular dissection of two distinct actions of melatonin on the suprachiasmatic circadian clock[END_REF]. Furthermore, 4P-PDOT and 4P-ADOT, that specifically displace 2-iodo-melatonin binding from MT2 binding sites expressed in transfected cells [START_REF] Dubocovich | Melatonin receptor antagonists that differentiate between the human Mel1a and Mel1b recombinant subtypes are used to assess the pharmacological profile of the rabbit retina ML1 presynaptic heteroreceptor[END_REF], do not displace 2-iodo-melatonin binding in SCN slices [START_REF] Dubocovich | Selective MT2 melatonin receptor antagonists block melatonin-mediated phase advances of circadian rhythms[END_REF]. Therefore, to date, MT1 mRNA is the only mRNA that can be correlated to the well known 2-iodo-melatonin binding sites. Several investigations have been conducted to characterize by in situ hybridization the tissues expressing the MT1 mRNA. However, only tissues previously known to express strong 2-iodo-melatonin binding were shown to present detectable MT1 antisense riboprobe hybridization, such as PT, SCN, paraventricular nuclei of the thalamus or cerebellum [START_REF] Reppert | Cloning and characterization of a mammalian melatonin receptor that mediates reproductive and circadian responses[END_REF][START_REF] Mazzucchelli | The melatonin receptor in the human brain: cloning experiments and distribution studies[END_REF]. In the present study, we decided to reevaluate the inventory of the central and peripheral rat tissues that express MT1 by a highly sensitive RT-PCR method.

MATERIALS AND METHODS:

Animals Male young adult Wistar rats (200-250 g, Faculté de Médecine, Strasbourg, France) were kept from birth in a 12 h light/12 h dark regime (12L/12D). A dim red light remained on throughout the 24 h period. Food and water were provided ad libitum.

cDNA isolation, RT-PCR :

All the animals were sacrificed during the second part of the light period (5 hours before lights off). Central and peripheral tissues were immediately dissected and frozen on dry ice. Total RNA was extracted using RNAble (Eurobio, France) according to the protocol of Chomczynski and Sacchi [START_REF] Chomczynski | Single-step method of RNA isolation by acid guanidinium thiocyanate-phenol-chloroform extraction[END_REF]. One µg of total RNA was reverse-transcribed in a reaction volume of 10 µl, containing 2.5 µM of oligo dT, 250 µM of each dNTP (Eurobio, France), 200 units of M-MLV reverse transcriptase (Promega, France), 50 mM Tris-HCL (pH 7.5), 75 mM KCl, 3 mM MgCl2, 10 mM of dithiothreitol and 10 units of RNAsin ribonuclease inhibitor (Promega). The samples were incubated for 60 min at 37 o C and then at 95°C for 5 min.

Typically, PCR reactions were carried out in a total volume of 25 µl. The PCR mixtures contained: 1 µl of cDNA, 100 µM of each dNTP, 10 mM Tris-HCL (pH 9.0), 50 mM KCL, 0.1% Triton X-100, 1.5 mM MgCl2, 0.5 µl S 35 dATP (12.5 mCi/ml), (NEN, France) and 0.625 units of Taq DNA Polymerase (Eurobio, France). In order to control for variations in cDNA synthesis and/or RNA preparations between samples, β-actin cDNA were simultaneously amplified with MT1 cDNA. Therefore, both MT1 and β-actin primers were mixed together in each PCR reaction tubes. The MT1 primers were selected to allow the amplification of a 375 bp specific product. The β actin primers were selected to frame an intron domain in order to generate either a long PCR product in case on genomic DNA contamination, or a 301 bp long PCR amplification of β-actin cDNA. Several primer concentrations have been tested and finally all PCR reactions contained : 1 µM of both MT1 upstream and downstream oligonucleotide primers : TGCTACATTTGCCACAGTCTCA and CCAGCACAGGGCAAAAAGTA, respectively, and 0.1µM of both upstream and downstream ß actin oligonucleotide primers : TGAACCCTAAGGCCAGT and CGCACGATTTCCCTCTCAGC, respectively. The PCR conditions were: 94°C, 3min and 32 cycles at 94 o C, 1 min; 58 o C, 1 min; 72 o C, 1 min, followed by a 20 min extension at 72 o C. Products were run on a 6% polyacrylamide and 7M-urea gel at 50°C. The gel was washed with acid acetic 10% and methanol 10%, then rinsed twice in water and dried 2 hours at 80°C. The amplified radio labeled products were revealed by autoradiography. Amplified bands were semi-quantitative analyzed by using the computerized analysis system Biocom-program RAG 200 [START_REF] Schuster | Photic regulation of mt1 melatonin receptors in the Siberian hamster pars tuberalis and suprachiasmatic nuclei: involvement of the circadian clock and intergeniculate leaflet[END_REF][START_REF] Gauer | Daily rhythms of melatonin binding sites in the rat pars tuberalis and suprachiasmatic nuclei; evidence for a regulation of melatonin receptors by melatonin itself[END_REF]. Data normalization of MT1 cDNA abundance was carried out by calculating the product density ratios of MT1/β-actin signals.

The whole procedure was repeated 6 times from 3 independent lots of tissues. An independent serie was treated with Rnase free Dnase (Promega, France) prior cDNA synthesis to clearly demonstrate that observed amplifications derived from real tissular MT1 mRNA expression and not from genomic DNA. All the results were expressed as mean ± SEM.

RESULTS and DISCUSSION :

The typical amplification pattern that we observed on autoradiograms clearly showed that only two bands were observed after the 32 amplification cycles, the 301 bp and the 375 bp products resulting from β actin and MT1 cDNA amplification, respectively. The absence of smears, background, or even surnumerous amplification products attest that the primer specificity, the annealing temperature, the number of cycles and salt concentrations (especially Mg 2+ ) were suitably selected (Fig. 1).

Experiments have been carried out to check that the selected primer concentrations and annealing temperature would guarantee that PCR experiments were conducted and stopped in the exponential phase of the amplification.

The β-actin amplification pattern appeared to be highly reproducible between samples and independent experiments, attesting that cDNA MT1 amplification level variations that we observed between the different tested tissues were not due to cDNA synthesis variations. Therefore, the absence of MT1 amplification product at the end of the PCR amplification clearly signifies that MT1 mRNA level in the starting tissue samples is either absent or below the detection limit. On the other hand, PCR positive results could signify genomic DNA contamination in the extracted total RNA. β-actin primers were selected to span an intron boundary in order to allow detection of genomic DNA contamination [START_REF] El Husseini | PCR detection of the rat brain basic fibroblast growth factor (bFGF) mRNA containing a unique 3' untranslated region[END_REF]. The band that would witness DNA contamination, was never observed in any tissue sample. However, a complete set of total RNA was treated with RNAse free DNAse prior cDNA synthesis. The obtained amplification pattern was similar for DNAse treated and untreated RNAs. Direct PCR product sequencing reactions were carried out to check the identity of the predicted 375 bp MT1 amplification product. Therefore, positive or negative results, we observed on PAGE gels, really attest the presence or the absence of MT1 mRNA expression in the starting tissue samples, respectively.

Quantitative estimation of MT1 PCR products was normalized by calculating the ratio between MT1/β-actin PCR amplification level (Fig. 2). There is, however, no precise evidence attesting that an housekeeping gene like β-actin has the same mRNA expression level in all the tissues. Actually, our repetitive experiments showed that slight β-actin mRNA expression tissular variations do exist, but that they are however strongly weaker than MT1 mRNA expression inter-tissular variations. Therefore, even if normalizing quantitative estimation of MT1 PCR by calculating the MT1/β-actin ratios does not provide an absolute quantification of inter-tissular MT1 mRNA expression level, we nevertheless used this tool to classify the different tissues in 4 different categories : tissues with strong, middle, low or no MT1 mRNA expression (Fig. 2). Surprisingly, only very few tissues have no detectable MT1 mRNA expression. Actually, among all the tested tissues, only enthorhynal and occipital cortex never exhibited any MT1 cDNA amplification signal.

Most of the tissues presenting a very high amplification signal, (listed in Fig. 2), were previously reported to bind specifically the 2-iodo-melatonin ligand. Particularly, the PT of the adenohypophysis, which is the tissue with the highest MT1 mRNA expression is also known to be the tissue with the highest melatonin receptor density. Data accumulated so far brought evidences that these receptors are implicated in mediating some of the photoperiodic effects on the neuroendocrine system like the photoperiodic control of prolactin secretion [START_REF] Morgan | The pars tuberalis: the missing link in the photoperiodic regulation of prolactin secretion?[END_REF]. The caudal artery is also known to be a clear melatonin signal target as melatonin is known to potentiate adrenergic vasoconstrictor effect on the tail artery [START_REF] Viswanathan | Expression of melatonin receptors in arteries involved in thermoregulation[END_REF][START_REF] Krause | Melatonin receptors mediate potentiation of contractile responses to adrenergic nerve stimulation in rat caudal artery[END_REF][START_REF] Bucher | Vasoconstrictor effects of various melatonin analogs on the rat tail artery in the presence of phenylephrine[END_REF]. In colonic tissues, the MT1 mRNA expression correlates with the high melatonin binding site density reported in both the mucosa, and the blood vessels of rodent and human colon sections [START_REF] Poon | Autoradiographic localization of 2[125I]iodomelatonin binding sites in the gastrointestinal tract of mammals including humans and birds[END_REF]. At the present time, it is difficult to correlate these melatonin binding sites to any particular function of the melatonin signal. Indeed, data have suggested that melatonin could inhibit the spontaneous motility of rat intestine tissue through an indirect inhibition of the contractile response to serotonin [START_REF] Harlow | Effect of melatonin on the force of spontaneous contractions of in vitro rat small and large intestine[END_REF][START_REF] Bubenik | The effect of serotonin, N-acetylserotonin, and melatonin on spontaneous contractions of isolated rat intestine[END_REF]. However, the large melatonin concentrations requested to get these effects suggest that they do not result from MT1 picomolar melatonin receptor stimulation [START_REF] Poon | Autoradiographic localization of 2[125I]iodomelatonin binding sites in the gastrointestinal tract of mammals including humans and birds[END_REF]. More surprinsingly, lung and heart expressed also very high MT1 mRNA. Indeed, until now, only bird lungs and hearts had been reported to putatively bind melatonin [START_REF] Pang | Melatonin receptors in peripheral tissues: a new area of melatonin research[END_REF][START_REF] Pang | Differential effects of short photoperiod on 2-[125I]iodomelatonin binding in the testis and brain of quail[END_REF]. The shared characteristics between colon, heart and lung is their high content in blood vessels, suggesting that MT1 mRNA expression positive signal in these tissues could be interpreted as vascular melatonin receptors.

The pineal gland presents also a quite strong PCR signal. This endocrine organ is known to present a high melatonin receptor density at birth [START_REF] Zitouni | Influence of maternal melatonin on melatonin receptors in rat offspring[END_REF], but the receptor density decreases rapidly afterwards, and in adult any signal is detected by classical autoradiographic methodology [START_REF] Zitouni | Influence of maternal melatonin on melatonin receptors in rat offspring[END_REF]. This situation is quite similar with melatonin receptor expression in the SCN of the Syrian hamster, where a severe loss of receptor is observed in adulthood, even if the MT1 mRNA is still largely expressed [START_REF] Gauer | Cloning experiments and developmental expression of both melatonin receptor Mel1A mRNA and melatonin binding sites in the Syrian hamster suprachiasmatic nuclei[END_REF]. Among all the other tissues exhibiting either middle or low MT1 cDNA amplification signal are previously described binding melatonin tissues such as pituitary, area postrema, cerebellum, olfactive bulbs, retina, spinal cord, SCN, spleen, prostate, epididymis, testis, thymus, harderian and adrenal glands [START_REF] Masson-Pevet | An attempt to correlate brain areas containing melatonin-binding sites with rhythmic functions: a study in five hibernator species[END_REF][START_REF] Pang | Neuroendocrinology of melatonin in reproduction : recent developments[END_REF]. However, some new putative targets of the melatonin signal which are detected in this study have to be added to this non exhaustive list such as liver, pancreas, salivary and seminal glands, bone narrow, fatty tissues, kidney, squeletic muscle and thyroid. Actually, at a first glance, the most intriguing results is that the SCN, one of the most studied melatonin target, does not belong to the tissues with the highest MT1 cDNA amplification. However, among all the tested tissues, the SCN belongs to the few whose dissection into a pure preparation is rendered very difficult due to the minute size of these cerebral nuclei. Therefore, the SCN specific cDNA concentration submitted to the MT1 primers is decreased by non SCN mRNA in the total RNA extracted from explants, which could explain the obtained results.

The main question raised by our data is the significance of the expression of MT1 mRNA in almost all the rat tissues tested. Melatonin is known to play a major role in regulating several circadian and seasonal rhythms. In the SCN and PT, well known melatonin responsive tissues, precise actions of melatonin have been established, even if the underlying mechanisms are not yet elucidated. The same holds true for the well characterized action of melatonin in potentiating the norepinephrine induced blood vessel vasoconstriction. In all the other tissues, many role have been proposed, but the real physiological function of the melatonin hormone remains speculative. In some way, until now, each new discovered melatonin target had to be presented with its particular role in decoding the melatonin message. Actually, new discoveries in biological rhythmicity should invite us to investigate the melatonin message function on different basis. Our data clearly suggest that the vast majority of the mammalian tissues might be responsive to melatonin. The recently discovered clock genes, which are the molecular basis of the SCN circadian oscillation [START_REF] Shearman | Interacting molecular loops in the mammalian circadian clock[END_REF], are not actually restricted to the circadian biological clock. They are also expressed in a vast variety of tissues such as heart, liver, retina, PT, fibroblast, etc [START_REF] Balsalobre | Clock genes in mammalian peripheral tissues[END_REF]. However, their specificity in the SCN is to generate circadian oscillations which persists for long times in vitro [START_REF] Shearman | Interacting molecular loops in the mammalian circadian clock[END_REF] while in other tissues, they generally oscillate on a circadian basis for a only few days [START_REF] Balsalobre | A serum shock induces circadian gene expression in mammalian tissue culture cells[END_REF]. It thus seems that intra-tissular oscillation loops in these tissues requires an external zeitgeber to both maintain the oscillation, and synchronize it with the whole body time, since most biological parameters appear to be synchronized together. However, most of them apparently do not receive any SCN out-coming neuronal message suggesting a endocrine synchronizing message. The melatonin hormone could be one internal zeitgeber that distributes the endogenous circadian rhythmicity originating in the SCN. As its secretion is nocturnally restricted melatonin could distribute the day/night information to all tissues that putatively express its receptors. Our present data support the idea that the number of both central and peripheric structures which receive this temporal message is extremely large.

In a large majority of the tissues studied by PCR, the attempt to detect MT1 mRNA by in situ hybridization failed. From the general opinion, the MT1 mRNA expression level is very low, making difficult to differentiate specific signal from background noise. The apparent discrepancy between the high melatonin receptor density measured for example in the rodent PT and the apparently low MT1 mRNA expression level [START_REF] Schuster | Photic regulation of mt1 melatonin receptors in the Siberian hamster pars tuberalis and suprachiasmatic nuclei: involvement of the circadian clock and intergeniculate leaflet[END_REF][START_REF] Poirel | MT1 melatonin receptor mRNA expression exhibits a circadian variation in the rat suprachiasmatic nuclei[END_REF] could be explained by a long half-life of the melatonin receptor. Pilot experiments using antisense oligonucleotide invalidation of MT1 melatonin receptor were carried out in the SCN of the Siberian hamster which lacks the MT2 melatonin receptor subtype [START_REF] Weaver | Nature's knockout: the Mel1b receptor is not necessary for reproductive and circadian responses to melatonin in Siberian hamsters[END_REF]. The SCN melatonin binding sites expressed in the Siberian hamster are then only due to MT1 mRNA translation. We observed that daily injections of oligonucleotides for 7 days decreased the melatonin receptor density in all the injected animals, but only by 30 %, strongly suggesting that the melatonin receptor is a very stable protein (Fig. 3). Continuous renewal of the protein might therefore only require a quite low MT1 mRNA expression level which would explain the encountered difficulties to detect MT1 mRNA by in situ hybridization. More than ever, new specific pharmacological tools are required to investigate widely the physiological function of the pineal hormone melatonin, and especially to demonstrate its possible implication as a master endocrine internal synchronizer of biological rhythmic functions.

ACKNOLEDGEMENTS :

The authors wish to thank D. Bonn for his excellent technical assistance. This study was partly supported by IRIS (Courbevoie, France), the Simone and Cino Del Duca Foundation (France), "la Région Alsace" (France) and by a grant from La Fondation pour la Recherche Médicale (France).

FIGURES LEGEND :

Figure 1 : Electrophoretic mobility of MT1 and β-actin S 35 -dATP radiolabeled RT-PCR fragments amplified from cDNA transcribed from total RNA extracted from different rat tissues. The cDNA were synthesized with reverse transcriptase using oligo dT. PCR parameters and gel electrophoresis conditions are described in the text. Tissues by alphabetical order and their corresponding electrophoresis pattern number : Adenohypophysis: 1 ; Adrenal glands : 35; Area postrema :

  Poirel Ph. D, Cathy Cailotto MS, Dominique Streicher BS, Paul Pévet Ph. D, Mireille Masson-Pévet Ph. D and François Gauer Ph. D Laboratoire de Neurobiologie des Rythmes, CNRS-UMR 7518, Université Louis Pasteur, Strasbourg, France. Please send correspondence to: Prof. François Gauer, "Laboratoire de Neurobiologie des Rythmes" ; CNRS-UMR 7518, 05 28, E-mail: gauer@neurochem.u-strasbg.fr Running headline : Melatonin receptor tissular localization

Figure 3 :

 3 Figure 3 :In antisense experiments, injection cannula (final external and internal diameters were 150 and 53 µm respectively) were implanted using standard stereotaxic procedures close by the SCN. Distribution of the antisense oligonucleotides was estimated after injecting fluorescein labeled oligonucleotides (Molecular research Laboratories, LCC, USA). Melatonin receptor density was measured in the SCN of Siberian hamster daily injected for 6 days with 1.5 µl of either a saline solution (control), a saline solutions containing 0.75 nmol of a MT1 antisense (TGCTGCGAGGCATTGAGCAG) or a mismatch (TGGTGCGACGCTATGAGCAG) oligonucleotide, respectively (Molecular research Laboratories, LCC, USA). Animals were sacrificed on day 7 at 2.00 pm. Autoradiographic procedures were conducted as previously described[START_REF] Chomczynski | Single-step method of RNA isolation by acid guanidinium thiocyanate-phenol-chloroform extraction[END_REF][START_REF] Pang | Melatonin receptors in peripheral tissues: a new area of melatonin research[END_REF]. The results are expressed in fmol/mg of protein. Values are the mean ± SEM of 7 animals.

19 21 23 25 27 29 31 33 35 37 39 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 Figure 2 :

  2; Blood : 32; Bone Marrow : 39; Brown fat : 17; Caudal artery : 3; Cerebellum : 5; Colon : 7 ; Duodenum : 10 ; Entorhinal cortex : 8 ; Epididymis : 11; Harderian gland : 18; Heart : 6 ; Ileon : 19 ; Kidney : 30 ; Liver : 13 ; Lung : 27 ; Muscle : 21 ; Neurohypophysis : 22 ; Occipital cortex : 9 ; Olfactory bulb : 4 ; Pancreas : 23 ; Pars tuberalis : 24 ; Pineal gland : 26 ; Prostate : 27 ; Retina : 31 ; SCN : 33 ; Seminal glands : 34 ; Skin : 25 ; Spinal cord : 20 ; Spleen : 29 ; Stomach : 12 ; Superior cervical ganglia : 14 ; Salivary glands : 15 ; Testis : 36 ; Thymus : 37 ; Thyroid : 38 ; White fat : 16. Expression of the MT1 melatonin receptor subtype mRNA in rat various tissues. The values are expressed as a ratio between the optical density of the MT1 and β-actin PCR amplification signal. Values are the mean ± SEM of 6 independent PCR amplifications. The number below the axis bar refers to tissue names listed below. Pars tuberalis ; 10: Pineal gland ; 11: Retina ; 12: Spinal cord ; 13: Suprachiasmatic nuclei ; 14: Blood ; 15: Caudal artery ; 16: Heart ; 17: Lung ; 18: Spleen ; 19: Colon ; 20: Duodenum ; 21: Ileon ; 22: Liver ; 23: Pancreas ; 24: Salivary glands ; 25: Stomach ; 26: Prostate ; 27: Seminal glands ; 28: Epididymis ; 29: Testis ; 30: Bone Marrow ; 31: Brown fat ; 32: Harderian gland ; 33: Kidney ; 34: Skin ; 35: Squeletic muscle ; 36: Adrenal glands ; 37: Thymus ; 38: Thyroid ; 39: White fat.

	ß-actin MT1	1 low amplification level were arbitrary define with a ratio value above 2, between 1 and 2, or 3 5 7 9 11 13 15 below 1 respectively. 1: Adenohypophysis; 2: Area postrema; 3: Cerebellum ; 4: Superior cervical ganglia ; 5: Entorhinal cortex ; 6: Neurohypophysis ; 7: Occipital cortex ; 8: Olfactory 17 Tissues with a high, middle or bulb ; 9: