Guillaume Godefroy
email: guillaume.godefroy@univ-grenoble-alpes.fr

Bastien Arnal

Emmanuel Bossy

Compensating for visibility artefacts in photoacoustic imaging with a deep learning approach providing prediction uncertainties

Keywords: Photoacoustic imaging, Deep learning, Visibility artefacts, Monte Carlo dropout, Bayesian neural network Photoacoustic imaging, Deep learning, Visibility artefacts, Monte Carlo dropout, Bayesian neural network

come

Compensating for visibility artefacts in photoacoustic

imaging with a deep learning approach providing prediction uncertainties

Introduction

Photoacoustic (PA) imaging is an emerging biomedical modality based on the generation of acoustic waves by light absorption. This modality is promising, as it enables imaging at large depths with high spatial and temporal resolution, and can provide images of the optical absorption [1] with specific molecular contrast which can be enhanced by spectroscopy.

In conventional PA imaging, a short nanosecond laser pulse is sent into the medium and the emitted ultrasonic waves are collected by a conventional ultrasound (US) probe. At the US propagation time scale, the object illumination is quasi instantaneous as the speed of light is several orders of magnitude higher than the speed of sound, resulting in the emission of strongly coherent acoustics waves [START_REF] Guo | On the speckle-free nature of photoacoustic tomography[END_REF]. These waves interfere constructively or destructively depending on the structure of the object, often leading to two well-known artefacts on the reconstructed image: the limited bandwidth and the limited view artefacts [START_REF] Deán-Ben | On the link between the speckle free nature of optoacoustics and visibility of structures in limited-view tomography[END_REF]. With a resonant detection bandwidth, when an object larger than the acoustic central wavelength of the transducer is illuminated, the strong low frequency component of the PA signals is filtered out. With a limited view (limited detection aperture) , for a structure elongated alongthe axis of the probe, the PA waves interfere constructively perpendicularly to the probe but mostly destructively throughout the elongation. As a result, very few signals are collected in this case by linear or matrix array probes due to their limited angular view. Both type of artefacts will further be referred to as the visibility problem in this paper.

The limited view problem has been addressed in several studies. The most intuitive approach is to either rotate the object relatively to the probe [START_REF] Kruger | Thermoacoustic computed tomography using a conventional linear transducer array[END_REF] (or vice versa [START_REF] Yang | Fast full-view photoacoustic imaging by combined scanning with a linear transducer array[END_REF]) or use ring shaped transducer arrays [START_REF] Xia | Whole-body ring-shaped confocal photoacoustic computed tomography of small animals in vivo[END_REF] [START_REF] Deán-Ben | On the link between the speckle free nature of optoacoustics and visibility of structures in limited-view tomography[END_REF] in order to cover all angles. However, a clinical implementation would benefit from a handheld real-time system as currently used in ultrasound imaging. Other approaches rely on the introduction of a spatial modulation of optical absorption of the sample, either using injection of sparse absorbing particles [START_REF] Dean-Ben | Localization optoacoustic tomography[END_REF], by a modulation of physical properties [START_REF] Wang | Ultrasonic-heating-encoded photoacoustic tomography with virtually augmented detection view[END_REF], or by computing statistical properties of the PA signal generated by fluctuating sources in the medium [START_REF] Chaigne | Superresolution photoacoustic imaging via flow-induced absorption fluctuations[END_REF] [START_REF] Vilov | Photoacoustic fluctuation imaging: theory and application to blood flow imaging[END_REF]. Fluctuations of the PA signals can be produced by random optical speckle pattern illuminations [START_REF] Chaigne | Superresolution photoacoustic imaging via flow-induced absorption fluctuations[END_REF] or by flowing red blood cells [START_REF] Vilov | Photoacoustic fluctuation imaging: theory and application to blood flow imaging[END_REF], naturally present in the blood vessels. Nevertheless, these methods require long acquisition times in order to get significant statistical properties, and therefore have a poor temporal resolution.

In this work, a deep learning approach is proposed to overcome the visibility problem and improve the image quality in a real-time single shot configuration. A neural network can be viewed as an algorithm composed of many parameters, called weights, designed to compute input data into a desired form of output [START_REF] Bengio | Deep learning[END_REF]. This algorithm is trained over multiple examples to obtain the best representation of the studied phenomenon. After training, the network transforms an input raw image into an output image that is expected to resemble the (unknown) ground truth. The training set consists of multiple raw data/ground truth pairs that will be used to optimize the weights of the network.

Convolutional neural networks (CNN)are amongst the most popular category of deep learning algorithms (DLA) [START_REF] Lecun | Gradientbased learning applied to document recognition[END_REF] , and have reached state of the art performances in several imaging problems including segmentation [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF], classification [START_REF] Russakovsky | Imagenet large scale visual recognition challenge[END_REF], artefacts removal [START_REF] Dong | Compression artifacts reduction by a deep convolutional network[END_REF] or denoising [START_REF] Xie | Image denoising and inpainting with deep neural networks[END_REF]. CNN have been introduced recently in biomedical imaging, showing impressive results in various tasks [START_REF] Akkus | Deep learning for brain mri segmentation: state of the art and future directions[END_REF][18] [START_REF] Liu | Deep learning in medical ultrasound analysis: a review[END_REF]. Over the past two years, a few groups started investigating deep learning applied to PA imaging for several purposes including direct reconstruction of the initial pressure [START_REF] Waibel | Reconstruction of initial pressure from limited view photoacoustic images using deep learning[END_REF], handling artefacts coming from sparse data [START_REF] Hauptmann | Model-based learning for accelerated, limited-view 3-d photoacoustic tomography[END_REF][22] [START_REF] Davoudi | Deep learning optoacoustic tomography with sparse data[END_REF], reflection artefacts removal [START_REF] Allman | Deep neural networks to remove photoacoustic reflection artifacts in ex vivo and in vivo tissue[END_REF][25], point source localization [START_REF] Reiter | A machine learning approach to identifying point source locations in photoacoustic data[END_REF][27] and quantitative measurements [START_REF] Cai | End-to-end deep neural network for optical inversion in quantitative photoacoustic imaging[END_REF] [START_REF] Kirchner | Context encoding enables machine learning-based quantitative photoacoustics[END_REF]. The correction of the limited bandwidth problem was also investigated on very simple objects [START_REF] Gutta | Deep neural network-based bandwidth enhancement of photoacoustic data[END_REF]. Some of these studies [START_REF] Guan | Fully dense unet for 2d sparse photoacoustic tomography artifact removal[END_REF][22] [START_REF] Davoudi | Deep learning optoacoustic tomography with sparse data[END_REF] showed that deep learning can also reduce the limited view artefacts although results were either numerical or obtained with non-conventional imaging devices. A linear array was used in experiments [START_REF] Kim | Deeplearning image reconstruction for real-time photoacoustic system[END_REF] but a ground truth was missing to assess the success of the approach. Finally, in most of the cited studies, experimental results were predicted from models trained only on simulation data.

In this work, we focus on the correction of the whole visibility problem, induced both by the limited view and limited bandwidth of a conventional linear US probe. The originality of our approach resides in the design of a dedicated model object and a method to create an experimental training dataset. The method is used to assess the capacity of a neural network to remove these artefacts on experimental images that were not used during the training. In this study, a ground truth is known also for the test set, which consisted of some of those unseen images. Thus, evaluation of the quality of the reconstruction can be performed. We point out that this study is not designed to produce quantitative PA images, as our ground truth does not directly represent the optical absorption, and focuses on providing morphological images. As a consequence, it can currently not be applied for quantitative imaging including multispectral investigations. Moreover, our study is limited to a given class of objects, and an investigation of the ability of the network to generalize to other classes of objects was out of the scope of this work. A preliminary discussion on generalization is provided in the supplementary materials.

