N Bianchi

L Balucani

P Podio

S Caselli

C Codella
email: codella@arcetri.astro.it

C Ceccarelli

E Bianchi

N Balucani

L Podio

P Caselli

S Feng

B Lefloch

A López-Sepulcre

R Neri

S Spezzano

M De Simone

Seeds of Life in Space (SOLIS) V. Methanol and acetaldehyde in the protostellar jet-driven shocks L1157-B0 and B1

Keywords: Stars, formation -ISM, jets and outflows -ISM, molecules -ISM, individual objects, L1157

Seeds of Life in Space (SOLIS). V. Methanol and acetaldehyde in the protostellar jet-driven shocks

L1157-B0 and B1.

Introduction

One of the open questions in astrochemistry refers to how chemistry richness develops and evolves in star forming regions. More specifically, the molecular complexity during the evolutionary stages that lead a molecular cloud to form a solar-type star and its planetary system has yet to be elucidated. Particularly relevant are the abundances of the so-called interstellar complex organic molecules (iCOMs: molecules with at least six atoms: Herbst & van Dishoeck 2009;Ceccarelli et al. 2017), as they may be considered as small blocks from which complex pre-biotic species can build up (e.g. Caselli & Ceccarelli 2012;Belloche et al. 2014). The present consensus is that after an initial cold phase in which icy mantles enriched with hydrogenated species form on the dust grains (e.g. Tielens & Hagen 1980;Watanabe & Kouchi 2002;Rimola et al. 2014), the iCOMs synthesis can either occur on the grain mantle surfaces initiated by energetic processing (e.g. Garrod & Herbst 2006;Garrod et al. 2008;Ruaud et al. 2015;Vasyunin et al. 2017) or in the gas phase initiated by the sublimation of the mantles themselves (e.g. Charnley et al. 1992;Balucani et al. 2015;Codella et al. 2017;Skouteris et al. 2018). In both cases, the composition of the grain mantles therefore plays a paramount role in the gaseous, observable abundance of iCOMs (see e.g. Burkhardt et al. 2019) -in the former possible scenario, because the mantle species are directly injected from the mantles into the gas; and in the latter possible scenario, because they feed gas reactions leading to iCOMs.

In this context, protostellar shocks are particularly precious because they can provide not only the observed iCOM abun-1 dances but also the approximate times at which they are formed. Indeed, at the passage of the shock, molecules present in the grain mantles are injected into the gas phase either because of sputtering (gas-grain collisions) or shattering (grain-grain) processes (e.g. Flower & Pineau Des Forêts 1994;Caselli et al. 1997;Schilke et al. 1997;Gusdorf et al. 2008). Once in the gas phase, the injected mantle species undergo reactions that destroy and/or form iCOMs (unless multiple shocks occur at the very same place, which is unlikely): therefore, the gaseous iCOM abundances are intrinsically time dependent. Fortunately, the chemical evolution time, which is a few thousand years at most, is of the order of the age of young protostellar shocks (e.g. Herbst & van Dishoeck 2009;Ceccarelli et al. 2017). Thus, sufficiently high spatial resolution observations towards these objects allow us to decipher the different time-dependent abundances of iCOMs and consequently provide extremely strong constraints on the formation and destruction of iCOMs. For example, this method was applied to constrain the formation route of formamide in the jet-driven protostellar shock L1157-B1 (Codella et al. 2017: hereinafter CC17). Previous observations have shown that this latter shock is extremely rich in iCOMs (Arce et al. 2008;Lefloch et al. 2017), suggesting that it may be an ideal case to look for others with the same method as that used by Codella et al. In this article, we continue the systematic study of the formation and destruction routes for iCOMs in the L1157-B1 shock started in CC17. Here we focus on the acetaldehyde (CH 3 CHO) and methanol (CH 3 OH). As in CC17, we make use of the observations from the IRAM/NOEMA (NOrthern Extended Millimeter Array 1 ) Large Project SOLIS (Seeds Of Life In Space: Ceccarelli et al. 2017), whose goal is to investigate the formation and destruction of iCOMs during the early stages of solar-type star formation. The article is organised as follows: in §2 we describe the source and summarise the previous studies of its astrochemistry. In §3 we describe the observations, in §4 we report the results of the analysis of these observations, in §5 we discuss the impact of those results, and in §6 we present our conclusions.

L1157-B1: an astrochemical laboratory

The chemically rich molecular outflow driven by the L1157-mm Class 0 protostar (d = 352 pc; Zucker et al. 2019) has been extensively used to investigate how the gas chemistry is altered by the injection of the species frozen onto the dust mantles. L1157mm drives an episodic and precessing jet (Gueth et al. 1996;Podio et al. 2016), which produced several cavities well observed using both single-dish antennas and interferometers (e.g. Gueth et al. 1998;Lefloch et al. 2012). In particular, the southern blueshifted outflow is associated with several shocked regions produced by impacts between jets and cavities. called B0, B1, and B2 (see also Fig. 1). The brightest shocked region, B1, is one of the targets of the IRAM NOEMA SOLIS Large Program. In particular, B1 consists of a series of shocks caused by different episodes of ejection hitting the cavity wall. These shocks produced a clumpy structure (identified using several molecular species: see e.g. Tafalla & Bachiller 1995;Codella et al. 2009;Lefloch et al. 2012;Benedettini et al. 2013), which can be classified as follows: (i) a northern arch-like structure composed of the clumps B1a-e-f-b (see Fig. 1), plus (ii) two clumps, labelled B1c and B1d, which are the oldest ones (kinematical age ∼ 1400 yr, from Podio et al. 2016, corrected for the new distance of 352 pc 1 http://iram-institute.org/EN/noema-project.php by Zucker et al. 2019), being the farthest away from the source, close to the B1 apex.

Within the context of the study of iCOMs, the advantages of observing L1157-B1 are that: (i) the observations are not polluted by the protostar emission, as it is located at a projected distance of 0.11 pc, and (ii) several previous studies have shown that grain-mantle species have been injected, due to shock-induced sputtering and shattering, into the gas phase (e.g. Bachiller et al. 2001;Arce et al. 2008;Codella et al. 2010;Fontani et al. 2014;Lefloch et al. 2017).

In the first SOLIS paper on L1157-B1 (CC17), we compared the spatial distribution of the acetaldehyde (CH 3 CHO) line emission with that from formamide (NH 2 CHO), finding a striking chemical stratification. This implies a NH 2 CHO/CH 3 COH abundance ratio which varies within the B1 structure: specifically, the formamide line emission is associated with older shocked gas with respect to acetaldehyde. Thanks to the comparison of the observations with astrochemical model predictions, we provided evidence that formamide is formed in L1157-B1 by gasphase chemical reactions. However, we were unable to provide any indication on the formation route of acetaldehyde.

Observations

The L1157-B0/B1 shocks were observed at 3 mm with the IRAM NOEMA array during two tracks in October 2016 (7 antennas), and three tracks in January 2017 (8 antennas) using both the C and D configurations. The shortest and longest baselines are 24 m and 304 m, respectively. These configurations allowed us to recover emission at scales (with a full efficiency) up to ∼ 14 . The primary beam was 64 . The phase centre is α(J2000) = 20 h 39 m 10. s 2, δ(J2000) = +68 • 01 10. 5. The CH 3 OH and CH 3 CHO lines listed in Table 1 were observed using the WideX backends, with a total spectral band of ∼ 4 GHz, and a spectral resolution of 1.95 MHz (∼ 6 km s -1 ). Several CH 3 OH lines were also observed with 80 MHz backends with a spectral resolution of 156 kHz (∼ 0.48 km s -1 ). Calibration was carried out following standard procedures using GILDAS-CLIC2 . The bandpass was calibrated on 3C454.3, while the absolute flux was fixed by observing MWC349, 2200+420, 2013+370, and 2007+659, the latter also being used to set the gains in phase and amplitude. The final uncertainty on the absolute flux scale is ≤ 15%. The phase rms was ≤ 50 • , the typical precipitable water vapour (PWV) was from 2 mm to 20 mm, and the system temperatures were ∼ 50-100 K (D) and ∼ 50-250 K (C). The rms noise in the 1.95 MHz channels was 4-20 mJy beam -1 , depending on the frequency (see Table 1). Images were produced using natural weighting, and restored with a clean beam of 2. 97 × 2. 26 (PA=-155 • ).

Results

We imaged six methanol emission lines towards the L1157 blueshifted outflow, covering upper level excitations (E u ) from 7 to 28 K (see Table 1). Figure 1 (Top) shows, as an example, the imaging of the CH 3 OH(2 1,2 -1 1,1 ) A line emission on top of the CO (1-0) image (Gueth et al. 1996), which accurately outlines the B0, B1, and B2 cavities opened by the precessing jet driven by L1157-mm, which is located at ∆α = -25 and ∆δ = +63. 5 with respect to the centre of the Fig. 1 map. Figure A.1 shows the images of the whole of the CH 3 OH line emission reported in Table 1. In addition, we detected and imaged eight acetaldehyde lines (E u = 14-23 K), which are also reported in Table 1. Figure 1 (Bottom) shows the CH 3 CHO(5 0,4 -4 0,4 ) E map, while Fig . A.2 is for all the other CH 3 CHO images.

Using the ASAI (Astrochemical Surveys At IRAM 3 ; Lefloch et al. 2017Lefloch et al. , 2018) ) unbiased spectral survey obtained with the IRAM 30m antenna, we evaluated the missing flux of the present NOEMA dataset. Comparison between the ASAI and NOEMA spectra indicates that a large fraction (about 60%) of the CH 3 OH and CH 3 CHO emission is filtered out, indicating that their emission is extended across structures larger than ∼14 . Indeed, this matches the goal of the present project, which is the analysis of the chemical content of the small structures associated with multiple shocked regions, and without any contribution due to more extended emission. In §4.3 and §5.1 the present results will be compared with those derived by the analysis based on singledish observations, showing there is no bias introduced by the missing flux.

CH 3 OH and CH 3 CHO spatial distribution

The present NOEMA CH 3 OH data improve the images and spectra based on IRAM-PdBI observations of the methanol spectral pattern at 97.94 GHz performed in 1997 using five antennas, and with a synthesised beam of ∼ 5 , as well as those of the 2 1,1 -1 1,0 A line, with a spectral resolution of 0.5 km s -1 (Benedettini et al. 2007). The eight CH 3 CHO images also allow us to perform a deeper analysis with respect to the spatial distributions obtained by Codella et al. (2015), who stacked the CH 3 CHO (7 0,7 -6 0,6 ) E and A lines with an angular resolution of ∼ 2. 5.

The present results are consistent with the previous observations that show that the CH 3 OH and CH 3 CHO spatial distributions are in good agreement; both trace the B0 eastern molecular cavity as well as the northern arch-like structure containing B1ae-f-b. In particular, the peak emission of both species is clearly associated with the B1b clump. In other words, although emission is observed also towards B1c, the bulk of the CH 3 OH and CH 3 CHO emission does not come from the southern B1 apex (Codella et al. 2017). Instead, they preferentially trace the region associated with the youngest B1 shock occurring next to the B1a position, where HDCO, a selective tracer of dust mantle release, has been observed (Fontani et al. 2014), as well as high-velocity SiO emission (Gueth et al. 1998;Podio et al. 2016), locating the likely current sputtering of the dust refractory cores due to the shock. These findings are not surprising for methanol, which is efficiently produced on the grain surfaces by CO hydrogenation. On the other hand, they indicate that acetaldehyde is either formed on the mantles as well or that it is quickly (in less than ∼ 1000 yr) produced in the gas phase using simpler mantle products (as suggested by Codella et al. 2015). This topic is further discussed in §5.2.

Line spectra

