Romane Le Gal

Karin I Öberg

Bertrand Lefloch

Charlotte Vastel

Cécile Favre

Eleonora Bianchi

Typeset using L A T E X twocolumn style in AASTeX63

Keywords: YSOs: individual (IRAS 04181+2654A, IRAS 04169+2702, IRAS 04166+2706, IRAS 04016+2610, IRAS 04365+2535, IRAS 04295+2251, IRAS 04302+2247) -

There is mounting evidence that the composition and structure of planetary systems are intimately linked to their birth environments. During the past decade, several spectral surveys probed the chemistry of the earliest stages of star formation and of late planet-forming disks. However, very little is known about the chemistry of intermediate protostellar stages, i.e. Class I Young Stellar Objects (YSOs), where planet formation may have already begun. We present here the first results of a 3mm spectral survey performed with the IRAM-30m telescope to investigate the chemistry of a sample of seven Class I YSOs located in the Taurus star-forming region. These sources were selected to embrace the wide diversity identified for low-mass protostellar envelope and disk systems. We present detections and upper limits of thirteen small (N atoms ≤ 3) C, N, O, and S carriers -namely CO, HCO + , HCN, HNC, CN, N 2 H + , C 2 H, CS, SO, HCS + , C 2 S, SO 2 , OCS -and some of their D, 13 C, 15 N, 18 O, 17 O, and 34 S isotopologues. Together, these species provide constraints on gas-phase C/N/O ratios, D-and 15 N-fractionation, source temperature and UV exposure, as well as the overall S-chemistry. We find substantial evidence of chemical differentiation among our source sample, some of which can be traced back to Class I physical parameters, such as the disk-to-envelope mass ratio (proxy for Class I evolutionary stage), the source luminosity, and the UV-field strength. Overall, these first results allow us to start investigating the astrochemistry of Class I objects, however, interferometric observations are needed to differentiate envelope versus disk chemistry.

INTRODUCTION

Planets have long been supposed to form by material accretion in mature planet-forming disks, a.k.a. Class II disks, i.e. when disk parental envelopes have been cleared away (Lada & Wilking 1984;Andrews & Williams 2005;[START_REF] Dullemond | Protostars and Planets V[END_REF]). This evolutionary stage was thus identified as the birthplace of planets, containing both the chemical and mass reservoirs for planet formation. However, during the past decade, evidence started to challenge this assumption.

First, Class II disk masses are found to be, on average, well below the so-called minimum mass solar nebula (MMSN, Hayashi 1981;Weidenschilling 1977) of 0.01-0.1 M required to form planetary systems such as our own, and more generally below the typical masses of exoplanetary systems (e.g., Barenfeld et al. 2016;Ansdell et al. 2016;Law et al. 2017;Tazzari et al. 2017;Long et al. 2017;Eisner et al. 2018;Cazzoletti et al. 2019). Although Class II disk gas and dust masses could be underestimated due to, for example, poorly-understood CO depletion and/or optical depth effects (Yu et al. 2017a,b;Zhang et al. 2019;Andrews 2020), if Class II disk masses are effectively well below the MMSN and, more generally, below the typical masses of exoplane-tary systems, planet formation would need either 1) to start at earlier evolutionary stages, i.e. in the Class I stage, where the circumstellar envelope starts to fade away as matter accretes onto the nascent protoplanetary disk and the central protostar (e.g., Harsono et al. 2018) ; or 2) to still be fed by a regular replenishment of the protoplanetary disk from the parental molecular cloud (e.g., Kuffmeier et al. 2017;Manara et al. 2018;[START_REF] Dullemond | Protostars and Planets V[END_REF]Kuffmeier et al. 2020).

Second, recent ALMA observations revealed that most Class II disks and even a few Class I disks (e.g., HL Tau, L1551 NE, and GY 91, Yen et al. 2017;Takakuwa et al. 2017;Sheehan & Eisner 2018) possess sub-structures such as rings, clumps, spirals and large-scale crescents (e.g., ALMA Partnership et al. 2015;Andrews et al. 2018;Huang et al. 2018;Long et al. 2018). These substructures are somehow puzzling since they are at the same time supposed to (i) trace planets in the making (i.e., by-products of planet formation) and (ii) facilitate/initiate planet formation. Even though their primary origin remains to be elucidated, these substructures likely indicate that planet formation starts at earlier stages of disk evolution.

Thus, to get the full picture of planet genesis and chemical origins of planet compositions, it is crucial to study earlier evolutionary stages starting with Class I Young Stellar Objects (YSOs), which are still embedded in their parental envelope (e.g., [START_REF] White | Protostars and Planets V[END_REF]), i.e. protostellar envelope and disk systems. In the Class I phase, the disk chemistry is likely influenced both by the infall of fresh envelope material (e.g., Yen et al. 2013;Sakai et al. 2014;Jacobsen et al. 2019) and by the increasing importance of UV and X-ray emissions from the star (e.g., Feigelson & Montmerle 1999;Igea & Glassgold 1999;Bergin et al. 2003;Jonkheid et al. 2004;[START_REF] Gorti | [END_REF]Woitke et al. 2010Woitke et al. , 2016;;Dutrey et al. 2014). It is unclear how much of the resulting chemical composition is inherited from their natal molecular cloud versus reset by chemical reprocessing. To address this question, we carried out a spectral line survey of a small sample of Class I systems to explore the chemical nature of these objects.

Here, we present IRAM-30m observations of a sample of seven low-mass protostellar systems located in the Taurus star-forming region and selected to encompass a wide range of Class I physical properties. We aim to explore the chemical content of Class I YSOs and the relationship between molecular emission patterns and the physical properties of these objects. The source sample and observations are described in Section 2. Our main findings are reported in Section 3 where, after presenting the different line detections ( § 3.1) and line flux variation across the source sample ( § 3.2), we focus on line flux ratios proposed to probe C/N/O gas-phase ratios, deuterium fractionation, temperature, UV-field strength, and chemical evolution ( § 3.3). In Section 4, we discuss the chemical differentiation found across our sample, and how it compares to other Class I studies, and to Class 0 and Class II chemistry. We present our conclusions in Section 5.

OBSERVATIONS

Source Sample

We selected seven Class I sources with wellcharacterized physical properties (e.g., Sheehan & Eisner 2017), located in the nearby star-forming region of Taurus (d < 200 pc, Galli et al. 2019, see Table 1). The targeted sources were chosen to maximize the diversity of physical properties (disk-to-envelope ratio, luminosity and disk size) in order to elucidate their relationship with the disk chemical properties. The coordinates and properties of each source are listed in Table 1; the sources span luminosities ranging from 0.3 to 7 L , envelope masses from 0.02 to 1.2 M and masses of their embedded nascent disks from 0.01 to 0.1 M . We provide a brief description of our selected sources in the following.

IRAS 04302+2247

IRAS 04302+2247 (hereafter I-04302), the so-called Butterfly star (Lucas & Roche 1997), is a well-known YSO with an edge-on circumstellar disk located near the L1536b dark cloud (e.g., Wolf et al. 2003;Connelley & Greene 2010). Of the sources in our sample, this one has the largest disk-to-envelope mass ratio (≈ 6.7), presenting a massive disk, with M disk ≈ 0.114M and a 240 au radius, and is close to fully dispersing its natal envelope (M envelope ≈ 0.017M , Sheehan & Eisner 2017). This indicates the advanced evolutionary stage of this YSO with respect to the six other sources of our sample (see Table 1).

IRAS 04295+2251

IRAS 04295+2251 (hereafter I-04295), also known as L1536 IRS, is classified as a Class I object according to its bolometric temperature of T bol = 270 K (Young et al. 2003). Although Chiang & Goldreich (1999) suggested that it could be an edge-on disk, later on observations rather suggest an inclined disk (Furlan et al. 2008;Eisner 2012;Sheehan & Eisner 2017). Although its disk-toenvelope mass ratio is ≈ 0.5, one order of magnitude lower than for I-04302, it is, however, the second largest of our sample. Thus, contrary to I-04302, this indicates that the disk is still well embedded in its natal envelope, with M disk ≈ 0.02M and a 130 au radius. (a) Dec. (a) T bol e From Beck (2007) f From Galli et al. (2019), using the distance of the molecular clouds where the sources are located (see Section 2.1).

L (b) M (b) Env. M (b) Disk M (b) Disk / R (b) Env. R (b) Disk VLSR Dist. (J2000) (J2000) (K) (L ) (M ) (M ) MEnv. (
g Identified as an outlier in Galli et al. (2019), so we kept here the former distance used for the Taurus star-forming regions.

Note-[1] Deduced from the rotational transition N = 1 -0 of HCN from this work, [2] Deduced from the rotational transition N = 1 -0 of CN from this work, [3] from Santiago-García et al. (2009) and in agreement with the one deduced from the rotational transition N = 1 -0 of HCN from this work.

IRAS 04365+2535

IRAS 04365+2535 (hereafter I-04365), also known as TMC1A, is located at the center of the L1534 dark cloud (Benson & Myers 1989). I-04365 shows an inclined Keplerian rotating disk (e.g., Harsono et al. 2014;Aso et al. 2015;Sheehan & Eisner 2017), with a disk-to-envelope mass ratio of ≈ 0.4, third largest of our sample, with M disk ≈ 0.04M and a 140 au disk radius. The disk is driving a wind (Bjerkeli et al. 2016) and large-scale outflow [START_REF] Ceccarelli | Protostars and Planets VI[END_REF]Aso et al. 2015) probably both at the origin of a forming cavity in the surrounding envelope.

IRAS 04016+2610

IRAS 04016+2610 (hereafter I-04016), also known as L1489 IRS, is a Class I system (T bol = 226 K, Green et al. 2013) likely located behind a foreground starless cloud (e.g., Hogerheijde & Sandell 2000;Brinch et al. 2007a). This source shows evidences of a disk in Keplerian rotation misaligned from its surrounding envelope (Hogerheijde 2001;Brinch et al. 2007b) and a bipolar outflow (Myers et al. 1988;Hogerheijde et al. 1998;Yen et al. 2014). Although this source has the largest disk of our sample, it also has one of the smallest disk masses (M disk = 0.009M , see Table 1), with a disk-to-envelope mass ratio ≈ 0.4 indicative of the embedded nature of this source.

