Cyrielle Garcia

Cécile Duby

Daniel D Catheline

Pablo Gutierrez Toral

Laurence Bernard

Philippe P Legrand

Vincent V Rioux
email: vincent.rioux@agrocampus-ouest.fr

Synthesis of the suspected trans-11,cis-13 conjugated linoleic acid isomer in ruminant mammary tissue by FADS3-catalyzed ∆13-desaturation of vaccenic acid

Keywords: vaccenic acid, trans-11, cis-13 conjugated linoleic acid, fatty acid desaturase 3, stearoyl-CoA desaturase

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

INTRODUCTION

Conjugated linoleic acid isomers consist of a group of octadecadienoic acids (18:2) with double bond pairs differing by position and geometrical configuration [START_REF] Palmquist | Biosynthesis of conjugated linoleic acid in ruminants and humans[END_REF]. In ruminant milk, rumenic acid (cis-9,trans-11 CLA) is by far (70-90% of total CLA isomers) the most predominant isomer [START_REF] Lock | Modifying milk fat composition of dairy cows to enhance fatty acids beneficial to human health[END_REF]. Rumenic acid derives partly from ruminal biohydrogenation of linoleic acid (18:2n-6) and mainly from endogenous host Δ9-desaturation of vaccenic acid (trans-11-18:1 or VA; [START_REF] Griinari | Conjugated linoleic acid is synthesized endogenously in lactating dairy cows by Delta(9)-desaturase[END_REF], itself produced as an intermediate from linoleic and α-linolenic acid (18:3n-3) ruminal hydrogenation [START_REF] Harfoot | Lipid metabolism in the rumen[END_REF][START_REF] Jenkins | Board-invited review: Recent advances in biohydrogenation of unsaturated fatty acids within the rumen microbial ecosystem[END_REF].

In studies quantifying milk CLA, the trans-7,cis-9 [START_REF] Corl | Trans-7,cis-9 CLA is synthesized endogenously by delta9-desaturase in dairy cows[END_REF][START_REF] Piperova | Duodenal and milk trans octadecenoic acid and conjugated linoleic acid (CLA) isomers indicate that postabsorptive synthesis is the predominant source of cis-9-containing CLA in lactating dairy cows[END_REF]Shingfield et al., 2008a) and the trans-11,cis-13 [START_REF] Sehat | Improved separation of conjugated fatty acid methyl esters by silver ion-high-performance liquid chromatography[END_REF][START_REF] Kraft | Differences in CLA isomer distribution of cow's milk lipids[END_REF][START_REF] Luna | Conjugated linoleic acid in processed cheeses during the manufacturing stages[END_REF] are described as the second most prevalent isomers. The trans-7,cis-9 CLA is thought to be exclusively synthesized by the action of tissue Δ9-desaturase on trans-7-18:1, which is produced in the rumen [START_REF] Corl | Trans-7,cis-9 CLA is synthesized endogenously by delta9-desaturase in dairy cows[END_REF]. By contrast, the trans-11,cis-13 CLA is for now suspected to be exclusively of rumen origin [START_REF] Lock | Modifying milk fat composition of dairy cows to enhance fatty acids beneficial to human health[END_REF]. Ruminal biohydrogenation (Figure 1) can indeed convert dietary α-linolenic acid (cis-9,cis-12,cis-15-18:3, or 18:3n-3) to cis-9,trans-11,cis-15-conjugated triene (Harfoot and Hazlewood, 1997;[START_REF] Jenkins | Board-invited review: Recent advances in biohydrogenation of unsaturated fatty acids within the rumen microbial ecosystem[END_REF], which is subsequently metabolized to trans-11,cis-15-18:2 by incomplete bacteria hydrogenation. The last step corresponding to the isomerization of trans-11,cis-15-18:2 to trans-11,cis-13 CLA has been shown to be catalyzed in vitro by specific Butyrivibrio fibrisolvens strains [START_REF] Fukuda | Evaluation and characterization of bacterial metabolic dynamics with a novel profiling technique, real-time metabolotyping[END_REF] but still remains speculative in vivo. In addition to milk, the minor presence of the trans-11,cis-13 CLA isomer has been described in the digesta (Shingfield et al., 2008b;[START_REF] Toral | Effect of fish oil and sunflower oil on rumen fermentation characteristics and fatty acid composition of digesta in ewes fed a high concentrate diet[END_REF][START_REF] Shen | Trans-18:1 and CLA isomers in rumen and duodenal digesta of bulls fed n-3 and n-6 PUFA-based diets[END_REF] and subcutaneous adipose tissue [START_REF] Dugan | Comparing subcutaneous adipose tissue in beef and muskox with emphasis on trans 18:1 and conjugated linoleic acids[END_REF][START_REF] Petri | Subcutaneous adipose Fatty Acid profiles and related rumen bacterial populations of steers fed red clover or grass hay diets containing flax or sunflower-seed[END_REF] of ruminants.

So far, this specific trans-11,cis-13 CLA isomer has not drawn too much attention, for the following reasons. First, its separation from other CLA isomers present in biological samples (dairy products) remains a considerable analytical challenge [START_REF] De La Fuente | Chromatographic techniques to determine conjugated linoleic acid idomers[END_REF][START_REF] Christie | A practical guide to the isolation, analysis and identification of conjugated linoleic acid[END_REF]. Second, no standard of the trans-11,cis-13 CLA is commercially available [START_REF] Sehat | Identification of conjugated linoleic acid isomers in cheese by gas chromatography, silver ion high performance liquid chromatography and mass spectral reconstructed ion profiles. Comparison of chromatographic elution sequences[END_REF][START_REF] Sehat | Improved separation of conjugated fatty acid methyl esters by silver ion-high-performance liquid chromatography[END_REF], whereas a commercial mixture of cis-11,trans-13, cis-11,cis-13, and trans-11,trans-13 CLA is available [START_REF] Rioux | Trans-vaccenate is Delta13-desaturated by FADS3 in rodents[END_REF]. Third, great varia-tions of its amount (0.07-8.00% of total isomers) are reported [START_REF] Lock | Modifying milk fat composition of dairy cows to enhance fatty acids beneficial to human health[END_REF]. Finally, no function is known for this CLA isomer although it might be physiologically important like other isomers such as cis-9,trans-11 or trans-10,cis-12 [START_REF] Wahle | Conjugated linoleic acids: Are they beneficial or detrimental to health?[END_REF]).

Yet, in rodents, the fatty acid desaturase 3 (FADS3) gene was recently shown to code for an enzyme able to catalyze the Δ13-desaturation of VA [START_REF] Rioux | Trans-vaccenate is Delta13-desaturated by FADS3 in rodents[END_REF]. Structural characterization of the Δ11,13 CLA produced by this FADS3-catalyzed desaturation strongly suggested that it was the trans-11,cis-13 CLA isomer. Because the FADS3 gene has been described in the genomic sequence of the FADS (fatty acid desaturase) gene cluster in cattle [START_REF] Park | Novel fatty acid desaturase 3 (FADS3) transcripts generated by alternative splicing[END_REF][START_REF] Taniguchi | Genomic structural analysis of porcine fatty acid desaturase cluster on chromosome 2[END_REF], we hypothesized in the present study that, by analogy with the dual origin of rumenic acid, an alternative and direct host FADS3-catalyzed Δ13desaturation of VA (Figure 1) may therefore co-exist with incomplete ruminal biohydrogenation of 18:3n-3 to explain the final presence of the trans-11,cis-13 CLA in ruminant milk fat.

MATERIALS AND METHODS

Bovine Mammary Epithelial Cell Culture and Incubation with VA

Bovine MAC-T and BME-UV mammary epithelial cells were kindly provided by F. Dessauge (UMR Pegase, INRA Saint-Gilles, France). The BME-UV cells were grown in 5:3:2 (vol/vol/vol) mixture of Dulbecco's modified Eagle medium-F12/RPMI-1640/ NCTC35 supplemented with 10% fetal bovine calf serum, 1% lactose, 0.5% lactalbumin hydrolysate, 1.2 mM l-glutathione reduced, 10 μg/mL l-ascorbic acid, 5 μg/mL of holo-transferrin, 1 μg/mL of hydrocortisone and insulin, 0.5 μg/mL of progesterone, and 1% penicillin/streptomycin. The MAC-T cells were grown in Dulbecco's modified Eagle medium supplemented with 10% fetal bovine calf serum, 5 μg/mL insulin, and 1% penicillin/streptomycin.