Figure 2 shows the spectra extracted from the CH 3 OH and CH 3 CHO emission peaks at the B1b position. In addition to methanol and acetaldehyde, in B1b we observe emission at 99.3 GHz due to another iCOM, dimethyl ether (CH 3 OCH 3 ). In particular, the detected line is due to the 4 1,4 -3 0,3 transitions of the EA, AE, EE, and AA isomers, blended at the ∼ 2 MHz spectral resolution (see Fig. 2). Finally, we detect several S-bearing species (CS, SO, and OCS isotopologues), which will be anal-3 www.oan.es/asai Fig. 1. The L1157 southern blueshifted lobe in CO (1-0) (white contours; Gueth et al. 1996). The precessing jet ejected by the central object L1157-mm (outside the frame and toward the northwest) excavated several clumpy cavities, named B0, B1, and B2, respectively. The maps are centred at α(J2000) = 20 h 39 m 10. s 2, δ(J2000) = +68 • 01 10. 5 (∆α = +25 and ∆δ = -63. 5 from the L1157-mm protostar). Top panel: Map (in colour scale) of the emission due to the CH 3 OH(2 1,2 -1 1,1 ) A transition (integrated over the whole velocity range). The first contour and step are 3σ (32 mJy beam -1 km s -1 ). For the CO image, the first contour and step are 6σ (1σ = 0.5 Jy beam -1 km s -1 ) and 4σ, respectively. The dashed circle shows the primary beam of the CH 3 OH image (64 ). Yellow labels are for the B0 and B1 clumps (black triangles) previously identified using molecular tracers (see text, and e.g. Benedettini et al. 2007;Codella et al. 2009Codella et al. , 2017)). The white and magenta ellipses depict the synthesised beams of the CO (3. 65 × 2. 96, PA=+88 • ), and CH 3 OH (2. 97 × 2. 26, PA=-155 • ) observations, respectively. Bottom panel: Same as top panel but for CH 3 CHO (5 0,4 -4 0,4 ) E. The first contour and step are 3σ (7 mJy beam -1 km s -1 ). The synthesised beams of the CH 3 CHO image are similar to those of the CH 3 OH map.

Table 1. CH 3 OH and CH 3 CHO emission lines observed using NOEMA and the WideX backend towards the L1157-B1b peak (see Fig. 2). 2 0,2 -1 0,1 A c 96.74138 7 1.6 -5.6 9 57.2(1.5) d 2 0,2 -1 0, 1.50(0.12) a Frequencies and spectroscopic parameters were extracted from the Jet Propulsion Laboratory molecular database (Pickett et al. 1998) for CH 3 CHO and CH 3 OCHO. Upper level energies refer to the ground state of each symmetry. b Integrated over the whole velocity emission range from the WideX spectra extracted at the CH 3 CHO and CH 3 OH emission peak (B1b; see text). In case of non-detections we report the 3σ limit. We note that the errors on the integrated areas I int are propagated using the r.m.s. and the spectral resolution of the spectra. In addition, the uncertainty on the absolute flux scale due to calibration is ≤ 15% (see Sect. 3). c The three lines at 96.74 GHz have been deblended using the 0.5 km s -1 spectral resolution (see Sect. 3.2 and Fig. 4).

Transition a ν 0 a E u a S µ 2 a log(A i j ) a rms I int b T rot N tot (GHz) (K) (D 2 ) (mK) (K km s -1 ) (K) (cm -2 ) CH 3 OH 10.0(1.1) 1.3(0.3) × 10 16 2 1,2 -1 1,1 A 95.91431 21 1.2 -5.6 14 9.05(0.15) 2 -1,2 -1 -1,1 E c 96.73936 13 1.2 -4.6 9 40.5(1.5) d
ysed in a forthcoming paper on sulphur chemistry (Feng et al., in preparation).

For each detected line, in addition to the spectroscopic parameters, Table 1 reports the intensity integrated over the whole velocity emission range (down to velocities blueshifted by 15 km s -1 ). The 2 MHz spectral resolution of the WideX backend does not allow us to investigate the iCOMs kinematics. However, we can confirm that all the spectra are blueshifted, in agreement with all the molecular lines so far observed towards L1157-B1 (e.g. Bachiller et al. 2001;Lefloch et al. 2017). In addition, the CH 3 OH spectral pattern at 96.7 GHz (see Table 1) has been observed with a 0.48 MHz (0.16 km s -1 ) spectral resolution. Figure 3 reports the spectra as observed towards the B0e and B1b clumps: four transitions of the 2 k,k -1 k,k pattern are reported, with the corresponding frequencies marked by small vertical red lines. The spectra are centred at 96741.38 MHz, which is the frequency of the 2 0,2 -1 0,1 A line. The methanol lines extracted from the B0e position have a peak velocity of ∼ +0.5 km s -1 (V LSR = +2.6 km s -1 ), and are about 3-4 km s -1 in width (FWHM), with evidence of bluer wings. The B1b spectrum is definitely more complex showing in addition a secondary peak at about -4 km s -1 . This high-velocity peak was tentatively detected by Benedettini et al. (2013) using the CH 3 OH(3 k,k -2 k,k ) emission at 2mm. In the present dataset, a clear Gaussian-like secondary peak is emerging, definitely offset in velocity from the mean peak. These findings confirm that the northern part of the B1 structure is associated with young shocked gas, here traced by methanol molecules injected into the gas phase at high velocities.

Figure 4 shows the spectral fit of the 2 -1,2 -1 -1,1 E, 2 0,2 -1 0,1 A, and 2 0,2 -1 0,1 E triplets as observed towards B1b. The velocity scale is based on the rest frequency of the brighter line, 2 0,2 -1 0,1 A, at 96741.38 MHz. In order to perform the fit, we assumed that the three low-excitation methanol transitions are associated with the same spectral profile. As a consequence, the rest frequencies of the transitions give us the expected velocity shifts among the three lines (see the red small vertical segments in Fig. 4). In addition, both 2 0,2 -1 0,1 E and A lines (blue and magenta) show a secondary peak (3.0 km s -1 broad) offset by 5 km s -1 with respect the main one (which has a line width of 4.6 km s -1 ).

In conclusion, we fixed the position and the line width of the secondary peak expected for the 2 -1,2 -1 -1,1 E line (grey). The overall fit is satisfactory, with a residual with an rms of 480 mK.

Large Velocity Gradient analysis and rotational diagrams: column densities and abundances

Given the lack of CH 3 CHO collisional rates in the literature, we only analysed the CH 3 OH emission with the non-local thermodinamic equilibrium (NLTE) large velocity gradient (LVG) approach using the model described in Ceccarelli et al. (2003). We used the CH 3 OH-H 2 collisional coefficients computed by Rabli & Flower (2010) and provided by the BASECOL database (Dubernet et al. 2013). In the computations, we assumed a H 2 ortho-to-para ratio of one (Nisini et al. 2010). We run grids of models varying the kinetic temperature (T kin ) from 50 to 300 K, the volume density (n H 2 ) from 10 4 to 10 8 cm -3 , and the methanol column density (N(CH 3 OH)) from 10 14 to 10 18 cm -2 . The best fit is obtained by the minimum χ 2 in the three parameters. The errors on the observed line intensities have been obtained by propagating the spectral rms with the uncertainties due to calibration (15%). Figure 5 shows: (i) the χ 2 contour plot in the T kin -n H 2 plane obtained with the best fit of the methanol column density N(CH 3 OH) = 5× 10 15 cm -2 (3-8 × 10 15 cm -2 , including uncertainties), (ii) the χ 2 -N(CH 3 OH) distribution, where it is possible to identify the minimum χ 2 value (2.4), and (iii) the ratio between the measured and predicted intensity of the observed lines. Overall, the best fit is very good, with a reduced χ 2 equal to 1.2, which for two degrees of freedom corresponds to a probability of 30% to exceeding χ 2 (i.e. approximately the standard 1 σ). The methanol line opacities are predicted to be less than Fig. 2. Upper and Lower panels: IRAM-WideX spectra (in F ν scale) extracted at the L1157-B1b blueshifted position (see Fig. 1 andtext). The spectral resolution is 1.95 MHz (see Sect. 3). A zoom-in in intensity is shown so that the weakest lines can be seen. The transitions producing the emission lines are labelled in blue, with the corresponding frequencies marked by vertical blue lines. The detected CH 3 OH and CH 3 COH lines (analysed here) are reported in Table 1. 0.1. The volume density is well constrained, n H 2 = 4× 10 5 cm -3 (2-10 × 10 5 cm -3 ), in agreement with the high densities (≥ 10 5 cm -3 ) inferred toward the molecular cavities by Gómez-Ruiz et al. ( 2015) from a CS multiline analysis. The kinetic temperature is also well constrained, T kin = 90 K (70-130 K), again in agreement with the temperatures of the B1 cavity as derived by Lefloch et al. (2012) using a CO multiline single-dish analysis. We note that a recent LVG analysis of methanol emission as observed using single-dish towards a large number of protostellar jet-driven shocks (Holdship et al. 2019) shows volume densities in the 10 5 -10 6 cm -3 range and kinetic temperatures between 20 and 60 K, values consistent with what is found in L1157-B1.

The analysis described immediately above allows us to constrain the density and temperature of the methanol-emitting gas. However, in order to derive the relative abundances of methanol and acetaldehyde, we carried out a rotational diagram (RD) analysis for both molecules, where a LTE population and optically thin lines are assumed (in agreement with the NLTE LVG analysis of methanol described above). Under these assumptions, for a given molecule, the relative population distribution of all the energy levels is described by a Boltzmann temperature, that is the Fig. 3. CH 3 OH emission spectra at 96.7 GHz (in T B scale) extracted at the B0e and B1b positions (see Fig. 1 andtext). The spectral resolution is 0.48 MHz, corresponding to 0.16 km s -1 . The transitions producing the emission lines (see Table 1) are reported, with the corresponding frequencies marked by small vertical red lines. The spectra are centred at the frequency of the 2 0,2 -1 0,1 A transition: 96741.38 MHz. The dashed vertical line stands for the systemic velocity: +2.6 km s -1 (e.g. Bachiller et al. 2001). rotational temperature T rot . As demonstrated by Lefloch et al. (2102), for example, by analysing CO emission from extended outflowing gas in L1157, the column density obtained with RD gives results that are in relatively good agreement with a NLTE LVG analysis, usually within a factor two, as is also the case for our analysis.

Figure 6 shows the RD of CH 3 OH, which provides a column density of (1.3±0.3) ×10 16 cm -2 , close to what was found using the LVG method. The rotational temperature is (10.0±1.1) K, in agreement with what was found using previous interferometric and single dish measurements of CH 3 OH lines in the same excitation range (12 K; Bachiller et al. 1995;Benedettini et al. 2007;Codella et al. 2010). Thanks to the many CH 3 CHO detected lines (covering an upper level energy E u range similar to that traced by the methanol lines) we can carry out the RD analysis also for acetaldehyde, obtaining, as for CH 3 OH, a low rotational temperature, (8±1) K, and a total column density N(CH 3 CHO) = (7±3) ×10 13 cm -2 . Lefloch et al. (2017) report for CH 3 CHO a rotational temperature of 17 K, but this was measured including transitions with E u up to 94 K. Figure 7 shows the spatial distribution of the CH 3 CHO (upper panel) and CH 3 OH (lower panel) rotational temperature in the B0 and B1 regions of the L1157 outflow. We note that the CH 3 CHO rotational temperature map covers a smaller portion of the L1157 region due to the weakness of the acetaldehyde lines with E u = 23 K. Considering the uncertainties (B0e and B1e: ∼ 6 K; B1b: ∼ 3 K) there is no sig-Fig. 4. Spectral fit of the 2 -1,2 -1 -1,1 E (grey lines), 2 0,2 -1 0,1 A (magenta), and 2 0,2 -1 0,1 E (blue) triplets as observed towards B1b using the NOEMA narrow-band spectrometer. The spectra are centred at the frequency of the 2 0,2 -1 0,1 A transition: 96741.38 MHz. The transitions producing the emission lines (see Table 1) are reported, with the corresponding frequencies marked by small vertical red lines. The dashed vertical line is for the systemic velocity: +2.6 km s -1 (e.g. Bachiller et al. 2001). Each transition shows a main peak at +0.5 km s -1 plus a secondary peak further blueshifted by ∼ 5 km s -1 . nificant difference between the rotational temperature of the two species or spatial trends throughout the L1157-B1 structure.