Despite the impressive performances of DLA to reconstruct PA images, errors can be made by the algorithm which may misinterpret the data. This is one of the main limitations of neural network approaches in the medical field: the lack of confidence in the results. In this work, we estimate the uncertainty in our prediction through a Bayesian machine learning framework. We followed the approach proposed by Ghahramani and Gal [START_REF] Gal | Dropout as a bayesian approximation: Representing model uncertainty in deep learning[END_REF], referred to as Monte Carlo dropout (MC dropout), which has been recently applied for phase imaging [START_REF] Xue | Reliable deep-learning-based phase imaging with uncertainty quantification[END_REF]. Uncertainty estimation using a Bayesian framework has already been studied in PA imaging to reduce artefacts induced by approximation of reconstruction parameters [START_REF] Sahlström | Modeling of errors due to uncertainties in ultrasound sensor locations in photoacoustic tomography[END_REF] [START_REF] Tick | Modelling of errors due to speed of sound variations in photoacoustic tomography using a bayesian framework[END_REF]. A specific deep learning approach has also been developed to estimate uncertainty of the optical parameter estimation in quantitative PA imaging [START_REF] Gröhl | Confidence estimation for machine learning-based quantitative photoacoustics[END_REF]. Our case is different, as we do not want to take into account approximation of some parameters, but estimate the uncertainty linked to the deep learning process. To do so, our CNN is converted into a Bayesian neural network to introduce randomness in the prediction process, which makes the prediction no longer deterministic: the model will predict different outputs for the same input. Then, for a given input, several outputs are generated and are interpreted as samples of a probabilistic distribution, from which parameters can be estimated, such as the mean value and the confidence measure. The uncertainty estimation provides positions of invented and poorly reconstructed structures. This estimation is very useful for real-time navigation as a feedback for the user, who may eventually choose to display only the reliable parts of the images. However, we point out that this estimation does not provide information for structures missing in the reconstruction, and is therefore not a measurement of the fidelity.

We also study the DLA performance over different input data types. Usually, a conventionally reconstructed image is used. This prior reconstruction is obtained by applying delays and summation (DAS) on the Hilbert transform of the radiofrequency (RF) signals.

This operation produces a complex image whose modulus is computed to form the demodulated beamformed (dmBF) image, which is the one displayed for the end user. While the input of the DLA for PA image reconstruction is usually the dmBF image, we choose to train our network with the modulated beamformed (mBF) image. The mBF image (Fig. 1.a) is obtained by applying DAS directly on the real-valued RF time signals. Consequently, the mBF image is modulated by the impulse response of the transducer, resulting in axial oscillations. While the mBF image represents the object less faithfully than the dmBF image (because of the probe oscillations), we show that it carries more information that the DLA can exploit. The two approaches are compared in Supplementary Materials.

Finally, we investigate the design of the training sets. Indeed, processing experimental data with CNN trained solely on simulated data seems to produce poor reconstruction [START_REF] Davoudi | Deep learning optoacoustic tomography with sparse data[END_REF] which we confirm here (Supplementary Materials), while constructing a large experimental dataset is complex and time consuming. We varied the relative sizes of the combined experimental and simulated datasets and observed its impact on the reconstruction performances.

Material and methods

Conventional reconstruction methods