The ability of the bovine FADS3 protein to Δ13desaturate VA was investigated by incubating the cells with VA (Sigma, Saint-Quentin-Fallavier, France) and [1-13 C]-VA (chemical purity 95%, isotopic enrichment 99%, synthesized by the Commissariat à l'Energie Atomique, CEA, Gif-sur-Yvette, France) [START_REF] Loreau | Large-scale preparation of (9Z,12E)-[1-(13)C]-octadeca-9,12-dienoic acid[END_REF]. Incubations with cis-11-18:1, cis-9-18:1, and cis-9,trans-11 CLA were also realized. Fatty acid (FA) albuminic complexes were prepared as described [START_REF] Rioux | Myristic acid, unlike palmitic acid, is rapidly metabolized in cultured rat hepatocytes[END_REF]. The final FA concentration of the incubation medium was 200 μM (unless indicated otherwise) and albumin concentration was 150 μM, leading to a FA:albumin molar ratio of 1.3:1. The incubation was initiated by replacing the culture medium with the FAcontaining medium. Incubation was then carried out for 24 h at 37°C in a 5% CO 2 atmosphere.

Fatty Acid Analysis

Cellular lipids were extracted with hexane-isopropanol (3:2, vol/vol) after acidification with 1 mL HCl 3 M, as described [START_REF] Rioux | Influence of the cis-9, cis-12 and cis-15 double bond position in octadecenoic acid (18:1) isomers on the rat FADS2-catalyzed Δ6-desaturation[END_REF]. After saponification for 15 min at 70°C with 1 mL of 0.5 M NaOH in methanol, FA were methylated with BF 3 (14% in methanol) at 70°C for 10 min. The FAME were extracted with pentane and analyzed by GC-MS using an Agilent Technologies 7890N (Bios Analytic, Santa Clara, CA) with a bonded silica capillary column (60 m × 0.25 mm; BPX 70; SGE). Helium was used as carrier gas (average velocity 24 cm/s). The column temperature program started at 150°C, gradually increased at 2°C/min to 220°C, and held at 220°C for 10 min. Mass spectra were recorded with an Agilent Technologies 5975C inert MSD with triple axis detector.

The mass spectrometer was operated under electron impact ionization conditions (electron energy 70 eV, source temperature 230°C). Data were obtained in scan mode with a mass range of m/z 50-550 atomic mass units (amu). Identification of the FAME was based upon retention times (R t ) obtained for methyl ester of authentic standards, when available. The commercial FA used as standards were rumenic acid (cis-9,trans-11 CLA), a mixture of trans13,trans13 CLA (Matreya,Pleasant Gap,PA), and pure trans-11,trans13 CLA (Larodan, Malmö, Sweden). The National Institute of Standards and Technology database (version 2.0) was also sometimes used to identify unknown FA.

Identification of CLA After Their Conversion

To confirm the identity of CLA isomers, FAME were first pre-concentrated by separation on preparative silver ion thin layer chromatography using silicagel H impregnated with 10% (wt/wt) AgNO 3 . Separation was performed according to the number and configuration of double bonds, using a mixture of 90/10 (vol/vol) hexane/diethylether for development. At the end of the chromatographic runs, the plates were sprayed with a solution of 10 mg per 100 mL of primuline in acetone (Sigma, Saint-Quentin-Fallavier, France) and visualized under UV light. Bands corresponding to CLA according to a standard of methylated rumenic acid, were scraped off and transferred into a test tube. Collected CLA were then converted to 4,4-dimethyloxazoline (DMOX) derivatives [START_REF] Guillou | Conversion of hexadecanoic acid to hexadecenoic acid by rat Delta 6-desaturase[END_REF][START_REF] Rioux | Trans-vaccenate is Delta13-desaturated by FADS3 in rodents[END_REF][START_REF] Rioux | Influence of the cis-9, cis-12 and cis-15 double bond position in octadecenoic acid (18:1) isomers on the rat FADS2-catalyzed Δ6-desaturation[END_REF] and analyzed by GC-MS in electron impact ionization mode. Briefly, a mild 2-step reaction was used for preparing DMOX. The FAME in dichloromethane (1 mL) were converted to corresponding amides by incubation at 45°C for 30 min with 200 μL of 2-amino-2-methyl-1-propanol and 100 μL of potassium tert-butoxide in tetrahydrofuran (2:1 vol/vol). The resulting 2-(methylpropanol) amides were isolated by partition between hexane/diethylether (50/50 vol/vol) and water, and then converted to DMOX by treatment with trifluoroacetic anhydride under mild conditions (50°C for 45 min). The column, gas vector, and detector were similar to those used for the analysis of FAME (see above).

Identification of CLA After Their Reaction

To further confirm the identity of minor CLA isomers, thin layer chromatography pre-concentrated CLA were also converted to Diels-Alder adducts [START_REF] Berdeaux | Large-scale synthesis of methyl cis-9, trans-11-octadecadi-enoate from methyl ricinoleate[END_REF][START_REF] Dobson | Identification of conjugated fatty acids by gas chromatography-mass spectrometry of 4-methyl-1,2,4-triazoline-3,5dione adducts[END_REF] by reaction with 4-methyl-1,2,4triazoline-3,5-dione (MTAD). The CLA in dichloromethane were mixed with MTAD (1 mg/mL, Santa Cruz Biotechnology, Santa Cruz, CA) at 0°C for 10 s. The reaction was stopped immediately by addition of a 5-fold excess of 1,3-hexadiene (Santa Cruz Biotechnology). The MTAD adducts were analyzed by GC-MS using the same BPX70 column as for FAME analysis but with a temperature gradient starting at 150°C, then gradually increasing at 10°C/min to 260°C, and holding at 260°C for 90 min.

RNA Isolation and Quantitative Reverse-Transcription PCR

Total RNA was extracted from cells with Trizol (Invitrogen). The RNA was converted into cDNA with the IScript cDNA synthesis kit (Bio-Rad, Hercules, CA). The RNA was quantified by Real-Time PCR with the Taqman Universal PCR Master Mix (Applied Biosystem) containing 40 ng of retrotranscripted RNA, 0.5 μM primers, and 0.25 μM Taqman probe. Primers and 5 -FAM/TAMRA-3 Taqman probes specific to caprine and bovine desaturase genes are described in Table 1. Amplification was performed by a CFX 96 real-time PCR Detection System (Bio-Rad) over 40 cycles of 95°C for 15 s and 60°C for 1 min [START_REF] Ramakers | Assumption-free analysis of quantitative real-time polymerase chain reaction (PCR) data[END_REF]. The 18S gene expression was quantified as an endogenous control using the 18S RNA Control kit Yakima-Yellow Eclipse Dark Quencher (Eurogentec, Liège, Belgium) and relative gene expression was determined from the cycle thresholds (Ct) using the 2 -ΔCt method.

Lactating Goat Tissue Samples

Samples of the mammary secretory tissue and liver were collected from 8 multiparous nonpregnant lactating Alpine goats. The experimental design of this study has been previously published [START_REF] Bernard | Effects of fish oil and starch added to a diet containing sunflower-seed oil on dairy goat performance, milk fatty acid composition and in vivo delta9desaturation of [13C]vaccenic acid[END_REF][START_REF] Toral | Effects of fish oil and additional starch on tissue fatty acid profile and lipogenic gene mRNA abundance in lactating goats fed a diet containing sunflower-seed oil[END_REF]. The FA composition of mammary tissue and liver samples was re-evaluated in the present study using the GC-MS methodology described above, with special focus on the suspected trans-11,cis-13 CLA. The quantitative reverse-transcription PCR was carried out on mammary tissue and liver mRNA samples using specific primers and probes for goat FADS1, FADS2, FADS3, and stearoyl-CoA desaturase 1 (SCD1; Table 1).

Expression of Results and Statistical Analysis

All data values shown represent the mean ± standard deviation or standard error of the mean for the number of replicates (n) indicated in each figure legend. Student's t-test was used to determine the statistical significance between groups, where indicated.

RESULTS AND DISCUSSION

Presence of the FADS3 Gene Within the FADS Cluster in Ruminant Species

The FADS3 gene was first described in humans, within the genomic structure of the FADS gene cluster, which also includes FADS1 (Δ5-desaturase) and FADS2 (Δ6-desaturase) genes and is located in chromosome 11 [START_REF] Marquardt | cDNA cloning, genomic structure, and chromosomal localization of three members of the human fatty acid desaturase family[END_REF]. Since then, FADS3 orthologs similarly clustered with FADS1 and FADS2 have been identified in rat chromosome 1 [START_REF] Pedrono | The fatty acid desaturase 3 gene encodes for different FADS3 protein isoforms in mammalian tissues[END_REF], pig chromosome 2 [START_REF] Taniguchi | Genomic structural analysis of porcine fatty acid desaturase cluster on chromosome 2[END_REF], baboon chromosome 14 [START_REF] Park | Novel fatty acid desaturase 3 (FADS3) transcripts generated by alternative splicing[END_REF], mouse chromosome 19 [START_REF] Nakamura | Structure, function, and dietary regulation of delta6, delta5, and delta9 desaturases[END_REF], and cattle chromosome 29 [START_REF] Park | Novel fatty acid desaturase 3 (FADS3) transcripts generated by alternative splicing[END_REF][START_REF] Taniguchi | Genomic structural analysis of porcine fatty acid desaturase cluster on chromosome 2[END_REF]. In the present study, genomic database searches in other ruminant species revealed the presence of the FADS1/ FADS2/FADS3 gene cluster with similar order and direction in caprine (chromosome 29, NC_022321.1, Capra hircus) and ovine (chromosome 21, NC_019478.1, Ovis aries) species, FADS3 being located tail-to-tail with FADS2 [START_REF] Marquardt | cDNA cloning, genomic structure, and chromosomal localization of three members of the human fatty acid desaturase family[END_REF]. In ruminant species, the exon/intron organization of the FADS3 gene (12 exons) was also similar to that of other mammals, Table 1. Primers and probes used for quantification of the mRNA levels of fatty acid desaturase 3 (FADS3) and stearoyl-CoA desaturase 1 (SCD1) genes by real-time reverse-transcription PCR

Gene

Forward primer 5 -3 Reverse primer 5 -3 5 -FAM/TAMRA-3 probe

Bovine SCD1 TTGTCCACTTTCTCCTGCTG GTAGCCATCACTGCCTCTGA TGGACAGGCAAGCCTCTTCTGTG Bovine FADS3 GAAAATCTGGCCTACATGCTGG GGCGTAGAAGCTGGCGG TGCATGCAGTGGACGGACTTGCTC Goat FADS1 CCGCCTCCCTCTGCAAG GGCATTGCCCCCAAGCT AGGGTTATGGCACCAAAGCAGACGG Goat FADS2 GCACGGACAAGTGGCTGG TGCCGGCTGGACCGTT No probe (Sybr Green) Goat FADS3 GAAAATCTGGCCTACATGCTGG GGCGTAGAAGCTGGCGG TGCATGCAGTGGACGGACTTGCTC Goat SCD1
TTTCTTCTCTCACGTGGGTTGG AGCGTAGCACCCTTTTCTCTGGA TGCGCAAACACCCAGCT although the first 2 exons are lacking in the goat draft genome sequence of the latest assembly (NC_022321.1).

Concerning the FADS3 nucleotide sequence and predicted AA sequence (Figure 2), the alignment of goat (XM_013975910.1), sheep (XM_012146315.1), cattle (NM_001083691.2), human (NM_021727.4), and rat (NM_173137.1) sequences exhibited a very high degree of identity (>85%).

Identification and Quantification of the Suspected Trans-11,Cis-13 CLA in Cow Milk

We have previously used GC-MS analysis of FAME and DMOX derivatives to accurately discriminate the suspected trans-11,cis-13 CLA produced by recombinant rat FADS3 [START_REF] Rioux | Trans-vaccenate is Delta13-desaturated by FADS3 in rodents[END_REF] from other CLA isomers and more specifically from the commercially available cis-11,trans-13, cis-11,cis-13, and trans-11,trans-13 CLA isomers. To conclusively identify and quantify the presence of the trans-11,cis-13 CLA isomer in ruminant milk, we analyzed the FAME profile from the European certified reference material anhydrous milk fat (BCR164 obtained from Sigma), which corresponds to a pooled sample of cow milks. Figure 3A shows that rumenic acid (cis-9,trans-11 CLA), probably co-eluting with trans-7,cis-9 CLA and trans-8,cis-10 CLA in our experimental conditions (R t 27.5 min, m/z 294 amu), represented 1.1 ± 0.2% of total FA in this sample. A well-isolated small peak (0.06 ± 0.01% of total FA) with a molecular FAME ion at m/z 294 amu and with the exact same retention time (R t 28.2 min) as the highly probable trans-11,cis-13 CLA produced by recombinant rat FADS3 when incubated with VA [START_REF] Rioux | Trans-vaccenate is Delta13-desaturated by FADS3 in rodents[END_REF], was also detected in this sample (Figure 3A).

This peak was also compared with commercially available standards of Δ11,13 CLA (cis,trans, cis,cis, and trans,trans but not the trans,cis; Figure 3A). By comparing the GC retention time of the desaturated product of VA by FADS3 with cow milk fat and these standard (Figure 3A), we first excluded the cis,cis (+0.6 min) and trans,trans (+0.9 min) isomers with retention times significantly higher than the unknown Δ11,13 CLA, whereas the cis,trans isomer had a slightly lower retention time (-0.2 min). This result suggested that the cis,trans isomer could also be discarded. Therefore, the GC retention times offer several arguments to reasonably conclude that this FA is the trans-11,cis-13 CLA isomer.

In addition, after conversion to DMOX derivatives, the mass spectrum of this CLA present in cow milk fat (Figure 3B) showed a pattern completely similar to the suspected trans-11,cis-13 CLA already characterized in cells expressing rat FADS3 and incubated with VA [START_REF] Rioux | Trans-vaccenate is Delta13-desaturated by FADS3 in rodents[END_REF]. Indeed, in addition to the DMOX characteristic intense ions at m/z 113 and 126 amu and molecular ion at m/z 333 amu, a gap of 12 amu was detected between ions containing 10 and 11 carbons (224 and 236 amu) and between ions containing 12 and 13 carbons (250 and 262 amu), showing that this FA is indeed a Δ11,13 CLA (Figure 3B). More specifically, the ion at m/z 250 amu of this DMOX derivative was significantly less abundant than the ion at m/z 262 amu (Figure 3B), which has already been used to specifically differentiate the suspected trans-11,cis-13 CLA from the commercial cis-11,trans-13 CLA isomer (Figure 3C) with the nearest R t of 28.0 min (Figure 3A; [START_REF] Sehat | Identification of conjugated linoleic acid isomers in cheese by gas chromatography, silver ion high performance liquid chromatography and mass spectral reconstructed ion profiles. Comparison of chromatographic elution sequences[END_REF][START_REF] Sehat | Improved separation of conjugated fatty acid methyl esters by silver ion-high-performance liquid chromatography[END_REF][START_REF] Rioux | Trans-vaccenate is Delta13-desaturated by FADS3 in rodents[END_REF]. Therefore, these additional results strongly suggest that this milk FA is the trans-11,cis-13 CLA, representing 0.06% of total FA in this sample.

Expression of FADS3 mRNA and Synthesis of Suspected Trans-11,Cis-13 CLA

To investigate the potential synthesis of the suspected trans-11,cis-13 CLA by the FADS3-catalyzed Δ13-desaturation of VA in ruminant mammary tissue, we first quantified the relative mRNA expression of FADS3 and SCD1 in bovine mammary MAC-T and BME-UV epithelial cells. Figures 4A and4B show that FADS3 and SCD1 mRNA were expressed in these 2 bovine mammary cell models, FADS3 mRNA being surprisingly 10 times more abundant than SCD1 mRNA. When incubated with 200 μM VA, FADS3 mRNA levels were not modified, whereas SCD1 mRNA levels decreased, although not significantly (Figures 4A and4B). This observation is consistent with the described SCD1 gene expression modulation through the associate gene regulatory protein SREBP-1 downregulation [START_REF] Keating | Characterization and regulation of the bovine stearoyl-CoA desaturase gene promoter[END_REF][START_REF] Minville-Walz | Distinct regulation of stearoyl-CoA desaturase 1 gene expression by cis and trans C18:1 fatty acids in human aortic smooth muscle cells[END_REF][START_REF] Jacobs | Effects of short-and long-chain fatty acids on the expression of stearoyl-CoA desaturase and other lipogenic genes in bovine mammary epithelial cells[END_REF].

The MAC-T and BME-UV mammary epithelial cells were then incubated for 24 h with media supplemented with VA or [1-13 C]-VA, and a control medium without FA. Rumenic acid produced by Δ9-desaturation of VA was easily detected (1.88 ± 0.44% of FA in BME-UV cells and 2.22 ± 1.05% of FA in MAC-T cells, n = 3) after incubation with VA (Figures 4B and4D). A smaller peak (0.49 ± 0.24% of FA in BME-UV cells and 0.82 ± 0.25% of FA in MAC-T cells, n = 3) with R t 28.2 min and FAME molecular ion at m/z 294 amu (incubation with 12 C-VA) or m/z 295 amu (incubation with 13 C-VA), similar to the suspected trans-11,cis-13 CLA previously characterized in milk fat sample (Figure 3A), was also specifically detected (Figures 4C Figure 2. Alignment of the AA sequences of goat, sheep, cow, rat, and human fatty acid desaturase 3 (FADS3) deduced from their available nucleotide sequences. The alignment of goat (Capra, XM_013975910.1), sheep (Ovis, XM_012146315.1), cattle (Bos, NM_001083691.2), human (Homo, NM_021727.4), and rat (Rattus NM_173137.1) sequences exhibited a very high degree of identity (>85%) (http://multalin.toulouse. inra.fr/multalin/). -11,trans-13; cis-11,cis-13; and trans-11,trans-13 CLA as FAME standards. (B) Electron-impact mass spectrum of the 4,4-dimethyloxazoline (DMOX) derivative of the suspected trans-11,cis-13 CLA from the European certified reference material anhydrous milk fat (BCR164). (C) Electron-impact mass spectrum of the DMOX derivative of the commercial cis-11,trans-13 CLA. and 4D). Indeed, when incubated mammary epithelial cells were incubated with cis-9-18:1, cis-11-18:1 or cis-9,trans-11-18:2, this specific FA was not detected (data not shown).

To clearly identify this CLA isomer obtained from incubation of VA with bovine mammary epithelial cells, it was converted both to MTAD adduct and DMOX derivative (Figure 5), followed by GC-MS analysis. The mass spectrum of the MTAD adduct (Figure 5A) of the suspected trans-11,cis-13 CLA contained a molecular ion at m/z 407 amu and 3 very characteristic ions at m/z 350, 318, and 222 amu, identifying this CLA as a Δ11,13 isomer, whereas rumenic acid MTAD adduct (Figure 5B) showed a similar molecular ion at m/z 407 amu but 3 other characteristic ions at m/z 322, 290, and 250 amu [START_REF] Berdeaux | Large-scale synthesis of methyl cis-9, trans-11-octadecadi-enoate from methyl ricinoleate[END_REF]. Figure 5C shows again that the DMOX mass spectrum of the suspected trans-11,cis-13 CLA was characterized by 2 gaps between ions containing 10 and 11 carbons (224 and 236 amu) and between ions containing 12 and 13 carbons (250 and 262 amu), locating the double bonds. More specifically, its ion at m/z 250 amu was less abundant than its ion at m/z 262 amu, as already shown in milk fat sample (Figure 3B). Again, these results offer several arguments to reasonably conclude that this FA is the trans-11,cis-13 CLA isomer.

To conclude, incubation of mammary epithelial cells with VA resulted in a Δ9-desaturation index [assessed as RA/(VA+RA+trans-11,cis-13 CLA)] of 4.4 and 5.7%, and a Δ13-desaturation index [assessed as trans-11,cis-13 CLA/(VA+RA+trans-11,cis-13 CLA)] of 1.1 and 2% in BME-UV and MAC-T cells, respectively. (224 and 236 amu) and between ions containing 12 and 13 carbons (250 and 262 amu), locating the double bonds in position 11-12 and 13-14. Specifically, its ion at m/z 250 amu was less abundant than its ion at m/z 262 amu, as indicated by the arrow, and as shown in the milk fat sample (Figure 3B), which strongly suggests that this fatty acid is the trans-11,cis-13 CLA isomer.

FADS3 mRNA and Suspected Trans-11,Cis-13 CLA Concomitant Presence in Goat Mammary Tissue and Absence in Goat Liver

To analyze the potential concomitant presence of FADS3 mRNA and of the suspected trans-11,cis-13 CLA in vivo in another ruminant model, we took advantage of tissue samples collected from a study carried out with lactating goats, which has already been published [START_REF] Bernard | Effects of fish oil and starch added to a diet containing sunflower-seed oil on dairy goat performance, milk fatty acid composition and in vivo delta9desaturation of [13C]vaccenic acid[END_REF][START_REF] Toral | Effects of fish oil and additional starch on tissue fatty acid profile and lipogenic gene mRNA abundance in lactating goats fed a diet containing sunflower-seed oil[END_REF]. We re-analyzed the mammary tissue and liver FA compositions, paying particular attention to the suspected trans-11,cis-13 CLA and additionally quantified FADS3 mRNA. Figure 6A shows that the suspected trans-11,cis-13 CLA (R t 28.2 min) was detected in goat mammary tissue, representing 0.15 ± 0.05% of total FA (n = 8), whereas its trace amounts were not quantifiable in the liver (Figure 6B). Figure 6C additionally shows that FADS3 and SCD1 mRNA levels were significantly higher in the mammary secretory tissue than in the liver of goats. Therefore, the higher presence of the suspected trans-11,cis-13 CLA in goat mammary tissue is consistent with FADS3 mRNA being expressed in this tissue, compared with very low levels of both suspected trans-11,cis-13 CLA and FADS3 mRNA in goat liver. However, in contrast to the results obtained with MAC-T and BME-UV epithelial cells (Figure 4A and 4B), FADS3 mRNA was about 5,000 times less abundant than SCD1 mRNA in goat mammary tissue (Figure 6C). This result suggested first that FADS3 may be less expressed in goats than in cows. Of course, FADS3 mRNA levels may also be different in mammary tissue samples collected in vivo and in mammary cultured cells. As previously shown [START_REF] Toral | Effects of fish oil and additional starch on tissue fatty acid profile and lipogenic gene mRNA abundance in lactating goats fed a diet containing sunflower-seed oil[END_REF], FADS1 and FADS2 mRNA levels were conversely significantly higher in the liver than in the mammary tissue (Figure 6C).

In a second approach to get further information regarding the effect of dietary PUFA on the amount of the suspected trans-11,cis-13 CLA and FADS3 mRNA in goat mammary tissue, we took into consideration that the goats were fed during 21 d with 2 experimental diets containing exactly the same amounts of 18:3n-3 but different amounts of other PUFA [START_REF] Toral | Effects of fish oil and additional starch on tissue fatty acid profile and lipogenic gene mRNA abundance in lactating goats fed a diet containing sunflower-seed oil[END_REF]. Half of the goats (n = 4) were fed 90 g of sunflower oil per day (sunflower oil diet including 27.3 g/ kg of DM of 18:2n-6 and 6.2 g/kg of 18:3n-3). Another group of n = 4 goats was fed with a combination of 60 g of sunflower oil and 30 g of fish oil (sunflower and fish oil diet including 21.6 g/kg of DM of 18:2n-6, 1.4 g/kg of 20:5n-3, 1.0 g/kg of 22:6n-3 and again 6.2 g/ kg of 18:3n-3) plus 59 g/kg of DM of starch from rolled barley. Considering now the effect of the 2 diets, Figure 7A shows that the suspected trans-11,cis-13 CLA was significantly higher (0.19 ± 0.04% vs. 0.12 ± 0.03%) in the mammary tissue of goats fed the sunflower and fish oil diet compared with goats fed the sunflower oil diet. In addition, a trend toward the increase of FADS3 mRNA was observed (Figure 7B), even if the difference was not statistically significant (P = 0.13) due to the low number of animals. Conversely, no effect of the diets has already been shown on trans-11-18:1, cis-9,trans-11 CLA (Figure 7A), and on SCD1 mRNA level [START_REF] Toral | Effects of fish oil and additional starch on tissue fatty acid profile and lipogenic gene mRNA abundance in lactating goats fed a diet containing sunflower-seed oil[END_REF], as also re-analyzed in the present study (Figure 7B). Because the 2 diets contained the same amounts of 18:3n-3, we can therefore hypothesize that its ruminal biohydrogenation potentially leading to trans-11,cis-13 CLA production may have been similar and can speculate that mammary tissue FADS3 would therefore have been involved in the elevated trans-11,cis-13 CLA level. Interestingly, the expression of the FADS3 transcript was also upregulated by dietary long-chain PUFA in baboon liver, whereas other desaturase (FADS1 and FADS2) mRNA levels were downregulated [START_REF] Reardon | Dietary longchain polyunsaturated fatty acids upregulate expression of FADS3 transcripts[END_REF]. These results are also in agreement with the positive correlation between trans-11,cis-13 CLA and trans-11-18:1 levels described in milk [START_REF] Lerch | Rapeseed or linseed in grass-based diets: effects on conjugated linoleic and conjugated linolenic acid isomers in milk fat from Holstein cows over 2 consecutive lactations[END_REF].

CONCLUSIONS

A vaccenate Δ13-desaturase activity has recently been attributed to the FADS3 enzyme in rodents [START_REF] Rioux | Trans-vaccenate is Delta13-desaturated by FADS3 in rodents[END_REF]. We hypothesized in the present study that this activity would be common to all mammals, explaining at least partly the presence of the trans-11,cis-13 CLA in ruminant milk. Indeed, this study provides several lines of evidence that, by analogy with rumenic acid, the trans-11,cis-13 CLA may originate both from ruminal biohydrogenation and from direct FADS3-catalyzed Δ13-desaturation of VA in mammary tissue.

Figure 1 .

 1 Figure 1. Suspected metabolic pathways leading to the presence of suspected trans-11,cis-13 CLA in ruminant milk. Ruminal biohydrogenation can convert dietary α-linolenic acid (cis-9,cis-12,cis-15-18:3) to cis-9,trans-11,cis-15-conjugated triene, which is subsequently metabolized to trans-11,cis-15-18:2 by incomplete bacteria hydrogenation. The step corresponding to the isomerization of trans-11,cis-15-18:2 to trans-11,cis-13 CLA has been shown to be catalyzed in vitro by specific Butyrivibrio fibrisolvens strains but still remains speculative in vivo. As recently shown in rodents, this fatty acid may also come directly from vaccenate Δ13-desaturation catalyzed by fatty acid desaturase 3 (FADS3).

Figure 3 .

 3 Figure 3. Identification and quantification of the suspected trans-11,cis-13 CLA in cow milk standard sample. (A) The GC-MS profile of the European certified reference material anhydrous milk fat (BCR164) compared with FAME extracted from COS-7 cells (African green monkey fibroblasts) expressing recombinant rat fatty acid desaturase 3 (FADS3) and incubated with vaccenic acid and a commercial mixture of cis-11,trans-13; cis-11,cis-13; and trans-11,trans-13 CLA as FAME standards. (B) Electron-impact mass spectrum of the 4,4-dimethyloxazoline (DMOX) derivative of the suspected trans-11,cis-13 CLA from the European certified reference material anhydrous milk fat (BCR164). (C) Electron-impact mass spectrum of the DMOX derivative of the commercial cis-11,trans-13 CLA.

Figure 4 .

 4 Figure 4. Characterization of bovine BME-UV and MAC-T mammary epithelial cells incubated with vaccenic acid (VA). Cells were incubated for 24 h with albumin-bound VA (200 μM) or with no fatty acids (FA; 0 μM VA). The stearoyl-CoA desaturase 1 (SCD1) and fatty acid desaturase 3 (FADS3) mRNA levels were quantified in BME-UV (A) and MAC-T (B) cells. Data are presented as means ± SEM (n = 3). The FAME extracted from bovine BME-UV (C) and MAC-T (D) cells were identified and quantified by GC-MS. The identity of each important FA is indicated above its respective peak. Panels C and D are representative examples of n = 3 distinct experiments.

Figure 5 .

 5 Figure 5. Structural characterization of the CLA obtained after incubation of vaccenic acid (VA) with BME-UV epithelial cells. Electronimpact mass spectrum of the 4-methyl-1,2,4-triazoline-3,5-dione (MTAD) adduct of the suspected trans-11,cis-13 CLA (A) contains a molecular ion at m/z 407 atomic mass units (amu) and 3 very characteristic ions at m/z 350, 318, and 222 amu, identifying this CLA as a Δ11,13 isomer. Electron-impact mass spectrum of the MTAD adduct of cis-9,trans-11 CLA (B) also shows a similar molecular ion at m/z 407 amu but 3 other characteristic ions at m/z 322, 290, and 250 amu. The 4,4-dimethyloxazoline mass spectrum of the suspected trans-11,cis-13 CLA (C) was characterized by 2 gaps (boxed areas) between ions containing 10 and 11 carbons (224 and 236 amu) and between ions containing 12 and 13 carbons (250 and 262 amu), locating the double bonds in position 11-12 and 13-14. Specifically, its ion at m/z 250 amu was less abundant than its ion at m/z 262 amu, as indicated by the arrow, and as shown in the milk fat sample (Figure3B), which strongly suggests that this fatty acid is the trans-11,cis-13 CLA isomer.

Figure 7 .

 7 Figure 7. Effect of sunflower-seed oil (SO) and sunflower-seed oil and fish oil (2:1) plus starch from rolled barley (SFO) diets containing exactly the same level of dietary 18:3n-3 on the amount of trans-11-18:1, cis-9,trans-11 CLA and suspected trans-11,cis-13 CLA and the levels of fatty acid desaturase 3 (FADS3) and SCD1 mRNA in goat mammary tissue. (A) Quantification (% of FA) of trans-11-18:1, cis-9,trans-11 CLA, and suspected trans-11,cis-13 CLA after GC-MS analysis of FAME extracted from goat mammary tissue. (B) Relative expression of FADS3 and SCD1 mRNA levels in goat mammary tissue. Data are presented as means ± SD (n = 4).

Journal of Dairy Science Vol. 100 No. 1, 2017

ACKNOWLEDGMENTS

The authors thank Frédéric Dessauge (UMR Pegase, INRA Saint-Gilles, France) for the MAC-T and BME-UV cell lines. The authors also gratefully acknowledge O. Loreau (CEA, Saclay and Cadarache, France) for the supply of [1-13 C]-vaccenic acid and the Transqual project (Clermont-Ferrand, France; ANR-05-PNRA-No.5.E.24). This work was supported by the Groupe Lipides et Nutrition (GLN, Paris, France).


Eva Fernandez-Tendero

Arnaud Day

Sandrine Legros

Anouck Habrant

Simon Hawkins

Brigitte Chabbert
email: brigitte.chabbert@inra.fr

Changes in hemp secondary fiber production related to technical fiber variability revealed by light microscopy and attenuated total reflectance Fourier transform infrared spectroscopy

Interest in hemp (Cannabis sativa L.) is increasing due to the development of a new range of industrial applications based on bast fibers. However the variability of bast fiber yield and quality represents an important barrier to further exploitation. Primary and secondary fiber content was examined in two commercial hemp varieties (Fedora 17, Santhica 27) grown under contrasted sowing density and irrigation conditions. Both growing conditions and hemp varieties impact stem tissue architecture with a large effect on the proportion of secondary fibers but not primary fibers. Attenuated total reflectance infrared spectroscopy allowed the discrimination of manually-isolated native primary fibers and secondary fibers but did not reveal any clustering according to growing conditions and variety. Infrared data were confirmed by wet chemistry analyses that revealed slight but significant differences between primary and secondary fiber cell wall composition. Infrared spectroscopy of technical fibers obtained after mechanical defibering revealed differences with native primary, but not secondary fibers and also discriminated samples obtained from plants grown under different conditions. Altogether the results suggested that the observed variability of hemp technical fibers could be partially explained by i) differences in secondary fiber production and ii) differential behavior during mechanical defibering resulting in unequal separation of primary and secondary fibers.

Introduction

Hemp fibers have been traditionally used for centuries in the textile and paper industries [START_REF] Ranalli | Current status and future scenarios of hemp breeding[END_REF]. Over the last few years, these natural fibers have also been increasingly used to replace synthetic fibers for the production of more environmentally friendly materials [START_REF] Ranalli | Current status and future scenarios of hemp breeding[END_REF][START_REF] Faruk | Biocomposites reinforced with natural fibers: 2000-2010[END_REF]. In addition to their lower environmental impact, the use of natural fibers offers other advantages including lower density and good mechanical properties [START_REF] Joshi | Genomics of cellulose biosynthesis in poplars[END_REF]. However, variability in the industrial quality represents an important barrier to the large-scale exploitation of these fibers. Previous studies have suggested that differences in the mechanical properties of fibers could be related to plant growth and the production of biological material [START_REF] Thygesen | The effects of growth conditions and of processing into yarn on dislocations in hemp fibres[END_REF][START_REF] Abot | Effects of cultural conditions on the hemp (Cannabis sativa) phloem fibres: Biological development and mechanical properties[END_REF] as well as biological and/or physicochemical processing used for defibering and fiber purification [START_REF] Faruk | Biocomposites reinforced with natural fibers: 2000-2010[END_REF][START_REF] Liu | Effect of harvest time and field retting duration on the chemical composition, morphology and mechanical properties of hemp fibers[END_REF].

Hemp fibers are located in the outer stem tissues of the plant and are therefore referred to as bast fibers. Hemp bast fibers are in fact composed of two fiber types with different origins, primary bast fibers are sclerenchyma cells derived from the procambium whereas secondary bast fibers are produced by the vascular cambium [START_REF] Gorshkova | Plant fiber formation: state of the art, recent and expected progress, and open questions[END_REF]. These two fiber types show a number of differences in cell morphology and structure that could impact greatly on their suitability for different industrial uses. For example, the primary fibers with a greater cell length and higher crystalline cellulose content levels are preferred to secondary fibers for composite reinforcement [START_REF] Bonatti | Histochemical and supramolecular studies in determining quality of hemp fibres for textile applications[END_REF][START_REF] Cro ˆnier | Structure and chemical composition of bast fibers isolated from developing hemp stem[END_REF]. In practice, the separation of the two fiber types during the defibering process is difficult and secondary fibers remain attached to the primary fibers thereby affecting the overall quality of the harvested material [START_REF] Mediavilla | Influence of the growth stage of industrial hemp on the yield formation in relation to certain fibre quality traits[END_REF][START_REF] Pickering | Interfacial modification of hemp fiber reinforced composites using fungal and alkali treatment[END_REF][START_REF] Placet | Investigation of the internal structure of hemp fibres using optical coherence tomography and Focused Ion Beam transverse cutting[END_REF]. The amount of secondary fibers varies according to the position in the plant stem (greater at the stem base) and increases with plant age [START_REF] Liu | Effect of harvest time and field retting duration on the chemical composition, morphology and mechanical properties of hemp fibers[END_REF][START_REF] Cro ˆnier | Structure and chemical composition of bast fibers isolated from developing hemp stem[END_REF][START_REF] Amaducci | Influence of agronomic factors on yield and quality of hemp (Cannabis sativa L.) fibre and implication for an innovative production system[END_REF]. The quantity of secondary fibers produced by hemp plants also depends upon the genotype (variety and sex) of the plant and is modified by contrasted growing conditions (environment) [START_REF] Mediavilla | Influence of the growth stage of industrial hemp on the yield formation in relation to certain fibre quality traits[END_REF][START_REF] Sankari | Comparison of bast fibre yield and mechanical fibre properties of hemp (Cannabis sativa L.) cultivars[END_REF][START_REF] Scha ¨fer | Effect of sowing date and plant density on the cell morphology of hemp (Cannabis sativa L.)[END_REF][START_REF] Westerhuis | Sowing density and harvest time affect fibre content in hemp (Cannabis sativa) through their effects on stem weight[END_REF]. Changes in growing conditions also affect the size and wall thickness of primary fibers [START_REF] Amaducci | Influence of agronomic factors on yield and quality of hemp (Cannabis sativa L.) fibre and implication for an innovative production system[END_REF][START_REF] Scha ¨fer | Effect of sowing date and plant density on the cell morphology of hemp (Cannabis sativa L.)[END_REF] that undergo significant chemical and structural changes from the vegetative to the seed maturity stage [START_REF] Liu | Effect of harvest time and field retting duration on the chemical composition, morphology and mechanical properties of hemp fibers[END_REF][START_REF] Cro ˆnier | Structure and chemical composition of bast fibers isolated from developing hemp stem[END_REF][START_REF] Blake | In situ analysis of cell wall polymers associated with phloem fibre cells in stems of hemp, Cannabis sativa L[END_REF].

In addition to the fiber variability in planta, bast fiber quality and yield can be greatly modified by biological and/or physicochemical processing (retting, defibering, and all hydrothermal, chemical and enzymatic treatments that aim at minimizing damage during the mechanical separation of the fibers) [START_REF] Faruk | Biocomposites reinforced with natural fibers: 2000-2010[END_REF][START_REF] Liu | Effect of harvest time and field retting duration on the chemical composition, morphology and mechanical properties of hemp fibers[END_REF][START_REF] Ouajai | Morphology and structure of hemp fibre after bioscouring[END_REF][START_REF] Islam | Influence of hygrothermal ageing on the physico-mechanical properties of alkali treated industrial hemp fibre reinforced polylactic acid composites[END_REF][START_REF] Kostic | Effect of hemicelluloses and lignin on the sorption and electric properties of hemp fibers[END_REF][START_REF] Ha ¨nninen | Mechanical processing of bast fibres: The occurrence of damage and its effect on fibre structure[END_REF]. In this context a reliable method for rapidly analyzing bast fibers would greatly facilitate optimization of downstream industrial processing and applications. Infrared spectroscopy has been extensively used for the characterization of plant cell walls because of its rapid sampling and high sensitivity [START_REF] Abidi | Changes in the cell wall and cellulose content of developing cotton fibers investigated by FTIR spectroscopy[END_REF][START_REF] Lupoi | Assessment of lignocellulosic biomass using analytical spectroscopy: an evolution to high-throughput techniques[END_REF] and therefore represents a potentially interesting technique for characterizing hemp bast fibers. In hemp, Fourier transform infrared (FTIR) spectroscopy has been previously used to study the effect of enzymatic and chemical treatments, fiber dislocation and chemical composition [START_REF] Ouajai | Morphology and structure of hemp fibre after bioscouring[END_REF][START_REF] Troedec | Influence of various chemical treatments on the composition and structure of hemp fibres[END_REF][START_REF] Dai | Investigation of the dislocation of natural fibres by Fourier-transform infrared spectroscopy[END_REF][START_REF] Kabir | Mechanical properties of chemically-treated hemp fibre reinforced sandwich composites[END_REF]. Another important parameter, cellulose crystallinity, can also be assessed by FTIR spectroscopy, although this parameter is usually better assessed by X-ray diffraction (XRD) approaches [START_REF] Troedec | Influence of various chemical treatments on the composition and structure of hemp fibres[END_REF]. The development of the attenuated total reflectance technique (ATR) has also proved to be a useful approach by reducing sampling time and allowing non-destructive analysis of samples in their original environment [START_REF] Abidi | Changes in the cell wall and cellulose content of developing cotton fibers investigated by FTIR spectroscopy[END_REF].

The aim of this work was to evaluate the main sources of variability in fibers obtained from two hemp cultivars grown under different growing conditions using attenuated total reflectance-based infrared spectroscopy as a rapid and non-destructive tool.

Material and methods

Plant material

Hemp plants (Cannabis sativa L.) were grown at Marigny-le-Cha ˆtel, France (48˚24'N, 3˚44'E) by professional hemp producers. Field trials were conducted by the French Hemp Technical Institute 'Terres Inovia' (http://www.terresinovia.fr) that possesses all legal authorization and expertise necessary for hemp field studies. Two monoecïous fiber varieties (Fedora 17: F17 and Santhica 27: S27) were sown early April and cultivated under 3 different growing conditions on a chalky soil: Condition 1 (standard conditions): sowing density = 50 kg.ha -1 , no irrigation; Condition 2 (irrigation): sowing density = 50 kg.ha -1 and irrigation corresponding to 50% additional water supply calculated on the base of average rainfall in this region (326 mm, MeteoFrance station Troyes-Barberey) and provided every 10 days during 3 months before flowering stage; Condition 3 (100 kg.ha -1 ): sowing density = 100 kg.ha -1 , no irrigation. For all conditions, 120 U N.ha -1 ammonium nitrate was added under solid form at crop planting (one unit nitrogen correspond to 1 kg nitrogen provided by ammonium nitrate). Average rainfall was 319 mm for the cultivation year. Plants were harvested at the end of the flowering stage (17 th August and 22 nd August for F17 and S27, respectively).

Microscopy

For microscopic observations of the stem basal regions (between 15 and 45 cm above ground level) were isolated from 4 individual plants for each growing condition and variety. Small discs (2.5 cm thick) were then removed and rapidly frozen in liquid nitrogen and stored at-20˚C. Frozen samples were gently thawed in deionized water at 4˚C before hand-cutting thick transverse sections. The average thicknesses of bark (outer part of the stem) and wood (xylem, inner part of the stem) obtained were calculated from measurements made on thick transverse sections examined under a stereomicroscope (Stemi 2000-C, Zeiss, 30 measurements per section). The surface areas of primary and secondary fibers were determined using UV fluorescence microscopy (Axioskop, Zeiss equipped with filters at 340 nm excitation and 430 nm emission). Square images (715 X 715 μm) per section (4 samples per growing condition) were extracted and images analyzed with image-processing software (Axiovision 4, Zeiss and Image J).

Fiber infrared analysis

Fibers were manually isolated from the basal stem region (between 15 and 45 cm above ground level) of individual plants. The bark of three plants per growing condition was separated from the inner core and fiber bundles were then peeled out and the epidermis discarded as previously described [START_REF] Cro ˆnier | Structure and chemical composition of bast fibers isolated from developing hemp stem[END_REF]. The secondary fibers were separated manually from primary ones under a stereomicrosope (Zeiss) without chemical treatment and native fiber bundles were then dried at room temperature. Three stem samples (1.5 kg) per variety/growing condition were also dried at 50˚C in a forced-air oven and then mechanically defibered with a roller mill machine (Depoortere-Fibres Recherche De ´veloppement) to obtain technical fibers.

Manually-isolated and technical fibers were then characterized without grinding by infrared spectroscopy. Spectra were recorded by attenuated total reflectance (ATR) using a Nicolet 6700 spectrometer (ThermoScientific) equipped with a SmartPerformer ATR accessory. Measurements were made in triplicate on 3 samples per growing condition for each variety Spectra average were obtained from 64 scans recorded from 4,000 to 400 cm -1 at a resolution of 4 cm -1 , and corrected for background absorbance by subtraction of the spectrum of the empty ATR crystal. FTIR spectra were baseline corrected and vector-normalized in the range of 1850-675 cm -1 .

Chemical analysis

Fiber bundles were ground with a ball crusher (Retsch, MM 2000) and extractives were removed using 80% ethanol at 40˚C (20 ml/200mg, 4 times), then dried. Complete acid hydrolysis of cell wall polysaccharides into sugar monomers was performed using 1 M H 2 SO 4 for 2 h at 100˚C. The released monosaccharides were separated by high performance anion-exchange chromatography (PA1 column, Dionex) and analyses performed in triplicate as previously described [START_REF] Beaugrand | Impact and efficiencey of GH10 and GH11 thermostable endoxylanases on wheat bran and alkali-extractable arabinoxylans[END_REF]. Lignin content was estimated by a spectrophotometric procedure using acetyl bromide and 15 mg of cell wall residue (CWR) as previously described [START_REF] Iiyama | Determination of lignin in herbaceous plants by an improved acetyl bromide procedure[END_REF].

Data analysis

Statistical analysis was performed with Statgraphics software. Equality of group variances was checked with a Barlett Test. The normality assumption was assessed by Kolmogrov Smirnov Test. Differences between treatments were evaluated using Analysis of Variance (ANOVA) and the HSD Tukey test used to separate averages at p < 0.005. When the normal curve of distribution was not verified or the number of samples was reduced, the non-parametric Kruskal Wallis test was applied. Simple regression analysis was used to show the correlation between certain variables. Multivariate statistical analysis of infrared spectra was performed by principal component analysis using the Unscrambler software (CAMO, Germany). ). For both varieties, the highest bark and wood thickness were obtained under irrigation conditions (p < 0.001) (Fig 1b). Overall the bark thickness shows a positive relation with the wood thickness and the different growing conditions did not significantly affect the overall bark/wood ratio for the two varieties (average ratio: F17 = 0.2, S27 = 0.3). Such an observation would suggest that changes in growing conditions equally affect bark and wood formation leading to significant changes in stem diameter, bark and wood thickness, but not the bark/wood ratio. One possible hypothesis to explain these results could be that the changes in growing conditions impact on the meristematic activity of the vascular cambium thereby modifying the number of derived cells differentiating into wood or secondary fibers in the bark. In agreement with this idea it is interesting to note that alterations in groundwater regime modify cambial activity and the number of cell derivatives in poplar, leading to changes in the rate of stem growth and circumference [START_REF] Arend | Seasonal change in the drought response of wood cell development in poplar[END_REF][START_REF] Berta | Type 3 metallothioneins respond to water deficit in leaf and in the cambial zone of white poplar (Populus alba)[END_REF]. Irrigation would favor individual plant growth giving rise to taller and thicker plants. In contrast, high seed sowing density has been reported to limit individual plant growth most likely as a result of inter-plant competition [START_REF] Westerhuis | Sowing density and harvest time affect fibre content in hemp (Cannabis sativa) through their effects on stem weight[END_REF][START_REF] Amaducci | Crop yield and quality parameters of four annual fibre crops (hemp, kenaf, maize and sorghum) in the North of Italy[END_REF]. Interestingly Schafer [START_REF] Scha ¨fer | Effect of sowing date and plant density on the cell morphology of hemp (Cannabis sativa L.)[END_REF] observed that this reduction was more important in a dry year when compared to a wet year.

Results and discussion

Growing conditions affect bast fiber production

Examination of the sections by UV fluorescence microscopy was then used to observe fiber distribution. UV autofluorescence which indicates the presence of lignin is essentially restricted to the middle lamella-primary cell wall region in both primary and secondary bast fibers, whereas autofluorescence is also detected in the secondary cell walls of xylem tissue [START_REF] Bonatti | Histochemical and supramolecular studies in determining quality of hemp fibres for textile applications[END_REF][START_REF] Cro ˆnier | Structure and chemical composition of bast fibers isolated from developing hemp stem[END_REF][START_REF] Liu | Characterization and biological depectinization of hemp fibers originating from different stem sections[END_REF] (Fig 2). Microscopic observations show that growing conditions impact mainly on ). In contrast to primary fibers, secondary fibers are present as several clearly defined bundle layers, the number of which increases with irrigation.

Outer limits of fiber bundles were identified and used for image analysis to determine the areas corresponding to primary and secondary fibers thereby providing an overall indication of fiber quantity (Fig 3a). Results show that despite variations in bark thickness, growing conditions had no significant effect on the total area (cell wall plus lumen, expressed as μm 2 ) of primary fibers (Fig 3b) and no significant differences were observed between the varieties. In contrast, the secondary fiber area significantly increased with additional irrigation (p<0.001) and was correlated with the bark thickness (Fig 3c) thereby supporting our hypothesis that cambial activity in hemp is dependent on growing conditions and variety. Such an observation is in agreement with previous studies suggesting that secondary fiber production decreases with increasing plant density [START_REF] Amaducci | Influence of agronomic factors on yield and quality of hemp (Cannabis sativa L.) fibre and implication for an innovative production system[END_REF][START_REF] Scha ¨fer | Effect of sowing date and plant density on the cell morphology of hemp (Cannabis sativa L.)[END_REF] and increases under rainy conditions [START_REF] Hoppner | Cultivation experiments with two fibre hemp varieties[END_REF]. The S27 variety contains more secondary fibers than F17 under all tested conditions (p < 0.001) indicating that plant genotype also affects secondary fiber development as previously reported [START_REF] Amaducci | Influence of agronomic factors on yield and quality of hemp (Cannabis sativa L.) fibre and implication for an innovative production system[END_REF][START_REF] Mediavilla | Agronomic characteristics of hemp varieties[END_REF].

ATR FTIR analysis of native fibers

Our results showed that the amount of secondary fibers produced depends upon both the genotype and the growing conditions thereby affecting the overall quality of the harvested material. To confirm that ATR-FTIR spectroscopy could be used as a potential tool for rapidly characterizing the proportions of primary and secondary fibers in hemp plants we firstly recorded spectra on manually-dissected primary and secondary fibers from both varieties grown under different conditions (S2 Table ). Since the material used for ATR does not need any prior preparation the quantity required for analyses is greatly reduced. Principal component (PC) analysis indicates partial separation of the two fiber types (Fig 4a). The PC1 axis explains 72% of the statistically significant variation between samples showing that these results can be used to provide information about fiber type. Secondary fibers are positively correlated with PC1 positive loadings (Fig 4b ) showing several peaks. The four peaks at 1161 (C-O-C asymmetric valence vibration), 1054 (C-O valence vibration), 993 (C-C, C-OH, C-H ring and side group vibrations) and 1112 cm -1 can be mainly attributed to cellulose and hemicellulose absorption [START_REF] Dai | Investigation of the dislocation of natural fibres by Fourier-transform infrared spectroscopy[END_REF][START_REF] Faix | Classification of lignins from different botanical origins by FT-IR spectroscopy[END_REF][START_REF] Schwanninger | Effects of short-time vibratory ball milling on the shape of FT-IR spectra of wood and cellulose[END_REF][START_REF] Akerholm | Interactions between wood polymers studied by dynamic FT-IR spectroscopy[END_REF]. Negative loading from PC1 showed two main bands at 1598 cm -1 and 1412 cm -1 that can be assigned to the asymmetric and symmetric vibrations, respectively, of the COO (carboxylate) structure of esterified pectin [START_REF] Santos | Aqueous extraction of pectin from sisal waste[END_REF]. Overall these results suggest that ATR FTIR is an effective tool for discriminating primary and secondary fibers based on the differential chemical composition of their cell walls.

When analyzed separately, no clear clustering of the spectra recorded on primary (Fig 5a) and secondary fibers (Fig 5b) could be observed according to the different growing conditions. At this stage of development (end of flowering stage) most of the structural components of the fiber cell wall have been deposited and our data therefore suggest that contrasted growing conditions do not strongly affect cell wall composition of both primary and secondary hemp fibers [START_REF] Liu | Effect of harvest time and field retting duration on the chemical composition, morphology and mechanical properties of hemp fibers[END_REF][START_REF] Cro ˆnier | Structure and chemical composition of bast fibers isolated from developing hemp stem[END_REF].

ATR FTIR analysis of technical fibers

Since ATR FTIR successfully discriminated native primary and secondary fibers manually-isolated from the hemp stem we decided to see whether this technique could also be used to characterize technical fibers obtained by mechanical defibering of plants grown under the same secondary fibers thereby suggesting that mechanical defibering can result in poor separation of these cell types.

Principal component analysis of just technical fibers (i.e. without native fibers) showed that spectra were separated along PC1 (Fig 8a). This component explained 85% variability and separates samples according to growing conditions, but not variety. No specific band was positively associated in the PC1 loading for samples obtained from plants grown under additional irrigation conditions whereas bands at 1161, 1112 and 1054 cm -1 were positively associated to technical fibers from high seed density samples (Fig 8b). These bands were also observed when native secondary fibers were compared to primary fibers (Fig 4b).

Altogether, these results would suggest that differences in secondary fiber content might partly explain the observed variability of the technical fibers obtained from hemp grown under different conditions. The absence of bands at 1161, 1112 and 1054 cm -1 in the PC loading of ATR FTIR spectra of technical fibers obtained from plants grown under additional irrigation conditions would suggest that few/no secondary fibers are present. These results are in clear contradiction with our microscopic observations showing that greater amounts of secondary fibers are present in plants grown under additional irrigation conditions as compared to plants grown under high density conditions (Fig 3b). One possible explanation for this unexpected result is that growing conditions may affect mechanical behavior during defibering. Secondary fibers might remain attached to primary fibers in technical fiber samples obtained from plants grown under high seed density conditions. In contrast, they become separated from primary fibers in samples obtained from plants grown with additional irrigation. This hypothesis is supported by measurements (Table 1) showing that the proportion of technical fibers as a percentage of straw yields is slightly lowest for plants grown under additional irrigation (both varieties). In this scenario, although additional irrigation results in the production of higher amounts of secondary fibers, these fibers are more easily separated from primary fibers leading to i) reduced total fiber yield and ii) absence of secondary fiber ATR FTIR bands in technical fiber samples produced by mechanical defibering. Further analyses of technical fibers should allow clarification of this point and, if confirmed, would suggest that the presence/absence of secondary fibers represent an important source of the observed variability between technical fibers samples. Cell-cell cohesion in plants depends upon the structure of the middle lamella and it would obviously be interesting to see whether the observed differences in primarysecondary fiber separation are related to differences in outer stem morphology (e.g. quantity of parenchyma cells between the two fiber types), or to modifications in the middle lamella structure.