IRAS 04166+2706

IRAS 04166+2706 (hereafter I-04166) is a YSO located in the molecular filament of the L1495/B213 Taurus (Hacar et al. 2013), harbouring highly collimated jet and bipolar outflow (Bontemps et al. 1996;Tafalla et al. 2004;Santiago-García et al. 2009). Claimed to be a Class 0 system by some (e.g., Tafalla et al. 2004Tafalla et al. , 2017;;Wang et al. 2014Wang et al. , 2019) ) and a Class I source by others (e.g., Young et al. 2003;Furlan et al. 2008;Eisner 2012;Sheehan & Eisner 2017), it is probably in transition between these two stages. This is in agreement with the fact that I-04166 is the source with the lowest luminosity (≈ 0.2L ) and bolometric temperature (T bol = 75 K) in our sample (see Table 1).

IRAS 04169+2702

IRAS 04169+2702 (hereafter I-04169), located in the molecular filament of the L1495/B213 Taurus (Hacar et al. 2013), is a Class I YSO (T bol =133 K, Young et al. 2003) also characterized by a bipolar outflow (Bontemps et al. 1996) found to be perpendicular to an elongated flattened envelope of about 2200 au × 1100 au in size area (Ohashi et al. 1997). A compact ∼ 40 au radius disk is embedded in this huge envelope (Sheehan & Eisner 2017), with a disk-to-envelope mass ratio ≈ 0.2.

IRAS 04181+2654A

IRAS 04181+2654A (hereafter I-04181), located in the molecular filament of the L1495/B213 Taurus (Bontemps et al. 1996;Davis et al. 2010), is the least characterized source of our sample due to its lack of submillimeter study. It is part of a wide binary system, with a deeper embedded B component lying 31" away (i.e., ∼ 4300 au away at the distance of Taurus, Furlan et al. 2008). The A component is closer to a dense core and shows two bipolar outflows (Davis et al. 2010).

Description of the observations

The observations were performed at the IRAM-30m telescope in a single observing run of 49.5h from June 19th to June 23rd, 2019. Each of the seven sources targeted for this study was observed with the spectral line Eight MIxer Receivers (EMIR) in band E090 (3mm), with an almost complete coverage of the 3mm band, covering 72.0 -79.9 GHz and 84.2 -115.5 GHz, using the wide mode of the Fast Fourier Transform Spectrometers (FTS) (Klein et al. 2012) at a spectral resolution of 200 kHz (i.e., ∼ 0.5 -0.8 km.s -1 ), with an average sideband rejection of -16dB (Carter et al. 2012). The observations were carried out in Wobbler Switching (WSW) mode, with a throw of 180", to ensure a flat baseline across the spectral bandwidth observed. Frequency Switching (FSW) mode was also used to check if the absorption features appearing on the strongest molecular spectra (mainly 12 CO and HCO + ) were due to plausible contamination from the OFF positions or probing different physical components of the targeted sources. More observational time would have been required to also check the absorption features appearing on fainter lines with this method, but these first FSW checks provide some clues on the plausible origins of the absorption features observed, as described in Section 3. Moreover, the wobbler OFF position is given as a single fixed relative offset in azimuth with respect to the target position. This results in a slightly different absolute position in the sky of the OFF position as the Earth rotates. Thus, contamination lines can be checked in WSW data by analyzing if such lines are moving or disappearing during the observation scans. We used this technique to explore the origins of a handful of lines observed in absorption instead of emission (see Section 3.1).

Primary pointing and focus were done on Uranus and secondary focus were done on quasar 0439+360 and on Venus. Pointing was checked every 1.5 hours on these same two nearby sources. The average pointing corrections were between 3" and 4", i.e. 5 to 10 times smaller than the IRAM beam at the observed frequencies, which varies from 34" at 72 GHz to 22" at 110 GHz. The averaged achieved rms is ≈ 5 mK.

The observational data were reduced and analysed using the CLASS/GILDAS software package [START_REF] Pety | SF2A-2005: Semaine de l'Astrophysique Francaise[END_REF][START_REF] Team | GILDAS: Grenoble Image and Line Data Analysis Software[END_REF] and the line identification and fitting was performed with the CASSIS 1 software.

RESULTS

In this section, we present the line detections of small (N atoms ≤ 3) molecules in our 3mm spectral survey of seven Class I YSOs and explore their relationship with one another and with source properties across the sample -detections of larger molecules will be presented in a future article.

Line Detections

Figure 1 presents the full 3mm spectra obtained toward each of the seven Class I sources targeted in this work, sorted by decreasing disk-to-envelope mass ratios. All seven sources present a large number of spectral lines, but the line intensity and richness decreases with disk-to-envelope mass ratio. We detected thirty molecules with no more than 3 atoms toward three or more sources in our sample, where our detection criterion is set to be ≥ 5σ for the integrated intensity. Table 2 lists the most intense line of each species, and the same lines are presented in Figures 2 to 5.

Figure 2 displays spectra for the transitions J = 1 -0 of all CO and HCO + isotopologues observed, with the exception of 12 CO, which will be discussed in a forthcoming study. 13 CO, C 18 O and C 17 O are detected toward all Class I sources, while 13 C 18 O is only seen in the four sources with the smallest disk-to-envelope mass ratios. Since absorption features are seen in 13 CO and C 18 O, which could be indicative of multiple velocity components along the line-of-sight or simply due to the optical thickness of the lines, in the remainder of this work we use C 17 O as a proxy for CO. HCO + , DCO + , and H 13 CO + are detected toward all sources, but the fainter isotopologue HC 18 O + is only detected in emission in 5/7 sources. Faint absorption features are detected in the remaining two disk-dominated sources but are likely due to contamination from the OFF positions, because these features only appear in a few scans of the WSW observations.

Nitrogen-bearing species are presented Fig. 3, with four isotopologues of HCN and HNC each (i.e., main, D, 13 C and 15 N isotopologues), and two of CN and N 2 H + each, (i.e., main and 13 C, and main and D isotopologues, respectively). The main isotopologues of these four Figure 1. Full 3mm spectral line survey obtained with the IRAM-30m (project ID:014-19), toward each of the seven Class I objects targeted in this work, and sorted by decreasing disk-to-envelope mass ratios (from top to bottom). The gap from ∼ 80 to ∼ 84 GHz was not covered by our observations. Ten of the major species detected in this work are presented in red.

molecules are detected in all seven sources, while the rarer isotopologues are mainly seen in the low disk-toenvelope mass ratio sources. In addition, HN 13 C is seen in absorption toward I-04302, but we suspect this to be a foreground cloud contamination. We should also note that multiple velocity components are seen in N 2 D + toward I-04181 and in N 2 H + toward I-04302. These data will need higher resolution observations to be further explored. In the remainder of the present study, we focus on the narrow, central emission or absorption feature toward each source from which first general overview results of the chemistry of Class I objects can be straightforwardly derived.

Figure 4 presents the spectra for the sole 3-atom hydrocarbon observed, C 2 H J = 1 -0. It is detected toward all sources, but is observed in absorption toward I-04302. Why some lines are seen in absorption toward this source is unclear. Contrary to the HC 18 O + 1 -0 line, checking the WSW observation scans one by one suggests that it is less likely that these absorption fea-tures are due to cloud contamination. We speculate that the line could originate from absorption in a cold envelope against optically thick dust in the edge-on disk of this system, but further observations are required to check their origin.

Figure 5 presents nine S-bearing species of interest, three isotopologues of CS (2-1) (including the main, 34 S and 13 C isotopologues), two isotopologues of SO (2 3 -1 2 ) (including the main and 34 S isotopologues), HCS + (2-1), CCS (7 8 -6 7 ), OCS (9-8), and SO 2 (3 1,3 -2 0,2 ). Absorption, self-absorption, or multiple velocity components are apparent for CS toward I-04166, and in both SO isotopologues toward I-04181, I-04016 and I-04365. In these cases, we integrated the emission of the central, narrow component, as indicated above (see Table 3).

Integrated line Fluxes Across Our Source Sample

All line profiles were fitted by a single Gaussian to extract the integrated line intensities. A line is considered detected when S/N ≥ 5σ, and tentatively detected when 5σ > S/N ≥ 3σ. The upper limits on the integrated line intensities were derived assuming a Gaussian shape and a line width of ∼ 0.5 km.s -1 , similar to the typical line width detected for the other molecules in these sources. Table 3 presents the integrated line intensities derived for all species presented in this work, except for the CO isotopologues, where we only focused on C 17 O as mentioned earlier. Figure 6 presents the integrated line intensities of all twenty-seven species versus disk-to-envelope mass ratio. For most species, line intensities vary by 1-2 orders of magnitude across the source sample, even within the central narrow range of disk-toenvelope mass ratios varying from ∼ 0.2 to ∼ 0.5, i.e. excluding I-04181 and I-04302 (see Table 1).

A majority of lines are the weakest in the diskdominated sources, but otherwise there is a range of trends. 33% (9/27) of the targeted lines are the brightest toward the source with smallest disk-to-envelope mass ratio I-04181, namely C 17 O, HC 18 O + , H 13 CN, HC 15 N, C 2 H, C 34 S, HCS + , C 2 S, 13 CS. Most of these lines are typically associated with dense gas (e.g., Bergin et al. 1997;van Dishoeck & Blake 1998;Tafalla et al. 2002;Magalhães et al. 2018;Artur de la Villarmois et al. 2019), and their overabundance in I-04181 may simply reflect more material in this envelope-dominated line-ofsight. 26% (7/27) of the lines are instead brightest toward the source with the second smallest disk-to-envelope mass ratio, I-04169, namely DCO + , H 13 CO + , CS, SO, 34 SO, SO 2 , OCS. Note that this group contains all of the oxygen containing S-species, which is probably suggesting that shocked regions play a non-negligible role. An equal number of lines are brightest toward the source with the third smallest disk-to-envelope mass ratio, I-04166, namely HNC, DCN, DNC, HN 13 C, H 15 NC, N 2 H + , N 2 D + . These are all species associated with cold gas, which is consistent with the relatively low luminosity and low bolometric temperature of this source (see Table 1), likely resulting in a large cold gas reservoir in the surrounding disk and envelope. Note that the only deuterated species that is not the brightest in this source is DCO + , which may reflect different excitation conditions or chemistry for DCO + compared to the other deuterated species, such as its relationship with the depletion of CO (Caselli 2002;Öberg et al. 2015).