For comparison purposes, conventional DAS images (dmBF) and L2 deconvolution images are provided. DAS is fast and robust whereas deconvolution methods are more computational and more complex to implement since the knowledge of the point spread function of the system and regularization are necessary. Here, image deconvolution is achieved using a least-square minimization approach with a L2-regularization penalty term. It was performed by a fast iterative shrinkage thresholding algorithm (FISTA) [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. The inversion is defined as:

X = argmin X 1 2 Y -AX 2 2 + α 2 X 2 (1)
X is the expected reconstructed object, X the object at each iteration, Y represents the RF signals and A is the propagation matrix, containing the imaging system response at each point of the reconstruction grid. α is the regularization parameter, tuned heuristically by visually comparing the reconstructed image with the ground truth.

Creation of the experimental dataset

A model of leaf skeleton was chosen as imaging sample (see Fig. 1.a). This model has been used in previous studies as it provides a branching structure qualitatively similar to that of a vascular network, and produces conventional images with similar visibility artefacts [START_REF] Hojman | Photoacoustic imaging beyond the acoustic diffraction-limit with dynamic speckle illumination and sparse joint support recovery[END_REF]. To obtain a sufficient photoacoustic signal, the leaves veins were tainted with black ink and the limbs were dissolved by chemical treatment. The smallest veins of the leaves are finally manually cut out to remove unresolvable details. Each pair of the dataset consisted of a mBF image (input of the network) and the corresponding photograph (ground truth) of a 5.12×5.12 mm 2 patch of the considered leaf.

As shown in Fig. 1.a, the leaf is maintained inside a horizontal plane of an agarose gel, which stands within a tank filled with degassed and deionized water. Through a side window composed of a frame tightening a Mylar membrane, an ultrasonic transducer array (15.6 MHz central frequency, L22-8v, Verasonics, USA), connected to a multi-channel acquisition system (Vantage 256 High Frequency, Verasonics, USA) is coupled to the water tank with echographic gel. Thus, the leaf is in the imaging plane of the probe. It is illuminated from the top with 5 ns laser pulses at 10 Hz repetition rate (λ = 532 nm), produced by a frequency-doubled Nd:YAG laser (Surrelite, Continuum, USA). For each laser shot, PA signals are acquired and mBF images are reconstructed using DAS assuming a homogeneous medium with a speed of sound of 1500 m.s -1 , neglecting the presence of the agarose gel. To obtain several independent samples from a same leaf, the leaf is mechanically translated respectively to the probe and the light source.

We define our "ground truth" images as photographs of the leaf taken with a CMOS camera (X-E2, Fujifilm, Japan). These photographs are converted to grey scale (8 bits) and pixels below a threshold are set to 0 to suppress background noise. A registration between the PA image and the corresponding photograph is needed. The magnitude of the transformation to apply to co-register the two images (decomposed as rotation, translation and scaling) were found automatically by maximizing the correlation coefficient between the PA image (the reference) and the transformed photograph. 593 pairs of images from two leaves are obtained, split between the experimental training set and the experimental validation set with respectively 500 and 93 pairs. The validation set is used during the training to assess the optimization process is over. An experimental test set of 15 pairs is then constituted from two different leaves. It is used to evaluate the performance of our approach.

Creation of the simulation dataset

The two photographs used as ground truth for the experimental training set were also used to create the simulation dataset. Data augmentation is applied on those photographs to increase the dataset size by applying rotations, mirror transformations, horizontal or vertical shears, and center expansions or compressions. Then, 1×1 cm images are extracted to compute their PA signals.

The method used to simulate PA imaging is described in our previous work [START_REF] Vilov | Super-resolution photoacoustic and ultrasound imaging with sparse arrays[END_REF]. In brief, the imaging system response is experimentally measured for a single source at one spatial location and the synthesis of the RF signals coming from a whole object is obtained by summing the contributions of each pixel of the object. The object is then reconstructed with DAS as for the experimental data. The medium is assumed to be homogeneous with a constant speed of sound of 1500 m.s -1 . DAS is then applied to reconstruct an mBF image of 5.12×5.12 mm 2 area and the photograph is cropped to match the dimensions. Propagating PA waves from a larger area (1×1 cm 2 ) than the one viewed by the network (5.12×5.12 mm 2 ) enable to take into account the presence of surrounding structures which can affect reconstruction during experiments. A series of 1400 pairs of images are obtained. Around five days are needed to compute the dataset on Matlab with a desktop computer.

Deep learning framework

UNET [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] Tensorflow and Keras. Dropout [START_REF] Srivastava | Dropout: a simple way to prevent neural networks from overfitting[END_REF] and batch normalization [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF] layers are added to limit over-fitting of the model. The last layer contains only one filter instead of two in the original version, as the expected output is a single image. The last activation function is also suppressed as the prediction is no longer a binary image. It is worth mentioning that several modifications supposed to improve the result including skip connections between input and output [START_REF] Akkus | Deep learning for brain mri segmentation: state of the art and future directions[END_REF], residual blocks [START_REF] Diakogiannis | Resuneta: a deep learning framework for semantic segmentation of remotely sensed data[END_REF] and fully densely connected blocks [START_REF] Guan | Fully dense unet for 2d sparse photoacoustic tomography artifact removal[END_REF] have been tested without significant improvement of the prediction. The cost function is the classical mean squared error, and an Adam optimizer is used with a learning rate of 5.10 -4 and momentum of 8 [START_REF] Kingma | A method for stochastic optimization[END_REF] with batch sizes of 8 images. An early stopping approach based on the validation loss was chosen to limit under-and over-fitting [START_REF] Prechelt | Early stopping-but when?[END_REF]. The prediction phase must be random to model uncertainty. In the MC dropout approach, noise is injected in the model by activation of the dropout layers (dropout rate of 50%) both during training and prediction. In this study, 20 inferences are generated from forward passes through the model with a different dropout mask. The different resulting predictions allow to further estimate the distribution mean and its standard deviation which gives a map of uncertainty (see Fig. 1.b). The training and the evaluation of the network, composed by around 30,000,000 neurons, are performed on a NVIDIA Quadro P2000. Around 50 minutes are needed for the training on the simulation set and 30 minutes on the experimental set.

Quantitative assessment of the network performance

As mentioned previously, the same leaves are used to create the simulation and the experimental dataset. It means that from the same object (i.e. an area of the leaf), we will be able to obtain either the experimental RF signals or the simulated one. For comparison purpose, the reconstructed objects shown in the figures are the same for both simulations and experiments. A third example is used for the MC dropout results, the estimated uncertainties of the two previous examples being described in the supplementary materials. All images are normalized by their maximum and represented with a colorbar from 0 to 1, since quantitative information is not expected.

To evaluate the accuracy of a trained neural network, the normalized 2D cross-correlation [START_REF] Yoo | Fast normalized cross-correlation, Circuits[END_REF] (NCC) and the scaled and shifted structured similarity index (sSSIM) are computed between each output and ground truth. The first one uses local sum to normalize the cross-correlation for feature matching. SSIM [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF] is a widely-used metric to evaluate the perceived quality of an image. It is computed over several small windows of the image, quantifying the structure, contrast and intensity similarities. The sSSIM [START_REF] Schwab | Real-time photoacoustic projection imaging using deep learning[END_REF] is used for obtaining a scaled and unbiased score which was not disadvantaging for the other reconstruction methods. For an overall performance estimation of the network, the mean and standard deviation values among the test set are presented. On each result, we computed an uncertainty map with the MC dropout method, which only require a single-shot acquisition. For comparison purposes, as an alternative way to assess the local variability in the reconstructions, we also computed an uncertainty map from the standard deviations of the reconstructions of 50 acquisitions of the same sample. In this case, the variability comes from the experimental noise while the CNN remain deterministic. We also computed the absolute truth error of the reconstruction defined as the difference between the ground truth and the predicted image.

Results

Simulation results

Fig. 2 shows a comparison between the reconstructed image from simulated data provided by the DLA (d,h) (trained with simulated data), the dmBF image (DAS, b,f), a L2-regularized deconvolution (c,g) and the photograph of the object (ground truth, a,e). L2 minimization and DAS clearly do not provide vertical structures, i.e. the structures elongated in the axis of the probe. Veins having inclinations beyond the detection aperture (typically beyond 45 degrees) are missing due to the limited view problem. The inside of the thicker vessel is also missing and the thickness of the thinner vessel is underestimated for DAS reconstruction and overestimated for the L2 minimization (arrows). In contrast, the deep learning reconstruction yields an almost artefact-free reconstruction with errors located only on the smallest appendages resulting from the manual cutting, and on few vertical structures which are not completely recovered (stars on the images). The performances of the different reconstruction methods and their standard deviation, evaluated over the 15 pairs of the simulation test set with the metric described in 2.5 are shown in Tab. 1. These numbers clearly confirm the qualitative visual impression: when the DLA is used, the NCC and sSSIM are about three times higher compared to the simple DAS. Scores for the deconvolution method, not shown here, are on the same order than that of the DAS.

Experimental results

Fig. 3 shows a comparison between the reconstructed image from the experimental data provided by the DLA (d,h) (trained with experimental data), the conventional DAS reconstruction (b,f), a L2-regularized deconvolution (c,g) and the photograph of the object (ground truth, a,e). Similarly, the DAS approach and the L2-minimization both fail to recover vertical structures as well as to provide a good rendering of the vein thicknesses by filling the inside of the thicker ones. In contrast, the DLA trained on the experimental data yields a reconstruction with most of the vertical structures recovered and a correct thickness of the veins (arrows). A few structures are again not recovered, and some mistakes occurred especially for the reconstruction of vertical veins (stars). The orientation is not always perfectly respected, and in some places, some veins appear when there are none in reality. Quantitative performances are shown in Tab. 1 where as for simulations, we observe a large improvement for the deep learning approach comparing to the DAS, although lower than simulation results. Both the sSSIM (0.76) and the NCC (0.8) are significantly enhanced. It may also be noted that the DAS performs better on experimental data than on simulation data.

Uncertainty estimation

Results of the MC dropout procedure are presented in Fig. 4. A low standard deviation indicates a good robustness of the technique: the prediction remains stable over several realizations. Areas with high value are similar in the estimated uncertainty map (Fig. 4.e), the absolute error map (Fig. 4.d) and in the map of standard deviation from experimental noise (Fig. 4.f). Most of the errors in the prediction are captured like the small vessel at the bottom, that was not fully recovered, or the central one which was reconstructed but in a curved shape instead of a straight one (arrows).

Impact of the pretraining and the size of the training set on the performance

In this part, the efficacy of a pretraining session with a simulation dataset is investigated as a means to improve the general performance and for reducing the size of the experimental training set. The uncertainty prediction was not studied in this configuration. We increase the training set size with unseen part of the two leaves used for the testset. The DLA was trained with experimental datasets of different sizes, from 10 to 550 pairs (the entire dataset). For each size, the training is repeated 30 times with, for each of them, a training and validation set composed of different pairs randomly chosen. This is needed to limit the influence of individual pairs on the training set size (for example, it is likely that 10 examples very close to the test set will provide a better prediction than 20 very distant ones), especially for small set sizes.

The displayed sSSIM values are therefore an average over all the test sets from the 30 different realizations. To evaluate pretraining, we repeated this procedure with weights initialized by those obtained at the end of a training session on a simulation dataset composed of 1400 pairs. The results are shown in Fig. 5. As expected, the performance increases with the experimental dataset size. Below 200 pairs, errors remain present and the veins thickness is not always faithfully represented. From 200 pairs, the image quality seems visually stable, although the sSSIM value still slowly increases. With pretraining (blue curve) convergence is faster. When the full dataset is used, pretraining only slightly increases the performance of the network (sSSIM improved from 0.76 to 0.78). For a smaller size such as 50 pairs, the score improvement is better (from 0.63 to 0.72). A reconstructed image comparable to the one obtained with the total experimental training set is almost reached from this experimental dataset constituted of 50 pairs. In this situation, a pretraining session therefore enables to decrease the size of the experimental training set by a factor 4.

Discussion

The algorithm trained with simulated data is able to produce images that are free from visibility artefacts, when applied to simulated test data. When trained with simulated data, the algorithm applied to experimental data however fails to provide images of good quality, as illustrated in the Supplementary Materials (see fig. S3). When both training and prediction are made with experimental data, while a few errors may remain, most vessels are well recovered. A fundamental difference between simulation and experiment is the nature of the ground truth. In simulation, the ground truth exactly represent the distribution of absorbed light. In the experiments, the ground truth is a photographic 2D projection of a three dimensional absorption distribution. It then represents the integration of the optical absorption through the sample thickness (which may vary among leaf veins). Consequently, the photograph is not a quantitative representation of the sample absorption, and cannot therefore not be a quantitative ground truth. Thus, our method is not supposed to provide a quantitative reconstruction of the absorbed light, as the network is forced to learn from a 2D representation of a 3D object of finite and varying thickness. Given the nature of our ground truth, the purpose here was to demonstrate that the morphology of the sample could be retrieved free of visibility artefact. Quantitative ground truth would be required to provide quantitative information on the absorption coefficient, as needed for instance for spectroscopic applications.

The use of the mBF image as input of the network improves the performance both quantitatively and qualitatively compared to the dmBF image input (see supplementary materials, Fig. S2). In fact, the mBF image, although affected by the oscillatory impulse response to the ultrasound probe, turns out to carry more information to be captured by the network.

The results show that errors are often located at the edges of the reconstructed image. Indeed, in these areas, less information about the surrounding structures is available. One way to limit these artefacts could be to reconstruct on a larger area and crop the edges. Most of the errors remaining on experimental images are located where the manual cutting was performed resulting in small appendices. These structures, which do not belong to the initial object, turned out to mislead the network which seeks to elongate them to join all the veins together. It is reasonable to think that the number of errors would have been lower on a more regular object. More broadly speaking, these results can be enhanced by improving the quality of the training set.

Nevertheless, potential errors in the reconstructed image, especially the invented structures are problematic for end users (clinicians, biologists...). The MC dropout approach proposed in this article helps locating most of them. Importantly, the estimated uncertainty remains low where reconstruction is correct, leading to a clear distinction of the suspicious areas: false alarms, which could mislead the clinician, are rare. It is worth mentioning that if the method helps locating the invented structures or incorrect reconstruction, (true positive), it is less efficient when capturing missed structure (false negative), as illustrated in the supplementary materials, Fig. S4. This is understandable, as these errors are mostly related to a lack of information in the data. The map of the standard deviation over experimental realizations was computed to compare our result with the uncertainty map based on making the CNN stochastic. The areas showing high values in the two maps are co-located, showing that the MC dropout method with a single-shot acquisition provides an information similar to that resulting from noise induced variability. This feature is promising in the context of moving tissue imaging and real time navigation. We explored other methods that can provide uncertainty maps, such as Deep ensemble [START_REF] Lakshminarayanan | Simple and scalable predictive uncertainty estimation using deep ensembles[END_REF] and Dropout ensemble [START_REF] Xue | Reliable deep-learning-based phase imaging with uncertainty quantification[END_REF]. In our study, the estimation of the mean image turned out to be superior for MC dropout according to the metrics presented previously. This difference could be explained by the required modification of the loss function for the other methods, involving a decrease of the overall performance.

To obtain the result of our study, our model was trained on an experimental dataset. However in clinical context, large experimental datasets may be complex to build. Using only simulated data to train a model and reconstruct experimental images would be ideal as simulation data can easily by produced. In our study, this approach turned out to produce unsatisfactory results, as we illustrate in supplementary materials (Fig. S3): predictions on experimental data provided by a DLA trained on either simulated data or experimental one are shown, and show that many artefacts remain on the experimental images reconstructed from a DLA trained with simulation results. These results are also in agreement with observations made by Davoudi et al. [START_REF] Davoudi | Deep learning optoacoustic tomography with sparse data[END_REF]. Although simulations are supposed to correctly model experiments, the nature of the ground truth for both training is very different for simulations and experiments in our case, which is likely the reason why training from experimental data provides much better results as compared to training with simulation data. Nonetheless, as shown in the previous section, the incorporation of simulated data through a transfer learning approach allows reducing significantly the size of the experimental dataset. The algorithm only needs to update its parameters with the difference between the simulations and the experiments, which is easier than learning the overall procedure. In the medical field, such a pretraining session could be useful for reducing the number of patients necessary to create a training set.

While our objective was limited to a proof-of-concept demonstration, several challenges must be taken into account in order to apply our method in a practical context.

The ground truth used for our proof of concept has to be replaced in practice by a ground truth that can be measured in a realistic environment (including deep inside tissue). Such a ground truth measurement could come from another imaging modality such as X-ray computer tomography (X-ray CT) or magnetic resonance imaging (MRI), which can accurately retrieve morphological information. They would however not provide a quantitative ground truth for the optical absorption, which would prevent the use of the method for spectroscopic approaches. A quantitative estimate of the ground truth can be obtained with a more sophisticated PA imaging device [START_REF] Davoudi | Deep learning optoacoustic tomography with sparse data[END_REF], or any method providing a visibility artefact free image proportional to the optical absorption. The training could be done with such a device free of visibility artefacts, to train a DLA to be applied on a simpler device. Quantitative reconstruction has been obtained with such an approach in the context of simulations [START_REF] Bench | Towards accurate quantitative photoacoustic imaging: learning vascular blood oxygen saturation in 3d[END_REF] The influence of noise on RF signals should also be studied to assess the validity of our approach in a noisier environment. In our work, the signal to noise ratio (SNR) on RF signals is about 60. This value however represents the SNR of signals produced by horizontal structures, while this work mainly focuses on the reconstruction of vessels affected by the visibility problem, for which the signal is almost nonexistent. In addition, the background of our DAS images is polluted by clutter, an artefact located around the object originating from the lack of information for the reconstruction. The amplitude of the clutter is often higher than the one from vertical structures. The model situation considered in this work thus remains significantly challenging.

Finally, the quality of the prediction is strongly influenced by the class of the object to reconstruct. The relative homogeneity of the studied dataset is one of the reason the DLA performs well: while very good predictions were made for leaves that were never seen by the network during the training, the unseen leaves were for the same species than the leaves used for the training. This is however a quite inherent limitation of deep learning approaches. In the supplementary information, we provide a preliminary investigation on the ability of our approach to generalize prediction to objects that do not belong to the class used for training: we tried to reconstruct the vessel structure available in the k-wave package (http://www.k-wave.org/) with the algorithm trained on the simulated dataset constructed with leaves. The predicted image, showed in supplementary materials (Fig. S5), is well reconstructed without artefact, suggesting there is no, or low overfitting in our approach, and that the network may generalize well. In a more general context, the capacity of a network to generalize is crucial and must be investigated for each particular situation.

Aside from increasing the quality of the reconstructed image, DLA offers other interesting features. Only 10 ms is needed to make a prediction using a regular graphic card which is much lower than the reconstruction time for the deconvolution method. Real-time reconstruction during user navigation could be achieved. Besides, once trained, the network does not need any parameters to be set by the user, unlike for deconvolution approach where the regularization parameter has to be chosen carefully and in a rather subjective way.

Conclusion

The possibility of removing visibility artefacts with a neural network has been demonstrated both in simulations and experiments on a model class of complex objects. Vertical parts of objects and the inside of large structures, missing on conventional reconstruction approaches, are recovered. These qualitative assessments are confirmed by quantitative metrics, which are far better for the DLA approach compared to conventional reconstruction methods. However, some errors may still be present in the reconstructed images, such as invented or poorly reconstructed structures as well as missing structures, although their number might be reducible by improving the experimental protocol. A MC dropout approach was proposed and successfully applied to identify invented and poorly reconstructed structures: high values on the generated uncertainty map are in agreement with high values on the true error map. Besides, it has been shown that pretraining the network with simulated data enables to reduce the size of the experimental training set by a factor of 4 while maintaining a similar quality of reconstruction.

Network architecture

Unet is a well known network architecture first developed for segmentation task. Our implementation is shown in Fig. S1. It is a convolutional neural network composed of two paths: the contracting and expanding path. The first one, called the encoder, is a traditional stack of convolutional and pooling layers where the network extracts more and more complex features. The second one, called the decoder, is the symmetric expanding path where pooling operations are replaced by upsampling operators to recover at the output the resolution of the input. Context information is propagated from the encoder to the decoder through skip connections to provide local information to the global information while upsampling (black arrows). The weights are initialized with samples from a truncated normal distribution centered on 0 with standard deviation depending of the number of units in the weight tensor. Dropout layers are added to this architecture. Dropout is a popular regularization technique to limit overfitting. A certain set of neurons, chosen randomly, are disabled at each training step. This prevents units from co-adapting too much and forces the network to learn more robust features. Batch normalization normalizes the output of the previous activation layer by subtracting the batch mean and dividing by the batch standard deviation. It helps to speed up the learning and also reduces overfitting by adding some noise, similarly as dropout.

mBF or dmBF image as input of the network

The input of the network is obtained from the delay and sum algorithm (DAS) applied to time signals. When applied to real time signals, DAS provides a mBF image. When applied to complex signals obtained with a Hilbert transform, DAS provides a complex image whose modulus is the dmBF image. The mBF image and the dmBF image are the two types of input that we consider here. One DLA was trained for each type of input, the corresponding predictions are displayed in Fig. S2. The algorithm performs better on the mBF image, leading to a scaled and shifted structured similarity index (sSSIM) of 0.76 instead of 0.72. The prediction from the dmBF image suffers from more artefacts (arrows) and the DLA fails to recover the true vessels thicknesses, which are over estimated. The mBF image, despite being more different from the true physical structure of the object than the dmBF image, and thus from the ground truth, carries more information to be captured by the network.

Reconstruction of experimental data with DLA trained on simulation data

Predictions from neural networks trained with simulation datasets or experimental datasets are presented in Fig. S3. Although the DLA trained on simulation data still manages to find several vertical structures that are not visible on the DAS image, the predicted image is polluted by a lot of artefacts. Clearly, experimental data are necessary to train efficiently the model. However, as shown in the main text, pretraining the network on a simulation dataset allows reducing the size of the experimental training set. object is clearly different, the algorithm turns out to perform extremely well on this specific example, as visibility artefacts are removed, suggesting no or low overfitting. However, this is just one example, obtained from simulated data, and further studies are required to investigate the generalization limits of the network, especially for experimental data.

Figure 1 :

 1 Figure 1: a, Creation of the experimental training set. A linear probe is coupled to a water tank containing the leaf, through a window composed of a tight Mylar membrane. The leaf is in the imaging plane of the probe. The laser beam is shined from the top and the RF signals are acquired. A photograph of the leaf was previously taken. The mBF PA image of the ROI is reconstructed and the photograph is processed to extract the same area. b, Uncertainty prediction: Several images are generated using the same input. The mean and the standard deviation (std) of these samples are estimated pixel by pixel. The prediction is unstable in the marked area, resulting in a high std.

Figure 2 :

 2 Photograph ( Ground truth)

Figure 3 :

 3 Figure 3: Results on the experimental test set built from experimentally acquired RF signals, two examples. a, e, Ground truth: photograph of the object. b, f, dmBF image, delay and sum. c, g, L2-regularized deconvolution. d, h, Prediction of the deep learning algorithm.

Figure 4 :

 4 Figure 4: Uncertainty estimation on predicted experimental image. a, Ground truth: photograph of the object. b, e, Mean image and standard deviation of the object computed over 20 inferences generated from a unique acquisition with the deep learning algorithm, dropout activated. c, f, Mean image and standard deviation of the object computed over prediction generated from several RF signals acquisition (experimental variability) with the deep learning algorithm, dropout disabled. d, Absolute error between the ground truth (a) and the mean (b).

  pretrained with simulated data STD over the trainingset, pretraining Network not pretrained STD over the trainingset, no pretraining Ground truth sSSIM Size of the experimental trainingset

Figure 5 :

 5 Figure 5: Performance of a deep learning algorithm trained on an experimental dataset only (red) versus a deep learning algorithm pretrained on 1400 simulated data then trained on the same experimental dataset (blue), for different number of experimental data.

Figure S1 :

 S1 Figure S1: Network architecture.

Figure S2 :

 S2 Figure S2: Deep learning algorithm (DLA) prediction on experimental data for dmBF PA image or mBF PA image as input, two examples. a, f, dmBF image. b, g, mBF image. c, h, Prediction with DLA trained on dmBF data. d, i, Prediction on experimental data with DLA trained on mBF image. e, j, Ground truth.

Figure S4 :

 S4 Figure S4: Uncertainty estimation on experimental mBF image. a, g, Ground truth: photograph of the object. b, h, Mean image of the object computed over 20 inferences generated from a unique acquisition with the deep learning algorithm, dropout activated. e, k, Corresponding STD. c, i, Mean image of the object computed over prediction generated from several RF signals acquired at different time with the deep learning algorithm, dropout disabled. f, l, Corresponding STD d, j, Absolute error between the ground truth (a,g) and the mean (b,h).

  is a widely used CNN in the medical field. A slightly modified architecture, presented in supplementary materials is implemented with the open source libraries

							Training set
						Image	Ground truth
		Photograph			processing
				RF data
						DAS	DLA input
						reconstruction
	Test set Unet input	Trained CNN N prediction	1	2	…	N	Mean	Std
		dropout
		activated

Table 1 :

 1 Quantitative measurement of reconstruction quality with the normalized 2D cross-correlation and the scaled and shifted structured similarity index .02 0.89 ± .01 0.44 ± .06 0.80 ± .03 sSSIM 0.29 ± .02 0.87 ± .01 0.38 ± .05 0.76 ± .03

		Simulation		Experiment
		DAS	DLA	DAS	DLA
	NCC	0.31 ±

6. Funding Information This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 681514-COHERENCE).

DLA prediction, experimental test data

Uncertainty estimation

Uncertainty estimation of the two previous examples are presented in Fig. S4. Similarly to the example in the main text, the standard deviation (std) map helps to locate errors in the reconstruction such as invented structures and incorrect orientation or position (arrows). One can notice that although the value of the std is higher at location of some missing veins (star), some of them are not even displayed on the std (circle). The lack of information in the data may be more important for these structures, misleading the DLA.

Generalization ability

This work was mostly dedicated to demonstrate the capacity of a neural network to correct visibility artefacts on a specific class of object. However the possibility or reconstructing objects from another class is an open question. To illustrate this problematic on a specific example, we used the network trained on leafs (simulated signals) to reconstruct the PA image of a vessel-like structure (test image from the k-wave package, http://www.kwave.org/). The result is presented in Fig. S5


Océane C Salles

Glenn R Almany

Michael L Berumen

Geoffrey P Jones

Pablo Saenz-Agudelo

Maya Srinivasan

Simon R Thorrold

Benoit Pujol
email: benoit.pujol@univ-perp.fr

Serge Planes

Saenz-Agudelo

Strong habitat and weak genetic effects shape the lifetime reproductive success in a wild clownfish population

Keywords: multi-generational pedigree, additive genetic variation, maternal effects, environmental 32 effects, adaptation, selection, heritability, evolvability : 2. Number of Tables: 1. Number of boxes: 0 40

à la diffusion de documents scientifiques 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.

 [START_REF] Clutton-Brock | Reproductive success[END_REF]. Several factors including different environmental, parental, 61 and additive genetic effects can influence LRS [START_REF] Hendry | The contemporary evolution of fitness. 524[END_REF]). If LRS is exclusively a 62 phenotypic response to the conditions experienced by individuals, populations will be severely 63 impacted by rapid environmental change and there is no prospect of adaptive microevolution.

64 Maternally transmitted responses to environmental conditions can mediate the effect of those changes 65 on the next generation [START_REF] Mousseau | The adaptative significance of maternal effects[END_REF]. However, it is the additive genetic variation in LRS 66 that determines the rate of adaptation of a population to the environmental demand [START_REF] Fisher | The genetical theory of natural selection[END_REF]). To 67 date, few studies have distinguished these relative contributions to LRS over multiple generations in 68 wild populations [START_REF] Hendry | The contemporary evolution of fitness. 524[END_REF]. This situation is changing as long-term, individual-based 69 ecological studies with multi-generational pedigrees provide the necessary longitudinal information 70 to quantify the different components of LRS [START_REF] Pemberton | Wild pedigrees: the way forward[END_REF]. Such studies are imperative as we 71 seek to understand the ability of species to withstand or adapt to accelerating climate change (Munday 72 et al. 2017).