Overall our results indicate that the main difference in technical fiber quality appears to be based upon the relative proportions of primary and secondary fibers. In order to complete this survey of hemp bast fiber content we also determined the structure of manually dissected primary and secondary fibers by wet chemistry. For this we used the S27 variety that contained significantly higher (and sufficient for analysis) amounts of secondary fibers (Fig 3b). The results showed that lignin content in secondary fibers was significantly higher (5.5% on average) than in primary fibers (3.6% on average), (p < 0.001) (Table 2). Sulfuric acid hydrolysis also revealed that secondary fibers polysaccharides contained significantly more xylose (p < 0.001). Albeit weak these variations in chemical composition between primary and secondary fibers suggest that cell wall architecture differs slightly between the two types of fibers. Except for the proportion of mannose in secondary fibers (significant decrease under high seed density conditions), wet chemistry did not reveal any large variations in cell wall composition between plants grown under different growth conditions in agreement with infrared analysis.

Conclusion

Variability in technical fiber quality is generally considered as one of the major disadvantages of natural (plant) fibers and represents an important obstacle to their widespread use in composite materials. A rapid and reliable method of evaluating technical fiber quality would therefore represent a valuable help to industry. In this paper we have shown that ATR FTIR spectroscopy represents an efficient and powerful technique for evaluating variability in technical fibers obtained from different sources. Both growing conditions and genetic differences (variety effect) impact on hemp stem tissue architecture with a major effect on i) secondary (but not primary) fiber production and ii) primary and secondary fiber separation during mechanical defibering. Our results also indicate that a major source of variability in technical fibers is related to the presence of secondary fibers. This information represents an important contribution to our understanding of the biological reasons underlining the variability of natural fibers and should lead to the improvement of production and defibering protocols necessary for a better control of fiber quality. The use of ATR infrared analysis also provides chemical information about the surfaces of technical fibers. Such information is of great Table 2. Chemical composition of primary and secondary fibers isolated from hemp harvested at the end of flowering stage (variety S27). Lignin and total sugar contents as percent of cell wall residue; monosaccharides as percent of total sugars. Means not sharing a common letter are significantly different (p<0.05).