HCN and HCO + lines are the brightest toward the most luminous source, I-04016, with increased width likely reflecting their optically thickness, as confirmed by our line optical depth study presented in Appendix A. Finally, CN and 13 CN are the brightest toward I-04365, perhaps indicating that this lower disk-to-envelope mass ratio source is more influenced by stellar UV illumination, since CN can be enhanced through photodissociation of HCN driven by Lyα radiation from the central star (e.g., Bergin et al. 2003). Note that none of the targeted lines of this work are found to be the brightest in the two sources with the highest disk-to-envelope ratios, I-04295 and I-04302.

For each line presented here, we calculated the Spearman's rank correlation coefficient, ρ s , between the integrated line intensity and the disk-to-envelope mass ratio across our Class I source sample, as shown in Fig. 6. This correlation method measures monotonic re- lationships, which is a reasonable approach for our case study where there is no reason to expect strictly linear trends. Moreover, this method is less sensitive to outliers than other methods such as the Pearson correlation method. Nine species display statistically significant anti-correlations with disk-to-envelope mass ratio, i.e. with |ρ s | 0.8 and p-values ≤ 0.05 (where the p-values are used as proxies to assess the statistical significances, e.g., Wasserstein & Lazar 2016): HCS + , H 13 CO + , C 2 S, H 13 CN, HC 18 O + , DCO + , DNC, 13 CS, and DCN.

Integrated Line Flux Ratios

In the following sub-section, we explore the relationship between subsets of species that have been proposed to trace C/N/O gas-phase ratios, deuterium fractiona- b Upper limits are given for the non-detection cases.

c "..." indicates the lines seen in absorption that we suspect to be contamination.

d In case of significant absorption features within emission ones, only lower limits are given. ‡ For species with hyperfine structure, when several hyperfine components are detected we integrated the intensity for all of the components. ‡ ‡13 CN, is the only species with an hyperfine structure for which the integrated intensity is derived only from the main component. tion, temperature, UV-field strength, and evolutionary stage. We note that in the optically thin limit of local thermodynamic equilibrium (LTE), the integrated line flux ratio of two co-spatial molecules presenting similar upper state energies is roughly proportional to their column density ratio. This assumption is, however, only reasonable for a subset of the observed species, since many of the brightest lines are optically thick (see Appendix A). CO, C 2 H, and HCN have been proposed to trace the gas-phase C/N/O ratios (e.g., Du et al. 2015;Bergin et al. 2016;Cleeves et al. 2018;Bergner et al. 2019). Figure 7 shows the integrated 1 -0 line flux ratios of HCN/C 17 O, C 2 H/C 17 O, C 2 H/HCN, and H 13 CN/C 17 O, as functions of the disk-to-envelope mass ratio. There is no clear monotonic trend for any of the ratios, suggestive of a lack of consistent change in C/N/O ratios with disk-to-envelope mass ratio, a proxy for Class I evolutionary stage. We note, however, that this is preliminary since HCN is found to be highly optically thick in 4/7 sources and marginally thick in the three remaining sources (see Appendix A). We expect that this is also the case for C 2 H in most of the sources. Moreover, different sources likely present different excitation environments. The tentative negative correlation between C 2 H/HCN and disk-to-envelope mass ratio likely reflects a combined changing excitation and chemical environment as Class I disks evolve and requires spatially resolved observations to interpret further. If the trend is confirmed with spatially resolved optically thin tracers, part of the chemical explanation could suggest a decreasing C/N fraction due to C sequestration into grain mantles. However, other interpretations would also be possible, since C 2 H has been shown to be very sensitive to UV environment in disks (Bergin et al. 2016;Bergner et al. 2019;Miotello et al. 2019). For instance, an increase of atomic C and C + in the gas phase of the sources where radiation can propagate more easily could also explain this tentative trend.

Deuterium Fractionation

Deuterium enrichment of molecules is intimately linked to their environments (e.g., Lubowich et al. 2000;Dalgarno & Lepp 1984;Roberts et al. 2003;[START_REF] White | Protostars and Planets V[END_REF]Mumma & Charnley 2011;Huang et al. 2017;Aikawa et al. 2018). Constraints on deuterium fractionation are thus commonly used as a diagnostic of chemical inheritance versus chemical reset along the evolutionary path from molecular clouds to planetary systems (we refer the reader to the reviews from Caselli & Cecca-relli 2012;[START_REF] Ceccarelli | Protostars and Planets VI[END_REF], for more exhaustive overviews).

Figure 8 shows the line ratios N 2 D + /N 2 H + , DCO + /HCO + , DCN/HCN, DNC/HNC, DCO + /H 13 CO + , DCN/H 13 CN, and DNC/HN 13 C across our source sample. All ratios indicate highly deuterated gas, but this is in some part an opacity effect. Using the resolved hyperfine patterns of the HCN 1 -0 line and the HCO + 1 -0 isotopologues, we find that both, HCN and HCO + 1 -0, are optically thick toward most of the targeted sources (see Appendix A), and we suspect this is also the case for the HNC 1 -0 line. However, N 2 H + , H 13 CN, and HCO + are found to be optically thin in all sources where they are detected (see Appendix A) and can be used to directly trace the D/H ratio.

Focusing on the optically thin line ratios, we see between a factor of two and five difference for the different species. Assuming that the 12 C/ 13 C ratio remains unchanged among our source sample, it is striking how constant the D/H ratio is across the Class I sample. Furthermore, there is no evident trend with D/H ratio and disk-to-envelope mass ratio. Additional observations are required to further this preliminary D-fractionation study.

Comparison with Class II disk observations are represented in Fig. 8 for the 3 -2 ratios of DCO + /H 13 CO + and DCN/H 13 CN derived in a sample of six Class II disks by Huang et al. (2017), and for the 3 -2 line ratio N 2 D + /N 2 H + derived toward the disks surrounding the T Tauri star AS 209 (Huang & Öberg 2015) and the Herbig-Ae star HD 163296 (Salinas et al. 2017). Despite the facts that we are using a different set of transitions and that the DCN/H 13 CN ratio spanning over one order magnitude in Class II disks is loosely constrained, this preliminary comparison already shows that the DCO + /H 13 CO + and DCN/H 13 CN ratios are in the same order of magnitude from Class I to Class II sources. Especially, in Class II disks the DCO + /H 13 CO + only spans over a factor ∼ 3 and the Class I DCO + /H 13 CO + ratios are all within a factor of 2 of this range. However, additional observations of these ratios in both Class I and Class II sources are required to further this comparison. The N 2 D + /N 2 H + 1 -0 line ratios derived in our Class I source sample are lower by factors of a few to almost an order of magnitude compared to the N 2 D + /N 2 H + 3-2 ratios derived in Class II disks. However, there is a detection bias since N 2 D + has been difficult to detect in Class II disks, so any N 2 D + /N 2 H + ratios reported in Class II disks are going to necessarily be high. In un-shielded gas, such as disk atmospheres, selective photo-dissociation of 14 N 15 N, whose opacity is lower than that of N 2 , may result in 15 N fractionation (e.g., Clayton 2002;Liang et al. 2007;Heays et al. 2014). HCN/HC 15 N is a proposed probe of this process, and could potentially be linked to the UV-field strength of its environments. Chemical 15 N fractionation of HCN has also long been supposed to be favored at low temperature (< 20 K) due to the positive difference in zeropoint energies between the 14 N and 15 N isotopologues, through:

HC 14 NH + + 15 N --→ HC 15 NH + + 14 N, (1)
followed by the dissociative recombination with electron HC 15 NH + + e ---→ HC 15 N + H (Terzieva & Herbst 2000). But, subsequent work demonstrated that the chemical fractionation of Eq (1) is unlikely to happen at low temperatures due to an important activation barrier (Roueff et al. 2015), ruling out a temperature-driven chemical fractionation of HCN. However, Wirström & Charnley (2018) pointed out that the updated reaction rates from Roueff et al. (2015) should be taken with caution, as they were only calculated for one collision angle. Further investigations of the low temperature 15 Nfractionation pathway are still required to completely rule out the temperature-driven chemical differentiation of HCN.

Figure 9 shows the HCN/HC 15 N ratios derived toward the five sources of our sample where H 13 CN and HC 15 N are at least tentatively detected, using the integrated line flux ratio of H 13 CN/HC 15 N J = 1 -0 and assuming a constant 13 C/ 12 C isotopic ratio of 68 ± 15 [START_REF] Milam | [END_REF]Manfroid et al. 2009;Asplund et al. 2009), and how they compare with HCN/HC 15 N ratios from other kinds of sources. The Class I disks span the ratio between the TW Hya disk (Hily-Blant et al. 2019) and prestellar cores (Magalhães et al. 2018). The HCN/HC 15 N Class I measurements are higher by factors of ∼ 1.5 -4 than the ratios measured in a sample of five Class II disks (spanning from 80 to 160, Guzmán et al. 2017) 2017)). However, the ratio may be somewhat overestimated in I-04016 and I-04169 where HC 15 N was only tentatively detected (see Table 3). Moreover, these are all source-averaged ratios. New higher-resolution measurements of the HCN/HC 15 N toward the TW Hya disk found that the isotopic HCN ratio varies from 121 ± 11 to 339 ± 28 from inner to outer radii (Hily-Blant et al. 2019). These new results suggest that the isotopic fractionation is occurring within the Class II disk stage resulting in 15 N enrichment from outer to inner disk regions where comets are supposed to originate. This ratio may then also be a good tracer of how exposed the disk is to UV radiation in the younger Class I disks, but high-resolution Class I disk observations are needed to confirm this.

CN/HCN: Probe Of The UV-field Strength?