Finally, we computed the abundances of the various species using the values of their respective column densities divided by the H 2 column density, which was derived from the CO column density in the B1 molecular cavity (N CO = 2 × 10 17 cm -2 : Lefloch et al. 2012) assuming the standard [CO]/[H] = 5 × 10 -5 value. In this way, we obtain methanol and acetaldehyde abundances with respect to H equal to 6.5 × 10 -6 (1.2-12 × 10 -6 , considered the uncertainties) and 3.5 × 10 -8 (2-5 × 10 -6 ), respectively, in good agreement with that derived by CC17 for acetaldehyde (1-3 × 10 -8 ) using a coarser angular resolution (5 -6 ).

Discussion

5.1. The CH 3 OH to CH 3 CHO abundance ratio Codella et al. 2017 showed that the spatial distribution of the measured abundance ratio of iCOMs can lead to very stringent constraints on how they are formed and destroyed. In particular, these latter authors showed that the formamide to acetaldehyde abundance ratio towards the L1157-B1 shock demonstrates a gas-phase formation route of formamide (e.g. Skouteris et al. 2017). Further analysis of the abundance ratio of acetaldehyde with respect to methanol is therefore worthwhile, which can be safely assumed to form on dust mantles and then injected into the gas phase. Based on the interferometric maps showing that both species are populating the same gas portion, the ratio between their column densities can be used to derive their abundance ra- tio. The ratio of CH 3 OH to CH 3 CHO is about 1.5 × 10 2 and 1.9 × 10 2 as measured towards the B0e and B1b peaks, respectively, and it varies by less than a factor of three throughout the whole B0-B1 structure, as shown in Fig. 8. Considering that the uncertainty on these measurements is about 30%, there is no evidence for significant variation of this latter ratio in the observed region. Furthermore, these measurements are consistent with those de-Fig. 6. Rotational diagrams for CH 3 OH (upper panel), and CH 3 CHO (lower panel) derived using the emission lines observed towards B1b (see Table 1 and Fig. 3). The parameters N u , g u , and E up are the column density, the degeneracy, and the energy (with respect to the ground state of each symmetry), respectively, of the upper level. The derived values of the rotational temperature are reported in the panels. rived using the ASAI IRAM 30m single-dish survey (1-3 × 10 2 ; Codella et al. 2015, Lefloch et al. 2017). We note that Holdship et al. ( 2019) measured CH 3 OH/CH 3 CHO abundance ratios towards a sample of shocks using the IRAM 30m antenna, finding higher values (by a factor of six on average): 0.4-1.3 × 10 3 ; highspatial-resolution surveys are needed to confirm these values.

It is interesting to check whether the CH 3 OH/CH 3 CHO abundance ratios derived for L1157-B0 and L1157-B1 are different with respect to measurements towards known Class 0 hot corinos, which rather trace the inner 100 au from the protostar where the temperature is high enough (at least 100 K) to thermally evaporate the dust-mantle products into the gas phase, or where grain mantles are sputtered by shocks (see below). If we only include in the comparison those sources imaged with high-spatial-resolution interferometric observations (to be sure of a proper comparison between the CH 3 OH and CH 3 CHO spatial distributions), the sample size is limited (see Table 2): HH212-mm (Lee et al. 2017, 2019, Bianchi et al. 2017, Codella et al. 2018), IRAS16293-2422B (Jørgensen et al. 2016(Jørgensen et al. , 2018)), NGC1333-IRAS4A2 (López-Sepulcre et al. 2017), Barnard 1b-S (Marcelino et al. 2018), and L483 (Jacobsen et al. 2018). Table 2 shows that the CH 3 OH/CH 3 CHO ratio of the hot-corinos varies from ∼ 90 to ∼ 500, that is, by less than a factor of six. The Fig. 8. Map of the ratio between the CH 3 OH and CH 3 CHO column density of the B0 and B1 region of the L1157 outflow (derived from the NOEMA images, where the emission is at least 3σ), overlaid on the CO (black contours; Gueth et al. 1996) emission map. Uncertainties are 30% (i.e. 30-100 in absolute values). Symbols are as in Fig. 1.

Table 2. Comparison of the CH 3 OH/CH 3 CHO abundance ratio as derived towards L1157-B1 with those measured towards Class 0 hot corinos using interferometric observations. Object CH 3 OH/CH 3 CHO References L1157-B0e 1.5 × 10 2 1 L1157-B1b 1.9 × 10 2 1 HH212-mm 0.9-3 × 10 2 2,3,4,5 IRAS16293-2422B 1 × 10 2 6,7 NGC1333-IRAS4A2 3-6 × 10 2 8,9 Barnard 1b-S 5 × 10 2 10 L483 2 × 10 2 11 References: 1. Present work; 2. Lee et al. (2017) L1157 values fall inside this range, suggesting that the chemistry at work in the protostellar shocks, at least in the CH 3 OH and CH 3 CHO context, is similar to that ruling in the inner 100 au of the protostellar region. Of course, this is speculation based on one shocked region, and verification is required: (i) by increasing similar measurements towards different outflow shocks, and (ii) by investigating whether the chemical enrichment around a protostar is ruled by thermal evaporation or sputtering due to accretion shocks (see e.g. Sakai et al. 2014ab, Lee et al. 2017, 2019;Codella et al. 2018).

Chemical modelling

The overlap between the spatial distribution of the CH 3 OH and CH 3 CHO emission suggests that either acetaldehyde or a par-ent molecule was previously formed on the grain mantles, such as methanol for example, and was injected into the gas because of the passage of the shock. In addition, given that the CH 3 OH/CH 3 CHO abundance ratio is similar in the B0 and B1 knots despite them having a different kinematical age of less than about 1000 years, if CH 3 CHO is synthesised in the gas phases, the process must be relatively fast. We note that recently, Burkhardt et al. (2019) modelled the shocked chemistry in L1157-B1 finding that CH 3 CHO has a post-shock chemistry due to sputtering, and additionally that a significant amount of its formation is due to gas-phase processes before 10 3 yr. The goal of this section is to explore whether gas-phase reactions are indeed fast enough to be able to reproduce the observations.

To this end, we used the same astrochemical model adopted in CC17 to analyse the formamide and acetaldehyde emission towards L1157-B1. The full details of the model can be found in CC17. Briefly, we used the time-dependent model gas-phase code MyNahoon, and adopted a two-step procedure: (1) We first compute the chemical composition of a standard molecular cloud at 10 K with a H 2 density of 2 × 10 4 cm -3 and initial elemental abundances as in Wakelam & Agundez (2013) (their Table 3). ( 2) We then successively increase the gas temperature, density, and gaseous abundance of grain-mantle species to simulate the shock passage. The gas temperature and H 2 density are taken as those found by the NLTE LVG modelling, namely 90 K and 4×10 5 cm -3 , respectively. For the abundances of the species injected into the gas phase at the moment of the passage of the shock, when possible, we assumed abundances (with respect to H nuclei) similar to those measured by IR observations of the interstellar dust ices (Boogert et al. 2015): CO 2 (3 × 10 -5 ) and H 2 O (2 × 10 -4 ; see also Busquet et al. 2014). Otherwise, we used the values constrained by previous studies on the L1157-B1 chemistry: OCS (2 × 10 -6 ; Podio et al. 2014), NH 3 (2 × 10 -5 ; Tafalla & Bachiller 1995). We varied the injected CH 3 CH 2 abundance in the gas-phase acetaldehyde formation model so as to obtain the observed acetaldehyde abundance. We found it necessary to inject CH 3 CH 2 /H = 8 × 10 -8 . Similarly, in the grain-surface acetaldehyde formation model, we modified the acetaldehyde abundances to fit the measured value in B1b (see § 4.3), namely we assumed CH 3 CHO/H = 3.5 × 10 -8 . Finally, in both models, methanol is assumed to be a grain-surface product and the abundance injected in the gas phase is equal to that observed, namely 6.5 × 10 -6 (see § 4.3).

As in CC2017, we run both the case where acetaldehyde is directly injected into the gas phase from the grain mantles and the case in which it is formed in the gas phase, mainly by the reaction of ethyl radical (CH 3 CH 2 ) with atomic oxygen: CH 3 CH 2 + O → CH 3 CHO + H, as suggested by Charnley et al. (2004). With respect to the chemical network used by CC2017, we included the recent theoretical study by Gao et al. (2018) on the destruction of methanol by the reaction with OH.

Figure 9 reports the model predictions as a function of time after the shock passage compared with the values measured towards B1b. As written above, if acetaldehyde is directly sputtered from the grain mantles, in order to fit the observations, its injected abundance must be equal to 3.5 × 10 -8 . On the other hand, if acetaldehyde is formed via gas-phase reactions, the injected ethyl radical abundance must be 8 × 10 -8 . As shown in the figure, the gas phase reactions are relatively fast at synthesising acetaldehyde: within 300 years the acetaldehyde abundance reaches its peak, after which it slowly decreases because of the destruction from ions in the gas phase, on a timescale larger than about 2000 years. Therefore, it is impossible to elucidate the formation route of acetaldehyde in L1157-B1: a similar study to-Fig. 9. Model predictions of methanol and acetaldehyde after a shock passage. Upper panel: CH 3 OH/CH 3 CHO abundance ratio: The solid line refers to a model where acetaldehyde is synthesised in the gas phase, whereas the dashed line refers to predictions assuming that CH 3 OH is injected into the gas phase directly from the grain mantles. The hatched green zone is defined by the value (including uncertainty) measured towards the B1b peak, with a kinematical age of 1100 yr (Podio et al. 2016, revised according to the new distance of 352 pc by Zucker et al. 2019). Lower panel: Acetaldehyde (CH 3 CHO, red) and methanol (CH 3 OH, blue) abundances with respect to H nuclei as a function of time from the passage of the shock. The hatched blue and red regions show the measurements towards B1b (see text).

wards an even younger shock might be able to shed light on the question. Alternatively, a less dense shock region, or one that was less irradiated by cosmic rays could be good site for better constraining the acetaldehyde formation route, as the chemical evolution would be slowed down in those cases.

As a final remark, we would like to highlight the fact that, although relatively simple, the modelling described above and used in this work catches the basic processes at work after the passage of a shock and has proven to be relatively reliable for reproducing the observations (see e.g. Podio et al. 2014 andC2017). The most drastic shortcoming of our model compared to more accurate and/or sophisticated shock models is its lack of a higher temperature (up to 1000 K) over a short (a few hundred years at most) period (see e.g. Viti et al. 2011). This could be important in the presence of species whose abundances depend on reactions with activation barriers that make them inefficient at low (≤100 K) temperatures but efficient at higher ones. In the specific case modelled here, we verified that no similar reactions are present and that a high-temperature (1000 K) phase does not impact our results.

Conclusions

In the context of the IRAM NOEMA SOLIS Large Program we imaged six methanol and eight acetaldehyde emission lines toward the B0 and B1 shocks along the L1157 blueshifted outflow.

The measured abundances relative to H nuclei are 6.5 × 10 -6 (methanol) and 3.5 × 10 -8 (acetaldehyde). The CH 3 OH and CH 3 CHO spatial distributions overlap well, tracing the earliest shocked regions. Comparison with astrochemical model predictions shows that CH 3 CHO is either formed on the dust mantles or is quickly formed in the gas phase using simpler mantle products. Our modelling demonstrates that two species emitting lines in the same region is not enough to guarantee that they are grain-surface products: gas-phase reactions triggered by different species from the grain mantles could also display the same spatial behaviour. Detailed studies are required in order to discriminate between grain-surface and gas-phase chemistry as the major process leading to the observed species and to decipher the timescale on which this formation occurs.

