Vegan diet in young children remodels metabolism and challenges the statuses of essential nutrients

Abstract Vegan diets are gaining popularity, also in families with young children. However, the effects of strict plant‐based diets on metabolism and micronutrient status of children are unknown. We recruited 40 Finnish children with a median age 3.5 years—vegans, vegetarians, or omnivores from same daycare centers—for a cross‐sectional study. They enjoyed nutritionist‐planned vegan or omnivore meals in daycare, and the full diets were analyzed with questionnaires and food records. Detailed analysis of serum metabolomics and biomarkers indicated vitamin A insufficiency and border‐line sufficient vitamin D in all vegan participants. Their serum total, HDL and LDL cholesterol, essential amino acid, and docosahexaenoic n‐3 fatty acid (DHA) levels were markedly low and primary bile acid biosynthesis, and phospholipid balance was distinct from omnivores. Possible combination of low vitamin A and DHA status raise concern for their visual health. Our evidence indicates that (i) vitamin A and D status of vegan children requires special attention; (ii) dietary recommendations for children cannot be extrapolated from adult vegan studies; and (iii) longitudinal studies on infant‐onset vegan diets are warranted.


Online questionnaires
We emailed the parents links to four short online questionnaires inquiring background information of the parents and the participating child as well as the food consumption of the responding parent(s) and the participating child. These included questions on the type of fat and salt usually used in cooking and whether salt is usually added when boiling porridge, rice, pasta and potatoes. A question whether the child has been given dietary supplements during the past month, and if yes, what is the product name, dosage and frequency of use, was also included. Information on the nutrient values of dietary supplements was gathered into an Excel table from producers' internet pages and package labels.

Food record
We sent a four-day food record to parents of each participant and assigned the dates (three weekdays and one day on the weekend) for filling in the record. The instruction page of the food records advised parents to record all foods and beverages that their child consumed during the recording days, except for what they consumed at preschool. An example page was also included. We provided the families with a validated [S1] Children's Food Picture Book[S2], designed to assist in portion size estimation. The parents were instructed to list all the ingredients of composite dishes and to estimate portion sizes using the picture book, weighing, household measures such as teaspoons or tablespoons, or package labels. For packed food products, the exact brand and product name was required.
We gave the day care personnel a separate pre-coded food record for recording food consumption at the day care centre on the dates matching the home food record. We instructed the personnel orally, and the food record included written instructions. Breakfast, lunch, afternoon snack and possible additional snacks each had predefined sections. Different food groups, such as main courses, side dishes (potatoes, pasta, rice), and salad at lunch each had predetermined rows. We also gave the personnel the Children's Food Picture Book and they could also estimate the amounts in household measures.
The food data were recorded using AivoDiet dietary software, which included the Fineli Food Composition Database Release 16 (2013) of the National Institute for Health and Welfare. We had previously updated this as described elsewhere [S3]. We also added new food items, especially vegan products, to the database when necessary. We entered the nutrient values of fortified food products as per information on package labels. The database includes recipes for typical Finnish mixed dishes. For each individual meal, the research assistant used a suitable recipe from the database, modified an existing recipe, or created a new recipe according to the parents' reports. The salt content of home dishes was based on the recipes in the database and we also made use of the information on the type of fat and the type and use of salt in the household gained form the online questionnaires. For food eaten at day care, we used recipes provided by the city of Helsinki's early childhood education and care food services. The food services also provided information on the products (e.g. milk, non-dairy milk substitutes, bread, fat spread, yoghurts) used at each day care centre.
After data entry, we checked for outlying values of food consumption in grams and outlying energy and nutrient intakes. After extracting the data from the software, each food code (food item or mixed dish) appearing in the data set was assigned to a food group and nutrient retention factors [S4] were applied using a single factor per nutrient per food group. The food records were kept between May and October 2017. Two omnivores and three vegetarians had missing or incomplete recordings for one to three days and we omitted these days from the analysis. We included all four days of data in the analysis for the rest of the participants. To calculate total intake from food and supplements, we combined the intake data from dietary supplements with the food record data in R.

Proportion of retinoids in vitamin A intake
Values for retinol and carotenoids were not available in our food composition database and thus the intake of different forms of vitamin A could not be directly calculated. Therefore, we calculated the proportion of retinoids in vitamin A manually by disaggregating the food consumption data into ingredients, multiplying the weight of each food containing retinoid naturally or by fortification with its retinoid content per gram. To find the retinoid values, we used product details available online and food composition databases. The proportion of retinoids in the total RAE intake was calculated as follows: retinoids in uncooked ingredients (µg) for the diet group ∑ RAE in uncooked ingredients (µg) for the diet group * 100%

Supplement 2: Extraction protocols for untargeted metabolomics
The following protocol was used for extraction of metabolome and lipidome for untargeted metabolomics: Frozen serum samples were thawed on ice, 20ml of serum was mixed to a microtube with 180ml 80% methanol for metabolomics or 10ml of serum with 190ml 100% isopropanol for lipidomics. Solutions were vortexed for 15 seconds, incubated at 4 • C for 1 hour and centrifuged at 14000g for 15 minutes to achieve supernatant that was transferred to another microtube and frozen to -20 • C until analysis.