The CN/HCN ratio is known to be highly correlated with UV-field strength of the environment it probes (e.g., Schilke et al. 1992;Fuente et al. 2003;Riaz et al. 2018;Bublitz et al. 2019). In protostellar sources in particular, it has been shown to probe the UV flux from the accretion regions of the embedded protostar and its forming circumstellar disk (e.g., Stäuber et al. 2007). In these regions where the stellar FUV flux is higher than the interstellar UV radiation field, the CN/HCN ratio can be enhanced through photodissociation of HCN driven by Lyα radiation from the central star (Bergin et al. 2003), as mentioned earlier ( § 3.2). CN/HCN line ratios 10 have been reported in Class II T Tauri (DM Tau; [START_REF] Dutrey | [END_REF] and Herbig Ae disks (MWC 480 and HD 163296;Thi et al. 2004;[START_REF] Chapillon | [END_REF]. In Class II disks the disk averaged CN/HCN does not vary much from one disk to another ( Öberg et al. 2011), however, but it does vary within disks (Guzmán et al. 2015).

Figure 10 shows the integrated line flux ratio of 13 CN/H 13 CN J = 1 -0 of our Class I source sample, which varies by almost one order of magnitude across our source sample. The source with the lowest 13 CN/H 13 CN ratios, is the one that we previously identified as the coldest of our sample due to the high intensity of deuterated species and to its high HNC/HCN ratio. The highest 13 CN/H 13 CN ratio is found in I-04365.

Based on the 13 CN/H 13 CN J = 1 -0 line ratio as a proxy for the CN/HCN abundance ratio and assum- ing that the latter is a robust UV-field strength probe, I-04365 should be the most exposed sources to UV radiation. However, the H 13 CN/HC 15 N ratio, also proposed as a probe for the UV-field strength (see Section 3.3.3), is found to be the lowest in I-04016. I-04016 and I-04365 are the two most luminous sources of our sample (see Table 1). Thus, while the results found with the 13 CN/H 13 CN and H 13 CN/HC 15 N ratios tend to be consistent, the differences found might indicate that either one of these ratios also depends on other parameters; or that each of these ratios trace different layers or components of the Class I source probed. Additional observations are required to test these scenarios.

HNC/HCN: A Probe Of Temperature?

The HNC/HCN J = 1 -0 line ratio is suggested to be a direct probe of the gas kinetic temperature, such that a higher ratio probes colder material (e.g., Baan et al. 2008;Hacar et al. 2019). In Fig. 11 we show the HNC/HCN and HN 13 C/H 13 CN J = 1 -0 line ratios. When comparing the top and bottom panels of this figure, we see that the order of the two ratios from one source to another is generally replicated. The lowest ratios are found in the sources with the highest disk-to-envelope mass ratios, suggesting they are overall warmest. Among the envelope dominated sources, there is no trend with disk-to-envelope mass ratio. However, this is not surprising since source luminosity should also play a large role for setting the temperature. In fact, we do see the highest HNC/HCN ratio toward the lowest luminosity source, I-04166, which is also the coldest source of our sample (see Table 1). Interestingly, the H 13 CN/HC 15 N is also found to be the highest in this source, highlighting the fact that the latter is indeed related to the luminosity of the source and that 15 N fractionation of HCN is less efficient at low temperature (see Section 3.3.3).

Sulfur-bearing Species Ratios

Sulfur-bearing species ratios have long been thought to be powerful diagnostics to probe both the chemistry and the physics of their environments. For instance, the CS/SO ratio has been proposed as a probe for the elemental C/O ratio (e.g Swade 1989;Bergin et al. 1997;Fuente et al. 2019). Sulfur-bearing ratios are also thought to be good chemical clocks (e.g., SO/SO 2 , CS/OCS, SO/OCS), because of their rapid chemical evolution with gas temperature and density (e.g., Charnley 1997;Hatchell et al. 1998;Viti et al. 2004;Herpin et al. 2009;Wakelam et al. 2011;Vidal & Wakelam 2018). This idea came up from hot core/corino theory, with the hypothesis that the bulk of the sulfur reservoir is locked onto grain mantles in previous evolutionary stage and then released in the gas phase during the hot core/corino formation under the form of H 2 S, SO and atomic S. This is supposed to initiate an active S-chemistry, where SO 2 is a direct product of SO, through the radiative association of O and SO and the neutral-neutral reaction OH + SO --→ SO 2 + H, making the SO/SO 2 ratio an ideal chemical clock for such astrophysical objects. In the gas phase, OCS is mainly formed through the radiative association of S and CO, and one of the formation pathways for CS is the reaction SO + C, justifying that CS/OCS and SO/OCS were also claimed as additional chemical evolutionary clocks. Finally, finding out whether sulfur is reduced or oxidized in planet-forming environments is interesting for origins of life chemistry (Patel et al. 2015).

Figure 12 shows the S-bearing integrated line flux ratios of CS/SO, CS/OCS, CS/C 2 S, SO/OCS, and SO/SO 2 with corresponding C 34 S ratios as functions of increasing disk-to-envelope mass ratio across the Class I sample. These ratios are also found to vary by several orders of magnitude from one source to another, but only CS/C 2 S (and C 34 S/C 2 S) show a clear positive correlation with increasing disk-to-envelope mass ratio, respectively. This is a priori counter-intuitive because, for instance, C 2 S can be produced from the reaction CS + CH --→ C 2 S + H in the gas phase (see the KIDA database2 ). Hence, one would rather expect an increase in carbon chains with evolution. However, it should be noted that the chemistry of C 2 S has remained poorly studied experimentally and theoretically so far.

DISCUSSION

Chemical Differentiation Within Our Sample

In Section 3, we explored specific line emission patterns within our sample. In this sub-section we investigate global patterns of chemical differentiation. Fig. 13 shows a cluster correlation matrix of all sources and line fluxes presented in this work, based on Spearman's correlation coefficients. This correlation matrix method produces a table containing the correlation coefficients in between each species, allowing us to simultaneously assess their interdependence. The clustering is then used to arrange and group the species by correlation coefficients, as illustrated by the dendrogram in the figure .  This global statistical representation highlights several clusters of correlated molecules with two or more sub-clusters each. At the bottom right of the plot lies the biggest cluster consisting of N 2 D + , N 2 H + , DCN, HNC, DCO + , H 15 NC, DNC, H 13 CN, C 2 S, and HN 13 C, which are readily identified with dense and cold gas, since these species are known to have high critical densities (i.e., tracing dense gas, Shirley 2015), and were observed in their 1-0 ground-state transition in our survey (i.e., tracing cold gas). Furthermore, deuterated species and especially N 2 D + are characteristics of cold gas for chemical reasons. Bright HNC emission is also associated with cold gas, while HCN isotopologues are more of dense gas tracers. But, the fact that they show up in the same cluster may suggest that sources that are cold are also more likely to have large reservoirs of dense gas. In either case, the first differentiating factor among Class I sources affecting their temperature is whether or not they have a large cool gas reservoir.

At the top left of the figure lies the second biggest cluster, overall less strongly correlated than the previous one, but composed of several highly correlated subclusters. It contains most of the S-bearing species but OCS and C 2 S. With its abundance in oxygenated Sbearing species, this molecular cluster is probably characteristic of shock induced chemistry tracers, which appear to be closely correlated with CS and isotopologues, C 17 O, and with ions such as H 13 CO + and HC 18 O + . The coupling of these last species with shock tracers could either be due to an evolutionary effect, or to the fact that sources with shocks may also have more ionized gas and/or more CO in the gas-phase. C 17 O (our proxy of CO in the gas phase, see Section 3.1) is indeed found to be the brightest in I-04169, where most S-bearing species and H 13 CO + are also the brightest.

Finally, two smaller clusters lie in the center of the matrix: (i) a first one made of (i) HCO + , OCS, and HCN; and (ii) another one made of CN, 13 CN, and C 2 H. Interestingly, the first molecular cluster is however barely correlated with the fainter isotopologues of HCO + and HCN, which suggests that the correlations of these molecules are dominated by excitation and opacity effects rather than chemical relationships. The second cluster contains UV-driven dense gas tracers, with CN and C 2 H being commonly observed in UV-enhanced environments (e.g., Kastner et al. 2015;Bergin et al. 2016;Cazzoletti et al. 2019;van Terwisga et al. 2019;Bergner et al. 2019;Miotello et al. 2019). SO 2 emission was found to be the brightest toward the most luminous sources, whereas the opposite was observed for C 17 O. The fact that these two molecules are not tracing similar environments is not surprising, since these two molecules are known to trace different physical conditions, where the former is often associated with outflows and shocked gas (Jørgensen et al. 2004;[START_REF] Pety | SF2A-2005: Semaine de l'Astrophysique Francaise[END_REF], whereas the latter is rather tracing denser and colder gas (Caselli et al. 1999;Tafalla et al. 2002;Redman et al. 2002). However, Artur de la Villarmois et al. (2019) propose that in Class I sources, the ratio of these two molecules could provide an interesting evo- The redder the color is, the more positive the correlation is and, similarly, the bluer the color is, the higher the anti-correlation is. Note that correlation calculations include upper limits. lution diagnostic, assuming that the luminosity of the source is linked to the evolutionary stage of the source.

SO

Figure . 14 represents the integrated line flux ratios SO 2 /C 17 O across our Class I sample as function of the luminosity of the sources. First, we can see that in our sample the luminosity and disk-to-envelope mass ratio are not correlated and luminosity alone is thus not likely a good evolutionary indicator here. For the sources where the SO 2 line we targeted here is detected (i.e., 4/7 in our sample), there is a noticeable increase in the SO 2 /C 17 O ratio with luminosity, but the highest luminosity sources present informative upper limits on this ratio which dramatically breaks the trend.

Class 0 versus Class I versus Class II chemistry

As opposed to Class I systems, many more spectral surveys have been performed to study the chemistry of the earliest stages of star formation, i.e. Class 0 objects -e.g., TIMASSS (Caux et al. 2011), PILS (Jørgensen et al. 2016), ASAI (Lefloch et al. 2018), SOLIS (Ceccarelli et al. 2017) -and of late planet-forming disks, i.e. Class II objects -DISCS ( Öberg et al. 2010, 2011), CID (Guilloteau et al. 2016). The question that is still pending is to determine if the chemical composition of Class I objects is more Class 0-like, believed to be associated with CH 3 OH and long carbon chains (e.g., Bergner et al. 2017;Law et al. 2018) or Class II-like, dominated by small molecules and an O-poor organic chemistry (e.g., Öberg et al. 2011;Guilloteau et al. 2016). Although the non-detection of larger molecules in Class II objects could be due to a detection bias, as the column densities typically decrease from Class 0 to Class II stages, non-oxygenated complex molecules, such as HC 3 N and CH 3 CN, are commonly observed in Class II disks (e.g., Bergner et al. 2018) as compared to O-containing complex molecules, such as CH 3 OH, therefore highly favoring that disk organic chemistry is indeed O-poor (e.g., Le Gal et al. 2019a).

In the present work, we focused on small (N atoms ≤ 3) molecules, with the detections of a wide variety of C, N, O, and S carriers -namely CO, HCO + , HCN, HNC, CN, N 2 H + , C 2 H, CS, SO, HCS + , C 2 S, SO 2 , OCS -and some of their D, 13 C, 15 N, 18 O, 17 O, and 34 S isotopologues. This set of detections is thus already providing interesting information regarding the fact that Class I systems are richer than Class II objects in small molecules, since all the Class II disk species detected so far are also present in Class I stage and that, in addition, oxygenbearing species are also detected in Class I stage, such as SO, SO 2 , and OCS. In particular, we note that SO is systematically detected in all the sources of the Class I sample targeted here, contrary to Class II disks where it has only been reported toward a few young disks presenting active accretion signs (Fuente et al. 2010;Guilloteau et al. 2016;Pacheco-Vázquez et al. 2016;Booth et al. 2018).

However, interferometric observations are needed to distinguish which source components are traced by each of these species. This would be, in particular, beneficial to better understand the lack of O-bearing species in Class II disks, such as the lack of oxygenated S-bearing species reported so far (e.g., Semenov et al. 2018;Le Gal et al. 2019b). As for the comparison with Class 0 objects, complex organic detections from the same Class I survey will be presented in a forthcoming study and should thus allow to identify if long carbon chain and/or complex organic molecules are still present in Class I stage and in which quantities and proportions as compared to small organics.

SUMMARY AND CONCLUSION

We presented a new single-dish 3mm spectral line survey carried out with the IRAM-30m telescope toward a sample of seven Class I YSOs. Our main findings are summarized below:

1. Our 3mm survey found that Class I YSOs are molecule-rich, considering the large number of spectral lines we found (see Fig. 1). We report the detections of twenty-seven spectral lines characterizing the main small (N atoms ≤ 3) molecules of our survey -namely CO, HCO + , HCN, HNC, CN, N 2 H + , C 2 H, CS, SO, HCS + , C 2 S, SO 2 , OCS, including some of their D, 13 C, 15 N, 18 O, 17 O, and 34 S isotopologues.

2. For most species, integrated line intensities vary by 1-2 orders of magnitude across the source sample, even within the central narrow range of diskto-envelope mass ratios varying from 0.2 to 0.5. There is a general negative trend in line intensity with increasing disk-to-envelope ratio, used as a proxy for the Class I evolutionary stage, but the more salient point is the large intensity variations found among sources with similar total masses.

3. Line flux ratios proposed to trace C/N/O gasphase ratios, do not show clear monotonic trends with disk-to-envelope mass ratio. A tentative negative correlation between C 2 H/HCN and the diskto-envelope mass ratio may indicate a decreasing C/N fraction due to C sequestration into grain mantles, although other scenarios, e.g. driven by UV-penetration, might also be possible.

4. Line flux ratios that probe deuterium fractionation are rather constant across the Class I sample, which indicates similar thermal histories and/or current thermal structure toward Class I objects. Higher angular and spectral resolution observations are required to corroborate these conclusions.

5. Comparing the HCN/HC 15 N measured in our Class I sample with ratios from other astrophysical environments, we found that the Class I ratios span between prestellar core and Class II disk values. This might either be indicative of the span in evolution of our sample, if the isotopic ratio is inherited; or of the span in physical properties of our sample, if the HCN/HC 15 N is reprocessed in situ, since this ratio is likely to be strongly altered with UV radiation in un-shielded environments where the visual extinction lies within 1 to 3 magnitudes (e.g., Heays et al. 2014).

6. The CN/HCN ratio, proposed as a UV-field strength diagnostic, varies by almost one order of magnitude across our source sample. Its highest values are reached in I-04365 suggesting that this source might be the most UV-illuminated, but detailed abundance modelling is needed to confirm.

7. The HNC/HCN ratio, proposed as a gas temperature proxy, is lower toward the sources with higher disk-to-envelope ratios, and the highest toward the source with the lowest luminosity and bolometric temperature, I-04166, where species characteristic of cold gas are the brightest.

8. We explored the use of ratios of S-bearing species -CS/SO, CS/OCS, CS/C 2 S, SO/SO 2 , and SO/OCS -as tracers of the C/O elemental ratio and chemical evolutionary stage. Only CS/C 2 S and SO/SO 2 show correlations with the disk-toenvelope mass ratio, which are found to be opposite. These opposite trends could indicate that these ratios probe different physical components of the Class I sources, but it could also be due to oxidized versus reduced S-chemistry. Higher angular resolution observations are required to further explored these results. 9. To illustrate the chemical differentiation within our sample, we built a cluster correlation matrix that highlights three main groups of correlated species: (i) dense cold gas tracers, (ii) shocked gas tracers, and (iii) dense ionized gas tracerslikely probing more UV-illuminated regions than the previous ones.

10. To compare our results to other Class I surveys, we focused on the SO 2 /C 17 O ratio proposed to be a Class I evolutionary tracer, assuming that the latter is correlated with source luminosity (Artur de la Villarmois et al. 2019). First, in our sample, the disk-to-envelope mass ratio that we proposed as a proxy for evolutionary stage is not correlated with source luminosity. Second, although there could be a tentative correlation in between this ratio and source luminosity within our sample, the fact that the ratios derived in the highest luminosity sources rely on upper limits breaks the trend.

11. Finally, we discuss the similarities and differences we found when comparing our Class I chemical results to Class 0 and Class II chemistry. Class I sources are found to be more chemically rich and diverse than Class II sources, the latter being known to harbor an O-poor organic chemistry. Class 0 objects are usually associated with more complex molecules and long carbon chains. However, in the present work we focused on small (N atoms ≤ 3) molecules. The composition in bigger and more complex organic molecules within the same sample of Class I sources will be investigated in a forthcoming study.

These results already provide important information on the chemical content and the relationship between Class I chemistry, evolutionary stage and other source properties. Since all sources are located in the same starforming region, our analysis avoids cloud-to-cloud background variations, but as a result, it may not be rep-resentative of all regions of the sky. Additional observations are required to extend this Class I sample and increase the statistics and demography on Class I YSO chemistry. Interferometric observations will be beneficial to elucidate the contributions of envelope and disk chemistry to the observed emission lines. Future work on complex organics in the same sample should also provide additional clues on the evolution and chemical content of these objects. [START_REF] Milam | [END_REF]Asplund et al. 2009;Manfroid et al. 2009) to be compared to the HCO + /HC 18 O + ratio and to the local ISM 16 O/ 18 O ratio (see Section A.1), represented in gray.  (Chin et al. 1996) and 12 C/ 13 C= 68 ± 15 in the local ISM [START_REF] Milam | [END_REF]Asplund et al. 2009;Manfroid et al. 2009) are represented in gray.

toward the sources of our sample where their hyperfine patterns are well-resolved (see Fig. 3), using LTE modeling through the CASSIS tool. The results are listed in Table 4 and shown for each molecule for one source, I-04166, in Fig. 17, for illustration. HCN is found to be highly optically thick in 4/7 sources, and marginally optically thick in the remaining three sources. H 13 CN and N 2 H + are found to be optically thin in all sources where they are detected. Artur de la Villarmois, E., Jørgensen, J. K., Kristensen, L. E., et al. 2019, A&A, 626, A71 Aso, Y., Ohashi, N., Saigo, K., et al. 2015, ApJ, 812, 27 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

Figure 2 .

 2 Figure2. Spectra of the transition J = 1 -0 of 13 CO, C 18 O, C 17 O, 13 C 18 O (upper panels), HCO + , DCO + , H 13 CO + , and HC 18 O + (lower panels) observed toward each of the seven Class I sources targeted in this work. Spectra are sorted by decreasing disk-to-envelope mass ratios from top to bottom in each panel, and by decreasing intensity in each row of panels.

Figure 3 .

 3 Figure3. Spectra of the transition 1 -0 of HCN, DCN, H 13 CN, HC 15 N (upper panels), HNC, DNC, HN 13 C, H 15 NC (middle panels), N2H + , N2D + , CN and 13 CN (lower panels) observed toward each of the seven Class I sources of our sample. Spectra are sorted by decreasing disk-to-envelope mass ratios from top to bottom in each panel, and by decreasing intensity in each row of panels.

Figure 4 .

 4 Figure 4. Spectra of the strongest component of C2H N = 1 -0 observed toward each of the seven targeted sources of our sample, sorted by decreasing disk-to-envelope mass ratios from top to bottom in each panel.

Figure 5 .

 5 Figure5. Spectra of the strongest component of CS, SO, HCS + , CCS, 13 CS, C 34 S, 34 SO, SO2, and OCS observed toward each of the seven Class I sources of our sample. Spectra are sorted by decreasing disk-to-envelope mass ratios from top to bottom in each panel.

Figure 6 .

 6 Figure 6. Integrated line intensities of the species of interest presented in this study and represented as function of the diskto-envelope mass ratio across our sample of Class I YSOs. Error bars are represented by the vertical line segments. Upper limits are represented by the downward triangles and lower limits by the upward triangles. Spearman correlation coefficients are indicated in the lower left corner of each panel, with corresponding p-value in brackets. These statistical parameters are described in § 3.2. Note that correlation calculations include upper limits.

  3.3.1. CO, C2H and HCN: Probing C/N/O Elemental Ratios?

Figure 7 .

 7 Figure 7. Integrated line flux ratios of C2H/C 17 O, HCN/C 17 O, C2H/HCN, and H 13 CN/C 17 O derived from the transition J = 1 -0 of each species toward each of the seven Class I sources of our sample, represented as function of the disk-to-envelope mass ratio. Error bars are indicated by the vertical line segments. Upper limits are represented by the downward triangles.

Figure 8 .

 8 Figure8. Integrated line flux ratios of N2D + /N2H + , DCO + /HCO + , DCN/HCN, DNC/HNC, DCO + /H 13 CO + , DCN/H 13 CN and DNC/HN 13 C 1 -0 derived for the transition J = 1 -0 in each of the seven Class I sources of our sample, sorted by decreasing disk-to-envelope mass ratios. The gray shaded regions represent the line flux ratios derived toward Class II disks fromHuang & Öberg (2015) andSalinas et al. (2017) for the N2D + /N2H + and fromHuang et al. (2017) for the DCO + /H 13 CO + and DCN/H 13 CN ratios. Upper limits are represented by the downward triangles and lower limits by the upward triangles. Error bars are indicated by the vertical line segments. Bockelée-Morvan et al. 2015; Shinnaka et al. 2016, see also the compilation made in Fig. A1 of Hily-Blant et al. (2017)). However, the ratio may be somewhat overestimated in I-04016 and I-04169 where HC 15 N was only tentatively detected (see Table3).Moreover, these are all source-averaged ratios. New higher-resolution measurements of the HCN/HC 15 N toward the TW Hya disk found that the isotopic HCN ratio varies from 121 ± 11 to 339 ± 28 from inner to outer radii(Hily-Blant et al. 2019). These new results suggest that the isotopic fractionation is occurring within the Class II disk stage resulting in 15 N enrichment from outer to inner disk regions where comets are supposed to originate. This ratio may then also be a good tracer of how exposed the disk is to UV radiation in the younger Class I disks, but high-resolution Class I disk observations are needed to confirm this.