73

74 Intergenerational responses to selection are a product of the interplay between evolutionary and 75 ecological mechanisms that ultimately shape inherited variation in fitness-related traits.

76 Environmentally driven mechanisms (e.g., phenotypic plasticity, genetic assimilation) can facilitate 77 [START_REF] Ghalambor | Adaptive versus non-adaptive phenotypic 514 plasticity and the potential for contemporary adaptation in new environments[END_REF][START_REF] Danchin | Epigenetically facilitated 493 mutational assimilation: epigenetics as a hub within the Inclusive Evolutionary Synthesis[END_REF] or constrain the microevolutionary response to selection 78 [START_REF] Pujol | The missing 579 response to selection in the wild[END_REF]. However, in the absence of genetic variation for LRS, these mechanisms will 79 likely have little effect on a negligible rate of adaptive evolution. In quantifying additive genetic 80 variation, it is important to distinguish between heritability and evolvability (Wheelwright et al. 81 2014). Narrow sense heritability is widely used as a measurement of the population potential to 82 respond to selection [START_REF] Mousseau | Natural selection and the heritability of fitness components[END_REF]. It is the additive genetic variance standardized by the S1). They have largely confirmed low additive genetic variation and 97 evolvability (~0.08) for LRS, which nevertheless reflects some evolutionary potential (Hendry et al. 98 2018). Whether marine fish and terrestrial species reveal the same patterns remains unknown to date.