Biomarkers for vitamin A status
We used two approaches for measuring vitamin A status in the participants.
Retinol-binding protein (RBP) has been used as a biomarker for vitamin A deficiency in large epidemiological studies. We used similar equipment and limits for insufficiency (<1·17mmol/l) and deficiency (<0·83mmol/l) as Engle-Stone et al. [S5] who reported sensitivity and specificity of 94·7% and 88·9% for RBP in predicting vitamin A deficiency in children.
In addition, we used a regression model of RBP, transthyretin and CRP described by Talsma et al. [S6] for vitamin A status assessment. For convenience, we made slight adjustments to the model that do not affect its specificity and sensitivity (Talsma et al. reported AUC of 0·98): We adjusted the constant of the model by -0·496 so that the decision limit of the model was 0 instead of the original 0·496. We also multiplied the estimated regression coefficient by -1 so that the deficiency was denoted by a negative value of the model instead of a positive value.

Permutation tests for assessing statistical significance
In our study, several aspects contradict the use of Student's t-test. We had a relatively small sample size, making it impossible to evaluate predict if the concentration of each measured metabolite is normally distributed in reference population. Furthermore, some measurements include multiple outlier suspects. With small sample size it is problematic to exclude these values from analyses as there is no credible reason to believe the outliers would result from measurement-related rather than physiological reasons. However, these outliers have a strong effect on estimates of variance with small sample size. This decreases the power of e.g. Student's t-test. In addition, there may be age-and sex-dependent variance in physiological concentrations of different metabolites, particularly in children. This effect should be considered while conducting statistical tests.
Therefore, we tailored a following age-and sex-adjusted non-parametric permutation test for our study, which is robust against outliers and does not rely on the unverifiable assumption of normality or equivariance. Recently, the use of such tests and probability index models have become more popular [S7,S8] Stratification and the measure of difference Let us assume two datasets, X and Y , representing observed concentrations of a metabolite in two groups, say omnivore and vegan groups. Sample sizes may be unequal. For each data point, we have an ID label, age and sex of the corresponding participant and measured metabolite concentration.
Age and sex information is used in stratification (dividing the participants to age-and sex-standardized groups). With small sample size we encounter a risk of empty age-sex strata. Thus, it is necessary to divide strata by age somewhat sparsely. In our study, as there was only 6 vegans, we were forced to use only two age classes: (1) under 4-year-olds and (2) at least 4-year-olds. This yields four strata.
To avoid unnecessary outlier removal, we will use a non-parametric probability index as a measure of difference between these datasets. Probability index P (X i > Y j ) is defined as the probability that randomly selected individual values from each dataset, X i and Y j , satisfy the inequality X i > Y j . With a small sample size and decent computing power it is possible to compute these inequalities, a brute force method of checking this inequality over all possible datapairs of observations (n X · n X operations) in the sample is eligible.
Forming null distribution by permutation First, we estimate a statistical null distribution of the probability index based on our data. [S9] For k rounds: 1.
Step 1: Randomly permute the group ID labels inside each stratum. This basically divides the participants randomly to two groups in an age-and sex-standardized manner while also maintaining the group sizes. 2.
Step 2: Calculate P (X i > Y j ) with the permuted group labels and calculate quadratic deviation from the null P (X i > Y j ) = 0.5 to be the test statistic. Save test statistics to a vector T = (T 1 , T 2 , . . . , T k ) that will eventually form the null distribution.
In our simulations, we used k = 47500 1 that took approximately 330 seconds to run with the datasets of 24 controls and 6 vegans and 6-year-old 1·90GHz processor.
Statistical test To use proper two-sided statistical test, we will use a squared test statistic. As probability index P (X i > Y j ) ∈ [0, 1] and comparing two identical datasets yields P (X i > Y j ) = 0.5, we used a test statistic After forming the null distribution or probability indexes, we transform it to test statistics and compare the test statistic from the original data, t real , to the null distribution to achieve approximated p-value: p ≈ P (t real > T ).
Finally, the value of the test statistic from the original data (with the real group labels), t real is used with its the null distribution. After forming the null distribution or probability indexes, we transform it to test statistics and compare the test statistic from the original data, t real , to the null distribution to estimate the approximated p-value:p