Figure 9 .

 9 Figure9. Integrated line flux ratio of HCN/HC 15 N J = 1-0 represented as function of the increasing disk-to-envelope mass ratio toward 5/7 of the Class I targeted sources where HC 15 N is at least tentatively detected. These ratios were derived from the H 13 CN/HC 15 N J = 1 -0 ratio, assuming the local ISM isotopic ratio of 12 C/ 13 C of 68 (see text § 3.3.3). The horizontal lines and rectangles represent the Solar (orange), Prestellar core (light blue), Earth (green) Class II disks (red), and Comets (dark blue) ratios for references. Error bars are indicated by the vertical line segments.

Figure 10 .

 10 Figure 10. Integrated line flux ratios of the 13 CN/H 13 CN J = 1 -0 as function of the increasing disk-to-envelope mass ratio of our Class I source sample. Upper limits are represented by the downward triangles. Error bars are indicated by the vertical line segments.

Figure 11 .

 11 Figure 11. Integrated line flux ratios of the HNC/HCN (top) and HN 13 C/H 13 CN (bottom) J = 1 -0 as function of the increasing disk-to-envelope mass ratio of our Class I source sample. Upper limits are represented by the downward triangles. Error bars are indicated by the vertical line segments.

Figure 12 .

 12 Figure12. Integrated line flux ratios of CS/SO, CS/OCS, CS/C2S, SO/SO2, and SO/OCS, with corresponding C 34 S ratios on the bottom row panels, derived toward our Class I source sample. These ratios are represented as function of the increasing disk-to-envelope mass ratio. Upper limits are indicated by the downward triangles and lower limits by the upward triangles. Error bars are indicated by the vertical line segments. Spearman correlation coefficients are indicated in the lower left corner of each panel, with corresponding p-value in brackets. These statistical parameters are described in § 3.2. Note that correlation calculations include upper limits.