For acetaldehyde in particular, at present, it is not possible to measure the abundance of one of the two parent species (ethyl radical) possibly synthesising acetaldehyde in the gas phase. The measurement of the second parent species, atomic oxygen, is challenging, as its fine structure line at 63 µm lies in a spectral window that is totally blocked by the terrestrial atmosphere. The SOFIA telescope could overcome this problem (Kristensen et al. 2017). On the other hand, although the frequencies (and spectroscopic parameters) of ethyl radicals are unknown at present, they would be a very welcome target of specific spectroscopic studies. It would then very likely be possible to search for ethyl radicals in regions where acetaldehyde is detected, and thus decipher whether or not the proposed gas-phase route is correct. Gueth et al. 1996) and NH 2 CHO (4 1,4 -3 1,3 ) (cyan contours; Codella et al. 2017). The precessing jet ejected by the central object L1157-mm (outside the frame and toward the northwest) excavated several clumpy cavities, named B0, B1, and B2, respectively. The maps are centred at α(J2000) = 20 h 39 m 10. s 2, δ(J2000) = +68 • 01 10. 5 (∆α = +25 and ∆δ = -63. 5 from the L1157-mm protostar). Map (in colour scale) of the sum of the CH 3 OH (2 -1,2 -1 -1,1 ) E , (2 0,2 -1 0,1 ) A, and (2 0,2 -1 0,1 ) E emission lines (labelled 2 k,k -1 k,k ), (2 1,2 -1 1,1 ) A, (2 1,1 -1 1,0 ) A, and (2 1,1 -1 1,0 ) E (integrated over the whole velocity range). For the CO image, the first contour and step are 6σ (1σ = 0.5 Jy beam -1 km s -1 ) and 4σ, respectively. The first contour and step of the NH 2 CHO map (cyan contours) correspond to 3σ (15 mJy beam -1 km s -1 ) and 1σ, respectively. The dashed circle shows the primary beam of the CH 3 OH images (64 ). The first contour and step are 3σ (56 mJy beam -1 km s -1 , 32 mJy beam -1 km s -1 , 45 mJy beam -1 km s -1 , and 23 mJy beam -1 km s -1 for 2 k,k -1 k,k , 2 1,2 -1 1,1 A, 2 1,1 -1 1,0 A, and (2 1,1 -1 1,0 ) E, respectively. Yellow labels are for the B0, B1, and B2 regions: the positions of the different clumps inside these regions are shown in Fig. 1 (see also text, and e.g. Codella et al. 2009Codella et al. , 2017)). The cyan, white, and magenta ellipses depict the synthesised beams of the NH 2 CHO (5. 79 Gueth et al. 1996) and NH 2 CHO (4 1,4 -3 1,3 ) (cyan contours; Codella et al. 2017). The precessing jet ejected by the central object L1157-mm (outside the frame and toward the northwest) excavated several clumpy cavities, named B0, B1, and B2, respectively. The maps are centred at α(J2000) = 20 h 39 m 10. s 2, δ(J2000) = +68 • 01 10. 5 (∆α = +25 and ∆δ = -63. 5 from L1157-mm). For the CO image, the first contour and step are 6σ (1σ = 0.5 Jy beam -1 km s -1 ) and 4σ, respectively. The first contour and step of the NH 2 CHO map (cyan contours) correspond to 3σ (15 mJy beam -1 km s -1 ) and 1σ, respectively. The dashed circle shows the primary beam of the CH 3 CHO images(64 ). Upper panels: Map (in colour scale) of line emission (integrated over the whole velocity range) due to CH 3 CHO (5 0,4 -4 0,4 ) E, (5 0,4 -4 0,4 ) A, (5 1,4 -4 1,3 ) E, and (5 1,4 -4 1,3 ) A, respectively. The first contour and step are 3σ: 7 mJy beam -1 km s -1 (5 0,4 -4 0,4 E, 5 0,4 -4 0,4 A), and 8 mJy beam -1 km s -1 (5 1,4 -4 1,3 E, and 5 1,4 -4 1,3 A). The cyan, white, and magenta ellipses depict the synthesised beams of the NH 2 CHO (5. 79 × 4. 81, PA = -94 • ), CO (3. 65 × 2. 96, PA=+88 • ), and CH 3 CHO (2. 97 × 2. 26, PA=-155 • ) observations, respectively). Yellow labels are for the B0, B1, and B2 regions: the positions of the different clumps inside these regions are shown in Fig. 1 (see also text, and e.g. Codella et al. 2009Codella et al. , 2017)). Lower panels: Same as upper panels but for the CH 3 CHO (5 2,4 -4 2,3 ) A, (5 2,4 -4 2,3 ) E, (5 2,3 -5 2,2 ) E, and (5 2,3 -4 2,2 ) A. The first contour and step are 3σ: 5 mJy beam -1 km s -1 for 5 2,4 -4 2,3 A,E, and 4 mJy beam -1 km s -1 for 5 2,3 -5 2,2 A,E.

Fig. 5 .

 5 Fig. 5. Upper panel: Density-temperature contour plot of χ 2 obtained considering the non-LTE LVG-model-predicted and observed intensity of all the A and E CH 3 OH emission lines detected (beam averaged) towards the B1b clump. The best fit is obtained with N(CH 3 OH) = 5 × 10 15 cm -2 , T kin = 90 K, and n H 2 = 4× 10 5 cm -3 . Middle panel: χ 2 -N(CH 3 OH) distribution, identifying the minimum χ 2 value (2.4). Lower panel: Ratio between observations and model predictions of the CH 3 OH A (circles) and E (stars) line intensities as a function of the upper level energy of the lines.

Fig. 7 .

 7 Fig. 7. CH 3 CHO (upper panel) and CH 3 OH (lower panel) rotational temperature (T rot ) map (in K scale) of the B0 and B1 region of the L1157 outflow (derived from the NOEMA images, where the emission is at least 3σ), overlaid on the CO (black contours; Gueth et al. 1996) emission map. Uncertainties on T rot are ∼ 3 K towards B1b, and ∼ 6 K for B1e and B0. Symbols are as in Fig. 1.

  ; 3. Lee et al. (2019); 4. Bianchi et al. (2017); 5. Codella et al. (2018); 6. Jørgensen et al. (2016); 7. Jørgensen et al. 2018; 8. Taquet et al. (2015); 9. López-Sepulcre et al. (2017); 10. Marcelino et al. (2018); 11. Jacobsen et al. (2018).

Fig

  Fig. A.1. The L1157 southern blueshifted lobe in CO (1-0) (white contours;Gueth et al. 1996) and NH 2 CHO (4 1,4 -3 1,3 ) (cyan contours;Codella et al. 2017). The precessing jet ejected by the central object L1157-mm (outside the frame and toward the northwest) excavated several clumpy cavities, named B0, B1, and B2, respectively. The maps are centred at α(J2000) = 20 h 39 m 10. s 2, δ(J2000) = +68 • 01 10. 5 (∆α = +25 and ∆δ = -63. 5 from the L1157-mm protostar). Map (in colour scale) of the sum of the CH 3 OH (2 -1,2 -1 -1,1 ) E , (2 0,2 -1 0,1 ) A, and (2 0,2 -1 0,1 ) E emission lines (labelled 2 k,k -1 k,k ), (2 1,2 -1 1,1 ) A, (2 1,1 -1 1,0 ) A, and (2 1,1 -1 1,0 ) E (integrated over the whole velocity range). For the CO image, the first contour and step are 6σ (1σ = 0.5 Jy beam -1 km s -1 ) and 4σ, respectively. The first contour and step of the NH 2 CHO map (cyan contours) correspond to 3σ (15 mJy beam -1 km s -1 ) and 1σ, respectively. The dashed circle shows the primary beam of the CH 3 OH images (64 ). The first contour and step are 3σ (56 mJy beam -1 km s -1 , 32 mJy beam -1 km s -1 , 45 mJy beam -1 km s -1 , and 23 mJy beam -1 km s -1 for 2 k,k -1 k,k , 2 1,2 -1 1,1 A, 2 1,1 -1 1,0 A, and (2 1,1 -1 1,0 ) E, respectively. Yellow labels are for the B0, B1, and B2 regions: the positions of the different clumps inside these regions are shown in Fig.1(see also text, and e.g.Codella et al. 2009Codella et al. , 2017)). The cyan, white, and magenta ellipses depict the synthesised beams of the NH 2 CHO (5. 79 × 4. 81, PA = -94 • ), CO (3. 65 × 2. 96, PA=+88 • ), and CH 3 OH (2. 97 × 2. 26, PA=-155 • ) observations, respectively.

  Fig. A.1. The L1157 southern blueshifted lobe in CO (1-0) (white contours;Gueth et al. 1996) and NH 2 CHO (4 1,4 -3 1,3 ) (cyan contours;Codella et al. 2017). The precessing jet ejected by the central object L1157-mm (outside the frame and toward the northwest) excavated several clumpy cavities, named B0, B1, and B2, respectively. The maps are centred at α(J2000) = 20 h 39 m 10. s 2, δ(J2000) = +68 • 01 10. 5 (∆α = +25 and ∆δ = -63. 5 from the L1157-mm protostar). Map (in colour scale) of the sum of the CH 3 OH (2 -1,2 -1 -1,1 ) E , (2 0,2 -1 0,1 ) A, and (2 0,2 -1 0,1 ) E emission lines (labelled 2 k,k -1 k,k ), (2 1,2 -1 1,1 ) A, (2 1,1 -1 1,0 ) A, and (2 1,1 -1 1,0 ) E (integrated over the whole velocity range). For the CO image, the first contour and step are 6σ (1σ = 0.5 Jy beam -1 km s -1 ) and 4σ, respectively. The first contour and step of the NH 2 CHO map (cyan contours) correspond to 3σ (15 mJy beam -1 km s -1 ) and 1σ, respectively. The dashed circle shows the primary beam of the CH 3 OH images (64 ). The first contour and step are 3σ (56 mJy beam -1 km s -1 , 32 mJy beam -1 km s -1 , 45 mJy beam -1 km s -1 , and 23 mJy beam -1 km s -1 for 2 k,k -1 k,k , 2 1,2 -1 1,1 A, 2 1,1 -1 1,0 A, and (2 1,1 -1 1,0 ) E, respectively. Yellow labels are for the B0, B1, and B2 regions: the positions of the different clumps inside these regions are shown in Fig.1(see also text, and e.g.Codella et al. 2009Codella et al. , 2017)). The cyan, white, and magenta ellipses depict the synthesised beams of the NH 2 CHO (5. 79 × 4. 81, PA = -94 • ), CO (3. 65 × 2. 96, PA=+88 • ), and CH 3 OH (2. 97 × 2. 26, PA=-155 • ) observations, respectively.

Fig

  Fig. A.2. The L1157 southern blueshifted lobe in CO (1-0) (white contours;Gueth et al. 1996) and NH 2 CHO (4 1,4 -3 1,3 ) (cyan contours;Codella et al. 2017). The precessing jet ejected by the central object L1157-mm (outside the frame and toward the northwest) excavated several clumpy cavities, named B0, B1, and B2, respectively. The maps are centred at α(J2000) = 20 h 39 m 10. s 2, δ(J2000) = +68 • 01 10. 5 (∆α = +25 and ∆δ = -63. 5 from L1157-mm). For the CO image, the first contour and step are 6σ (1σ = 0.5 Jy beam -1 km s -1 ) and 4σ, respectively. The first contour and step of the NH 2 CHO map (cyan contours) correspond to 3σ (15 mJy beam -1 km s -1 ) and 1σ, respectively. The dashed circle shows the primary beam of the CH 3 CHO images(64 ). Upper panels: Map (in colour scale) of line emission (integrated over the whole velocity range) due to CH 3 CHO (5 0,4 -4 0,4 ) E, (5 0,4 -4 0,4 ) A, (5 1,4 -4 1,3 ) E, and (5 1,4 -4 1,3 ) A, respectively. The first contour and step are 3σ: 7 mJy beam -1 km s -1 (5 0,4 -4 0,4 E, 5 0,4 -4 0,4 A), and 8 mJy beam -1 km s -1 (5 1,4 -4 1,3 E, and 5 1,4 -4 1,3 A). The cyan, white, and magenta ellipses depict the synthesised beams of the NH 2 CHO (5. 79 × 4. 81, PA = -94 • ), CO (3. 65 × 2. 96, PA=+88 • ), and CH 3 CHO (2. 97 × 2. 26, PA=-155 • ) observations, respectively). Yellow labels are for the B0, B1, and B2 regions: the positions of the different clumps inside these regions are shown in Fig.1(see also text, and e.g.Codella et al. 2009Codella et al. , 2017)). Lower panels: Same as upper panels but for the CH 3 CHO (5 2,4 -4 2,3 ) A, (5 2,4 -4 2,3 ) E, (5 2,3 -5 2,2 ) E, and (5 2,3 -4 2,2 ) A. The first contour and step are 3σ: 5 mJy beam -1 km s -1 for 5 2,4 -4 2,3 A,E, and 4 mJy beam -1 km s -1 for 5 2,3 -5 2,2 A,E.

