Inﬂuences of Rearing Environment on Behaviour and Welfare of Captive Chilean Flamingos: A Case Study on Foster-Reared and Parent-Reared Birds

: Behaviour signals the internal states that relate to an individual’s welfare and its development is inﬂuenced by the early social environment that an animal experiences. Husbandry practices can alter this early social environment, for example different rearing conditions (e.g., foster rearing by a surrogate parent of another species). Widespread implementation of altered rearing can lack empirical support and non-parent-reared animals may experience poorer welfare resulting from maternal deprivation. An opportunity presented itself to measure the effect of foster-rearing on Chilean ﬂamingo behaviour and social preferences at WWT Slimbridge Wetland Centre and compare ﬁndings to parent-reared conspeciﬁcs in the same time period. Data were collected from April to July 2019 at three timepoints during each observation day. Binomial generalized linear mixed models were used to assess the relationship between focal chicks’ rearing background with behaviour, zone usage, and ﬂock position whilst accounting for climatic factors and visitor numbers. The development of social preferences was assessed using social network analysis. Our results showed limited impacts on ﬂamingo behavioural development due to foster rearing. Foster-reared chicks spent less time feeding, were more likely to occupy the nesting area of the enclosure, and had fewer signiﬁcant preferred associations than parent-reared chicks, but preferred social bonds were as equally strong and durable for both foster-reared and parent-reared chicks. Our results have important welfare implications for the use of foster-rearing in captive environments; altered early social rearing environments through cross-fostering in Chilean ﬂamingos is associated with limited differences in behavioural and social development.