99 Until recently, quantifying LRS in marine organisms with a pelagic larval stage has been considered 100 impossible because of the difficulties in following the fate of offspring from one generation to the 101 next. However, there is increasing evidence of some degree of natal philopatry or self-recruitment in 102 local marine populations [START_REF] Jones | 531 Larval retention and connectivity among populations of corals and reef fishes: history, advances 532 and challenges[END_REF]. Genetic parentage analysis is making it possible to assign 103 a significant proportion of offspring to their parents (Planes et al. 2009) and construct 104 multigenerational pedigrees for the offspring that return to their natal population [START_REF] Salles | First 590 genealogy for a wild marine fish population reveals multi-generational philopatry[END_REF].

105 As in most quantitative genetic studies of LRS, it is impossible to measure the recruitment of marine 106 fish dispersing juveniles at other locations. The regional component of LRS, which would inform us 107 on fitness variation beyond the local scale is impossible to obtain. Measuring the local component of 111 For coral reef fishes, quantifying environmental and genetic components of LRS and assessing its 112 evolvability in wild populations is of great contemporary importance. Between 30 to 50% of the 113 world's coral reefs have been lost and those remaining are considered highly vulnerable (van 114 Hooidonk et al. 2016). The rapid loss of suitable habitat is widely acknowledged to be contributing 115 to a decline in reef fish populations and biodiversity [START_REF] Jones | Coral decline threatens fish 534 biodiversity in marine reserves[END_REF]). The potential for adaptation 116 of reef fish is uncertain and near future environmental conditions predicted under climate change 117 should have a dramatic effect on their reproductive success [START_REF] Munday | Predicting 565 evolutionary responses to climate change in the sea[END_REF]. To date, 118 environmental, maternal and additive genetic contributions to LRS in wild coral reef fish populations 119 have not been assessed. However, recent work establishing high levels of natal philopatry in some 120 coral reef fishes [START_REF] Jones | Coral reef fish larvae settle close to home[END_REF]) that allowed for resolving their pedigree across multiple 121 generations [START_REF] Salles | First 590 genealogy for a wild marine fish population reveals multi-generational philopatry[END_REF] opened the way for quantifying the local component of LRS for the 122 first time in coral reef fishes.