2

  Versus C 17 O In Class I Sources A recent study of twelve low-mass Class I sources located in the Ophiuchus molecular cloud shed light on the chemical differentiation between SO 2 and C 17 O across their source sample (Artur de la Villarmois et al. 2019).

Figure 14 .

 14 Figure 14. SO2/C 17 O integrated line flux ratio as function of the source luminosity across our Class I sample. Upper limits are represented by the downward triangles.

Figure 15 .

 15 Figure15. Integrated line flux ratios of HCO + /HC 18 O + and H 13 CO + /HC 18 O + as function of the increasing disk-to-envelope mass ratio across the five Class I sources where HC 18 O + is detected. Error bars are indicated by the vertical line segments. The H 13 CO + /HC 18 O + ratio is multiplied by the elemental 12 C/ 13 C= 68 ± 15 in the local ISM[START_REF] Milam | [END_REF] Asplund et al. 2009;Manfroid et al. 2009) to be compared to the HCO + /HC 18 O + ratio and to the local ISM 16 O/ 18 O ratio (see Section A.1), represented in gray.

Figure 16 .

 16 Figure16. Integrated line flux ratios of CS/C 34 S and CS/ 13 CS as function of the increasing disk-to-envelope mass ratio across the seven Class I sources targeted in this work. Lower limits are indicated by the upward triangles. Error bars are indicated by the vertical line segments. The elemental 32 S/ 34 S= 24.4 ± 5.0 in the vicinity of the Sun(Chin et al. 1996) and 12 C/ 13 C= 68 ± 15 in the local ISM[START_REF] Milam | [END_REF] Asplund et al. 2009;Manfroid et al. 2009) are represented in gray.

Figure 17 .

 17 Figure17. LTE modeling obtained for the hyperfine transitions of HCN, H 13 CN, and N2H + J = 1 -0 toward I-04166.

Table 1 .

 1 Characteristics of the targeted sources sorted by decreasing disk-to-envelope mass ratio

	Source	R.A.

  Right ascension (R.A.) and declination (Dec.) of each source are from the Gaia DR2 catalog (Gaia Collaboration 2018) b These values correspond to the best-fit parameters derived in Sheehan & Eisner (2017) with their MCMC model, where the uncertainties are the range containing 68% of the posterior distribution. Systematic errors could be larger.

								au)	(au)	(km/s)	(pc)
	IRAS 04302+2247	04:33:16.501 22:53:20.400 122 (c)	0.4	0.017 +0.006 -0.004	0.114 +0.019 -0.026	6.7	1086	244	5.5 [1]	161±3 (f )
	IRAS 04295+2251	04:32:32.055 22:57:26.670 270 (c)	0.3	0.037 +0.008 -0.006	0.018±0.001	0.49	1081	127	5.3 [1]	161±3 (f )
	IRAS 04365+2535	04:39:35.194 25:41:44.730 164 (d)	2.1	0.071 +0.035 -0.019	0.030 +0.002 -0.003	0.42	1829	143	6.6 [2]	140±4 (f )
	IRAS 04016+2610	04:04:43.071 26:18:56.390 226 (d)	7.0	0.023 +0.010 -0.004	0.009±0.001	0.39	1446	497	6.8 [2]	∼140 (g)
	IRAS 04166+2706	04:19:42.627 27:13:38.430	75 (c)	0.2	0.100±0.009 0.027±0.003	0.27	1209	180	6.7 [3]	160±3 (f )
	IRAS 04169+2702	04:19:58.449 27:09:57.070 133 (c)	0.8	0.055 +0.004 -0.005	0.012±0.001	0.22	672	39	6.8 [2]	160±3 (f )
	IRAS 04181+2654A 04:21:11.469 27:01:09.400 346 (e)	0.3	1.234 +0.688 -0.391	0.006±0.001	4.8e-3	> 20000	47	7.1 [1]	160±3 (f )
	a c From Young et al. (2003)
			d From Green et al. (2013)

Table 2 .

 2 Properties of all detected lines from CDMS(a)

	Molecule	Frequency Eup gup	A ij	Transition
		(MHz)	(K)		(s -1 )
	13 CO	110201.354	5.3	6	6.33e-08	-0
	C 18 O	109782.173	5.3	3	6.27e-08	-0
	C 17 O ( ‡)	112358.982	5.4	8	6.70e-08	1 -0 3
	13 C 18 O	104711.404	5.0	4	5.46e-08	1 -0 1
	HCO +	89188.525	4.3	3	4.19e-05	-0
	DCO +	72039.312	3.5	3	2.21e-05	-0
	H 13 CO +	86754.288	4.2	3	3.85e-05	-0
	HC 18 O +	85162.223	4.1	3	3.65e-05	-0
	HCN ( ‡)	88631.847	4.3	5	2.43e-05	1 -0 1
	DCN ( ‡)	72414.905	3.5	5	1.31e-05	1 -0 1
	H 13 CN ( ‡)	86339.921	4.1	9	2.23e-05	-0
	HC 15 N	86054.966	4.1	3	2.20e-05	-0
	HNC	90663.568	4.4	3	2.69e-05	-0
	DNC	76305.700	3.7	3	1.60e-05	-0
	HN 13 C	87090.825	4.2	3	2.38e-05	-0
	H 15 NC	88865.715	4.3	3	1.98e-05	-0
	CN ( ‡)	113490.970	5.4	6	1.19e-05 1 0 2 -0 0 1 2
	13 CN ( ‡)	108780.201	5.2	7	1.05e-05 1 2 2 -0 1 1 2
	N2H +( ‡)	93173.700	4.5	15	3.63e-05	1 -0 1
	N2D +( ‡)	77109.610	3.7	5	2.06e-05	1 -0 1
	C2H ( ‡)	87316.898	4.2	5	1.53e-06	1 2 -0 1 1
	CS	97980.953	7.1	5	1.68e-05	2 -1 0
	13 CS	92494.308	6.7	10	1.41e-05	2 -1 0
	C 34 S	96412.950	6.9	5	1.60e-05	2 -1 0
	SO	99299.870	9.2.	7	1.13e-05	2 -1 2
	34 SO	97715.317	9.1	7	1.07e-05	2 -1 2
	HCS +	85347.890	6.1	5	1.11e-05	-1
	CCS	93870.107	19.9	17	3.74e-05	7 -6 7
	SO2	104029.418	7.7	7	1.01e-05	3 1 -2 0 2
	OCS	109463.063 26.3	19	3.70e-06	-8

a https://cdms.astro.uni-koeln.de/cdms/portal/;