http://www.iram.fr/IRAMFR/GILDAS

Acknowledgements. We are very grateful to all the IRAM staff, whose dedication allowed us to carry out the SOLIS project. This work was supported by (i) the PRIN-INAF 2016 "The Cradle of Life -GENESIS-SKA (General Conditions in Early Planetary Systems for the rise of life with SKA)", (ii) the program PRIN-MIUR 2015 STARS in the CAOS -Simulation Tools for Astrochemical Reactivity and Spectroscopy in the Cyberinfrastructure for Astrochemical Organic Species (2015F59J3R, MIUR Ministero dell'Istruzione, dell'Università della Ricerca e della Scuola Normale Superiore), and (iii) the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme, for the Project "The Dawn of Organic Chemistry" (DOC), grant agreement No 741002.


Alexandra Houssaye
email: houssaye@mnhn.fr

Jocerand Prévoteau

What about limb long bone nutrient canal(s)? -A 3D investigation in mammals

Keywords: Mammal limb long bone nutrient canals Bone growth, femur, growth center, humerus, Mammalia, nutrient artery, nutrient foramen, shape

What about limb long bone nutrient canal(s)? -A 3D investigation in mammals

A. Houssaye, Jocerand

Introduction

Bone is a connective tissue whose matrix mineralizes. This hard tissue constitutes the skeleton, which provides support and protection and enables locomotion. Bones can notably be distinguished based on their position within the skeleton and their shape. Limb long bones occur in the stylopod (e.g. femur) and zeugopod (e.g. tibia) regions of the limbs. They are morphologically divided into a shaft bounded by both proximal and distal ends, called epiphyses (with intermediate metaphyses). The diaphysis of long bones is often (but not always) tubular, i.e. compact bone deposits surround an empty medullar cavity. Conversely, epiphyses are generally spongy [START_REF] Carter | Skeletal function and form: mechanobiology of skeletal development, aging, and regeneration[END_REF]. Long bones form through two distinct modes of ossification: periosteal, i.e., direct ossification in the periosteum (membrane covering the outer surface of the bone), and endochondral, i.e. following the resorption of a preformed cartilaginous structure [START_REF] Ross | Histology: a Text and Atlas[END_REF]. In long bones, after the formation of a cartilage model, periosteal bone deposits laterally and centrifugally. The cartilage inside the shaft starts calcifying and blood vessels invade the calcified cartilage. Then resorption occurs and a primary center of ossification develops, from which endochondral ossification starts proximally and distally centrifugally, with synchronous periosteal bone deposition (laterally) [START_REF] Carter | Skeletal function and form: mechanobiology of skeletal development, aging, and regeneration[END_REF]. The development of secondary centers of ossification in the epiphyses can develop, pending on the taxa (e.g. it occurs in mammals, lepidosaurs but not in turtles and crocodiles; [START_REF] Haines | Epiphyses and sesamoids[END_REF]. Nutrient foramina are openings through which blood vessels and peripheral nerves connect the marrow. If there can be epiphyseal and metaphyseal arteries, the nutrient arteries are located in the long bone diaphysis. The nutrient arteries are the major blood supply to long bones, especially during the early phases of growth and ossification [START_REF] Gümüsburun | A study of the nutrient foramina of lower limb long bones[END_REF]. The ossification center is considered to correspond to the intersection of the nutrient canal(s) with the central axis of the medullary area [START_REF] Digby | The measurement of diaphysial growth in proximal and distal directions[END_REF].

If growth can be fast symmetrical in the proximal and distal directions, and thus the growth center localized around mid-shaft, asymmetrical proximo-distal growth has been reported since the 18 th Century [START_REF] Ollier | Traité expérimental et clinique de la régénération des os et de la production artificielle du tissu osseux[END_REF][START_REF] Bisgard | Longitudinal growth of long bones[END_REF][START_REF] Payton | The growth in length of the long bones in the madder-fed pig[END_REF]; it is considered rather constant in position for a given bone [START_REF] Pereira | Nutrient foramina in the upper and lower limb long bones: morphometric study in bones of Southern Brazilian adults[END_REF] or not [START_REF] Campos | A study of the nutrient foramina in human long bones[END_REF]. Moreover several studies mentioned the occurrence of several nutrient canals [START_REF] Mysorekar | Diaphysial nutrient foramina in human long bones[END_REF][START_REF] Sendemir | Nutrient foramina in the shafts of lower limb long bones: situation and number[END_REF][START_REF] Gümüsburun | A study of the nutrient foramina of lower limb long bones[END_REF], or the absence of nutrient canal [START_REF] Longia | Study of diaphyseal nutrient foramina in human long bones[END_REF]; they did not see any significant relationship between this number and the bone length or the number of ossification centers [START_REF] Patake | Diaphysial nutrient foramina in human metacarpals and metatarsals[END_REF]. The number of nutrient canals seems to greatly vary intraspecifically for a given bone (Sendemir & Çimen, 1991;[START_REF] Pereira | Nutrient foramina in the upper and lower limb long bones: morphometric study in bones of Southern Brazilian adults[END_REF].

The nutrient arteries are supposed to enter the bone through a discrete foramen, the nutrient foramen, on the diaphysis. This nutrient is considered as originally indicating the position of the growth center; however, when the shaft extends asymmetrically in length, it becomes proximally or distally located relative to the growth center [START_REF] Hughes | The factors determining the direction of the canal for the nutrient artery in the long bones of mammals and birds[END_REF].

Considering the possible occurrence of several nutrient foramina and the asymmetrical proximo-distal growth, the relationship between nutrient foramen, nutrient canal, and growth center requires to be further investigated.

Beyond knowledge about growth mode, the nutrient canals are also supposed to bear physiological information. Some studies have indeed concluded that a greater blood pressure engenders a thickening and strengthening of the walls, as well as an increase in the circumference of vascular canals [START_REF] Seymour | Blood flow to long bones indicates activity metabolism in mammals, reptiles and dinosaurs[END_REF]. The nutrient canal could thus be larger in taxa showing a higher metabolic rate.

Nutrient arteries/foramen have essentially been analyzed in humans, and only to a much lesser extent in a few mammals (e.g., Ahn, 2013 on dogs). It has also been visualized in a few turtles [START_REF] Nakajima | Turtle humeral microanatomy and its relationship to lifestyle[END_REF]. Some studies have taken measurements of the nutrient foramen [START_REF] Seymour | Blood flow to long bones indicates activity metabolism in mammals, reptiles and dinosaurs[END_REF]; others have investigated the shape and orientation of the nutrient canals but only partially, for example using needles inserted through the nutrient foramen to get insights into the diameter (considering the diameter of the largest needle that could be inserted more than 1cm deep) and the orientation (penetration direction of the needle; [START_REF] Sim | Anatomy of the diaphyseal nutrient foramen in the long bones of the pectoral limb of German Shepherds[END_REF], but no study has analyzed the nutrient canal in three dimensions inside the bone.

The current study proposes to investigate the shape of the nutrient canals in stylopod bones of various mammals. It proposes the definition of qualitative and quantitative parameters in order to discuss the diversity in morphology, orientation, diameter… encountered, resorting to two different datasets in order to maximize differences within mammals and then analyze these variations within morphologically and phylogenetically closer taxa. This enables to discuss the growth mode of the stylopod bones in various mammals, the biological signal of some features of the nutrient canals, and the degree of interspecific and intraspecific variation observed.

Material and methods

Material

The material includes humeri and femora from two datasets (see Table 1) already available from previous studies (Houssaye et al. 2018;[START_REF] Houssaye | From land to water: evolutionary changes in long bone microanatomy of otters (Mammalia: Mustelidae)[END_REF], and representing adult specimens. The first one (hereafter referred to as QM) consists in 15 quadrupedal mammal species of various sizes, morphologies, locomotor modes, and they are widespread on mammal phylogeny (see Table 1). The aim of this sample is to show diversity, despite a limited sample size. The other sample (hereafter referred to as MU) is more restricted (in e.g., morphology, size and phylogenetic position). It focuses on mustelids, with 10 species (3 common with the previous sample) illustrating diverse degrees of adaptation to an aquatic lifestyle (terrestrial, semi-aquatic, aquatic).

Methods

Data acquisition

Bones were scanned using high-resolution computed tomography (GE_phoenix|X-ray v|tome|xs 240) at the Steinmann-Institut, University of Bonn (Germany) and at the AST-RX platform at the Museum national d'Histoire naturelle, Paris (France), with reconstructions performed using DATOX/RES software. Voxel size naturally varies pending on specimen size, from 19 to 246 µm. Image segmentation and visualization were performed from the reconstructed image data using Avizo 9.4 (VSG, Burlington, MA, USA).

The nutrient canals were segmented. The limits of the canals were difficult to clearly fix: at the outer surface, notably because of the obliquity of the canal; in the medullary area, because of bone resorption. Some canals were located more distally or proximally on the bone and not linked to the diaphyseal growth center; they thus appeared rather as metaphyseal or epiphyseal canals and were not taken into consideration.

Data analysis

The following parameters were defined in order to qualitatively and quantitatively characterize the nutrient canals: 1) number of nutrient canals in the bone (NNC); 2) mean diaphyseal radius (R), averaged for the whole diaphysis, based on the averaged sectional area (averaged CSA parameter in BoneJ), as √(averaged CSA/π); 3) the foraminal index (FI); this index was defined by [START_REF] Hughes | The factors determining the direction of the canal for the nutrient artery in the long bones of mammals and birds[END_REF] to estimate the relative position of the nutrient foramen on the bone; it is calculated as the ratio of the distance between the proximal extremity of the bone and the nutrient foramen (d) over the bone total length [START_REF] Kizilkanat | Location, number and clinical significance of nutrient foramina in human long bones[END_REF]. Because some epiphyses were incomplete (on the scans), the missing part was estimated based on comparisons with similar bones (e.g., on archeozootheque; https://www.archeozoo.org/archeozootheque/), and the points chosen to calculate the distances do not correspond to the most extremes points: in the humerus, the proximal point corresponds to the most distal point between the humeral head and the greater tubercle on the proximal surface; the distal point corresponds to the core of the trochlea; in the femur, the proximal point is the most distal point between the femoral head and the greater trochanter and the distal point is located at the core of the patella groove (Fig. 1); 4) canal diameter (CD) calculated based on the mid-canal sectional area (CA) obtained in Artec Studio Professional v12.1.5.1 software (Artec 3D, 2018) with the measuring tool "section", as √(CA/π). The parameter is calculated at the mid-length of what is preserved (not resorbed) of the canal walls, the section of the canal being usually rather constant except around its extremities; 5) canal diameter at the outer extremity (CDE), in order to take into consideration the variation of the canal diameter along its length, calculated as for CD but based on the most external section perpendicular to the canal; 6) canal obliquity (CO); this parameter had previously been analyzed based on the insertion of a needle in the canal and the measurement of the angle between the needle and the long bone longitudinal axis [START_REF] Sim | Anatomy of the diaphyseal nutrient foramen in the long bones of the pectoral limb of German Shepherds[END_REF]. Here since we have access to the 3D organization of the canal and sometimes observe orientation changes along the canal, we measured the general obliquity as the angle between a line linking the canal extremities as preserved and the horizontal to the bone longitudinal axis; 7) Canal cortical height (H) representing the relative height of the reconstructed canal (that could essentially be reconstructed in the cortex) divided by the bone length as measured for FI.

Parameters linked to canal length (and thus also volume) could not be obtained because of resorption occurring in the medullary area. Two additional parameters were defined in order to estimate the bone length and shape: L-bone length, as measured for FI; and LR-bone length (L) divided by the bone mean diaphyseal radius (R).

In addition to these quantitative parameters, some qualitative ones were also defined: We tested the phylogenetic signal of the continuous parameters in the data. For that we averaged the values obtained for all canals (when several canals occurred) and all specimens (when several specimens were available) for each species. Then we calculated the K-statistic following [START_REF] Blomberg | Testing for phylogenetic signal in comparative data: behavioral traits are more labile[END_REF] for each parameter and performed randomization tests. The K-statistic compares the observed phylogenetic signal in a trait with the signal under a Brownian motion model of trait evolution. A K-value > 1 implies more similarity between relatives than expected under Brownian motion; K < 1 highlights convergences. We used the phylogeny from [START_REF] Meredith | Impacts of the Cretaceous terrestrial revolution and KPg extinction on mammal diversification[END_REF] and from [START_REF] Slater | Integrating fossils with molecular phylogenies improves inference of trait evolution[END_REF] for the quadrupedal mammal (QM) and mustelid (MU) samples, respectively (Fig. S1).

We tested the influence of size on all parameters, first on all canals independently, then for each specimen, and thus grouping data for all canals for each specimen. In order to do so, we averaged all parameters but also calculated the sum of the canal diameters: CDsum and CDEsum. We performed linear regressions of each parameter to R.

In order to evaluate how these different parameters drive the variation observable in our samples, we conducted normalized PCAs [START_REF] David | Principal component analysis: a method for determining the essential dynamics of proteins, in: protein dynamics, methods in molecular biology[END_REF] for all datasets, those with all canals and those averaged by specimen.

These analyses and additional statistical ones (Student's t-test, Anova) were performed using statistical software R (R Core Team, 2014).

Results

The values obtained for the different parameters are listed in Tables S1-4.

Discrete parameters

The number of nutrient canals (NNC)

It varies from one to four in our samples. It varies among a single species (e.g. from 1 to 3 in the humerus of Mustela putorius) and between humeri and femora from the same species (1 in the humerus, 4 in the femur of Erinaceus europaeus). Correlation tests between the humerus and femur datasets averaged for each species indicate no correlation between the two bones.

Anovas revealed that the number of canals is always correlated with LR, and thus the bone proportions, and CDr (but not always with CDEr; Table 2). LR decreases when the bone is more robust for a given length. In less elongated bones, canals appear generally more

Direction

Canals are mainly directed distalo-proximally (from the core to the outer surface of the bone) in the humerus (Fig. 2A). There are a few exception, with proximo-distal orientations in Meles meles, one specimen of Mustela putorius, Pteronura brasiliensis and Lontra felina, and one of the two canals of the Mustela eversmanni specimen. From our sample, it thus only varies from the common condition in mustelids. In addition, a few canals are sub-horizontal, which usually corresponds to CO values lower than 15° (Fig. 2B). This occurs only in bones with several canals (except in the specimen of Lutra lutra). Conversely, in the femur, canals are mainly proximo-distally directed (Fig. 2C-D The orientation of the canals in the bone. It is extremely variable in the humeri and femora of the QM sample; they can extend either anteriorly, posteriorly, medially, laterally or in between. In the humerus of the mustelids, most canals are anteriorly directed, some medially and a few posteriorly directed. No canal is laterally oriented. In the femora, canals are almost always posteriorly oriented, only a few being medially oriented.

Continuous parameters

Phylogenetic signal

There is no phylogenetic signal in the parameters used to describe the canals for the QM dataset (Table 3). However, within mustelids, a phylogenetic signal is observed in some parameters (sometimes with a high K value) but not for the same parameters in the humerus as in the femur, except for CDE (Table 3). There is a significant phylogenetic signal for the size index R only for the Humerus QM (with no other parameter showing a significant

Influence of size

Linear regressions clearly show that the canal diameter is correlated with size, whereas its geometry is not, in all datasets (Table 4). For further analyses we used ratios for the canal diameter parameters, dividing their values by R (CDr,CDEr, CDrsum, CDErsum).

Values of the parameters

The continuous parameters are discussed based on datasets considering all canals (Fig. 3). The foraminal index (FI) is clearly higher in humeri than femora (t.tests; p<0.001; Fig. 3A). It is also slightly lower in the QM sample than in the MU one (t.tests; p<0.001; Fig. 3A). For the canal diameter parameters, no significant difference is observed between the various datasets (Fig. 3B-C). The canal obliquity (CO) and the canal cortical height (H) are much higher in otter femora (t.tests; p<0.001; Fig. 3D-E).

Principal component analyses

The PCAs enable not only to study the distribution of the different specimens in the morphospace, but also and above all in our study, to analyze the relative contribution of each parameter to the others and to compare results obtained based on the different datasets.

Quadrupedal mammals

The first two axes of the humerus PCA including all canals (representing 87.1% of the total variance; Fig. 4) show that several canals from the same specimen can be either close (as for Talpa; Te) or very distant (as for Choeropsis; Cl). The intrabone variation seems thus potentially very high. In this small sample, the intraspecific variation seems essentially driven by the canal geometry (H, CO, FI). About the variables, the two canal diameter parameters co-vary. H and FI vary antagonistically. When data are averaged for each specimen, the relationships between the variables remain similar.

The first two axes of the femur PCA including all canals (representing 78.6% of the total variance; Fig. 5A) also show a strong covariation between the two canal diameter parameters. The obliquity (CO) and canal cortical height (H) also strongly co-vary. As for the humerus, the intrabone variation varies between taxa (rather weak in Talpa [Te], high in CDErsum for the humerus when data are averaged for each specimen, this is not the case for the femur, where the first ones drive variability essentially along the first axis and the others along the second axis (Fig. 5B). The relationships between the variables are rather consistent between the two femur PCAs, with the averaged CDr and CDEr replacing the CDr and CDEr values. The addition of the CDrsum and CDErsum does not change much the relative distribution of the specimens in the morphospace.

Mustelids

The first two axes of the PCA for all canals of the humeri (representing 82.3% of the variance; Fig. 6A) clearly show that variation within a bone can be greater than between specimens (e.g. in Enhydra lutris [El], Lontra canadensis [Lc]). In this taxonomically more narrow sample, the relationships between the different variables is similar as in the QM sample, with a slightly lower co-variation between CDr and CDEr. When data are grouped by specimens, the CDEr and CDErsum variables separate from the CDr and CDrsum ones that strongly co-vary (Fig. 6B). The distribution of the specimens in the morphospace is this time different between the two samples (e.g.

Meles [Mm]).
The first two axes of PCA for all canals of the femora (representing 90.7% of the variance) tend to separate otters (lutrines) from other mustelids, with a limited overlap. H and FI strongly co-vary, slightly less with CO, whereas, at about right angle to those variables, CDr and CDEr also strongly co-vary. As for the humerus, variation within a bone can be higher than between bones of a same species and even between species (Fig. 7). When canals from single specimens are globally considered, all canal diameter variables co-vary and with the same relationship to the other variables as CDr and CDEr in the previous analysis. The distribution of the specimens is also similar.

Considering the divergence between lutrines and mustelines in the mustelid sample of femora and based on the strong differences in femur shape and proportions between these two groups [START_REF] Botton-Divet | Morphological analysis of long bones in semi-aquatic mustelids and their terrestrial relatives[END_REF], we analyzed the possible correlation between bone proportions and the various parameters, assuming a possible link with the canal geometry parameters. For that we performed linear regressions of all datasets with L and LR (Table 5).

These analyses reveal that there is no correlation between bone length (L) and the canal parameters (only two significant results for distinct parameters and datasets), but also no general correlation between bone proportions, i.e. notably between small and large bones versus elongated ones, and the canal parameters. Indeed only CDr shows a significant link but only for three of the four datasets. However, it is clear that the link is stronger for otter femora, all parameters except CO showing a significant correlation with LR in both datasets.

The canal features (except CO) are thus strongly linked with the bone proportions in the mustelid femora.

Correlation tests between the humerus and femur datasets averaged for each species indicate a significant correlation for the canal diameter parameters but not for the others.

Discussion

Our study was the first one, to our knowledge, to analyze in three dimensions both qualitatively and quantitatively the nutrient canals in limb bones thanks to the use of microtomography. We proposed various parameters to describe their features, and analyzed their values in two mammal samples in order to discuss their biological significance.

General comparisons

We did not observe a general phylogenetic signal in the quantitative parameters used to describe the canals. Indeed, a phylogenetic signal was observed for some parameters in the MU sample, but not for the same parameters in the two bones (except CDE). Moreover, no parameter was showing a significant phylogenetic signal in the QM sample. The occurrence of a phylogenetic signal seems thus very dependent on the sample. In all our analyses, intraspecific variation and even intrabone variation was often very high and sometimes much higher than interspecific variation. As a consequence, the nutrient canals appear to be poor indicators of biological features of a taxon or even organism. Nevertheless, some specific features could be identified, like the correlation between the canal features and the bone proportions in the MU femora (and not in other samples) notably distinguishing between otters and mustelines. This raises the potential interest in analyzing in detail the biological signal of nutrient foramina in taxonomically smaller taxa. The parameters CD and CDE were co-varying in all analyses, so that we would suggest the use of only one of these parameters in further studies. The other parameters analyzed appear suitable to describe the diversity of the canal shape. This study showed in limb bones a number of nutrient canals varying from one to four, but more generally between one and two, in its samples. Previous studies on human long bones described a majority of humeri with a single nutrient foramen and of femora with one or two nutrient foramina, some also showed up to nine et even no foramen at all [START_REF] Longia | Study of diaphyseal nutrient foramina in human long bones[END_REF][START_REF] Sendemir | Nutrient foramina in the shafts of lower limb long bones: situation and number[END_REF][START_REF] Nagel | The clinical significance of the nutrient artery[END_REF][START_REF] Gümüsburun | A study of the nutrient foramina of lower limb long bones[END_REF][START_REF] Kizilkanat | Location, number and clinical significance of nutrient foramina in human long bones[END_REF]).

Number and diameter of the nutrient canals

During our study, we could observe that some metaphyseal canals exit the bone in the diaphysis, so that some nutrient foramina on the diaphysis can belong to metaphyseal rather than diaphyseal nutrient canals.

In our samples, the number of canals varies also within a single species (e.g. from 1 to 3 in the humerus of Mustela putorius) and between humeri and femora from the same species. This parameter appears thus strongly variable intraspecifically. However, our study reveals a trend for nutrient canals to be generally more numerous and relatively thinner in less elongated bones. As for the canal diameter, the diameters taken at mid-length or at the periphery co-vary. These parameters are also the only one to co-vary between the humeri and the femora. This shows that the diameter of the canals is rather constrained. [START_REF] Seymour | Blood flow to long bones indicates activity metabolism in mammals, reptiles and dinosaurs[END_REF] suggested a link between blood vessel circumference and metabolic rates (larger blood vessels to service higher flow rates) and thus that the diameter of the nutrient canals would be linked to the physiology of the organisms. Our results confirm that the size of the canals appears constrained and rather homogeneous within a specimen.

Canal direction and orientation

In humeri, canals are essentially distalo-proximally directed, whereas it is the opposite in the femur. However foraminal indices are much higher in humeri than in femora, in accordance with studies on human bones (Pereire et al. 2011). It is because the growth center is located much more distally in the humerus than in the femur (in our samples). The nutrient canal seems thus generally oriented towards the mid-diaphysis. Unfortunately, because the centralmost part of the bones generally underwent resorption, the medialmost extension of the nutrient canals is not observable. As a consequence, the position of the growth center cannot precisely be ascertained. If in the case of straight canals, it would be tempting to extend the canal up to the core of the shaft, the occurrence of some bending in a few canals or of a slight change of orientation close to the growth center (e.g. femur of Enhydra) prevents the use of such an approach. When several canals occur, the intersection of the canals can also be ascertained as the position of the growth center. Despite the difficulty in determining the exact position of the growth center, the position of the medialmost extension of the canal clearly shows that the growth center is much distal (at about 70% of the bone length for both samples) in the humerus, and much proximal (at about 30% for QM and 40% for MU) in the femur. Growth is thus clearly not symmetrical in most taxa of our sample, being much more intense proximally in the humerus and distally in the femur. This confirms previous observations made on humans (e.g., [START_REF] Ollier | Traité expérimental et clinique de la régénération des os et de la production artificielle du tissu osseux[END_REF] but also preliminary observations made on fossil whales [START_REF] Houssaye | Transition of Eocene whales from land to sea: evidence from bone microstructure[END_REF]. However, a proximal growth center has been observed in the humerus of turtles [START_REF] Nakajima | Turtle humeral microanatomy and its relationship to lifestyle[END_REF]. If our result appears thus rather general in mammals, in cannot be generalized outside of Mammalia and further studies are required to investigate this question in other amniotes.

Canals are generally straight with only a few exceptions occurring. The relative extension of these canals along the diaphyseal cortex (H) is limited (<15% L with average values between 3 and 9% L). However this corresponds to minimal estimates since the position of the growth center cannot be ascertained precisely. Moreover, a few taxa display high values (e.g. 26% L in Choeropsis humerus, 15% in Rangifer and Lontra felina femur).

Our study also clearly reveals that horizontal or subhorizontal canals are rare. The nutrient foramen thus does not indicate the location of the growth center. It has been assumed that the nutrient foramen was originally located at the level of the growth center but that an asymmetrical growth could engender the foramen to move proximally or distally, towards the faster growing end [START_REF] Hughes | The factors determining the direction of the canal for the nutrient artery in the long bones of mammals and birds[END_REF][START_REF] Henderson | The position of the nutrient foramen in the growing tibia and femur of the rat[END_REF]. Our results confirm that the canals are essentially oriented towards the faster growing end.

There is a high variation in the orientation of the canals in the bone within mammals.

However, it was interesting to note preferential or absent orientations in the mustelid sample, that probably reflect morphological constraints during development that would be specific to these taxa, and thus that would imply stronger constraints at smaller taxonomical scales.

Similarly, the much higher canal obliquity and cortical height in mustelid femora than in the other samples could also reflect specific developmental constraints. As a consequence, if the biological information could be flooded in a large sample of diverse taxa, more specific patterns might be revealed at a smaller taxonomic scale. 3 and Table 1. A-all canals; B-all specimens. Abbreviations as in Fig. 3 and Table 1. Figure 6. Distribution of the mustelid specimens in the morphospace and contribution of the various parameters along the two first axes of the humerus PCA. A-all canals; B-all specimens. Yellow: Mustelinae; Blue: Lutrinae. Abbreviations as in Fig. 3 and Table 1. Mustelinae; Blue: Lutrinae. Abbreviations as in Fig. 3 Abbreviations as in Fig. 3 and Table 1. Blue: Lutrinae. Abbreviations as in Fig. 3 and Table 1. Abbreviations as in Fig. 3 and Table 1.

  a) direction: proximal to distal (PD) or distal to proximal (DP); b) shape; it defines the global shape of the canal: straight or angled; c) orientation: it describes the anatomical orientation of the canal in the bone.

  thinner. But, as shown by the intraspecific variation, this is a common but not absolute trend.

  ). Exceptions occur in Rangifer tarandus, Vulpes vulpes, and one of the canals of Erinaceus europaeus, Rupicapra rupicapra, and Mustela eversmanni. One horizontal canal is observed in one of the canals of Choeropsis liberiensis. Canals are essentially straight, except in Rangifer's humerus and femur and in several bones of mustelids where they are angled: the humerus of Meles, Martes, two M. putorius, one canal of M. eversmanni, one canal of one Enhydra lutris. Moreover, a humerus of Lontra felina shows a serpentiform canal. Conversely all canals in the mustelid femora are straight.

  this sample). The occurrence of a phylogenetic signal in the data appears thus clearly independent of size.

  Erinaceus[Ee] and Rupicapra [Rr]) and appears clearly affected by various parameters specimens. Whereas CDr and CDEr rather co-vary with CDrsum and

  , Bonn, Germany), S. Merker (Staatliches Museum für Naturkunde Stuttgart, Stuttgart, Germany), and L.E. Olson and A. Gunderson (University of Alaska Museum, Fairbanks, AK, USA) for the loan of the specimens. We thank the Steinmann-Institut (University of Bonn, Germany) for providing beamtime and support, and M. Garcia Sanz for performing scans and reconstructions at the AST-RX platform (UMS 2700, MNHN). We also thank A. Delapré for technical advices. A.H. acknowledges financial support from the A. v. Humboldt Foundation, the ANR-13-PDOC-0011, and the ERC-2016-STG-715300.