Supplement 3c: GSEA-based pathway analysis
Gene-Set Enrichment Analysis (GSEA) was originally introduced by Subramanian et al. [S10] for obtaining pathway-level knowledge from complex gene expression datasets. We used a GSEA-based approach with HMDB 3.0 pathway database as a framework for metabolome-level pathway analysis. The algorithm was carried out in R and the principles are discussed in the following.
Data preprocessing and statistical significance of metabolite-level change Original data from untargeted metabolomics can be seen in Table S3. For each ion found, there are usually multiple alternative metabolites that correspond to same mass per charge -ratio, thus indistinguishable by mass spectrometer. In total, our data had 10155 alternative metabolites contributing to data signal. We matched each of these metabolite ID's to HMDB 3.0 pathway database to discover the amount of metabolites for each pathway that was recognized at all in our data. Metabolites in one pathway have very similar mass per charge -ratios as many modifications include deleting or adding single ions. This poses an obvious problem of metabolites of a certain pathway appearing under the same mass per charge -peak in mass spectrometry. To avoid overestimation of amount of "found metabolites", we allowed counting each ion for each pathway only once. For example, out of 49 metabolites in bile acid biosynthesis pathway, total of 23 metabolites were found in the list of alternative metabolites. 10 of these were not included in assessing statistical significance, as they appeared under the same mass per charge -peak as others, yielding 13 found metabolites. However, these 10 discarded metabolites may contribute to statistical significance of another pathway, assuming they were the only metabolite with same mass per charge in that particular pathway. List of pathways in the database and corresponding data can be seen in Table S2.
These found metabolites were further divided to significantly changed metabolites. Significance criteria in our pathway analysis was that fold change for the metabolite had to be > 0.25 and uncorrected p-value had to be < 0.05. As the goal of pathway analysis is to recognize changes that correlate in one pathway, it is not necessary to use as tight significance levels as when analyzing single metabolites. Currently, there is no golden standard for pathway analysis significance levels, partly because there is no golden standard for pathway analysis method altogether.
Statistical significance of pathway-level change with Fisher's exact test We analyzed each pathway separately for positive changes, negative changes and all changes regardless of direction. In contrast to GSEA, we did not use permutation tests based on Enrichment Score for assessing pathway-level statistical significance. To save computation time, we preferred Fisher's exact test with null hypothesis that the proportion of significant changes to certain direction in a set of found metabolites in pathway i, is equal to the proportion of significantly changed ions to that direction in the set of all found ions 2 : P-values from Fisher's exact test were corrected for multiple hypothesis testing with Benjamini-Hochberg method.
GSEA-based Enrichment Score as a measure of pathway impact Enrichment Score (ES) describes how well changes in metabolite levels of a certain pathway correlate with each other. If a pathway has multiple significant changes to same direction, this pathway will have high ES. If significant changes are randomly scattered in terms of direction and magnitude, the pathway will have low ES. Thus, ES is affected by (1) the magnitude and (2) amount of changes in metabolite levels of a pathway as well as (3) consistency of changes. We calculated ES for each pathway as described by Subramanian et al. [S10] with the exception that we used fold change instead of correlation to sort metabolites according to their magnitude of change.
In addition, only the significant changes, i.e. the changes satisfying the conditions defined above, were counted for the Enrichment Score cumulative function.
Disadvantages of using Enrichment Score as a measure of pathway impact include overestimating the impact on pathways with high amount of metabolites.
Appendix Figure S1: Participant flow chart     Figure S3: Phospholipids detected by untargeted metabolomics. Phospholipids include two fatty acid tails, and the figure is grouped by sum of fatty acid length/double bonds, and phospholipid type. Notation of fatty acid is "sum of carbon atoms in the two fatty acids : sum of double bonds in the two fatty acids". Each bar represents group median and error bar is interquartile range. Note that almost all detected phospholipids are found in smaller amounts in vegan serum than in omnivores and vegetarians. However, this is most probably related to most serum phospholipids being located in lipoprotein particles, and thus correlating well with low cholesterol in vegans (see figure S6). PA = phosphatidic acid, PC = phosphatidylcholine, PE = phosphatidylethanolamine, PG = phosphatidylglycerol, PI = phosphatidylinositol, PS = phosphatidylserine.  Figure S4: Triglycerides detected by untargeted metabolomics. The figure is grouped by fatty acid length and molecule type. Notation of fatty acid is "sum of carbon atoms in the two fatty acids : sum of double bonds in the two fatty acids". Each bar represents group median and error bar is interquartile range. MCFA = medium-chain fatty acid (the triglyceride includes most probably three fatty acids with carbon atoms between 8-13 in each), SCFA = short-chain fatty acid (the triglyceride includes most probably three fatty acids with carbon atoms between 3-7 in each), (V)LCFA = (very) long-chain fatty acid (the triglyceride includes most probably three fatty acids with at least 14 carbon atoms in each).    Figure S6: Correlation between cholesterol and detected molecules with one, two or three fatty acids. Nearly all phospholipids, sphingolipids and (V)LCFA triglycerides show great correlation with serum cholesterol, while MCFA triglycerides and most lysophospholipids and free fatty acids ("one-chained" fatty acid molecules) vary more independently. Notation of fatty acid is "sum of carbon atoms in the two fatty acids : sum of double bonds in the two fatty acids". Some lollipop bars show multiple candy due to e.g. one lysophospholipid and one free fatty acid were detected with identical fatty acid type.
Appendix  Table S2: Results of the pathway analysis showing a list of pathways included in the analysis, total amount of metabolites in each pathway, amount of found metabolites by the untargeted metabolomics method and significant changes to positive (concentration in vegans > concentration in omnivores) or negative direction with corresponding p-values and Enrichment Score (defined in Supplement 3c).