(Müller et al. 2001(Müller et al. , 2005

) ‡ For species with hyperfine structure, properties are only given for the main component.

Table 3 .

 3 Integrated line intensities in mK km/s of the twenty-seven species of interest for this work.

	Molecule	Source	I-04302	I-04295	I-04365	I-04016	I-04166	I-04169	I-04181
	C 17 O ( ‡)		26.7 ± 2.3	33.0 ± 1.6	21.9 ± 1.6	31.6 ± 3.9	21.7 ± 1.6	76.6 ± 1.9	100.5 ± 4.2
	HCO +		177.8 ± 2.5	91.7 ± 2.8	165.6 ± 3.1	1228.1 ± 7.3	398.5 ± 5.2	346.6 ± 15.4	632.6 ± 5.9
	DCO +		6.4 ± 1.8 (a)	< 10.1 (b)	136.9 ± 2.7	109.8 ± 2.9	363.0 ± 3.3	576.4 ± 4.7	199.2 ± 5.3
	H 13 CO +		3.7 ± 1.0	25.5 ± 1.0	159.0 ± 1.0	200.2 ± 1.1	238.6 ± 1.1	387.8 ± 1.3	271.4 ± 1.3
	HC 18 O +		... (c)	...	20.1 ± 1.4	17.9 ± 1.3	28.7 ± 1.1	40.5 ± 1.3	42.3 ± 1.3
	HCN ( ‡)		77.7 ± 3.1	38.4 ± 2.5	311.7 ± 1.9	603.1 ± 5.2	492.7 ± 3.9	460.5 ± 6.0	364.8 ± 2.5
	DCN ( ‡)		< 11.5	< 11.5	30.1 ± 4.0	29.0 ± 3.4	147.5 ± 2.8	87.6 ± 4.3	63.8 ± 6.6
	H 13 CN ( ‡)		< 10.5	< 10.5	33.3 ± 3.5	27.2 ± 3.3	62.3 ± 1.9	49.3 ± 2.5	71.0 ± 3.3
	HC 15 N		< 10.6	< 4.7	7.0 ± 1.3	8.3 ± 1.8	11.2 ± 1.6	9.9 ± 2.9	20.0 ± 1.4
	HNC		23.9 ± 2.7	20.8 ± 2.2	558.8 ± 0.05	497.9 ± 1.8	969.1 ± 6.5	742.8 ± 2.1	348.9 ± 1.5
	DNC		< 4.6	< 4.6	93.2 ± 1.9	34.4 ± 1.9	342.4 ± 1.8	180.3 ± 2.7	198.6 ± 2.5
	HN 13 C		...	...	82.1 ± 1.6	19.1 ± 1.1	173.4 ± 1.4	75.8 ± 1.8	115.7 ± 1.3
	H 15 NC		< 4.7	< 4.7	13.5 ± 1.1	5.5 ± 1.9	34.4 ± 1.1	12.3 ± 1.0	29.0 ± 1.2
	CN		279.5 ± 22.5	259.2 ± 15.2 1682.0 ± 9.3 698.7 ± 15.8	780.0 ± 8.9	1354.2 ± 7.7	1178.7 ± 10.0
	13 CN		< 7.2	< 4.8	30.6 ± 4.1	11.8 ± 2.8	6.6 ± 1.7	16.1 ± 1.8	21.9 ± 6.4
	N 2 H +		< 5.9	17.0 ± 1.8	1044.0 ± 0.1	149.3 ± 3.3	1906.4 ± 0.2	997.2 ± 1.1	520.8 ± 4.5
	N 2 D +		< 4.6	< 8.0	40.2 ± 3.3	< 13.8	202.5 ± 2.7	105.2 ± 4.5	> 15.9 ± 3.1 (d)
	C 2 H		...	53.3 ± 4.3	922.2 ± 2.7	460.5 ± 2.7	589.8 ± 3.7	847.8 ± 4.6	1193.8 ± 5.9
	CS		228.5 ± 4.5	143.2 ± 2.7	180.7 ± 2.3	312.0 ± 3.5	> 98.6 ± 4.1	660.5 ± 1.3	335.6 ± 4.7
	C 34 S		15.5 ± 3.3	10.8 ± 2.3	14.5 ± 1.6	8.9 ± 2.2	15.2 ± 1.6	67.2 ± 2.0	104.9 ± 2.2
	13 CS		< 2.9	7.4 ± 1.6	4.4 ± 1.1	< 3.8	7.4 ± 1.7	20.5 ± 1.4	65.5 ± 1.7
	SO		37.0 ± 4.2	44.3 ± 2.4	> 16.2 ± 2.0 > 23.1 ± 5.3	68.1 ± 6.4	564.0 ± 1.8	> 129.0 ± 5.1
	34 SO		< 3.0	< 3.0	< 3.0	4.4 ± 1.3	< 3.0	78.6 ± 1.3	14.5 ± 1.6
	HCS +		< 2.9	5.9 ± 1.4	4.6 ± 1.2	7.4 ± 2.0	9.3 ± 1.1	10.9 ± 1.9	69.4 ± 1.8
	C 2 S		< 3.8	3.3 ± 1.0	15.3 ± 1.3	7.8 ± 1.3	22.9 ± 0.8	28.0 ± 1.2	181.6 ± 1.6
	SO 2		< 4.2	8.8 ± 1.8	< 2.4	< 2.4	< 2.4	54.4 ± 1.9	8.7 ± 2.8
	OCS		< 4.3	< 3.0	< 3.0	9.8 ± 2.2	...	23.3 ± 2.8	9.5 ± 2.5

a Tentative detection are indicated in italics, where our tentative detection criterion is set to be 5σ > S/N ≥ 3σ.

  and in comets (with an average ratio of 140 ± 3,Jehin et al. 2009;Mumma & Charnley 2011;

  Cluster correlation matrix plot for all the 27 species presented in this work, based on Spearman's correlation coefficients indicated in each box. This statistical parameter is described in § 3.2. The dendrogram illustrates the arrangement produced by the clustering analysis.

	1.00 0.91 0.75 0.48 0.54 0.46 0.50 0.50 0.31 0.39 0.93 0.43 0.25 0.49 0.36 0.57 0.25 -0.50 0.50 0.25 0.07 -0.07 0.29 0.00 0.11 -0.18 -0.04	CS
	0.91 1.00 0.77 0.49 0.81 0.77 0.56 0.47 0.58 0.61 0.89 0.57 0.47 0.39 0.39 0.45 0.32 -0.41 0.67 0.50 0.36 0.18 0.57 0.22 0.40 0.08 0.26	34 SO
	0.75 0.77 1.00 0.78 0.43 0.57 0.60 0.43 0.61 0.57 0.52 0.18 -0.11 0.31 0.00 0.21 0.11 -0.30 0.32 0.18 -0.04 -0.18 0.14 -0.39 -0.07 -0.36 -0.21	C 17 O
	0.48 0.49 0.78 1.00 0.26 0.33 0.78 0.48 0.56 0.70 0.28 -0.30 -0.41 -0.03 -0.11 -0.11 -0.04 -0.11 0.11 -0.04 -0.15 -0.34 0.15 -0.33 -0.07 -0.33 -0.15	SO 2
	0.54 0.81 0.43 0.26 1.00 0.93 0.80 0.50 0.77 0.71 0.75 0.57 0.64 0.54 0.61 0.39 0.46 0.10 0.89 0.85 0.81 0.67 0.93 0.64 0.85 0.61 0.75	H 13 CO +
	0.46 0.77 0.57 0.33 0.93 1.00 0.90 0.50 0.86 0.79 0.64 0.61 0.54 0.43 0.36 0.21 0.54 0.30 0.82 0.83 0.76 0.61 0.79 0.39 0.72 0.43 0.57	HCS +
	0.50 0.56 0.60 0.78 0.80 0.90 1.00 1.00 1.00 0.80 0.00 -0.10 -0.40 0.70 0.30 0.20 0.80 0.50 1.00 0.90 0.70 0.50 0.70 -0.10 0.60 0.10 0.30	HC 18 O +
	0.50 0.47 0.43 0.48 0.50 0.50 1.00 1.00 0.58 0.71 0.29 0.14 -0.07 0.77 0.39 0.36 0.75 0.50 0.71 0.59 0.49 0.38 0.46 0.14 0.45 0.14 0.25	C 34 S
	0.31 0.58 0.61 0.56 0.77 0.86 1.00 0.58 1.00 0.81 0.23 0.14 0.09 0.49 0.38 0.20 0.34 0.50 0.72 0.74 0.65 0.54 0.72 0.22 0.61 0.41 0.52	13 CS
	0.39 0.61 0.57 0.70 0.71 0.79 0.80 0.71 0.81 1.00 0.52 0.21 0.14 0.26 0.11 -0.11 0.54 0.10 0.61 0.56 0.52 0.31 0.64 0.21 0.56 0.18 0.39	SO
	0.93 0.89 0.52 0.28 0.75 0.64 0.00 0.29 0.23 0.52 1.00 0.81 0.84 0.05 0.23 0.12 0.52 -0.60 0.52 0.49 0.49 0.22 0.55 0.55 0.60 0.14 0.41	OCS
	0.43 0.57 0.18 -0.30 0.57 0.61 -0.10 0.14 0.14 0.21 0.81 1.00 0.86 0.20 0.14 0.21 0.61 -0.20 0.54 0.54 0.50 0.43 0.36 0.36 0.43 0.18 0.25	HCO +
	0.25 0.47 -0.11 -0.41 0.64 0.54 -0.40 -0.07 0.09 0.14 0.84 0.86 1.00 -0.09 0.32 0.18 0.36 -0.30 0.50 0.52 0.59 0.52 0.57 0.71 0.63 0.50 0.57	HCN
	0.49 0.39 0.31 -0.03 0.54 0.43 0.70 0.77 0.49 0.26 0.05 0.20 -0.09 1.00 0.83 0.89 0.60 0.40 0.83 0.77 0.60 0.66 0.54 0.20 0.49 0.49 0.43	C 2 H
	0.36 0.39 0.00 -0.11 0.61 0.36 0.30 0.39 0.38 0.11 0.23 0.14 0.32 0.83 1.00 0.86 0.18 0.10 0.75 0.67 0.63 0.67 0.71 0.68 0.67 0.75 0.71	CN
	0.57 0.45 0.21 -0.11 0.39 0.21 0.20 0.36 0.20 -0.11 0.12 0.21 0.18 0.89 0.86 1.00 0.14 -0.10 0.61 0.47 0.32 0.40 0.36 0.29 0.29 0.36 0.29	13 CN
	0.25 0.32 0.11 -0.04 0.46 0.54 0.80 0.75 0.34 0.54 0.52 0.61 0.36 0.60 0.18 0.14 1.00 0.60 0.68 0.68 0.65 0.58 0.39 0.29 0.54 0.21 0.29	HC 15 N
	-0.50 -0.41 -0.30 -0.11 0.10 0.30 0.50 0.50 0.50 0.10 -0.60 -0.20 -0.30 0.40 0.10 -0.10 0.60 1.00 0.50 0.80 0.90 1.00 0.40 0.30 0.70 0.70 0.60	HN 13 C
	0.50 0.67 0.32 0.11 0.89 0.82 1.00 0.71 0.72 0.61 0.52 0.54 0.50 0.83 0.75 0.61 0.68 0.50 1.00 0.95 0.88 0.81 0.86 0.61 0.85 0.64 0.71	C 2 S
	0.25 0.50 0.18 -0.04 0.85 0.83 0.90 0.59 0.74 0.56 0.49 0.54 0.52 0.77 0.67 0.47 0.68 0.80 0.95 1.00 0.96 0.93 0.85 0.63 0.89 0.74 0.77	H 13 CN
	0.07 0.36 -0.04 -0.15 0.81 0.76 0.70 0.49 0.65 0.52 0.49 0.50 0.59 0.60 0.63 0.32 0.65 0.90 0.88 0.96 1.00 0.96 0.88 0.77 0.96 0.85 0.88	DNC
	-0.07 0.18 -0.18 -0.34 0.67 0.61 0.50 0.38 0.54 0.31 0.22 0.43 0.52 0.66 0.67 0.40 0.58 1.00 0.81 0.93 0.96 1.00 0.77 0.74 0.89 0.88 0.85	H 15 NC
	0.29 0.57 0.14 0.15 0.93 0.79 0.70 0.46 0.72 0.64 0.55 0.36 0.57 0.54 0.71 0.36 0.39 0.40 0.86 0.85 0.88 0.77 1.00 0.82 0.95 0.82 0.93	DCO +
	0.00 0.22 -0.39 -0.33 0.64 0.39 -0.10 0.14 0.22 0.21 0.55 0.36 0.71 0.20 0.68 0.29 0.29 0.30 0.61 0.63 0.77 0.74 0.82 1.00 0.88 0.89 0.93	HNC
	0.11 0.40 -0.07 -0.07 0.85 0.72 0.60 0.45 0.61 0.56 0.60 0.43 0.63 0.49 0.67 0.29 0.54 0.70 0.85 0.89 0.96 0.89 0.95 0.88 1.00 0.88 0.95	DCN
	-0.18 0.08 -0.36 -0.33 0.61 0.43 0.10 0.14 0.41 0.18 0.14 0.18 0.50 0.49 0.75 0.36 0.21 0.70 0.64 0.74 0.85 0.88 0.82 0.89 0.88 1.00 0.96	N 2 H +
	-0.04 0.26 -0.21 -0.15 0.75 0.57 0.30 0.25 0.52 0.39 0.41 0.25 0.57 0.43 0.71 0.29 0.29 0.60 0.71 0.77 0.88 0.85 0.93 0.93 0.95 0.96 1.00	N 2 D +
	CS	34 SO	C 17 O	SO 2	H 13 CO +	HCS +	HC 18 O +	C 34 S	13 CS	SO	OCS	HCO +	HCN	C 2 H	CN	13 CN	HC 15 N	HN 13 C	C 2 S	H 13 CN	DNC	H 15 NC	DCO +	HNC	DCN	N 2 H +	N 2 D +
	Figure 13.

Table 4 .

 4 Opacity of the HCN, H 13 CN, and N2H + 1 -0 lines toward our Class I source sample. '...' indicates that the corresponding observations were insufficient (see Table3) for that species and source.

	Species				Opacity
		I-04302 I-04295 I-04365 I-04016 I-04166 I-04169 I-04181
	HCN	≈ 1	≈ 1	1	≈ 1	1	1	1
	H 13 CN	... (a)	...	1	1	1	1	1
	N2H +	...	1	1	1	1	1	1

a

http://cassis.irap.omp.eu

http://kida.obs.u-bordeaux1.fr/

ACKNOWLEDGMENTS

We thank the anonymous reviewer for a careful reading and constructive comments that helped us to improve the paper. This paper make use of IRAM 30M data, project ID: 014-19. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). RLG also thanks the IRAM staff and in particular the IRAM AOD A. Ritacco, and the two IRAM operators V. Peula and F. Damour for their support to run the telescope during the IRAM shift. J.H. acknowledges support from the National Science Foundation Graduate Research Fellowship under grant No. DGE-1144152. C.J.L. acknowledges funding from the National Science Foundation Graduate Research Fellowship under grant No. DGE-1745303.

Software:

CASSIS: http://cassis.irap.omp.eu, CLASS/GILDAS: http:www.iram.fr/IRAMFR/GILDAS [START_REF] Pety | SF2A-2005: Semaine de l'Astrophysique Francaise[END_REF][START_REF] Team | GILDAS: Grenoble Image and Line Data Analysis Software[END_REF], Pandas (McKinney 2010), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), SciPy (Virtanen et al. 2020), Seaborn (Waskom et al. 2018) .

APPENDIX

A. OPTICAL DEPTH EFFECTS Line optical depth can be assessed when several isotopologues of the same line are detected or when the line possesses a hyperfine pattern that is resolved by the observations. In this appendix, we investigate the line opacity of some of the molecular species we detected possessing one of these characteristics and for which puzzling results were found. For instance, the brightest lines of CS, HCN, and HCO + occur in different sources than their respective fainter isotopologues. Similarly, the brightest H 13 CO + and HC 18 O + lines also occur in different sources. Furthermore, it seems that the more abundant isotopologues H 13 CO + and CS are brighter in the source with the second smallest disk-to-envelope ratio, I-04169, whereas the rarer ones are brighter in the one with the smallest ratio, I-04181. All these results could be due to optical depth effects coming either from the lines themselves, or even from the dust continuum. Dust opacity estimate being out of the scope of this work, we focus here on line opacity study.

A.1. Isotopologue Ratios

To investigate the opacity of the HCO + 1 -0 and H 13 CO + 1 -0 lines, we derived their isotopologue line flux ratios with the HC 18 O + 1 -0 line, which is expected to be thin. As depicted in Fig. 15, across the seven Class I sources of our sample, HCO + /HC 18 O + line ratio does not exceed 60. Across the five sources of our sample where HC 18 O + is detected, the H 13 CO + /HC 18 O + ratio varies from 11.2 (in I-04016) to 6.4 (in I-04181) (see Table 3), which leads to an HCO + /HC 18 O + ratio varying from 762 (in I-04016) to 436 (in I-04181), assuming the local ISM 12 C/ 13 C ratio of 68 ± 15 [START_REF] Milam | [END_REF]Asplund et al. 2009;Manfroid et al. 2009). In the ISM, the 16 O/ 18 O ratio is found to be 557 ± 30 [START_REF] Wilson | [END_REF], close to the Solar System value of 530 in the Solar wind (McKeegan et al. 2011), 511 ± 10 in the photosphere of the Sun (Ayres et al. 2013), and ≈ 500 in comets (Jehin et al. 2009;Bockelée-Morvan et al. 2012) and meteorites (Lodders 2003). In the diffuse ISM, the 16 O/ 18 O ratio is found to be a bit higher, with a value of 672 ± 110 derived from the HCO + /HC 18 O + ratio (Lucas & Liszt 1998). Hence, these results confirm that HCO + 1 -0 is highly thick. On the contrary, the H 13 CO + 1 -0 line seems to be rather thin except in one source, I-04181. However, variation in the 12 C/ 13 C and/or 16 O/ 18 O isotopic ratios in Class I YSOs cannot be excluded. Higher spatial resolution observations are required to disentangle between these two scenarios, to determine whether the lines are tracing the same regions and physical components of the sources.