F

Figure 1 .

 1 Figure 1. 3D reconstructions showing how the foraminal index (FI) was calculated for the nutrient canal (in purple), based on bone length (L) and the distance between the proximal extremity and the nutrient foramen (d). A-Humerus of Rangifer tarandus STIPB M47. B-Femur of Mustela eversmannii MNHN 2005-668.

Figure 2 .

 2 Figure 2. 3D reconstructions showing the number and orientation of the nutrient canals. A-Two canals distalo-proximally oriented in the humerus of Dama dama STIPB M1; B-four horizontal canals in the humerus of Talpa europaea STIPB Unnumbered; C-two canals proximo-distally oriented and in the femur of Lontra felina MNHN 1884-874; and D-a single canal proximo-distally oriented in the femur of Felis silvestris UFGK Unnumbered.

Figure 3 .

 3 Figure 3. Boxplots illustrating the variation of the various parameters pending on the samples. FI: foraminal index; CDr: relative canal diameter; CDEr: relative canal diameter at the outer extremity; CO: canal obliquity; H: canal cortical height; R: bone mean diaphyseal radius.

Figure 4 .

 4 Figure 4. Distribution of the specimens in the morphospace and contribution of the various parameters along the two first axes of the humerus PCA for the quadrupedal mammal (all canals) sample. Abbreviations as in Fig. 3 andTable 1.