123 124 Here, we focus on the entire local population of the orange clownfish Amphiprion percula at Kimbe 125 Island, Papua New Guinea where each year ~half the juveniles successfully recruiting are progeny of 126 local breeding pairs [START_REF] Salles | First 590 genealogy for a wild marine fish population reveals multi-generational philopatry[END_REF]. We use multi-generational pedigrees of up to 5 generations 127 obtained from biennial DNA sampling over 10 years [START_REF] Salles | First 590 genealogy for a wild marine fish population reveals multi-generational philopatry[END_REF]) and apply a quantitative 128 genetic linear mixed model approach [START_REF] Kruuk | Introduction. Evolutionary dynamics of wild populations: the use 543 of long-term pedigree data[END_REF] to quantify the additive genetic, maternal 129 and environmental components of variation in LRS for the self-recruiting portion of the population. 142 with one of two host sea anemone species, Heteractis magnifica (Fig. 1b) and Stichodactyla gigantea 143 (Fig. 1c). We geographically located and tagged a total of 310 anemones (176 H. magnifica and 134 144 S. gigantea) that were occupied by A. percula on the entire reef surrounding the island. 145 146 These two anemone species are remarkably different in terms of the micro-habitat they provide, 147 including a wide range of shapes, sizes, depth distributions and surrounding substrata [START_REF] Dunn | The clownfish sea-anemones -Stichodactylidae (coelenterata, actiniaria) and 499 other sea-anemones symbiotic with pomacentrid fishes[END_REF].

148 Although we did not directly measure these variables (other than depth), the combination of host 149 anemone species (H. magnifica or S. gigantea) with one of the three geographical areas covering the 150 entire reef around the island (northern, western or eastern areas, Figure 1a) where it is located 151 describes a 'habitat' variable for each fish that encompasses a suite of biotic and abiotic 152 environmental conditions. These geographical areas correspond to the different coasts of the island 153 that reflect possible environmental effects of the geographic location. In total, the combination of the 154 two host anemone species with the three geographical areas allowed us to describe six different 155 habitats. [START_REF] Fricke | Mating system, resource defense and sex change in the anemonefish 508 Amphiprion akallopisos[END_REF]). The female is the largest, the male is the second

F o r R e v i e w O n l y
8 160 largest, and the non-breeders rank progressively lower in the hierarchy as they decrease in size. If the 161 single female adult of a group dies, then the male changes sex to female, and the largest non-breeder 162 from the anemone becomes sexually mature as male. Reproduction occurs year round, with females 163 laying several hundred eggs in a clutch near the pedal disk of the host anemone each lunar month.

164 The eggs hatch after ~7 days of paternal care into larvae that spend ~10 days in the pelagic 165 environment [START_REF] Berumen | Otolith 481 geochemistry does not reflect dispersal history of clownfish larvae[END_REF] before settling settling to an anemone, either at their natal location 166 (Kimbe Island) or elsewhere (Planes et al. 2009).