Similarly, in Fig. 16 we compare the CS/C 34 S, CS/ 13 CS, and C 34 S/ 13 CS, isotopologue line flux ratios with the elemental 32 S/ 34 S= 24.4 ± 5.0 in the vicinity of the Sun (Chin et al. 1996) and 12 C/ 13 C= 68 ± 15 in the local ISM [START_REF] Milam | [END_REF]Asplund et al. 2009;Manfroid et al. 2009). The CS/ 13 CS and CS/C 34 S isotopologue ratios are found to be lower than the elemental ratios in 5/7 and 6/7 sources, respectively, from a factor of two up to one order of magnitude in I-04181, the source with the smallest disk-to-envelope ratio. The isotopologue ratios are only found to be compatible with elemental ratios in I-04016, for both of them, and in I-04302, for the CS/ 13 CS. As for the C 34 S/ 13 CS ratios in all seven sources, they are all within a factor of two of the assumed elemental isotopic ratios, which is a very reasonable factor of uncertainty due to fractionation. Thus, it is very likely that both fainter isotopologues C 34 S and 13 CS are optically thin in our source sample. Hence, the fact that the brightest lines of CS appear in a different source than its fainter isotopologues is indeed likely due to optical depth effects.

A.2. Resolved hyperfine structure HCN, H 13 CN, and N 2 H + lines split into a resolved set of hyperfine (hf) components, whose relative opacities are given by their Einstein coefficient for spontaneous decay and level degeneracy. In the LTE approximation, assuming that for each species all of the hf lines have the same excitation temperature, they can be used to infer the total opacity τ of the hf multiplet. Each hf component k at rest frequency ν k has an opacity τ k = r k τ that scales accordingly with its relative intensity r k . Assuming Gaussian profiles, the total opacity can thus be written as follows:

where the sum is over all hf components, v 0 is the systemic source velocity, and v k the velocity shift of a particular hf component k. Considering these assumptions, we derived the line opacity of the HCN, H 13 CN, and N 2 H + 1 -0 lines