Figure 5 .

 5 Figure 5. Distribution of the specimens in the morphospace and contribution of the various parameters along the two first axes of the femur PCA for the quadrupedal mammal sample.

Figure 7 .

 7 Figure 7. Distribution of the mustelid specimens in the morphospace and contribution of the various parameters along the two first axes of the femur PCA on all canals. Yellow:

Figure 1 .Figure 2 .Figure 3 .

 123 Figure 1. 3D reconstructions showing how the foraminal index (FI) was calculated for the nutrient canal (in purple), based on bone length (L) and the distance between the proximal extremity and the nutrient foramen (d). A-Humerus of Rangifer tarandus STIPB M47. B-Femur of Mustela eversmannii MNHN 2005-668.

Figure 4 .

 4 Figure 4. Distribution of the specimens in the morphospace and contribution of the various parameters along the two first axes of the humerus PCA for the quadrupedal mammal (all canals) sample. Abbreviations as inFig.3and Table1.

Figure 5 .

 5 Figure 5. Distribution of the specimens in the morphospace and contribution of the various parameters along the two first axes of the femur PCA for the quadrupedal mammal sample. A-all canals; B-all specimens.Abbreviations as in Fig.3and Table1.

Figure 6 .

 6 Figure 6. Distribution of the mustelid specimens in the morphospace and contribution of the various parameters along the two first axes of the humerus PCA. A-all canals; B-all specimens. Yellow: Mustelinae;Blue: Lutrinae. Abbreviations as in Fig.3and Table1.

Figure 7 .

 7 Figure 7. Distribution of the mustelid specimens in the morphospace and contribution of the various parameters along the two first axes of the femur PCA on all canals. Yellow: Mustelinae; Blue: Lutrinae.Abbreviations as in Fig.3and Table1.

Table 1 .

 1 List of the material analyzed in this study. QM: dataset of the various quadrupedal mammal species; MU: dataset of the mustelids; H: humerus; F: femur. Abb: abbreviations,

	o r
	P
	e e r
	R e v
	i e w
	O n
	l
	y

Table 2 .

 2 P-values from the Anovas on the link between the number of nutrient canals (NNC) and the various continuous characters (for the datasets including all canals). Abbreviations as in Material and Methods. In bold when significant at 5%.

	D/P	R	FI	CDr	CDEr	CO	H	L	LR
	H QM p=0.39 p=0.47 p<0.001 p=0.008 p=0.004 p=0.68 p=0.08 p<0.001
	F QM p=0.77 p=0.98 p=0.02 p=0.70 p=0.07 p=0.20 p=0.31 p=0.01
	H MU p=0.02 p=0.44 p=0.04 p=0.10 p=0.85 p=0.67 p=0.009 p=0.03
	F MU p=0.14 p=0.43 p=0.003 p=0.02 p=0.83 p=0.27 p=0.49 p=0.04
	Table. 3. K-statistics with the associated p-values for the continuous parameters. In bold when
	significant at 5%.
	Dataset/Parameter	FI	CD		CDE	CO	H	R
	Humerus QM	K=0.44	K=0.65	K=0.45	K=0.32	K=0.46	K=1.40
			p=0.22	p=0.08	p=0.15	p=0.80	p=0.32	p<0.001
	Femur QM		K=0.46	K=0.42	K=0.40	K=0.41	K=0.36	K=0.48
			p=0.24	p=0.47	p=0.43	p=0.66	p=0.80	p=0.23
	Humerus MU	K=1.01	K=1.34	K=1.41	K=0.73	K=0.73	K=0.74
			p=0.049	p=0.003	p=0.002	p=0.08	p=0.11	p=0.10
	Femur MU		K=0.81	K=0.84	K=1.46	K=0.69	K=1.57	K=0.73
			p=0.07	p=0.10	p=0.03	p=0.10	p=0.005	p=0.13

Table

4

. Results of the linear regressions performed on the datasets (with R as a size estimate). AC: all canals; I: averaged by individual; In bold when significant at 5%.