Introduction
Behaviour is considered a relevant and sensitive signal that informs us of individuals' specific preferences, requirements, and internal states that correspond to their welfare and subjective wellbeing [1][2][3][4]. It is well established that postnatal environmental factors during development play a significant role in determining behavioural responses [5]. Amongst these factors is the early social rearing environment, which is often manipulated within zoological institutions for conservation purposes [6]. As a result, the postnatal rearing conditions for captive young tend to vary. Parent-rearing (i.e., by the natural parents of the offspring), hand-rearing (i.e., by a human caregiver), peer-rearing (i.e., reared in the presence of same age peers), foster-rearing, and cross foster-rearing are common practice, and it is important to understand their unique impact on ontogenetic behaviour. Foster-rearing and cross foster-foster rearing refers to "rearing non-maternal young by intraspecific or interspecific surrogate parents", respectively [7].
The behavioural alterations incurred through hand-rearing and peer-rearing practices suggests that secure, maternal bonds promote the development of relevant and adaptive behaviours [21]. As such, foster-rearing and cross foster-rearing, which may provide a suitable alternative biologically relevant social environments, presents themselves as useful husbandry practices in zoological institutions [22]. Cross-fostering has been used extensively in captivity with many species, to let experienced parents care for abandoned young or to give inexperienced parents a chance to learn how to raise young [23]. This has benefited conservation efforts in rodents [24], marsupials [25], and birds [26] to increase the productivity of breeding individuals in managed programmes.
Although its use is extensive, few reports have investigated the developmental impacts of foster rearing to the same extent as the hand-rearing and peer-rearing literature. Studies that use foster-rearing to assist with conservation efforts typically utilise survival analyses as a measure of success (e.g., [7,[27][28][29]). Although survival prerequisites the appropriateness of husbandry manipulations, it does not necessarily equate to the development of adaptive behaviours or positive welfare. Altered patterns of behaviour and sociality induced by disruptions of the early rearing environment are likely to go unnoticed by survival analyses (Beck 2002 as cited in [12]). A research bias toward primates and megafauna within this literature also hinders our understanding of the welfare needs across taxa [3].
Where behavioural analyses have been conducted on birds, cross-fostering has shown to result in misprinting (i.e., chicks show similarities and preferences to their foster species). Misimprinting in blue tits (Cyanistes caeruleus) and great tits (Parus major) is well documented for several behaviours, including mate choice [30], aggression [31], social dominance [32], alarm calls [33], bird song [34], and paternal behaviour [35]. Similar findings for communication and mate preference are demonstrated research on male zebra finches cross-fostered to Bengalese finches (Lonchura striata var. domestica) [36][37][38]. Consequently, where cross fostering is used as a husbandry tool, monitoring of the behavioural development throughout the juvenile stages into adulthood should be performed to evaluate the efficacy of this procedure.
Indications that the early rearing environment influences the development of social networks is evident in sandhill cranes (Antigone canadensis) who show association preferences for their foster species rather than conspecifics [23]. Foster-reared Hawaiian goose (nene) goslings (Branta sandvicensis) were less able to integrate and adjust into new environments compared to parent-reared conspecifics [39]. Foster-reared goslings displayed longer latencies to associate with flock members and performed behaviours indicative of an inferior antipredator strategy (e.g., reduced vigilance and predator avoidance) [39]. However, these foster-reared goslings only had visual and auditory contact with their foster parents due to a mesh wire barrier. Parent-reared goslings had more social opportunities and were housed in larger enclosures. Hence, full maternal care was precluded and environmental conditions were not consistent between rearing practices.
When interactions with parents and exploration opportunities have been provided, no differences in the activity patterns and associations can be seen, as highlighted by research comparing the development of captive killdeer (Charadrius vociferous) cross-fostered to wild spotted sandpipers (Actitis macularia) and wild parent-reared killdeer [40] (Charadrius vociferous); Nonetheless, the use of fostering and cross-fostering alike has expected and unexpected fitness outcomes [41]. As a result, authors have advised against, and limited, the use of fostering as a husbandry technique for reintroduction programmes; due to the potential for unsuitable coping strategies and mate choice both in the wild and in captivity (e.g., [23,40,42]). However, the use of foster rearing techniques can be widespread for particular species across and within zoos. To fully understand the long term of effects of such strategies, consistent methods and repeatable data collection techniques need to be applied to the young animals involved and provide a contemporary investigation of behavioural development.
Captive flamingos offer excellent opportunities for strong empirical behavioural research due to their large sample sizes, diverse behavioural repertoire, ability to identify ringed individuals, and their global captive presence permitting research replication [43]. As such, this study investigates how disruption of the early rearing environment affects the behavioural development of Chilean flamingo (Phoenicopterus chilensis) chick's crossfostered to Andean flamingos (Phoenicoparrus andinus) at WWT Slimbridge Wetland Centre. This unique opportunity to collect data on an understudied area of zoo management became apparent when a sudden and unexpected breeding event in the group of Andean flamingos occurred. Behavioural differences between Andean and Chilean flamingos are few and the two species potentially have a close evolutionary relationship [44]. Food selectivity varies between the two, and Chilean flamingo chicks show looser social bonds out of the breeding season [45][46][47][48]. Differences between these two close taxonomic relatives are most prominent in their visual characteristics as adults ( Figure 1). wing-coverts are more heavily marked red. Flight feathers are black, more obvious than in other species, as not covered by back plumes. Legs are yellow; bill is black, pale yellow at base sometimes with pink flush on culmen; eyes are red-brown, with a small area of red skin on the lores. (b) Chilean flamingo: Plumage of head, neck and body is pale salmon pink, deeper on lower neck and upper breast. Scapulars, long back feathers, and upper and lower coverts are red. Legs are gray with contrasting red-pink feet and joints. Bills have an extensive black tip. The base of the bill is paler, almost white and the eyes are yellow. Morphological descriptions are taken from Brown and King (2005 [49]).
There is potential that these behavioural and physiological differences may manifest into altered behavioural repertoires and social preferences in cross-fostered Chilean flamingo chicks. We aimed to use well established welfare measures of captive flamingo activity and enclosure usage (e.g., [48]), as well as social network analyses (e.g., [50,51]) to compare temporal patterns of behaviour and sociality of cross-fostered Chilean flamingo chicks with that of parent-reared Chilean flamingo chicks housed in the same enclosure.
Chilean flamingos are a commonly housed zoo flamingo, whereas the Andean flamingo is not (being held, as of 2021 at two institutions). Although husbandry enclosure factors were similar for these two study flocks, limited information on captive Andean flamingo breeding behaviour means that no prior assumptions on the behavioural development of cross fostered chicks were made. Our findings have the potential to inform flamingo husbandry regimes, with an aim of helping produce self-sustaining populations across many zoological institutions. Our results also provide a baseline of flamingo behavioural development that are useful for others wishing to evaluate breeding behaviour in captive flocks.

Study Species
The Chilean flamingo and Andean flamingo are both highly gregarious, diurnal Phoenicopteriformes that feed, reproduce, and communicate in similar ways. They live across overlapping geographical ranges in South America where their population growth is limited by anthropogenic threats and environmental changes [52,53].
The behaviour of 13 Chilean flamingo chicks was observed at WWT Slimbridge Wetland Centre, Gloucestershire, UK, from March to July 2019. Prior to the chicks hatching, six eggs were translocated from their Chilean flock to the Andean flamingo flock at WWT Slimbridge Wetland Centre on the 29 June 2018. Upon hatching at the beginning of August 2018, the six chicks continued to be reared by Andean flamingo foster-surrogates. The Andean flamingos resided in the Andean Flamingo Pen (1093 m 2 ) that also housed several species of wildfowl (Anseriformes). The Andean Flamingo Pen contained a nesting island (21 m 2 ) and a pool that encompassed c30% of the enclosure. The foster-reared chicks were then moved back into their Chilean flamingo flock on the 27 February 2019. This unique event gave rise to a situation where (once the cross-fostered chicks had been moved from the Andean Flamingo Pen to the South American Pen) six foster-reared (F = 4, M = 2) and seven parent-reared Chilean flamingo chicks (F = 4, M = 3) could be easily observed within the same Chilean flamingo flock (N = 137).
The flock of Chilean flamingos resided within the South America Pen, a walk-through enclosure (4921 m 2 ) containing several other species of captive wildfowl. The enclosure consisted of distinct features including a large pond, several islands, grassland, bushes, steppingstones, and public footpaths ( Figure 2). The flock had constant access to a large shelter with sloped ground and shallow water, where they are fed during winter. In the summer, they are fed outdoors in a shallow feeding pool near the front of the enclosure. The birds were fed twice a day on a complete flamingo pellet; once in the morning (c.08:30) and once in the afternoon (c.15:00), with water available ad libitum. Visitors could enter the enclosure between c.09:30 and c.17:30 and observe the flamingos from several accessible viewpoints. For the purposes of brevity, cross-fostered Chilean flamingo juveniles are simply referred to as "foster-reared chicks".

Data Collection
On each day of data collection, the flock was sampled from 10:00-11:00, 12:00-13:00, and 15:00-16:00 by one researcher. Photographic records were used to collect flamingo behavioural data, employing an instantaneous scan sampling approach with a photo being taken at five-minute intervals to record the activity, enclosure usage and social associations of the focal Chilean flamingo chicks. Foster-reared chicks were individually identified using leg rings. Parent-reared chicks were not provided with leg rings and were instead identified through differences in plumage with the rest of the flock and sexed based on differences in physical size with other chicks (see [54] for a detailed report on flamingo plumage development). Male and female parent-reared chicks were grouped and given artificial codes for analysis (PAM and PAF, respectively). 'South America' walk through enclosure at WWT Slimbridge split into zones for analysis. 1/purple = shelter; 2/red = back grasslands; 3/green = nesting area; 4/orange = bushlands; 5/blue = mid-waters; 6/pink = island; 7/dark green = front waters; 8/white = public viewing areas.
Photographs were taken with a Canon EOS 750D Digital SLR camera and later a Casio Exilim EX-FH20 Digital camera due to technical issues. Flamingo states behaviour and enclosure usage were scored from the photos and social preferences were deduced from photos in the same as per Rose, Brereton, and Croft (2018, [48]) and Rose and Croft (2020 [51]). Observations were made from public viewing areas to ensure flamingos were habituated with the presence of an observer. The total number of visitors visible in the 'South America' enclosure were recorded every 15 min within each hour observation period. Total number of visitors entering WWT Slimbridge on each day of data collection were also provided by the wetland centre. Temperature, humidity, daily sunshine, and rainfall were recorded from worldweatheronline.com [55]. Supplemental data for social network analysis were collected by a second researcher, adopting the same procedure. A total of 75 h and 25 min of data were recorded.

Activity
Photographic data were assessed to record the behavioural states of individual flamingo chicks using a previously established ethogram adapted from Rose, Brereton, and Croft (2018 [48]; Table 1). Counts of behavioural states were used to produce time-activity budgets for individual foster-reared chicks. Activity budgets for the average female and male parent-reared chick were also produced. Chick raises and displays scapular feathers ("chrysthanthemuming"); or postures bill towards a nearby bird, neck swaying side-to-side with directed shouting vocalisations (hooking); or jabs another bird with its bill (jousting)

Courtship
Chick stretches its wings out from its body then snaps them shut (wing salute); or stretches out both one wing to its side and one leg to its rear (wing-leg stretch); or chick is involved in a large, ritualised group display characterised by quick movements of the flock (marching) and exaggerated side-to-side movement of its head in an upright position when neck is erect (head-flagging)

Unknown
Chick performs a behaviour not described by the ethogram

Enclosure Usage
The location of focal chicks within predetermined enclosure zones, adapted from Rose, Brereton and Croft (2018 [54]), were recorded from photographic data. Counts of enclosure zone usage were inputted into the modified Spread of Participation Index (SPI [56] to determine use of available space, given as: where ƒo is the observed frequency of observations in a zone, ƒe the expected frequency of observations in a zone based on zone size and even use of the whole enclosure, N the total number of observations in all zones and ƒemin the expected frequency of observations in the smallest zone.

Flock Position
The positions held by flamingo chicks within their flock were also recorded. A focal chick was considered to hold a 'central' position in their flock if they had counts of ≥5 adults and/or subadults (i.e., non-chicks) within two neck lengths around them at multiple angles; 'outer' if they had ≥5 adult and/or subadults within two neck lengths to one side of them; 'peripheral' if they had 1-4 flamingos within two neck lengths; or 'alone' if there were 0 flamingos within two neck lengths. Chicks were included in the count of the peripheral location to avoid considering chicks who were only associating with other chicks as being alone.

Social Preferences
Partner preferences were determined by identifying the ring codes of the nearest neighbours to the focal chicks. Nearest neighbours were defined as flamingos within a twoneck length radius of the focal chick. Where ring codes were not available, artificial codes were given to individuals which contained information of their age group (foster-chick, parent-chick, juvenile, subadult, adult, or unknown) and sex (male, female or unknown) as determined by the researcher.

Activity, Enclosure Usage, and Flock Position
Behavioural count data were first analysed using RStudio version 1.3.1073 [57]. Several statistical methods of modelling count data were considered to address our hypotheses that disruption of the early rearing environment would influence the development of behaviour and social preferences in Chilean flamingo chicks.
First, we considered fitting Generalized Linear Mixed Model's (GLMM's) with a Poisson error distribution. Counts of each behaviour were summed for each day of observation, for each individual, and tabulated into a new dataset. Exploration of the data indicated that dependent variables were zero-inflated and overdispersed, thus violating the assumption of a Poisson distribution where the mean is assumed to equal the variance.
Using the glmmTMB package [58] we then attempted to fit zero-inflated GLMM's (ZIGLMM"s) with a negative binomial error distribution to the same dataset. ID was entered into the zero-inflation portion of the model, as excessive zeros were believed to derive from individuals not being observed. ZIGLMM's failed to converge.
Following recommended procedure [59], fitting the variation among individuals in the zero-inflation parameter as a random effect rather than a fixed effect did not resolve convergence issues. The models were likely to have been overparameterised in that the data did not contain enough information to reliably estimate the parameters [59].
To maximise the amount of data and increase power, each key state behaviour was treated as a dependent variable with a binomial response to indicate whether the behaviour was performed at the given observation point, rather than a sum of total counts observed in each day. Using the glmer function of the lme4 R package [60], we fitted full conditional (FC) GLMM's with a binomial error distribution. Models contained chick sex and rearing background as categorical predictors. The time and date of observation, temperature, humidity, sunshine, rainfall, daily number of total park visitors, and number of visitors in the Chilean flamingo enclosure were entered as continuous predictors. Chick ID was entered as a random intercept, and the time and date of observation were also entered as a random slope. This random effect parameter allowed for individual behavioural responses to vary over time. We used a logit link function and estimated our binomial GLMM's using maximum likelihood.
All models converged with singularity warnings indicated by a random-effects variance of nearly zero and estimates of correlations of exactly −1.00. Consensus on dealing with singularity issues is mixed and the extent of its impact has been debated. Barr et al. (2013;[61]) suggest to "fit the most complex model consistent with the experimental design, removing only terms required to allow a non-singular fit". A non-singular fit was allowed when the fixed-effect and random slope parameters of time were removed from the models. Since longitudinal analysis requires the inclusion of time parameters, we opted to keep them within our models and accept warnings of singularity as an indication of minimal variance in the random effect. Their inclusion has no effect on any of the estimated quantities [59].
Data for zone occupancy and flock occupancy were treated in the same manner, resulting in binomial GLMM's for individual key state behaviours, zones, and flock positions. The model investigating usage of the island zone did not converge. Models investigating the bush, back grassland, and front water zones, converged but gave the additional warnings indicating that the standard error estimates might be less accurate [62]. Variance Inflation Factors (VIF's) also indicated that the model investigating zone use in the front water and public viewing areas each contained four variables that were overdispersed (i.e.,VIF's ≥ 5). Due to the collinearity between independent variables indicated by high VIF's, these models were not discussed to avoid misleading interpretations. Variables in all other models were not overdispersed (i.e., VIF's < 5) and could be interpreted reliably. A total of 13 full conditional (FC) models were fitted, six investigating chicks' performance of key state behaviours, three investigating chicks' zone usage, and four investigating chicks' position within the flock.
Using an information theoretic approach, we compared each FC model to four other models to reveal the importance of rearing background and time as predictors of behaviour, zone occupancy, and flock position: NULL: A null model containing the intercept and the random slope parameter; 2.
RO: A rearing-only model containing rearing background and the random slope parameter; 3.
TO: A time-only model containing time and date of observation and the random slope parameter; 4.
FC-R: A full conditional model with rearing background omitted.
We ranked the relative support of each model using Akaike Criterion scores (AIC [63]) corrected for small sample size (AICc [64]) calculated using the MuMIn R package [65]. The coefficient of determination (R 2 ) adjusted for mixed models were calculated for each model using r.squaredGLMM function of the MuMIn package [66].
Hypothesis driven model comparisons were carried out using the anova function, revealing whether the addition of specific variables produced a significantly better model fit. Comparing RO to NULL reveals whether knowing information on individual chicks rearing background is useful for predicting behaviour, zone occupancy, and/or flock position. Comparing FC to FC-R reveals whether rearing background explains any unique variance not captured by all other identity, climatic, and visitor variables. Comparing TO to NULL indicates whether patterns of behaviour, zone use, and flock position varied over time.
We then investigated the relative relationship of each significant predictor to its outcome variable within each FC model. We did not necessarily investigate the best fitting models as we were interested in understanding the relationship of all predictors; not seeking to provide a parsimonious model from which to predict future behaviour. Predictoroutcome relationships were inferred through odds ratios; calculated by exponentiating the Beta coefficient estimates.

Social Preferences
Association measures were analysed using Socprog version 2.9 [67] and social network diagrams were constructed using NetDraw [68]. Association indices were defined using a Half-Weight Association Index to account for biases and error during sampling [69].
The Chilean flamingo chicks' network of associations from March to July 2019 were first assessed for social differentiation and standard errors obtained from bootstrapping with 10,000 replications. Social differentiation scores of greater than 0.5 indicate welldifferentiated societies [67]. Modularity analysis was run to determine whether the social network could be divided into smaller communities, with outputs of 0.3 or greater indicating suitability for division [67]. The cophenetic correlation coefficient was also examined to determine the suitability of cluster analysis, with coefficients greater than 0.8 indicating appropriateness of cluster analysis [67].
Permutation tests were run to reveal significant preferred and avoided associations within the network. The number of permutations were incrementally increased until the pvalues for the coefficient of variation (CV) stabilised. Final p-values stabilised at 1000 trials with 10,000 permutations. The number of observed significant dyads were compared to the number of significant dyads expected from a randomly associating network. The observed and expected CV's were also compared to identify long-term preferred and avoided associations. Rearing background was set as a class variable to enable betweenclass comparisons.
Multiple regression quadratic assignment procedure (MRQAP) tests were used to assess how well individual predictors explain the variation in the association matrix, whilst controlling for all other predictors in the regression model. Rearing background, sex, age, climatic factors, and visitor numbers were treated as predictors of the association matrix.
Socprog computed network measures (strength, eigenvector centrality, reach, clustering coefficient, and affinity) were also entered as predictors.
Lagged association rates were plotted against null association rates for both parentreared and foster-reared chicks to assess temporal patterns and longevity of dyadic associations. Quasi Akaike Information Criterion (QAIC [70]) were examined to identify and fit the best explanatory model to the lagged association rate.
Finally, Mantel Z tests were used to analyse the stability of the association matrix over time, by comparing matrices from the first month of sampling to the following months. Stability is inferred if association matrices are correlated. The same procedures were used to assess the stability of associations across behaviours, enclosure zones, flock positions, varying climates, and varying numbers of visitors.

Activity
Differences in the development of chick behaviour by month are shown in Figure 3. Chicks spent most of their time resting, 35%, preening, 24%, and feeding, 25%. The effect of rearing background, sex, climate, and visitors with the random slope, explained 1.8 to 12.4% of the variance in key state behaviours, with rearing background accounting for 0.3% to 4.1% of the variance. Model fit values are displayed in Appendix A. Average monthly time-activity budgets of individual foster-reared chicks, foster-reared chicks as a group (FOS), and grouped female and male parent-reared chicks (PAF and PAM, respectively). White = resting; grey = preening; blue = other; yellow = movement; green = feeding; red = alert. Standard error bars represent variation within parent-reared chicks. Chicks spent most of their time resting, preening, and feeding. Over time, feeding, movement, and vigilance decreased, whereas preening and resting increased. Foster-reared chicks spent less time feeding than parent-reared chicks, overall average time spent feeding 23% and 26%, respectively. Key state behaviour model comparison results are displayed in Table 2. FC models for each key state behaviour performed significantly better than null models, X 2 s (10, N = 8059) ≥ 22.01, p's ≤ 0.015. RO models performed significantly better than null models when explaining feeding and movement behaviour, X 2 s (1, N = 8059) ≥ 4.37, p's ≤ 0.037. FC models performed significantly better than FC-R models when explaining feeding behaviour, X 2 (2, N = 8059) = 6.21, p = 0.045. TO models performed significantly better than the null model when explaining preening and movement behaviour, X 2 s(1, N = 8059) ≥ 11.24, p's ≤ 0.001.

Enclosure Usage
Differences in the development of SPIs and zone usage by month are shown in Figure 4. SPI values were high and showed an increasing trend over time, indicating discriminative use of enclosure zones. Foster-reared and parent-reared chicks spent most of their time in the nesting area, 60% and 55%, respectively. The effect rearing background, sex, climate, and visitors with the random slope, explained 9.2% to 20.7% of the variance in usage of the nest, shelter, and mid-water areas, with rearing background accounting for 1.9% to 10.5% of the variance (Appendix C). Overall and monthly Spread of Participation Index (SPI) for individual foster-reared chicks, foster-reared chicks as a group (FOS), and grouped female and male parent-reared chicks (PAF and PAM, respectively). Grey = overall; blue = April; green = May; yellow = June; orange = July. Standard error bars represent variation within parent-reared chicks. (b) Average monthly time-zone usage budgets of individual foster-reared chicks, grouped foster-reared chicks (FOS), and grouped female and male parent-reared chicks (PAF and PAM, respectively). Dark grey = shelter; white = public viewing area's; light grey = nesting area; yellow = mid-waters; orange = island; blue = front water; green = back grasslands; red = bushland. Standard error bars represent variation within parent-reared chicks. Chicks spent most time in the nesting area. Over time, occupancy in the nesting area increased, and occupancy in the middle-water and shelter area decreased. Foster-reared chicks occupied the nesting area more than parent-reared chicks, overall average of 60% and 55%, respectively.
Zone usage model comparisons are displayed in Table 3. FC models performed significantly better than null models, X 2 s (10, N = 8059) ≥ 106.65, p's ≤ 0.001. RO models performed significantly better than null models when explaining usage of the nesting area, X 2 (1, N = 8059) ≥ 9.50, p's = 0.002. FC models did not perform significantly better than FC-R models when explaining occupancy in the nesting area, X 2 (2, N = 8059) ≤ 3.63, p's ≥ 0.0163. TO models performed significantly better than the null model when explaining usage of shelter and mid-water areas, X 2 (1, N = 8059) ≥ 8.03, p's ≤ 0.005. The influence of predictors on zone usage, indicated through binomial GLMMs, are displayed in Appendix D. Foster-reared chicks were 1.281 (95% CI: 1.099, 1.494) times more likely to occupy the nesting area than parent-reared chicks, β = 0.248, SE = 0.078, p = 0.002.

Flock Position
Differences in the position of chicks within the flock by month are shown in Figure 5. Chilean flamingo spent on average 92% of their time associating with other flamingos, spending most of their time in peripheral, 45%, and outer, 33%, positions to the flock. The effect of rearing background, sex, climate, and visitors with the random slope, explained 2.8% to 8.6% of the variance in flock position, with rearing background accounting for 0.6% to 2.1% of the variance (Appendix E).

Social Preferences
Networks of associations of parent-reared and foster-reared chicks are displayed in Figure 6. The overall mean association rate of both foster-reared (Mean Assoc. = 0.08, SD = 0.08) and parent-reared chicks (Mean Assoc. = 0.10, SD = 0.07) indicated that their associations within this network were weak; but large group sizes may distort this . Estimation of social differentiation using a likelihood method revealed that this network was well differentiated (CV = 1.139, SE = 0.072) as defined by Whitehead (2009) and the power to detect the true social network was moderate (P = 0.494, SE = 0.006). Modularity analysis indicated that the network could not be divided into separate communities (Modularity = 0.079 for 13 clusters). Furthermore, a cophenetic correlation coefficient of 0.778 did not justify the use of cluster analysis. Permutation tests revealed that foster-reared chicks held less significant dyads than expected (obs = 78, exp = 133.2) and parent-reared chicks held more significant dyads than expected (obs = 155, exp = 125.25). See Appendix G for significant avoided and preferred associations. The observed CV was significantly greater than the expected CV for both foster-reared and parent-reared chicks, indicating that Chilean flamingo chicks were preferentially associating or disassociating with other flamingos, CV observed = 1.908, CV expected = 1.886, p = 0.001; CV observed = 1.644, CV expected = 1.526, p < 0.001.
A Mantel test revealed that patterns of associations were significantly different between rearing classes, t = −4.632, p < 0.001. MRQAP tests also revealed that rearing background predicted the association matrix whilst controlling for sex, age, and network measures, r = −0.215, p < 0.001. Sex and age were not significant predictors of the association matrix, r = −0.0347, p = 0.066; r = −0.036, p = 0.262. Amongst the network measures, only eigenvector centrality was a significant predictor of the association matrix, r = −0.048, p = 0.009.
Differences between the plotted lagged and null association rates indicated that foster-reared and parent-reared chicks were both showing preferred associations over time ( Figure 7) and that the strength of a typical dyadic association is similar between the two groups of chicks. The best fitting model (i.e., model with lowest QAIC value) explaining the lagged association rate of foster-reared chicks to all other flamingos indicated two levels of casual acquaintances. The best fitting model explaining the lagged association rate of parent-reared chicks to all other flamingos indicated rapid dispersal and casual acquaintances. Other models also showed support (see Appendix H).
Results from Mantel tests are displayed in Table 5. The overall associations matrix from the first month of sampling consistently and accurately correlated with all other monthly association matrices, suggesting that the association patterns of chicks were consistent over the sampling period. Behavioural association matrices also accurately correlated with one another. However, the partial correlation coefficient reduced when looking at the associations when chicks were alert or performing other behaviours, indicating that associations varied across some behaviours. Similarly, association matrices observed in different enclosure zones accurately correlated but the strength of correlation was not consistent, indicating that associations varied across zones. Associations observed in public viewing areas did not correlate with associations elsewhere in the enclosure. Associations across central, outer, and peripheral locations relative to the flock strongly and accurately correlated; indicating that chicks moved around and within the flock with a consistent set of associates. As expected, the association matrix of when chicks were alone did not correlate with the association's chicks held when they were closer to the flock. Patterns of association across visitor and climatic variables strongly and accurately correlated with little variation.  Lagged and null association rates of a typical (a) foster-reared and (b) parent-reared Chilean flamingo chick dyad. Lagged association and null association rates are plotted using a moving average of 10,000 recorded associations. Best fitting models indicated two levels of casual acquaintances for foster-reared chicks, and rapid dispersal and casual acquaintances for parent-reared chicks. Differences between the null and lagged association rate indicate that chicks held stronger bonds than expected from random association.

Discussion
Our results show no major differences between the juvenile Chilean flamingos reared as chicks by Andean flamingo foster-surrogates when compared to juvenile Chilean flamingos reared by their own parents. Such zoo husbandry interventions can be successful and enhance the welfare of individuals involved (e.g., in this case enabling a full breeding cycle to experience by Andean flamingos) if the ecology and behaviour of all species involved is carefully considered, and the intervention and output carefully monitored. Comparison of full conditional (FC) against null models indicated that parent-reared and fostered chick behaviour, enclosure usage, and flock positioning were non-random. Rearing background explained a small proportion of variance in behaviour. Comparisons of rearing background models (RO) against null models suggested that the predictive use of rearing background was limited to feeding, movement, and occupancy in the nesting zone. Comparison of FC models to FC models without rearing background (FC-R) revealed that the proportion of variance in movement and nest occupancy explained by rearing background did not significantly improve model fit. Despite this, the fixed effect of rearing background on nesting occupancy within the FC model was significant. Although the variance in nest occupancy explained by rearing background was small, this variance was unique and not encapsulated by other predictors in the model. Conversely, rearing background did not influence movement; the initial variance explained by rearing background could be better explained by other predictors in the model. For feeding behaviour, rearing background significantly improved model fit and significantly influenced time spent feeding.
The Chilean flamingo chicks formed a highly variable and well differentiated network of associations, indicative of non-random, discriminative patterns of association (Figure 6a). The filtered network shows that all chicks held several preferred associations (Figure 6b). Although foster-reared chicks showed fewer significant preferred associations, differences the null and lagged association rates (Figure 7) show that these bonds were equally strong and stable as the bonds that parent-reared chicks formed. Both groups of chicks formed long-lasting preferential bonds and demonstrated several short-lived relationships.
Chicks showed an impressive ability to maintain their pattern of associations across the sampling period, behavioural states, enclosure zones, positions within the flock, varying climates, and varying numbers of visitors. Chosen associates significantly differed between foster-reared and parent-reared chicks. Rearing background predicted the association patterns, but the sex or age of a bird had no influence. Foster-reared chicks associated with other foster-reared chicks, but also formed relationships with other flamingos. Eigenvector centrality was also a predictor of the association matrix, as birds with more associates are more likely to be connected to the well connected.
Patterns of behaviour (Figure 3), enclosure usage (Figure 4), and positioning within the flock ( Figure 5) were also non-random. Chicks spent most of their time feeding, resting, and preening, occupying the nesting area, and occupying peripheral and outer locations to the flock. Foster-reared and parent-reared chicks' behavioural states, use of the enclosure, and positions within their flock over the sampling period showed complex relationships with temporal, identity, and environmental factors. However, early rearing experience did not explain a significant amount of variation in most behaviours. The effect of rearing background on Chilean flamingo chick behaviour was limited to feeding and use of the nesting zone. Independent of their sex, foster-reared chicks had greater odds of occupying the nesting zone than parent-reared chicks and were also less likely to be observed feeding.
Taken together, our results show that altering the early social rearing environment of captive Chilean flamingo chicks, through foster-rearing intervention with Andean flamingos, is associated with few behavioural and social differences relative to parent-reared conspecifics within the same flock. Our findings conflict with previous studies where cross-fostered birds display numerous behavioural similarities to their foster species [30,31,34] and develop altered social preferences [23,39]. Instead, both foster-reared and parent-reared Chilean flamingo chicks were able to express patterns of activity and association similar to their captive and wild counterparts [47,50,51,71].
Cross-fostering studies are typically used to assess the relative influence of genetics and the environment [72]. Experimental procedures utilise two species that are closely related enough that fostering is successful, but also dissimilar in observable aspects of behaviour so that misimprinting is easily recognisable (e.g., vocalisations [34]). The resulting behavioural differences between cross-fostered and parent-reared conspecifics are therefore the result of pre-existing behavioural differences between the two species as a whole. It follows, then, that the numerous similarities in various aspects of the Chilean flamingo chick's behaviour can be largely attributed to the ecological similarities between Andean and Chilean flamingos.
Previous social network analyses on the Chilean and Andean flocks involved in our study indicate that both show similar non-random patterns of assortment [47]. Individuals formed strong preferential bonds which were maintained in and outside of the breeding season [47]. The social dynamics of the foster flock were therefore a good representative of the biological flock. The early social rearing environment equipped foster-chicks with relevant experiences that facilitated their ability to form long-lasting bonds with other flamingos whilst avoiding others.
Both parent-reared and foster-reared chicks showed a preference for the nesting area ( Figure 4b). Over time, chicks were less likely to be found alone or on the periphery of the flock and were more likely to be associating with five or more older flamingos. Foster-reared chicks also occupied the nest area more than parent-reared chicks. Captive greater flamingos have also shown biased use of the nesting area when studied during the breeding season, as in this study [73]. Flamingos naturally group together and are more likely to congregate in the nesting zone during breeding season [49,74]. By assembling together into the popular nesting area, chicks can settle into the flock and benefit from the welfare advantages of flocking (e.g., foraging efficiency, predator detection and avoidance, access to mates [75][76][77]. Wild Andean and Chilean flamingos inhabit intersecting geographical regions in South America [51,52], utilising similar wetland resources. This ecology is represented in their enclosure design at this animal collection. The two enclosures at WWT Slimbridge (South American Pen and the Andean Flamingo Pen) incorporated biologically relevant features to both species providing, as best possible, opportunities for normal time-activity patterns to be performed. Pre-existing preferences developed during early rearing mean that when relocating, birds should prefer to settle in habitats similar to their natal foster-rearing site [78,79]. When later life environments are different to the natal rearing environment, cross-fostered young show negative altered behaviour patterns that reduce fitness [39]. Dispersing birds are under strong selection pressure to quickly settle into a habitat, and so make discriminative choices based on their natal experience [78,80]. Once settled, individuals are then predicted to increase activity in the breeding area [81]. For instance, after dispersal, cross-fostered pied flycatchers show preferences to their foster species habitat and choose to breed there later in life [82].
Andean and Chilean flamingos spend most of their time preening, resting, and feeding [48,70]; analogous to chicks behavioural patterns (Figure 3). Less time spent feeding by foster-reared chicks may initially seem concerning since reduced feeding is linked to deferred maturity and higher mortality in birds [83]. Wild juvenile Chilean flamingo chicks feed less than adults due to aggressive displacement [84]. However, the abundance of food appropriately dispersed across the Chilean flamingo enclosure and freedom of movement would minimise aggression between birds [85]. This is seen in other wild flamingo flocks; for instance, the superabundant supplies of food at Kamfers Dam in South Africa reduced the effects of any lost feeding time by allowing lesser flamingos to spread out and avoid aggressive interactions [75]. Foster-reared Chilean flamingo chicks, motivated by hunger, could freely move to non-competitive zones to feed. Instead, greater foraging efficiency requires individuals to spend less time feeding [84,86]. Generalist feeding birds crossfostered to specialist feeding birds have previously shown foraging patterns that align with their foster species and that are more efficient than parent-reared controls [35,87].
It is not well-known if feeding behaviours of the generalist feeding Chilean flamingo are shared by the specialist feeding Andean flamingo 49]. Andean and Chilean flamingos show different food preferences which enable them to coexist together in the wild [45,46]. Differences in bill structure [88] and an understanding that flamingos display distinct feeding behaviours that are relevant to their ecological niche (Rooth, 1976 as cited in 49]) suggest it is likely that feeding behaviour varies between Andean and Chilean flamingos. Understanding species-specific feeding behaviour of the understudied Andean flamingo would help clarify if distinct foraging niches exist and if these can be culturally transmitted into non-biological offspring.
Whilst we present these results to show the efficacy of a husbandry intervention, such information can benefit ex situ conservation efforts by encouraging increased reproductive output in populations that may need assistance for breeding. Conservation efforts may require animals being moved between groups to ensure genetic management of the population [89]. For captive flamingos, this is almost a certainty, since zoos can no longer remove flamingos from the wild and their current captive populations are not self-sustaining 49]. Foster-rearing presents a way to address this reproductive dilemma. Introducing reproductively viable individuals into foreign populations increases genetic diversity and reduces inbreeding [22]. Providing inexperienced individuals with paternal experience also increases future rearing success in flamingos and removing eggs from breeding pairs encourages further egg production 49]. The individuals that can benefit from foster-rearing are vast and understanding the effects are imperative to its success. The evidence-based approach is key here; assessing the ecology of the selected species for cross fostering and monitoring development (i.e., growth and behavioural) in young animals to ensure any impacts on adult characteristics and activity patterns are not detrimental or long-lasting.
Cases against the use of foster-rearing are valid, e.g., [23,40,42]. Although our small sample size restricts generalisability, we suggest that foster-rearing husbandry practices need to be implemented on a case-by-case basis to suit individual welfare needs. Rather than generalising knowledge across species, future fostering interventions should be informed by comprehensive habitat, behavioural, and social assessments of the specific animals involved; analogous to reintroduction programmes [26,[90][91][92]. Previous experiences place specific demands on individuals that must be factored into husbandry provision in captivity [89]. Further research assessing the reproductive behaviours of the foster-reared chicks, once they reach sexual maturity, would help reveal the propagative viability of these captive foster-reared flamingos. It would also be beneficial to study patterns of sociality over longer sampling periods since captive Chilean flamingos show looser association patterns outside of the breeding season than Andean flamingos [47].

Conclusions
Our study indicates that foster-reared Chilean flamingo chicks can develop and integrate well if pre-existing behavioural differences between flocks are minimal. We also show that success of integration is improved if individuals are reintroduced into social settings and habitats similar to their natal rearing environment. Our findings highlight the importance of assessing the behaviour of the specific individuals involved and the habitats in which they reside, prior to implementing fostering intervention. Understanding how to successfully apply foster-rearing husbandry has the power to strengthen the sustainability of captive species by increasing reproductive outputs without impacting on animal welfare.

Funding: This research received no external funding.
Institutional Review Board Statement: Ethical approval was obtained prior to data collection by the Exeter Psychology Ethics Committee, 12 February 2019 (Application ID: eCLESPsy000361 v1.1). This project's methods were also reviewed by the Animal Welfare & Ethics Reviewers at WWT.

Informed Consent Statement:
Not applicable for studies not involving humans.
Data Availability Statement: Data can be made available upon reasonable request from the corresponding author.
Acknowledgments: Thank you to Ruth Cromie and Michelle O'Brien for project review at WWT. Thank you to Mark Roberts and Phil Tovey for assisting with the banding of flamingo chicks and the movement of foster chicks from the Andean flamingo flock to the Chilean flamingo flock, and for providing information on the husbandry, management and foster rearing process. Thank you to three anonymous reviewers for their helpful and development feedback on the original manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. Models are ranked highest to lowest by Akaike weights within each key state behaviour with full conditional models highlighted for each behaviour. AICc = calculation driven by log likelihood but penalises models for every additional parameter used. ∆AICc = AICc score of the respective model minus that of the best model. AIC weights = the likelihood that the respective model is the best in the set. K = number of parameters in the model. Models are ranked highest to lowest by Akaike weights within each key state behaviour with full conditional models highlighted for each behaviour. AICc = calculation driven by log likelihood but penalises models for every additional parameter used. ∆AICc = AICc score of the respective model minus that of the best model. AIC weights = the likelihood that the respective model is the best in the set. K = number of parameters in the model. Models are ranked highest to lowest by Akaike weights within each key state behaviour with full conditional models highlighted for each behaviour. AICc = calculation driven by log likelihood but penalises models for every additional parameter used. ∆AICc = AICc score of the respective model minus that of the best model. AIC weights = the likelihood that the respective model is the best in the set. K = number of parameters in the model. DNC = did not converge.