Fibers

Growing importance for a better understanding of bast fiber cohesion in agrocomposites and is complementary to information about the mechanical properties of technical fibers.

  Stem morphological features including bark and wood thickness and bast fiber distribution were determined by microscopy observations of hemp stem transverse sections. Examination of basal stem regions indicates significant variation in the thicknesses of stem outer tissues (designated as bark) and stem inner tissues (wood layer, but not pith) in transverse cross-sections (Fig 1a, S1 Table

Fig 1 .

 1 Fig 1. Stereomicroscopy imaging of transverse cross section of the hemp stem. (a) showing wood and bark including primary fibers (f1), secondary fibers (f2); arrows indicate the dimensions measured to determine wood-/barkthickness (b) in the basal stem region of F17 and S27 grown under different conditions (standard, irrigation, high sowing density: 100 kg.ha -1 ). https://doi.org/10.1371/journal.pone.0179794.g001

Fig 2 .

 2 Fig 2. UV autofluorescence of cross-sections from the basal stem region of hemp. S27 (a,c,e) and F17 (b,d,f) grown under standard (a,b), high sowing density: 100 kg.ha -1 (c,d) irrigation (e,f). Primary fibers (f1), secondary fibers (f2). https://doi.org/10.1371/journal.pone.0179794.g002

Fig 3 .

 3 Fig 3. Image treatment of UV auto fluorescence of hemp stem cross-sections (a) used to determine areas of primary fibers (b) and secondary fibers (c) in F17 and S27 varieties grown under different conditions (standard, irrigation, high sowing density: 100 kg.ha -1 ) (expressed as μm 2 ) with regards to bark thickness (expressed as μm). https://doi.org/10.1371/journal.pone.0179794.g003

Fig 4 .Fig 5 .

 45 Fig 4. Principal component analysis of ATR spectra of hemp native primary fibers (f1) and secondary native fibers (f2) (a) obtained from F17 and S27 with all growing conditions. Score loading corresponding to PC1 (b). https://doi.org/10.1371/journal.pone.0179794.g004

Fig 6 .Fig 7 .

 67 Fig 6. Principal component analysis of ATR spectra of hemp fibers showing dispersion with regards to technical fibers (tech) and native primary fibers (f1) (a). Score loading corresponding to PC1 (b). https://doi.org/10.1371/journal.pone.0179794.g006

Fig 8 .

 8 Fig 8. Principal component analysis of ATR spectra of technical fibers showing dispersion with regards to growing conditions (a). Score loading corresponding to PC1 (b). https://doi.org/10.1371/journal.pone.0179794.g008

Table 1 . Fiber yields of F17 and S27 hemp varieties harvested at the end of flowering stage (as percent of straw dry matter).

 1

	Variety	Standard	Irrigation	100 kg.ha -1
	F17	27.7 ± 0.6	26.5 ± 0.8	28.1 ± 0.2
	S27	32.4 ± 0.7	31.3 ± 0.1	32.9 ± 0.1

https://doi.org/10.1371/journal.pone.0179794.t001

  ://doi.org/10.1371/journal.pone.0179794.t002

		conditions	Lignin	Total Sugars	Glc	Gal	Xyl	Man	GalU
	Primary	Standard	3.0 a	79.9 a	90.1 a	2.5 a	0.2 a	5.5 a	4.1 a
		Irrigation	3.0 a	74.7 a	90.0 a	2.7 a	0.3 a	5.4 a	3.9 a
		100 kg.ha -1	3.4 a	75.0 a	91.1 a	2.2 a	0.3 a	4.7 a	3.5 a
	Secondary	Standard	5.5 a	72.3 a	87.1 a	2.9 a	1.8 b	6.2 a	4.4 a
		Irrigation	5.5 a	70.1 a	86.0 a	3.1 a	2.8 b	5.7 a	4.4 a
		100 kg.ha -1	5.6 a	67.7 a	90.4 a	2.6 a	1.9 b	2.9 b	4.8 a

https

PLOS ONE | https://doi.org/10.1371/journal.pone.0179794June 22, 2017

Acknowledgments

The authors thank Nicolas Dasilva for mechanical defibering of hemp, and Antoine Portelette for technical assistance in chemical analysis.

This work was performed in the frame of the project MAPROFI (MAı ˆtrise des PROprie ´te ´s des fibres de chanvre (fibre / chènevotte), dans le cadre de valorisations en agro-composites base polymère ») funded by the Champagne-Ardenne and Nord -Pas de Calais Regions, the European Union (FEDER funds), the French Association for Research and Technology, and Troyes Champagne Metropole.

All relevant data are within the paper and its Supporting Information files.


Y Masuda

I Misztal

Andres Legarra

S Tsuruta

D A L Lourenco

B O Fragomeni

I Aguilar

Avoiding the direct inversion of the numerator relationship matrix for genotyped animals in single-step genomic best linear unbiased prediction solved with the preconditioned conjugate gradient

Keywords: computation, genomic selection, inversion, numerator relationship matrix, preconditioned conjugate gradient, sparse matrix

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

INTRODUCTION

Single-step genomic BLUP (ssGBLUP) is a unified approach for genomic evaluation to combine phenotypes, pedigree, and genotypes [START_REF] Aguilar | A unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[END_REF]. Mixed model equations (MME) in ssGBLUP require the inverse of a genomic relationship matrix (G -1 ; VanRaden, 2008) and the inverse of a numerator relationship matrix ( 1 22  A ) for genotyped animals. When the number of genotyped animals is limited, possibly less than 100,000, both inverses can be efficiently calculated [START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF][START_REF] Aguilar | Efficient computation of the genomic relationship matrix and other matrices used in single-step evaluation[END_REF][START_REF] Fragomeni | Use of genomic recursions in single-step genomic BLUP with a large number of genotypes[END_REF]. When more animals are genotyped, G-1 can be obtained efficiently using the "algorithm for proven and young" (APY; [START_REF] Misztal | Using recursion to compute the inverse of the genomic relationship matrix[END_REF][START_REF] Misztal | Inexpensive computation of the inverse of the genomic relationship matrix in populations with small effective population size[END_REF], which exploits a limited rank of G due to a small effective population size. [START_REF] Faux | Inversion of a part of the numerator relationship matrix using pedigree information[END_REF] developed an algorithm to create 1 22  A directly from a pedigree. However, when the number of generations or the number of ancestors is large, the matrix becomes dense, resulting in expensive computations and more memory requirements. [START_REF] Strandén | Comparison of some equivalent equations to solve single-step GBLUP[END_REF] showed that 1 22  A could be decomposed into a product of several sparse matrices. When MME are solved with the preconditioned conjugate gradient (PCG), the explicit 1 22



A is not needed because the computations require only a product of ) by a vector (q). The computation is required for solving mixed model equations in single-step genomic BLUP (ssGBLUP) with the preconditioned conjugate gradient (PCG). The inverse can be decomposed into sparse matrices that are blocks of the sparse inverse of a numerator relationship matrix (A -1 ) including genotyped animals and their ancestors. The elements of A -1 were rapidly calculated with the Henderson's rule and stored as sparse matrices in memory. Implementation of 1 22  A q was by a series of sparse matrix-vector multiplications. Diagonal elements of 1 22  A , which were required as preconditioners in PCG, were approximated with a Monte Carlo method using 1,000 samples. The efficient implementation of 1 22  A q was compared with explicit inversion of A 22 with 3 data sets including about 15,000, 81,000, and 570,000 genotyped animals selected from populations with 213,000, 8.2 million, and 10.7 million pedigree animals, respectively. The explicit inversion required 1.8 GB, 49 GB, and 2,415 GB (estimated) of memory, respectively, and 42 s, 56 min, and 13.5 d (estimated), respectively, for the computations. The efficient implementation required <1 MB, 2.9 GB, and 2.3 GB of memory, respectively, and <1 sec, 3 min, and 5 min, respectively, for setting up. Only <1 sec was required for the multiplication in each PCG iteration for any data sets. When the equations in ssG-BLUP are solved with the PCG algorithm, 1 22  A is no longer a limiting factor in the computations.

that is, only

1 22  A q is needed. With PCG, diagonals of 1 22
 A are needed as preconditioners [START_REF] Tsuruta | Use of the preconditioned conjugate gradient algorithm as a generic solver for mixed-model equations in animal breeding applications[END_REF]. The purposes of this study were to implement efficient computations of the product 1 22  A q and of the diagonal elements of 1 22  A and compare their computing costs with that of an explicit inversion.

MATERIALS AND METHODS

Indirect Multiplication of 1 22



A q

A numerator relationship matrix involving only genotyped animals and their ancestors is defined as

00 02 * 20 22 = é ù ê ú ê ú ë û A A A A A
and its inverse is defined as

00 02 * 1 20 22  = é ù ê ú ê ú ë û A A A A
A , in which the subscript "2" represents a group of genotyped animals and "0" represents a group of their ancestors. The inverse considering inbreeding coefficients can be easily calculated using the rapid rule from a pedigree [START_REF] Henderson | A simple method for computing the inverse of a numerator relationship matrix used in prediction of breeding values[END_REF][START_REF] Quaas | Computing the diagonal elements and inverse of a large numerator relationship matrix[END_REF]. The matrix A *-1 does not need either descendants of or the animals unrelated to the genotyped animals.

Using a property of the inverse (Searle, 1982, p. 260), the inverse of A 22 is a function of submatrices of A *-1 (Strandén and Mäntysaari, 2014):

( )

1 1 22 20 00 02 22   =  A A A A A .
A product of the inverse by a vector, q, is expressed as

( ) ( ) { } 1 1 22 20 00 02 22   é ù =  ê ú ë û A q A q A A A q
.

[1] Let 4 temporary vectors be v, w, x, and y; the above product can be computed by a sequence of matrix-vector multiplications: v = A 02 q, w = (A 00 ) -1 v, x = A 20 w, y = A 22 q, and A q 22 1 -= y -x. The matrices A 22 , A 20 , and A 02 are sparse and hence the matrixvector multiplications can be efficiently computed. We do not need explicitly to compute (A 00 ) -1 because the product w = (A 00 ) -1 v can be computed by solving the sparse equation A 00 w = v. The Cholesky factor of A 00 is typically required to solve the equation. The sparse matrices A 22 , A 20 , and A 02 and the Cholesky factor of A00 are calculated and stored in memory before the PCG iterations. The product 1 22  A q will be calculated in each round because q is a vector of the current solution or search direction, and it changes every round. Although the explicit computation of the diagonal elements is expensive, the diagonals can be approximated by a Monte Carlo method [START_REF] García-Cortés | Multiple trait estimation of variance components in animal models with different design matrices[END_REF][START_REF] Dong | Stochastic estimation with Z 2 noise[END_REF][START_REF] García-Cortés | A Monte Carlo algorithm for efficient large matrix inversion[END_REF] as ( ) ( ) ( )

Diagonals of

1 1 22 22 1 diag 1/ n i i i n   = » å A s A s 
, in which n is the number of samples,  is the direct product operator, and s i is a vector of random numbers containing either +1 or -1 with equal probability and the product ( 1

22 i  A s ) is computed with Eq. [1]. The computa- tion of diag( 1 22  A
) has to be performed only once before the PCG iterations. Small errors in the preconditioner affect only the convergence rate in PCG, not the solutions of MME. In our tests, 1,000 samples were always sufficient to provide the same convergence rate and the solutions compared with the exact values of diag( 1 22  A ).

Data

Three data sets were used for testing (Table 1). Data 1 consisted of 15,723 genotyped animals and 16,694 pedigree animals in a commercial broiler population (Lourenco et al., 2015a). Data 2 consisted of 80,993 genotyped animals and 375,946 pedigree animals in the U.S. Angus population (Lourenco et al., 2015b). Data 3 consisted of 569,404 genotyped animals and 1,436,112 pedigree animals in the U.S. Holstein population [START_REF] Masuda | Implementation of genomic recursions in single-step genomic BLUP for US Holsteins with a large number of genotyped animals[END_REF].

Implementation and Computations

The inverse of A was calculated using the rapid method by [START_REF] Quaas | Computing the diagonal elements and inverse of a large numerator relationship matrix[END_REF] and stored as a sparse matrix using the SPARSEM module [START_REF] Misztal | Complex models, larger data, simpler computing?[END_REF], as present in the blupf90 package [START_REF] Misztal | Manual for BLUPF90 family of programs[END_REF]; http:// nce.ads.uga.edu/wiki/BLUPmanual). The system of equations involving A 00 were solved with a sparse matrix package, YAMS [START_REF] Masuda | Application of supernodal sparse factorization and inversion to the estimation of (co)variance components by residual maximum likelihood[END_REF], which used optimized dense-matrix subroutines from the Intel Math Kernel Library (Intel Corporation, Santa Clara, CA). For comparison, A 22 was explicitly created and inverted as in [START_REF] Aguilar | Efficient computation of the genomic relationship matrix and other matrices used in single-step evaluation[END_REF]. The diagonals of 1 Genotyped animals and their ancestors.

1 22  A were calculated with the exact inverse and with the Monte Carlo method. All programs were written in Fortran 95 and compiled with the Intel Fortran Compiler 13.1 (Intel Corporation). Analyses were performed on a computer running Linux (x86_64) with Intel Xeon E5-2689 central processing unit (2.9 GHz) processor with 8 cores.

RESULTS AND DISCUSSION

Computing time and memory requirement for the explicit calculation of A 22 and its inverse are shown in Table 2. Memory requirements were 1.8 GB for Data 1 and 49 GB for Data 2. Computing time was 42 s for Data 1 and 56 min for Data 2. For Data 3, the memory requirement was too much to allocate for the available computer. Extrapolating from Data 2 assuming a quadratic cost for memory and cubic costs for computing time, the computing for Data 3 would require 2,415 GB and more than 13 d to finish the inversion.

Table 2 also presents the computing time and memory requirement to prepare the sparse components of 1 22  A . Total wall-clock time for the preparation of sparse components was short for all the data sets (less than 1 s in Data 1, 5.3 min in Data 2, and 4.4 min in Data 3). The required memory for the sparse components was much less than for the full inverse. Although Data 2 had fewer genotyped animals than Data 3, longer running time and more memory were required for the preparation. This was caused by more nonzero elements in the factor of A 00 , as almost all of them were newly created as "fill-in" during the factorization [START_REF] Masuda | Application of supernodal sparse factorization and inversion to the estimation of (co)variance components by residual maximum likelihood[END_REF]. This illustrates the fact that1 22  A can be relatively dense for the limited number of animals. Therefore, as shown in this study, the indirect computation of 1 22  A q takes advantage of sparsity in storage and computing time compared with the direct inversion of A 22 . The factor of A 00 as used with YAMS actually occupied more than 99% of the required memory, but this amount of memory is still negligible with current computers.

Wall-clock time for 1 22  A q with the sparse components for one round of the PCG iteration was also negligible (<1 s) as shown in Table 3. Solving the sparse equation A 00 w = v was inconsequential in the computation. One iteration of PCG in ssGBLUP with Data 3 took about 12 s [START_REF] Masuda | Implementation of genomic recursions in single-step genomic BLUP for US Holsteins with a large number of genotyped animals[END_REF].

In ssGBLUP runs, the use of the Monte Carlo approximation of diag( 1 22  A ) with 1,000 samples resulted in the same convergence rate and solutions as with the exact diagonals (results not provided). The approximation was less expensive to compute, as it required only 2.5 min for Data 2 and 2.2 min for Data 3 compared with 56 min and 13.4 d (estimated), respectively, to compute the exact diagonals. In summary, the computation of 1 22  A q removes computing limits from ssGBLUP as these computations take a small fraction of memory and computing time even for the largest pedigrees.

 = é ù ê ú ê ú ë û A A A A
A , in which the subscript "2" and "0" represent a group of genotyped animals and their ancestors, respectively.

2 Parallel processed with 8 computing cores.

3 Using a Monte Carlo method with 1,000 samples.

4 Extrapolated.

A

  by an arbitrary vector, say q; Technical note: Avoiding the direct inversion of the numerator relationship matrix for genotyped animals in single-step genomic best linear unbiased prediction solved with the preconditioned conjugate gradient 1 Y. Masuda,* 2 I. Misztal,* A. Legarra, † S. Tsuruta,* D. A. L. Lourenco,* B. O. Fragomeni,* and I. Aguilar ‡ *Department of Animal and Dairy Science, University of Georgia, Athens 30602; †INRA, UR631 SAGA, BP 52627, 31326 Castanet-Tolosan Cedex, France; and ‡Instituto Nacional de Investigación Agropecuaria, Canelones, Uruguay 90200ABSTRACT: This paper evaluates an efficient implementation to multiply the inverse of a numerator relationship matrix for genotyped animals (

A

  as a part of a preconditioner.

Table 1 .

 1 Description of data used in this study

	Item	Data 1	Data 2	Data 3
	Number of animals
	Genotyped	15,723	80,993	569,404
	In selected pedigree 1	16,694	375,946	1,436,112
	In whole pedigree	213,297	8,234,208	10,710,380
	Equations
	Order	219,226	8,526,614	21,985,710

Table 2 .

 2 Wall-clock time and storage memory required for the preparation of sparse components 1 of the inverse of a subset of numerator relationship matrix (A 22 ) for genotyped animals and the results from the computations with the direct inverse of A 22 for comparisons

	Item				Data 1	Data 2	Data 3
	Wall-clock time for preparation 2
	Indirect approach
	Setting up A 00 , A 20 , A 02 , and A 22 <1 s	<1 s	2 s
	Factorization of A 00			<1 s	167 s	125 s
	Computing diagonals 3 of	A	1 22 	<1 s	148 s	134 s
	Direct approach
	Setting up A 22			<1 s	199 s	2.7 h 4
	Inversion of A 22			42 s	56 min	13.4 d 4
	Required memory for storage
	Indirect approach			<1 MB	2.9 GB	2.3 GB
	Direct approach			1.8 GB 49.0 GB 2,415 GB 4
	1 Submatrices of the inverse of a numerator relationship matrix
	00	02
	* 1
	20	22

Table 3 .

 3 Wall-clock time for the indirect multiplication of the inverse of a subset of numerator relationship matrix (A 22 ) with a vector (q) in an iteration in the preconditioned conjugate gradient (PCG) and the results from the computations with the full inverse for comparisons

	Wall-clock time for	1 22  A q per PCG iteration 1 Data 1 Data 2 Data 3
	Indirect approach	<1 s	<1 s	<1 s
	Direct approach	<1 s	<1 s	59 s 2

Parallel processed with 8 computing cores.

Extrapolated.

1 This research was partially funded by the United States Department of Agriculture's National Institute of Food and Agriculture (Agriculture and Food Research Initiative competitive grant 2015-67015-22936).