Table 5 .

 5 Results of the linear regressions performed on the datasets with L and LR. In bold

	D/	FI		CDr		CDEr		CDrsum	CDErsum	CO		H
	P
	V	L	LR	L	LR	L	LR	L	LR	L	LR	L	LR	L	LR
	ar
	H	R=-	R=0.	R=-	R=0.	R=-	R=0.	-	-	-	-	R=0	R=0.	R=0	R=0.
	Q	0.20	15	0.16	72	0.06	45					.41	63	.36	08
	M	p=0.	p=0.	p=0.	p<0.	p=0.	p=0.					p=0.	p=0.	p=0.	p=0.
	A	40	53	50	001	79	047					08	003	12	73
	C
	F	R=0	R=-	R=-	R=0.	R=-	R=0.	-	-	-	-	R=-	R=0.	R=0	R=0.
	Q	.06	0.16	0.08	44	0.09	08					0.12	26	.27	06
	M	p=0.	p=0.	p=0.	p=0.	p=0.	p=0.					p=0.	p=0.	p=0.	p=0.
	A	79	46	71	03	69	70					56	22	20	78
	C
	H	R=0	R=0.	R=0	R=0.	R=0	R=0.	-	-	-	-	R=-	R=0.	R=-	R=-
	M	.19	20	.10	37	.07	31					0.05	13	0.07	0.21
	U	p=0.	p=0.	p=0.	p=0.	p=0.	p=0.					p=0.	p=0.	p=0.	p=0.
	A	40	37	65	09	75	17					83	56	75	36
	C
	F	R=-	R=-	R=-	R=0.	R=-	R=0.	-	-	-	-	R=-	R=-	R=-	R=-
	M	0.03	0.60	0.32	56	0.21	60					0.14	0.28	0.14	0.70
	U	p=0.	p<0.	p=0.	p=0.	p=0.	p<0.					p=0.	p=0.	p=0.	p<0.
	A	89	001	09	002	28	001					47	14	47	001
	C
	H	R=-	R=0.	R=-	R=0.	R=-	R=0.	R=-	R=0.	R=-	R=0.	R=0	R=0.	R=0	R=-
	Q	0.28	23	0.46	64	0.19	31	0.56	21	0.23	00	.18	46	.39	0.02
	M	p=0.	p=0.	p=0.	p=0.	p=0.	p=0.	p=0.	p=0.	p=0	p=1.	p=0.	p=0.	p=0.	p=0.
	I	33	43	10	01	52	27	04	48	.44	00	55	10	17	95
	F	R=0	R=-	R=-	R=0.	R=-	R=0.	R=-	R=-	R=-	R=-	R=0	R=0.	R=0	R=-
	Q	.25	0.24	0.12	39	0.12	06	0.23	0.06	0.19	0.32	.07	23	.57	0.15
	M	p=0.	p=0.	p=0.	p=0.	p=0.	p=0.	p=0.	p=0.	p=0	p=0.	p=0.	p=0.	p=0.	p=0.
	I	39	40	67	17	68	84	42	84	.52	27	82	43	03	61

when significant at 5%. D/P: Dataset/Parameter; Var: variable; other abbreviations as in Table

3

.

Table S1 .

 S1 and Table 1. Values of the parameters for the humerus QM dataset. Med: medial; Lat: lateral;

	566
	567	Ant: anterior; Post: posterior; DP: distalo-proximal; H: horizontal; PD: proximo-distal; S:
	568	straight; A: angled.
		Taxon	Direction Orientation Shape	(mm) NNC R	(mm) FI CD	(mm) CDE	(°C) CO	H
		Marmota Med	DP	S			1	3.93 0.67	0.64	1.40	61.7 0.02
		Tupaia	Lat	DP	S			1	1.70 0.75	0.21	0.25	37.6 0.01
		Erinaceus Post	DP	S			1	2.53 0.53	0.19	0.20	71.8 0.09
		Talpa	Post	DP	S			4	2.35 0.60	0.10	0.17	28.0 0.04
			Post	DP	S			4	2.35 0.60	0.04	0.04	27.6 0.04
			Post	H	S			4	2.35 0.62	0.06	0.11	3.0 0.01
		F o r Ant-lat H F Post DP Choeropsis Post Sus DP o r Med H	S S S S			4 1 12.49 0.57 2.35 0.68 2 16.32 0.34 2 16.32 0.75	0.09 0.91 1.02 0.81	0.09 2.51 1.67 0.47	3.1 0.03 60.5 0.09 69.5 0.26 13.1 0.01
		P e e r DP A P DP S Ant-med DP Ant Post S e Rupicapra Lat Rangifer Dama DP S e r Post DP S	1 14.06 0.74 2 13.42 0.61 2 13.42 0.69 2 8.76 0.66 2 8.76 0.43	0.95 0.65 0.31 0.43 0.56	2.63 0.98 0.46 0.28 0.53	34.0 0.01 63.0 0.05 50.0 0.03 70.8 0.06 80.1 0.11
	569 570	Felis Vulpes Meles Martes Mustela	Med Post Med Med Ant	DP DP PD DP PD	S S A A A	R e v i e w 1 4.90 0.67 1 5.19 0.52 R 1 7.02 0.73 e 1 2.72 0.76 v 1 2.43 0.76 i e w	0.51 0.44 0.77 0.44 0.30	0.64 0.46 0.54 0.75 0.43	64.7 0.06 80.8 0.13 49.2 0.03 33.7 0.01 57.8 0.03
	571									O n n O
										l l
										y y
										23

Table S2 .

 S2 Values of the parameters for the femur QM dataset. Abbreviations as in TableS1.

		Ant-med H	S	3 17.84 0.33	0.81	2.47	4.1 0.00
	Rangifer	Ant	DP	A	1 14.94 0.36	1.96	2.44	61.4 0.15
	Rupicapra Ant	PD	S	2	8.49 0.24	0.79	1.18	36.6 0.01
		Post	DP	S	2	8.49 0.61	0.51	0.41	88.3 0.13
	Felis	Post	PD	S	1	4.78 0.41	0.36	0.39	70.0 0.04
	Vulpes	Post	DP	S	1	5.45 0.30	0.56	0.60	62.5 0.03
	Meles	Post	PD	S	2	6.94 0.52	0.37	0.67	68.1 0.08
		Med	PD	S	2	6.94 0.41	0.32	0.60	41.0 0.02
	Martes	Post-lat PD	S	2	2.91 0.37	0.36	0.61	46.3 0.04
		Post-lat PD	S	2	2.91 0.37	0.18	0.41	44.2 0.03
	Mustela	Post	PD-H	S	1	2.65 0.34	0.21	0.14	18.8 0.01
	572 573 574	F o r
			P
			e e r
				R e v
					i e w
	Taxon Cavia	Direction Orientation Shape Post-med PD S	NNC R (mm) 1 2.74 0.30 FI CD (mm) 0.23 O n l CDE (mm) 0.39	CO (°C) 52.9 0.05 H
	Marmota Med Tupaia Med	PD PD	S S	1 1	4.00 0.33 1.27 0.29	0.35 0.23 y	0.87 0.23	70.4 0.06 63.9 0.02
	Erinaceus Ant-med PD	S	4	2.32 0.38	0.26	0.73	51.9 0.06
		Med-ant PD	S	4	2.32 0.32	0.12	0.19	44.8 0.05
		Ant	PD	S	4	2.32 0.30	0.11	0.07	36.8 0.02
		Ant	DP	S	4	2.32 0.26	0.09	0.08	24.0 0.01
	Talpa	Lat	PD	S	3	1.13 0.53	0.06	0.21	64.2 0.05
		Med-post	PD	S	3	1.13 0.54	0.08	0.14	68.4 0.04
		Ant	PD	S	3	1.13 0.52	0.06	0.09	70.0 0.03
	Sus	Post	PD	S	1 11.37 0.56	0.66	1.28	80.7 0.08
	Choeropsis Med-ant PD	S	3 17.84 0.42	0.72	1.72	47.1 0.04
		Ant	PD	S	3 17.84 0.36	1.28	2.16	21.0 0.01

Table S3 .

 S3 Values of the parameters for the humerus MU dataset. Abbreviations as in Table

	578
			F
			o r
	575	S1. SE: serpentiform. Taxon Direction Orientation Shape P Meles Med PD A Martes Med DP A e e r	NNC R (mm) 1 7.02 0.73 FI CD (mm) 0.77 1 2.72 0.76 0.44	CDE (mm) 0.54 49.2 0.03 CO (°C) H 0.75 33.7 0.01
		Mustela putorius Mustela eversmanni	Ant Post Ant Ant Post Ant Ant	PD DP DP H DP PD DP	A S S S R e 1 3 3 3 v A 1 A 2 S 2 i 2.43 0.76 0.96 0.59 0.96 0.73 0.96 0.72 1.22 0.62 1.15 0.75 1.15 0.71 e w	0.30 0.09 0.08 0.04 0.18 0.14 0.10		0.43 57.8 0.03 0.09 73.6 0.08 0.13 31.1 0.01 0.03 10.7 0.01 0.13 55.6 0.05 0.23 25.4 0.02 0.13 60.6 0.02
		Neovison Pteronura	Ant Ant Ant Ant	DP DP PD DP	S S A A	2 2 1 1	1.15 0.69 1.15 0.71 3.36 0.75 3.18 0.71 O n 0.11 0.15 0.25 0.23 l		0.14 45.8 0.03 0.16 34.8 0.01 0.31 17.1 0.02 0.27 47.5 0.01
		Lontra felina	Ant	PD	SE	1	2.04 0.69	0.23 y	0.20 53.1 0.06
			Med	DP	S	1	2.58 0.62	0.22		0.49 41.4 0.05
		Lontra canadensis	Ant-med DP	S	1	3.03 0.65	0.39		0.40 61.6 0.04
			Med-post	DP	S	2	2.81 0.46	0.09		0.09 49.6 0.14
			Med	DP	S	2	2.81 0.7	0.27		0.33 48.1 0.03
		Lutra	Med	H	S	1	2.78 0.73	0.28		0.45	7.9 0.00
		Enhydra	Ant	DP	S	1	3.67 0.72	0.34		0.63 22.5 0.02
			Ant	H	A	2	3.92 0.72	0.20		0.38	3.7 0.02
			Ant	DP	S	2	3.92 0.65	0.33		0.37 48.0 0.04
	576
	577
										25

Table S4 .

 S4 Values of the parameters for the femur MU dataset. Abbreviations as in TableS1.

	Taxon	Direction Orientation Shape	(mm) NNC R	(mm) FI CD	(mm) CDE	(°C) CO	H
	Meles	Post	PD	S	2	6.94 0.52	0.37	0.67 68.1 0.08
		Med	PD	S	2	6.94 0.41	0.32	0.60 41.0 0.02
	Martes	Post-lat PD	S	2	2.91 0.37	0.36	0.61 46.3 0.04
		Post-lat PD	S	2	2.91 0.37	0.18	0.41 44.2 0.03
	Mustela putorius	Post	PD-H	S	1	2.65 0.34	0.21	0.14 18.8 0.01
		Post-med PD	S	1	1.23 0.36	0.29	0.42 21.3 0.01
	Mustela eversmanni	Post Post-lat DP PD Post PD F o r	S S S	1 2 2	0.97 0.39 1.10 0.33 1.10 0.4	0.12 0.12 0.15	0.26 67.3 0.03 0.17 19.2 0.01 0.24 46.5 0.02
	Neovison	Post Post Post	PD PD PD P e S S S	1 1 1	1.12 0.4 1.09 0.45 1.11 0.47	0.26 0.11 0.09	0.49 66.3 0.03 0.11 71.0 0.06 0.11 75.2 0.08
	Pteronura	Post Post-med PD PD	S S e r	1 1	2.83 0.48 3.65 0.56	0.32 0.37	0.32 67.5 0.09 0.38 77.5 0.15
	Lontra felina	Post Post Post	PD PD PD	S S S R e 1 1 1	3.3 0.56 2.29 0.53 2.28 0.48	0.31 0.25 0.22	0.32 76.2 0.15 0.29 77.7 0.14 0.54 69.0 0.12
	Lontra canadensis	Post Med Post Post	PD PD PD PD	S S S S	2 2 v 1.89 0.5 1.89 0.55 2 2.76 0.52 i e 2 2.76 0.45 w	0.17 0.05 0.22 0.18	0.19 71.4 0.13 0.04 77.7 0.17 0.37 73.7 0.21 0.15 60.4 0.07
	Lutra Enhydra	Post-med PD Med PD Post PD Post PD Post-med PD Post PD	S S S S S S	3 3 3 1 2 2	2.64 0.58 2.64 0.54 2.64 0.43 O 0.16 0.16 0.08 2.66 0.45 0.24 4.42 0.42 0.40 n l 4.42 0.57 0.11 y 0.18 75.3 0.18 0.18 75.0 0.14 0.13 56.4 0.06 0.27 69.6 0.09 0.46 67.8 0.12 0.11 70.6 0.17
		Post	PD	S	2	4.48 0.58	0.30	0.31 54.4 0.21
		Med-post	PD	S	2	4.48 0.35	0.23	0.33 20.9 0.05
	579

Journal of Anatomy

Acknowledgments

We warmly thank G. Véron (Museum National d'Histoire Naturelle, Paris, France), the Steinmann-Institut (University of Bonn, Bonn, Germany), R. Hutterer (Zoologische

Supporting Information