167 168 Divers captured fish by using hand nets. Individuals were measured in situ using calipers, fin-clipped 169 (size > 35mm) or collected whole (size < 35mm) for genetic analysis and then released back on the 170 same anemone. Small pieces of fin tissue were preserved in 95% ethanol in 2-mL vials. The biggest 171 fish in each anemone was identified as the female, the second largest individual was assumed to be 172 the male, and all other individuals were classified as non-breeders. We extracted DNA from all 173 samples and genotyped them at 22 polymorphic microsatellite loci [START_REF] Bonin | 483 Characterisation and cross-amplification of microsatellite markers in four species of 484 anemonefish (Pomacentridae, Amphiprion spp.)[END_REF]. Then, we 174 identified the individuals sampled multiple times over the years by using the Excel macro GenAlex 175 v6.5 [START_REF] Peakall | GenAlex 6.5: Genetic analysis in Excel. Population genetic 569 software for teaching and research-an update[END_REF] to compare multilocus genotypes from 2003,2005,2007,2009,2011 176 and 2013. Individuals were in average sampled 2.88±0.04 times (mean±SE) over the six surveys (1% 177 of individuals persisted over the 10-yr period, [START_REF] Salles | First 590 genealogy for a wild marine fish population reveals multi-generational philopatry[END_REF]. The 2-yr sampling scheme 178 precluded calculating a precise measurement of the age of individuals, in particular for fish sampled 179 in 2003 during the first sampling period, which age was unknown. The total duration of this long term 180 survey did not allow us to obtain many replicated measurements within individuals before and after 181 sex change (n=41 individuals). Estimating sex dependent additive genetic variance is precluded in 182 this case because some effects cannot be disentangled as the clownfish only changes sex in one 183 direction (from male to female). This change is always associated with a change of sexual partner and 184 with an increase in female body size, which we expect to generate a confounding effect between a 188 Pedigree used for quantitative genetic analysis 189 Pedigree reconstruction was conducted by assigning juvenile fish to parental pairs on the basis of 190 their multilocus genotypes [START_REF] Salles | First 590 genealogy for a wild marine fish population reveals multi-generational philopatry[END_REF]. We used the software FaMoz [START_REF] Gerber | FaMoz: A software for parentage analysis using 510 dominant, codominant and uniparentally inherited markers[END_REF].

191 This approach is based on the calculation of LOD scores (Log of the odd-ratio comparison) for any 192 potential parentage relationship. It determined critical thresholds to accept or reject assignments by 193 simulating true and false parent-offspring pairs. Further details on parentage analyses and pedigree 194 reconstruction are given in Salles et al. (2016a). We kept assignments to known parental pairs, but 195 rejected assignments to single adults. In the context of overlapping generations, we used the year of 196 first sampling and the anemone of each parental couple as information to avoid possible false 197 assignments. As a result, sibship links could not be confused with parental links. The same individual 198 can be related to its offspring with either a paternal or maternal link because of sex changes. We 199 identified the mother and the father based on the size of the two parents on the year of first capture. [START_REF] Salles | First 590 genealogy for a wild marine fish population reveals multi-generational philopatry[END_REF]. For this study, we excluded from the original 204 pedigree the 1192 individuals that were removed from the habitat at the juvenile stage (size < 35mm, 205 10 to 458 days old). The final pedigree used for this study includes 1735 individuals from five 206 generations (Fig. S1). We used the R package 'pedantics' [START_REF] Morrissey | Pedantics: an r package for pedigree-based genetic 556 simulation and pedigree manipulation, characterization and viewing[END_REF] to assess the 207 power of the resolved pedigree to detect significant quantitative genetic parameters (Fig. S1).

209 LRS: the individual contribution to self-recruitment

F o r R e v i e w O n l y F o r R e v i e w O n l y
11 236 pedigree information to extract the additive genetic component. This approach is more powerful than 237 traditional analyses (e.g., parent-offspring regressions) because it takes into account every 238 relationship link in a pedigree. Maternal variance was modeled using the mother's identity as a 239 random effect, allowing maternal effects to include both genetic and environmental maternal effects.

240 Permanent environmental effects were modeled by including the identity of individuals as a random 241 effect. The LRS variance is the sum of six variance components: 271 transformed back on the data-scale are expected to be upward biased, their ratio is reliable, and hence 272 the estimators derived from their relative proportions such as h². It is necessary to point out two 273 specific aspects of this back transformation on the observed data scale. First, V R is estimated on the 274 basis of the additive over-dispersion term in the nonlinear model and its value cannot be interpreted 275 similarly to the usual residual variance term estimated by classical quantitative genetic linear mixed 276 models. Second, the sum of the variance components estimated on the data scale are not additive and 277 therefore not expected to sum up to the value of the phenotypic variance calculated directly on the 278 raw data. For the sake of clarity and comparison, we present the results on the latent scale and the 279 observed data scale. We calculated the 95% credibility intervals from the posterior distributions of 280 observed parameters for all the variance components and other estimates expressed on their basis by 281 using the 'HDInterval' package [START_REF] Meredith | HDInterval: highest (posterior) density intervals. R package 551 version 0[END_REF]. 14 314 average number of juveniles originating and recruiting in the population per generation. The 315 heritability estimate expressed on the observed data-scale was h²=0.013 (CI 95% 4.951×10 -5 to 316 1.227×10 -2 ) for LRS (Table 1). We can therefore estimate the maximum response (R) to selection (S), 317 expected in theory, in the presence of strong selection pressures acting on the Kimbe Island orange 318 clownfish population by using the Breeder's equation R=h²×S [START_REF] Lush | Very low levels of direct additive genetic variance in fitness and fitness components in 549[END_REF]. The low to negligible 319 value of the LRS heritability means that the maximum genetic change of the population average LRS 320 would never exceed ~0.020 offspring per generation. Similar results were obtained for DL (see 321 Supplementary Information for more details).

242 = + + + + + (

323 Weak maternal effects for LRS

324 Our analysis detected maternal variance for the LRS but it was also extremely small to the extent that 325 it might be considered as null (Table 1): V M =0.019 (CI 95% 2.966×10 -5 to 2.044×10 -2 ). It made very 326 little contribution to the total variance in LRS (m²=1.9%, expressed on observed data-scale, Fig. 2).

327 The habitat occupied by the mother (Natal Habitat) had a stronger effect on LRS than the mother 328 herself. The relative contribution of individuals to the population replenishment was indeed 329 influenced by the Natal Habitat to an extent of 19.1% for LRS. Similar results were obtained for DL 330 (see Supplementary Information for more details).

332 Discussion

Strong habitat and weak genetic effects on LRS

334 Our study revealed that biennial LRS in the Kimbe Island orange clownfish population quantified 335 over five generations was largely explained by ecological factors. Host anemone species, depth and 336 geographical location explained ~25%, and permanent environmental effects ~30% of its variation, 337 with only weak maternal (1.9%) and additive genetic effects (1.3%). Intrinsic biological 338 characteristics of anemone species (e.g., size, shape and toxicity) affecting the life-history traits of 339 their resident clownfish might explain strong habitat effects [START_REF] Chausson | Host anemone size as a determinant of social group 486 size and structure in the orange clownfish (Amphiprion percula)[END_REF]. In addition, the 355 Detrimental direct and indirect anthropogenic impacts on reef anemones habitats are already affecting 356 numerous clownfish species [START_REF] Saenz-Agudelo | Detrimental effects of host 588 anemone bleaching on anemonefish populations[END_REF]. S. gigantea anemones located in shallow 357 waters are likely to be disproportionately more impacted by increasing water temperatures and 358 irradiance [START_REF] Hobbs | 526 Taxonomic, spatial and temporal patterns of bleaching in anemones inhabited by 527 anemonefishes[END_REF]. If these habitats are differentially impacted, this will affect clownfish 359 contributions to the local replenishment of the population and compromise population persistence.

360

361 Low to negligible evolutionary potential 362 Our findings provide the first empirical support for a wild marine population to Fisher's fundamental 363 theorem of selection that additive genetic variance in fitness is depleted under selection and tends 364 towards zero in a population reaching evolutionary equilibrium [START_REF] Gustafsson | Lifetime reproductive success and heritability: Empirical support for Fisher's 521[END_REF]. Although normal 365 and expected, low heritability and evolvability in LRS is concerning giving the increasing rate of

F o r R e v i e w O n l y F o r R e v i e w O n l y
17 392 effects explained 1.9% of the LRS variance, which is quite small, even if it was more than additive 393 genetic effects. Maternal effects therefore increased slightly the low to negligible rate of LRS change 394 by adaptive evolution. One must note that this increase is nearly negligible. It is likely that this small 395 value reflects the genetic component of the maternal effect. The habitat of birth, which is also the 396 maternal habitat, probably encompassed a non-negligible part of the maternal environmental effects 397 [START_REF] Germain | Hidden responses to environmental variation: maternal effects 512 reveal species niche dimensions[END_REF]. There is growing awareness that maternal environmental effects can 398 contribute to adaptation in natural populations, especially when maternal and offspring environments 399 are positively correlated [START_REF] Shama | Bet hedging in a warming ocean: predictability of maternal environment shapes 593 offspring size variation in marine sticklebacks[END_REF]. It might even buy some time for adaptive evolution through 400 slow genetic change to occur [START_REF] Levis | Evaluating'Plasticity-First' Evolution in Nature: Key Criteria 545 and Empirical Approaches[END_REF].

402 Towards a wider sample of contemporary rates of adaptive evolution in the wild

403 In our study, LRS estimated local fish contribution to the local population replenishment. As with 404 most wild population pedigrees, this excludes the dispersal fitness because the amount of offspring 405 that successfully dispersed somewhere else is unknown [START_REF] Kruuk | Heritability of fitness in a wild mammal population[END_REF]. Its genetic variation 406 evaluates the rate of adaptive evolution inside the Kimbe island population. While there are no 407 comparable data from marine systems, 15 studies conducted on terrestrial vertebrates have also 408 estimated the additive genetic variation and the heritability of LRS (Table S1). The number of 419 Although useful insights were gained, one must note that the use of LRS as a surrogate for fitness is 420 not well suited for iteroparous species [START_REF] Grafen | On the uses of lifetime reproductive success[END_REF]. This is particularly true when population growth 421 varies between reproductive seasons [START_REF] Coulson | 490 Estimating individual contributions to population growth: evolutionary fitness in ecological 491 time[END_REF]. Furthermore, although quantitative genetics generalized linear mixed model analysis (for details see Table 1).

58

  Lifetime reproductive success (LRS) -the number of successful offspring an individual contributes 59 to the next generation -is a critical variable underpinning ecological and evolutionary responses to 60 the environment

  variance. Low heritability values can either reflect low additive genetic variance or 84 large environmental, or residual effects. Evolvability is the mean-standardized additive genetic 85 variance[START_REF] Houle | Comparing evolvability and variability of quantitative traits[END_REF]. Environmental or maternal effects do not affect it, which makes it a more 86 appropriate metric in the comparison of evolutionary potential between traits, populations and 87 species. When measured for a surrogate of fitness like LRS, it estimates the expected proportional 88 change per generation in population mean fitness given a unit strength of selection. Heritability for 89 LRS reveals whether its additive genetic variance represents a non-trivial proportion of its total 90 variance in the actual environmental context of a given wild population. Together, evolvability and 91 heritability inform us about how much environmental change a wild population can withstand on the 92 basis of its evolutionary potential.93 94The few (n=15) long-term, individual-based studies that have quantified additive genetic variation, 95 heritability and evolvability of LRS in wild populations have all focused on terrestrial species(Postma 96 2014;[START_REF] Hendry | The contemporary evolution of fitness. 524[END_REF] Table

  marine fish is nevertheless an opportunity as in any other species because it estimates the 109 relative contribution of local fish to the population self-recruitment and replenishment.

  110

130

  Habitat effects were quantified by examining this local LRS for individuals resident in two different 131 anemone species and from different geographic locations around the island. By integrating habitat 132 data with the pedigree information in a quantitative genetic generalized linear mixed model, we were 133 able to assess the relative contribution of additive genetic, maternal and habitat effects to local LRS.

134

  We also calculated the evolvability and heritability of local LRS to evaluate its evolutionary potential 135 to respond to selection at the scale of the Kimbe island population. 136 137 Material and methods 138 Study population and data collection 139 A natural population of orange clownfish (A. percula) living in the reef surrounding Kimbe Island 140 (Fig. 1a; 5°12'22.56'' S, 150°22'35.58'' E), West New Britain Province of Papua New Guinea, was 141 surveyed every second year from 2003 to 2013. Here, A. percula lives in a mutualistic association

156157

  Clownfish within one anemone live in group of typically three to five individuals in the Kimbe Island 158 population. The size-based dominance hierarchy in Amphiprion allows us to determine the 159 reproductive status of each individual

  and its genetic quality. We therefore did not consider sex in our model as a result of 186 data and analytical limitations.

  187

200

  The original population pedigree includes 2927 clownfish over five generations including 121 201 families, 987 paternal, 987 maternal, 1809 full-sib, 412 maternal half-sibs, 248 paternal half-sib, 135 202 maternal grandmothers, 135 maternal grandfathers, 278 paternal grandmothers, 278 paternal 203 grandfathers and 218 cousins

284

  Habitats dominantly shape the Lifetime Reproductive Success in the clownfish population 285 Biennial estimates of the Lifetime Reproductive Success (LRS) measured inside the area of the 286 Kimbe island population ranged from 0 to 13, with a phenotypic variance =1.31 and an average 𝑉 𝐿𝑅𝑆 287 value of 0.54±0.05 (mean±SE) offspring per individual for a two-year period. Because clownfish live

  of S. gigantea[START_REF] Nedosyko | Searching for a toxic key to 567 unlock the mystery of anemonefish and anemone symbiosis[END_REF]) might provide better protection against predators 341 of eggs attached to the host anemone, but this hypothesis remains to be tested. The geographical 342 location also appears to be important. Most successful individuals were close to the land in shallow 343 water on S. gigantea and in deeper lagoons for H. magnifica, which might promote greater local 344 retention of larvae. Mechanisms responsible for geographical differences in LRS around Kimbe 345 Island remain unknown (Berumen et al. 2010). Some breeders likely have a different reproductive 346 success beyond the sampled population, through dispersers rather than self-recruiters. Inside the 347 Kimbe island clownfish population, weak genetic effects on LRS indicate a low to negligible rate of 348 adaptive evolution in progress and raise concern about the ability of this population to adapt to rapid 349 climate change. 350 351 Susceptibility to habitat change 352 Habitat largely drove LRS variation. Individuals that happen to settle on particular anemones and 353 places do well. The dependence on habitat quality of LRS inside the Kimbe island population suggests 354 this species will be extremely susceptible to habitat degradation over ecological time scales.

  409 estimates of maternal effects on LRS variation are extremely rare(McFarlane et al. 2014 and 410 references therein). A majority of these studies similarly found low to negligible contributions of 411 additive genetic effects in the wild. Marine fish, mammals and birds revealed the same patterns. There 412 is little evidence of additive genetic variance for LRS in studies conducted in the wild. This situation 413 is less clear for maternal effects, partly because studies remain scarce. Very low but significant 414 additive and maternal genetic variations for fitness indicate that there was some genetic change over 415 the course of the survey but very limited potential for short-term adaptive evolution. 416 417 A wider use of LRS does not preclude caution in wild populations with a known pedigree have multiplied during the last thirty years.

422

  genetic approaches decomposing the genetic and environmental components of LRS are theoretically 423 correct, they might provide imprecise parameter estimates in the wild because temporal and spatial 424 environmental effects (e.g., condition dependent LRS, cohort effects) might not have been adequately 425 corrected for. Caution must therefore be taken when interpreting the additive genetic variance of LRS 426 as adaptive evolutionary rate or evolutionary potential in wild populations where it might in fact 427 reflect other ecological parameters. Our quantitative genetic model explained a large amount of 428 environmental variation in LRS in the Kimbe island clownfish population, which supports the 429 reliability of our genetic estimates. 430 431 Conclusion 432 The major outcome of this study is that the heterogeneity of the habitat of the Kimbe Island orange 433 clownfish had a profound influence on the individual contribution to the local population 434 replenishment over five generations. This finding implies that habitat ecology is crucial for this 435 clownfish population. In terms of future persistence, expected changes in habitat quality and 436 configuration over relatively short time scales might affect the ability of fish to self-recruit. This 437 ability harbored low to negligible additive genetic and maternal genetic variation. As a consequence, 438 this population potential for rapid evolutionary change of LRS by selection, and therefore its rate of 439 adaptive evolution, can be considered negligible in the current state of the population. This finding, 440 which is in line with other studies on the topic, stresses the importance of environmental mechanisms 441 (e.g., plasticity) that have the potential to enable rapid adaptive responses. From the perspective of 442 management, our results caution against hoping for local adaptive evolutionary responses and lend 443 support to focusing conservation efforts on maintaining habitat quality.

F

  of variation in the biennial estimate of the local Lifetime Reproductive 628 Success (LRS) of the Kimbe Island orange clownfish. 629 Distribution of the biennial estimate of the LRS (histograms). Variance components on both latent-630 scale and observed data-scale (pie charts) for the biennial estimate of the LRS explained by 631 Additive genetic (V A ), Maternal (V M ), Natal Habitat (V NH ), Resident Habitat (V RH ), Permanent 632 Environmental (V PE ) and Residual (V R ) variances. These proportions were calculated from the 633 values of the posterior modes of a quantitative genetics generalized linear mixed model analysis (for 634 details see Table 1

Figure 1 .Figure 2 .

 12 Figure 1. Variation of the total number of offspring orange clownfish produced on each anemone around Kimbe Island between 2003 and 2013.(a) The studied region was divided in three areas (northern, eastern and western areas). Colors correspond to the variation of the total number of juveniles locally self-recruited that were produced on each anemone (varying from 0 to 27) over a 10-year period. The expected value is interpolated from those around it (using default algorithms implemented in Origin software). Dots correspond respectively to the location of the two host anemones species: (b) Heteractis magnifica (black dots) and (c) Stichodactyla gigantea (white dots).Photos by Tane Sinclair-Taylor.

  1)𝑉 𝐿𝑅𝑆 𝑉 𝐴 𝑉 𝑀 𝑉 𝑁𝐻 𝑉 𝑅𝐻 𝑉 𝑃𝐸 𝑉 𝑅 243 Quantitative genetic models were computed as univariate GLMMs using the 'MCMCglmm' package standard narrow sense heritability estimates for LRS (h²) by applying the basic formula 263 (h² = / ), and similarly maternal effects by estimating the proportion of total phenotypic variance 𝑉 𝐴 𝑉 𝑃 264 explained by the maternal variance (m² = / ). Evolvability (I A ) of LRS is the additive genetic 𝑉 𝑀 𝑉 𝑃 265 variance divided by the squared mean of the LRS[START_REF] Houle | Comparing evolvability and variability of quantitative traits[END_REF]. Analyses assumed a Poisson 266 distribution and provided parameter estimates on a statistically convenient latent scale for non-267 Gaussian traits. We therefore back-transformed these estimates onto the observed data scale to 268 improve our inferences by using the 'QGlmm' package(de Villemereuil et al. 2016). We used the

	262 we obtained 269 function 'QGparams' to estimate additive components such as	𝑉 𝐴	and h², and 'QGicc' to estimate
	F r o 270 broader sense components such as 𝑉 𝑀	and m²,	, 𝑉 𝑁𝐻 𝑉 𝑅𝐻 𝑉 𝑃𝐸 ,	and . Although parameter estimates 𝑉 𝑅
	R
	e v i e
		w
			O n l
			y
			12

Table S1 .

 S1 Additional Supporting Information may be found in the online version 3 of this preprint on the Zenodo 472 repository: https://doi.org/10.5281/zenodo.3476529. Previous estimates of fitness heritability and maternal effects on fitness in 15 wild Pedigree data and power analysis for the Kimbe Island orange clownfish population. Wheelwright, N.T., Keller, L.K.[START_REF] Postma | Four decades of estimating heritabilities in wild vertebrate populations: Improved 576 methods, more data, better estimates?[END_REF]. The effect of trait type and strength of selection 600 on heritability and evolvability in an island bird population. Evolution 68-11, 3325-3336.601 Wilson, A.J., Réale, D., Clements, M.N., Morrissey, M.M.,Postma, E., Walling, C.A. et al. (2010).

	470 Supplementary material 471 473 The Additional Supporting Information file contains: 495 de Villemereuil, P., Schielzeth, H., Nakagawa, S. & Morrissey, M. (2016). General methods for 496 evolutionary quantitative genetic inference from generalized mixed models. Genetics, 204, An ecologist's guide to the animal model. J. Anim. Ecol., 79, 13-26. (CI) (1.71 to 3.65) (0.65 to 8.05) (CI) (6.827×10 -4 to 0.057) 599 602 h² 0.019 0.013
	603 474 475 populations. 604
	476 Figure S1. 477 Supplementary methods, results and discussion on De-lifing measures.
	F F F F F F F F
	o r o r o r o r o r o r o r o r o r
	R R R R R R R R R
	e v i e e v i e e v i e e v i e e v i e e v i e e v i e e v i e e v i e
	w w w w w w w w w
	O n l O n l O n l O n l O n l O n l O n l O n l O n l
	y y y y y y y y y
	20 21 25 27

Ecology Letters

McFarlane, S.E., Gorrell, J.C., Coltman, D.W., Humphries, M.M., Boutin, S. & McAdam, A.G.

LRS