Using Geometric Morphometrics to Quantify Variation of Shape and Magnitude of the Pattern of Milk Production of Dairy Cattle

Show more

1. Introduction

Milk produced by dairy species is important worldwide in providing food, employment, and income. Driven mainly by the income factor, researchers, following the requirement of dairy industries at some countries, have developed the mathematical representation of milk production during lactation. Mathematical models are the unparalleled data analysis tool generating knowledge for the progress of the dairy industries, and the past 90 years have seen their development and evolution under different setting (for details of the different mathematical functions see the review by [1] - [4] ).

Those mathematical models can describe and predict, with minimum error, the milk yields of lactation records, commonly of 305-day lactation length and two different phases, ascending and descending. They have been also used to investigate other aspects of the cow’s lactation including evaluating its variation by the influence of factors and variables related to management, environment, and physiology [5] - [11] ; detecting different shapes of lactation curves [12] [13] ; modeling extended lactations [14] - [17] ; and quantifying the shape variability of lactation models [18] [19] . Milk yield (or magnitude) and the shape (or pattern), which are determined by the physiological process of the cow’s mammary gland, are intrinsic features of any lactation curve. While mathematical models can successfully quantify the milk yields of lactation records, they fall short to appropriately describe shape, the projected surface on a plane by graphical representation. First, the parameters representing shape, often, have no direct inference to it or are associated to milk yield (magnitude) and lactation stage (time), which strictly speaking are not units of shape (linear m or areal m^{2}). Second, it is constructed fixed, which neglects important effects about the management, environment, and physiology of cow on shape.

The current approach to describe shape of milk production can be compared to a designer that is designing clothes without the person’s measures of body, and instead he is using the person’s weight and age. Working with those data would be of more confusion than usefulness because the subject is described as a blob, of known weight and age, but unknown shape. For the purpose of the tailor, use of body shape measures is the appropriate. As for the development of dairy cattle industry, regarding any lactation pattern as any other biological shape and being able to appropriately measure it, may not only complement to the analysis by mathematical models, but also reveal profound information about the biological process represented through the shape. The well- known shape of lactation curve is the phenotypical expression of the physiological process of milk secretion [4] : the ascending phase, from calving to peak, is the result of increase in number of epithelial cells and their proliferation, and the descending one, from peak to dry off, decrease in cell numbers by apoptosis [20] .

The appropriate analysis of shapes in this field has been hindered by difficulty and lack of a standard method [19] . However, the geometric morphometric (GM) methods, which provide analyses for quantifying, testing, and visualizing shape variation and its covariation with biotic and abiotic variables [21] - [23] , has already proved to contribute to increase the scientific rigor in the description of shape on Biology, Medicine, and Engineering [24] [25] . GM, a relatively new methodological achievement by scientists on Biology, analyzes the relative position of anatomical landmarks and sets of point representing outlines and surfaces to quantify size (magnitude) and shape [25] . The GM can be a powerful tool of analysis for measuring variation and co-variation of milk production patterns of data on individual lactations, farm, and whole industry. The objective of this study was to apply the landmark-based GM to quantify variation of shape and magnitude on the milk production pattern of dairy cattle. We used free software and a dataset from two leading dairy industries of distinct farming system and geographical location to present the application of the new methodology. The selected samples corresponded to monthly time series data of aggregated production from 2007 to 2015 of New Zealand and United States.

2. Theoretical Assumptions on Depiction and Data of Milk Production Shape

2.1. Depiction of Milk Production

Milk production data exhibit two intrinsical features, magnitude and shape; the former refers to the quantity expressed in units of mass (kg) or volume (L) and the later states the relationship between milk quantity and time (projected surface on a plane by graphical depiction). To measure the lactation curve shape, phenotypical expression of this physiological process, we have developed a graphical depiction of the data. Line graph has been essential in visualizing data for fitting lactation models (Figure 1(a)). Our intuition tells us that the pairwise relationship between milk yield and time could be represented in a circular way, since they are both cycles. A radar graph (Figure 1(b)), hereafter called orbital graph, can be alternatively used for depicting the milk production data. This choice of graphing presents the benefits of displaying a closed curved required by the landmark based GM analysis, ordering the time at milk production to the counterclockwise rotation of the earth spinning on its tilted axis as it orbits the sun to produce the seasons, and putting it into position to be transformed to a polar

Figure 1. Depiction of a hypothetical example of the lactation curve of cattle (in m^{3}) by line graph (left) and orbital graph (right). Lactation length is 12-month with 2 months dry off of cow; calving is assumed to be in January. Production is order to the yearly cycling, starting in January and following a counterclockwise rotation.

graph for further analysis. The orbital graph consisted of six axes to account for the twelve months of the year. The monthly data were plotted into the graph and the order of the months, calendar basis, was rotated counterclockwise to relate the time cows produced milk to the season.

The milk production data are converted to cubic meter to use its meter-side-length that is in accordance with the shape unit (linear or areal). For instance, the quantity of milk 9 m^{3} in January can be thought, for the objective of providing a linear measure, as a rectangular solid of 9 m length, 1 m height, and 1 m width; then the 9 m length of the rectangular solid provides the linear measure that is in accordance with units of shape. Each of the 12 linear measures in a year comprises the two-dimensional (2D) shape of milk production, which can be now treated as any other biological form and a group of shapes of milk production can be quantified with the appropriate analysis.

2.2. Data

Different types of data on milk production are generated in dairy industries. They are commonly found of complete individual and aggregated lactations, single farm and region, and whole industry. Such data can be recorded daily, weekly or monthly. For instance, the test-day-model system usually records milk quantity at month interval during the entire lactation of a cow. Such continuous records can be grouped into lactation cycles; then each lactation cycle can be viewed as an individual. The lactation curve represents an individual (a cow). This framework of representing an individual by a curve was first introduced by Kirpatrick and Lofsvold [26] as the infinite-dimensional approach to quantify growth trajectories, morphological shapes, and reaction norms in quantitative genetic. In growth trajectories, an individual is represented by measures of size at different ages, in which the size of the individual for each different age is treated as a different character (but correlated), which overall has been found to provide better quantifications of growth curves [26] - [29] . Similarly, we treated the lactation curve as an individual, the quantity of milk of a cow for each different month is viewed as separate (but correlated) character.

For instance, the milk quantity of a cow at January (time) can be viewed as a different character from the milk quantity of the same cow at February. This framework can be also applied to monthly data on milk production of a farm, region, and industry, and consideration should be given to the different factors influencing each data type.

3. Materials and Methods

3.1. Data Collection and Description of the Selected Dairy Industries

The complete lactation of cow has been the metric unit to develop the current analyses (lactation models) of data in the field. Those lactation models are the data analysis tool and provide the theory, concept, and idea from the research findings that are then applied to production records of farm or whole industry (e.g., [30] ). Magnitude and shape of the pattern of milk production from different type of data such as individual and aggregated complete lactations as well as farm, region, and whole industry can be used to be analyzed with the new methodology. Since data on cows’ lactation were not available to us and for the purpose of demonstrating the application of the new methodological framework to analyze the pattern of milk production (regardless of the data type) under different setting, we used as a study case New Zealand and United States dairy industries. The industries from these two countries were selected because are leading producers at different latitudes and have developed distinguishable farming systems (pasture-based and confinement) that derive different milk supply (seasonal and uniform).Monthly data on aggregated production, 2007-2015, were collected from the Dairy Companies Association of New Zealand [31] and the United States Department of Agriculture [32] . The selected periods were the only monthly data available for New Zealand. The data for the United States included the 23 major producer States, which comprised 92% of the country’s total production, and that of New Zealand included the 17 regions which comprised 100% of the country’s total production.

The weather conditions were typical climatically for the 9 years under study (temperature above the average and rainfall below the average); exceptions were 2008 and 2013 for New Zealand and 2011 for the United States that reported the incident of droughts. In New Zealand, the main effect of drought on milk production is through decreasing the pasture growth, which is main feed supply. In the United States, while drought also reduces pasture, corn, and soybean growths, the main effect on milk production is the heat stress on dairy cows. In addition, another effect is the lower reproduction rates, which causes problems in the subsequent production year.

The New Zealand and United States industries vary in size and structure, as shown in Table 1. Crossbred cows in New Zealand are managed grass-based, for which cows are scheduled to calve in a concentrated pattern, late-winter to early-spring [33] , to match the grass-growing season. This dairy industry presents peak production of milk from September to November and trough from May to August. High breed Holstein cows in the United States are mostly in house all-year-round and supplemented total mixed ration of stored forages. This management system is input intensive dependent that rewards with high yield of milk. The industry, though less pronounce, presents two calving patterns: spring and fall, in which a greater proportion of cows calving in the former season is reported and which overall result in some more milk in spring than in fall [34] [35] . The supply of milk derived from pasture and confinement operations are seasonal-based in New Zealand and almost uniform in United States.

3.2. Geometric Morphometrics (GM)

Morphometrics is overall defined in the biological-shape literature as the measurement of shape variation and its covariation with other variables [22] [24] . This relatively new method has undergone changes in the way it is applied. It has been employed in biology, medicine, and engineering [24] [25] to capture the shape of 2-dimen- sional (2D) or 3-dimensional (3D) structures images using specifically defined homologous landmarks [24] .

In the theory and measurement of morphometric, it is commonly distinguished

Table 1. Characteristics of dairy industries in New Zealand and United States.*

Calculated with data obtained from [32] [36] ; *The last five years were selected because it was believed to vary less than including the 9 years.

between the size (magnitude) and shape of a biological structure. Shape corresponds an outline with landmarks after differences in location, scale, and orientation has been removed [37] ; size refers to the magnitude in an area or volume of a structure, often measured through the centroid size: the square root of the summed squared distances between all landmarks and their centroid [38] ; and form includes both the shape and size of an object [24] .

GM relies on multivariate analysis such as General Procrustes Analysis, Principal Component Analysis, ANOVA, MANOVA, Regression, and other tests. Several freeware programs to conduct analysis on GM are available such as Morpho J [39] , Morphologika [40] , PAST [41] , geomorph from R [42] . We used Morpho J [39] because previous studies provide concise explanation through examples of how to use this program.

3.3. Landmark Identification

Landmark-based GM techniques begin with the collection of 2 or 3D coordinates (which are based on the x, y, and z Cartesian plane) of biologically definable landmarks [22] . Landmarks are special points associated with the contour of a biological shape [37] . In GM, the acquisition of landmarks is done by the aids of cameras, scanners or microscopes to obtain an image and locate landmarks. Further discussion of landlarms can be found in Zelditch [43] . For our case, the landmarks associated with the contour of the shape were already identified, the monthly quantities of milk production in a year in cubic meters unit (see Figure 1). A set of coordinates based on the x and y Cartesian plane were estimated for each of the 12 landmarks.

The coordinates were zero on the x-axis and the respective milk quantity on the y-axis for the January and July landmarks (see Figure 1), with positive sign for the former month and negative for the later one (it was to account for the orientation). Similarly, the coordinates were the respective milk quantity on x-axis and zero on the y-axis for the April and October landmarks, with negative sing for the former month and positive for the later one. The coordinates for landmarks on February, June, August, and December were calculated using Pythagorean Theorem. We knew the diagonal which was the monthly quantity of milk. From the diagonal, we drew a rectangle connecting the x- and y-axis (see dotted lines in Figure 2).

The calculation yielded the width and length of the rectangle which corresponded to the x and y coordinates, respectively. Similar calculation was done to obtain the coordinate values of the March, May, September, and November landmarks with the difference that the length of the rectangle represented the x-axis and the width the y-ones. The negative signs were also added to either the x- or y-axis depending on the direction on the plane.

3.4. General Procrustes Analysis (GPA)

This is the most widely used superimposition method in GM (others: full Procrustes superimposition, resistant-fit method, Bookstein registration, and sliding baseline registration [21] ), and once a configuration of x, y landmarks have been obtained, the first step is the application of GPA (also called Generalized Least Squares). It is applied to the landmarks configuration to separate the shape from size, and remove unwanted effects of differences in specimen position and rotation (see Figure 3) [38] , it superimposes landmark configurations using least square estimates for translation and rotation

Figure 2. Cartesian plane showing the rectangles drawn from the diagonal point on the February-and March-axis.

Figure 3. Three steps (translation, scaling, and rotation) Procrustest superimposition. The example presents four configurations of 12 landmarks each (New Zealand monthly pattern of milk production in a year). The raw landmarks configurations in (a) are translated to a common centroid (b). The centered configuration then are scaled to the same centroid size (c) and iteratively rotated until the summed squared distances between the landmarks and their corresponding sample average position is minimum (d).

parameters [22] .

Removal of location differences is done by centering configurations, for which the centroid of each configuration is calculated and used as the origin of the new coordinates system. The centroid of a configuration is the mean values of the x- and y- (and z-) coordinates for all k landmarks in the configuration. Removal of size differences between configurations is done by rescaling each configuration so that they share a common centroid size of 1. Centroid size is the square root of the sum of the squared distances between each landmark and the centroid of the form [21] [23] . Removal of orientation differences between two configurations is done by rotating one configuration (the target form) around its centroid until it show minimal offset in location of its landmarks relative to the other configuration (the reference form). The three steps process produces a new set of shape variables, Procrustes shape coordinates, which contain information only about the shape of the configurations, then suitable for statistical analysis. For details on the mathematical notations of the process of Procrustes super imposition see [44] .

3.5. Principal Component Analysis (PCA)

The Procrustest shape coordinates obtained by GPA is of a high dimensionality, which makes it difficult to visualize and interpret. PCA is a method that summarizes multivariate data by building linear combinations of the original variables that are uncorrelated and maximizes the sample total amount of variance explained [25] . It was employed on the variance-covariance matrix of the GPA shape coordinates to summarize variation in shape.

3.6. Visualization of Shape Variation

There are several options to visualized shape variation in a sample, such as thin plate spline, displace vector (called lollipop graph in MorphoJ), and wireframe. In our study, the last option was adopted because it was commonly observed in previous studies. Wireframe diagrams are lines that connect landmarks in to resemble a 2 or 3D shape.

3.7. Testing Variation within and between Dairy Industries

The previous methods (PCA and visualization) are largely exploratory and others statistical analysis are employed for testing hypotheses [24] . GM employs several techniques for this purpose such as Procrustes ANOVA, MANOVA, discriminant function analysis (DFA), canonical variate analysis (CVA).

Differences in the shape and size of the production pattern can occur within and between dairy industries. To test within country variations, we used one factor ANOVA for repeated measures, because the set of landmarks associated with the shape are generated from the monthly milk produce in a year by the same dairy herd. The calculation was done separately for each country using the freeware Real Statistic Resource Pack for Excel [45] . As to statistically compare the production patterns between the two industries, the DAF of MorphoJ is specially used when there are only two groups, providing a parametric t-square for the difference between means and if requested a permutation test.

4. Results

Followed from the depiction on orbital graph is the result of the geometric mophometrics assessment of shapes.

4.1. Orbital Graph

Given the farming system, geographical region, and the genetic merit and the lactation stage of the cows’ population, the New Zealand and United States dairy industries projected production patterns, hereafter named as, of cardioid and heart shapes, Figure 4 (NZ) and (US). Each shape represents the production of milk in a year; January is the starting month following a counterclockwise rotation to complete a year. The development of the cardioid shape of New Zealand may be easier to explain, since most of the farms follow a common rhythm of production: seasonal-based that begins in June and ends in May of the following year. Cows calve in a concentrated pattern during late-winter to early-spring, August-September, which ascends milk production from therein to a maximum in October followed by a smooth decrease to a minimum in Jun-July, which is caused by the dry off of the herd in previous months.

Figure 4. Monthly milk production of New Zealand (NZ) and the United States (US) (in 1000 m^{3}).

On the other hand, the development of the heart shape of the United States may be more difficult to explain, because it is the result of production by dairy farms distributed in seven geographical regions, among which, each one presents its own environmental conditions, adopts its own management practices, and develops its own pattern of production. Overall, a general pattern can be cast. Yet production in the United States is more uniform throughout the months, the seasonal pattern can still be seen. The greater proportion of cows calving in spring ascends milk production to a peak, from March-to-May, followed by two slightly-gradual decreases, first from June-to- August and second from September-to-October, and an increase of up to the June-to- July level, from November to December (caused perhaps by another proportion of cows calving during autumn, October onward, which may be practiced by some farms). The cleft of the heart shape (February), which is the minimum point in the year, may be the result by those cows calving in spring and that are dried off during December and January.

The influence of environment on milk production was exhibited. The peak production in both countries occurred in their respective grass growing season, observed on the bigger upper half of the March-September axis for New Zealand, and the lower one of the same axis for United States (Figure 4 (NZ) and (US)). The drought effect on the decreased milk production of New Zealand in 2008 and 2013 can be seen on the smaller left half of the January-July axis, compared to that of the other years; in the United States, the drought effect on the decreased milk production of 2012 can be seen in the overall smaller shape, compared to that of the other years. The two forms show how varies the milk production pattern, either widening or narrowing within and among the years. To make inference of the changes, in size and shape, over time, GM tools were best suited to summarize, visualize, and test the pattern of variations, which are presented in the proceeding sections.

4.2. General Procrustes Analysis

GPA superimposed the configuration of landmarks, then suitable for the statistical analysis commonly done in GM studies. The procedure produced an average shape with the scatter of landmark positions around the overall mean. Figure 5 shows the mean shape of New Zealand’s cardioid and United Sates’ heart shapes (the bigger dots represent the mean, and the smaller ones the scatter landmarks). The scatters of landmarks around the mean for the two shapes were not circularly positioned (Figure 5(a) and Figure 5(b)), which questioned the isotropic assumption: landmarks circularly distributed around the mean [44] [46] . Consequently, a test of significance that does not assume the isotropic variation would be more appropriate than the ANOVA for repeated measures.

4.3. Testing Variation of Size and Shape within and between Patterns

The assumption of sphericity variation at each landmark position was questionable because the scatters of landmark positions around the overall consensus were not circular (Figure 5). Therefore, the Hotelling’s Test was more appropriate and alternatively used

Figure 5. Mean landmark configuration for the New Zealand’s cardioid (a) and United States’ heart (b). The bigger dot indicates the mean (n = 9), while the smaller ones are the scatter values of the landmarks.

over the ANOVA for repeated measures. The results for size and shape variations are reported in Table 2. The test for one sample indicated that individual variation in size was statistically significant for the New Zealand’s cardioid and United States’ heart shapes. Conversely, the Hotelling’s Test indicated that individual variation in shape was statistically non-significant for the two countries (Table 2).

On the other hand, similar results to the significant of size within countries were obtained for their comparison (Table 3). The parametric and permutation test of the DAF indicated that size differences between the two countries were statistically significant. This could have been expected since annual milk production of the United States doubles that of New Zealand. Similarly, shape was statistically significant between the two countries. The production pattern is seasonal based in New Zealand, while it is mostly even throughout the year in the United States.

4.4. Principal Component Analysis (PCA)

We only present the summary for the pattern of individual variations, though the non-significance of shape, because for that of the cross comparison will produce an average shape from the cardioid and heart; this may not be appropriate since they correspond to two completely different forms of milk production. The pattern of variation showed by the PCs visualized the major changes in milk production over the 9 years under study; it also showed the specific time of the year that such changes occur. In the New Zealand, the shape variations were described by the first four PCs, which accounted for 93.1% of total variance in the cardioid shape. The analysis of variation of production pattern (the wireframe depiction in Figure 6 (NZ)) shows that the PC1 (black line and solid circle) explained 58.8% variation. Compared to the mean shape (gray line and open circle), the cardioid production shape of PC1 has a slightly narrower

Table 2. Variation within centroid size and shape of the New Zealand’s carioid and United States’ heart.

Table 3. Differences between the means of cardioid and heart shapes using discriminant functional analysis.

Figure 6. Patterns of the shape variations in New Zealand (NZ) and United States (US).The diagrams show the first four principal components (PCs) in New Zealand and the first two PCs in United States of the covariance matrix for variation among individual; that is, variation in the average (n = 9) of the individual. For each PC, the diagram shows the change from the mean shape (gray lines and open circles) to a configuration with an arbitrary PC score of +0.1 units of Procrustes distance (black line and solid circle). The percentages indicate the share of the total shape variation for which each PC accounts. This pattern of variation agrees well with the overall reported pattern of milk production, in which more milk is produce in spring, from March to June landmarks. The small differences show the uniform production not only throughout the years, but also among the months (observed in the more circular like of the shape).

right lobe, wider left lobe, and right shifted cleft. This pattern of variation agrees well with changes in the cows’ lactations, because majority of New Zealand’s cows calve in a concentrated pattern, in late-August to early-September, the overall shape of the national milk production could be similar to that of a standard lactation curve. The PC2 and PC3 were also associated with changes of the left lobe, summer to autumn production, and the PC4 featured variation in the right lobe.

Similarly, in the United States (Figure 6 (US)), the shape variations were described by the first four PCs, accounting 95.5% of total variance in the heart shape. PC1 variation was related to the cleft of the heart, compared with the mean shape (gray line and open circle).

5. Discussion

The application of GM as a new data analysis tool for dairy science was demonstrated using a dataset of monthly milk production at the industry level of two leading countries. The projected surface on a plane, representing milk yield over time by orbital graph, was the basis for the landmark based GM analysis. The monthly quantities of milk were the shapes’ landmarks and their respective coordinates were estimated using Pythagorean Theorem. Given the landmark (or a quantity of milk), the set of coordinates can be exactly calculated, if mathematical calculation are free of error, which overall is straightforward and convenient compared to that of GM that depends on special equipment and software to acquire image and locate landmarks. In addition, the counterclockwise rotation of the milk production data in orbital graph also provided a visual combination of the biological process of milk production coupled to its environment (seasons), which may provide understanding and guidance to appropriately analyze the interaction of these two processes. Rather than following the calendar year, the starting month can be adjusted to that of the enterprise since most, if not all, agricultural activities obey a seasonal pattern; the New Zealand milk production, for instance, follows a seasonal operation that begins in June prior to the months of grass growing season and ends in May of the followed year. The application of orbital graph to represent time series data is new, and its extent to lactation records has to be developed, especially in the number of axis representing lactations of different lengths.

GPA superimposed the landmark configurations, which influences the results in testing and visualizing shape variations. Different superimposition methods have different pros and cons, and those of GPA are well discussed by studies applying GM See [21] [22] [25] . The two most relevant to the current study are, first, in the rescaling process of GPA superimposition, the centroid size is not an immediately intuitive size measure, because the centroid size calculated for any subject will change if different landmarks configurations are used to summarize its shape [21] ; centroid size, however, is the most common measure of size used in GM studies [38] . Second, the variation of the shape may be related to only some of the landmarks, but after supper imposition the variation in shape is spread to all the landmarks configurations.

The individual variation of the production pattern was non-significant. This observation is reasonable about the genetic, management, and environment shared by the cows’ population throughout the years. The small differences may also provide evidence of the standardized production of milk by these two countries. Another factor to bear in mind is the type of data, national time series, which may also play a roll, since it is subject to methodological procedure in collection as well as the unequal number of days comprising each month. Differences between the production patterns were larger than within them, which may have been expected from Table 1, since we are comparing industries of distinct geographical region, farming system, and herd structure. It is important to mention that however similar or different these shapes are, one is not better than the other since they are the result of specific factors, above mentioned, so they should be considered as the best expression. A practical application of identifying the cardioid and heart shapes and their non-significant variation may be that they are validated to describe shape when modeling milk production as they significantly represent the pattern of milk production under given settings. As for that of the significant variation may be to discriminate groups, which raises caution as milk production performs different under different environmental and managerial conditions.

About 58.8% of shape variation in New Zealand was associated with the PC1 (Figure 6 (NZ)), which was interpreted to variation in the pattern of cows’ lactations. The narrowing part corresponds to the beginning of the aggregated lactations and the widening to its end. This pattern of variation may be explained by improvement on the genetic merit of cows, increase of day-in-milk, and cow productivity. Among the four groups of breeds presented in Table 1, Holstein-Friesian/Jersey crossbred and Holstein-Friesian groups comprised almost 80% of total herd, which growth compared to the other have both maintained steady in the last two years after a slightly increase of the former breed and a slightly decrease of the latter one; these two groups characterize by presenting longer lactations and higher productivity in the herd. Indeed, average days-in-milk increased by 7 days and daily milk per cow 1 unit in the last four years, compared with 2007-2010 [36] . The pattern of variation of PC1 (58.8% of total variation) may overall suggest an aggregated curve of a slightly lower peak (from calving to maximum) and a higher persistency (from maximum to end of lactation) for the 9 years under study. On the other hand, about 56.2% of shape variation in the United States was associated with PC1, which was interpreted to variation in the pattern of cows’ milk production. The sudden rise (PC1) and fall (PC2 and 3) of production at the cleft of the heart shape, correspond to the transition of the winter-spring production, which may be explained by the greater proportion of cows that are dried off during winter (which calved during the last spring) and that are freshen again in spring. Overall, the cleft of the cardioid and heart shapes marked the beginning (winter) and end (winter) of production in a year.

The application of GM and the results, however, should be seen as preliminary and exploratory. First, because it has to be considered whether the invariance procedure (translation, scaling, and rotation) done with Procrustes methods applied to the kind of data here used. For instance, the need for rotational invariance since the data might be already in a fixed orientation. The shape comparison in this study can benefit from using other GM methods such as Fourier analysis, which treat most of the landmarks as semi-landmarks. Regarding the results, because the statistical significance marked by the permutation test of DFA, though applicable to sample sizes very small [25] , may be reduced by small sample sizes as well as increase inaccuracies in estimating group means and variances. Second, because the number of observations (n = 9) is preferable to be greater than the dimensionality of the shape (number of landmark coordinates = 24), which is important to obtain a reliable estimate of the variance-covariance structure in the data [21] . Hence, a larger data set may confirm the significant difference between the patterns from New Zealand and United States as well as the non-significant difference of individual ones, since there must be some differences over the years.

6. Conclusions

We have developed the use of landmark-based geometric morphometric method to measure the form (magnitude and shape) projected on a plane by the graphical representation of milk production data. We demonstrated the new methodology in a small case study of milk production, monthly time series data, from New Zealand and United States dairy industries for their distinct farming system and geographical location. The monthly time series data were appropriate for implementing the methodology, which in turns demonstrated their use as morphological data. The closed curve depicted by orbital graph and regarded as the biological form of milk production was the basic for applying the landmark-based GM analysis.

The analysis revealed production patterns of cardioid shape in New Zealand and heart shape in United States. The data indicated that individual and group variations in size and shape of these patterns were significant, except for shape within country. Overall, GM method seems to be effective to quantify variation of shape and the magnitude of the milk production patterns. This may complement the analysis of milk prediction and reveal profound information of the biological process represented through the shape by controlling the different factors influencing it. The new methodological framework can be also applied to analyze the production pattern of individual and aggregated lactations and farm and region. The assumptions in applying the methodology and environmental, managerial, and physiological factors affecting each type of data may differ.

Acknowledgements

We profoundly thank to Edgardo Reyes Calderon and Wen-I Chan for their comments during the time we were conceiving this study and for their carefully review and suggestions on several version of this manuscript. We also thank to Meidiana Purnamasari and Thuline Phasha Zwane for their comments and suggestions on this manuscript as well as to Omkar Byadgi for his advice in preparing our response to the reviewers, to whom gratitude is also given for helping to craft the study.

References

[1] Beever, D.E., et al. (1991) A Review of Empirical and Mechanistic Model of Lactational Performance by the Dairy Cow. Livestock Production Science, 29, 115-130.

http://dx.doi.org/10.1016/0301-6226(91)90061-T

[2] Dongre, V.B., et al. (2011) A Brift Review on Lactation Curve Models for Predicting Milk Yield and Different Factors Affecting lactation Curve in Dairy Cattle. International Journal of Agriculture: Research and Review, 1, 6-15.

[3] Kawata, Y. (2011) Lactation Curves of Dairy Animals: An Interim Literature Review. Research Bulletin of Obihiro University, 32, 71-91.

[4] Macciotta, N.P.P., et al. (2011) The Mathematical Description of Lactation Curves in Dairy Cattle. Italian Journal of Animal Science, 10, 213-223.

http://dx.doi.org/10.4081/ijas.2011.e51

[5] Leclerc, H., et al. (2008) Environmental Effects on Lactation Curves Included in a Test Day Model Genetic Evaluation. Animal, 2, 344-353.

http://dx.doi.org/10.1017/S175173110700119X

[6] Rekik, B. and Ben Gara, A. (2004) Factors Affecting the Occurrence of Atypical Lactations for Holstein-Friesian Cows. Livestock Production Science, 87, 245-250.

http://dx.doi.org/10.1016/j.livprodsci.2003.09.023

[7] Tekerli, M., et al. (2000) Factors Affecting the Shape of Lactation Curves of Holstein Cows from the Balikesir Province of Turkey. Journal of Dairy Science, 83, 1381-1386.

http://dx.doi.org/10.3168/jds.S0022-0302(00)75006-5

[8] Wilmink, J.B.M. (1987) Adjustment of Test-Day Milk, Fat and Protein Yield for Age, Season and Stage of Lactation. Livestock Production Science, 16, 335-348.

http://dx.doi.org/10.1016/0301-6226(87)90003-0

[9] Silvestre, A.M., Petim-Batista, F. and Colaco, J. (2006) The Accuracy of Seven Mathematical Functions in Modeling Dairy Cattle Lactation Curves Based on Test-Day Records from Varying Sample Schemes. Journal of Dairy Science, 89, 1813-1821.

http://dx.doi.org/10.3168/jds.S0022-0302(06)72250-0

[10] Macciotta, N.P.P., et al. (2006) Factors Affecting Individual Lactation Curve Shape in Italian River Buffaloes. Livestock Science, 104, 33-37.

http://dx.doi.org/10.1016/j.livsci.2006.03.001

[11] Dematawewa, C.M.B. and Berger, P.J. (1998) Genertic and Phenotypic Parameters for 305-Day Yield, Fertility, and Survival in Holsteins. Journal of Dairy Science, 81, 2700-2709.

http://dx.doi.org/10.3168/jds.S0022-0302(98)75827-8

[12] Macciotta, N.P.P., Vicario, D. and Cappio-Borlino, A. (2005) Detection of Different Shapes of Lactation Curve for Milk Yield in Dairy Cattle by Empirical Mathematical Models. Journal of Dairy Science, 88, 1178-1191. http://dx.doi.org/10.3168/jds.S0022-0302(05)72784-3

[13] Olori, V.E., et al. (1999) Fit of Standard Models of the Lactation Curve to Weekly Records of Milk Production of Cows in a Single Herd. Livestock Production Science, 58, 55-63.

http://dx.doi.org/10.1016/S0301-6226(98)00194-8

[14] Cole, J.B., Null, D.J. and VanRaden, P.M. (2009) Best Prediction of Yields for Long Lactations. Journal of Dairy Science, 92, 1796-1810.

http://dx.doi.org/10.3168/jds.2007-0976

[15] Dematawewa, C.M.B., Pearson, R.E. and VanRaden, P.M. (2007) Modeling Extended Lactation of Holsteins. Journal of Dairy Science, 90, 3924-3936.

http://dx.doi.org/10.3168/jds.2006-790

[16] Grossman, M. and Koops, W.J. (2003) Modeling Extended Lactation Curves of Dairy Cattle: A Biological Basis for Multiphasic Approach. Journal of Dairy Science, 86, 988-998.

http://dx.doi.org/10.3168/jds.S0022-0302(03)73682-0

[17] Pollott, G.E. (2011) Short Communication: Do Holstein Lactations of Varied Lengths Have Different Characteristics? Journal of Dairy Science, 94, 6173-6180.

http://dx.doi.org/10.3168/jds.2011-4467

[18] Ehrlich, J.L. (2011) Quantifying Shape of Lactation Curves, and Benchmark Curves for Common Dairy Breeds and Parities. Bovine Practitioner, 45, 88-96.

[19] Ehrlich, J.L. (2013) Quantifying Inter-Group Variability in Lactation Curve Shape and Magnitude with the Milkbot Lactation Model. PeerJ, 1, e54.

http://dx.doi.org/10.7717/peerj.54

[20] Capuco, A.V., Wood, D.L., Baldwin, R., Mcleod, K. and Paape, M.J. (2001) Mammary Cell Number, Proliferation, and Apoptosis during a Bovine Lactation: Relation Milk Production and Effect of bST. Journal of Dairy Science, 84, 2177-2187.

http://dx.doi.org/10.3168/jds.S0022-0302(01)74664-4

[21] Webster, M. and Sheets, A.D. (2010) A Practical Introduction to Landmark-Based Geometric Morphometrics. Quantitative Methods in Paleobioogy, 16, 163-188.

[22] Adams, D.C., Rohlf, F.J. and Slice, D.E. (2004) Geometric Morphometrics: Ten Years of Progress Following the “Revolution”. Italian Journal of Zoology, 71, 5-16.

http://dx.doi.org/10.1080/11250000409356545

[23] Bookstein, F.L. (1991) Morphometric Tools for Landmark Data Geometric and Biology. Cambrige University, Cambrige.

[24] Cooke, S.B. and Terhune, C.E. (2015) Form, Function, and Geometric Morphometrics. The Anatomical Record, 298, 5-28.

http://dx.doi.org/10.1002/ar.23065

[25] Viscosi, V. and Cardini, A. (2011) Leaf Morphology, Taxonomy and Geometric Morphometrics: A Simple Protocol for Begginners. PLoS ONE, 6, e25630.

http://dx.doi.org/10.1371/journal.pone.0025630

[26] Kirkpatrick, A. and Lofsvold, D. (1989) The Evolution of Growth Trajectories and Other Quantitative Characters. Genome, 31, 778-783.

http://dx.doi.org/10.1139/g89-137

[27] Kirkpatrick, M. (1997) Genetic Improvement of Livestock Growth Using Infinite-Dimentional Analysis. Animal Biotechnology, 8, 55-61.

http://dx.doi.org/10.1080/10495399709525867

[28] Kirkpatrick, M., Hill, W.G. and Thomson, R. (1994) Estimating the Covariance Structure of Traits during Growth and Aging, Illustrated with Lactation in Dairy Cattle. Genetical Research, 64, 57-69. http://dx.doi.org/10.1017/S0016672300032559

[29] Kirkpatrick, M., Lofsvold, D. and Bulmer, M. (1990) Analysis of the Inheritance, Selection and Evolution of Growth Trajectores. Genetics, 124, 979-993.

[30] Murphy, M.D., O’Mahony, M.J., Shalloo, L., French, P. and Upton, J. (2014) Comparison of Modelling Technique for Milk-Production Forecasting. Journal of Dairy Science, 97, 3352- 3363. http://dx.doi.org/10.3168/jds.2013-7451

[31] Dairy Companies Association of New Zealand (DCANZ) (2015) New Zealand Milk Production 2007-2014. Dairy Companies Association of New Zealand.

[32] United Stated Department of Agriculture (USDA) (2015) Milk Production from 2007-2014. National Agricultural Statistics Service.

[33] Garcia, S.C. and Holmes, C.W. (1999) Effects of Time of Calving on the Productivity of Pastures-Based Dairy Systems: A Review. New Zealand Journal of Agricultural Research, 42, 347-362. http://dx.doi.org/10.1080/00288233.1999.9513384

[34] Kaiser, H.M., Oltenacu, P.A. and Smith, T.R. (1988) The Effect of Alternative Seasonal Price Differentials on Milk Production in New York. Northeastern Journal of Agricultural and Resource Economics, 17, 46-55.

[35] Sun, C.-H., Kaiser, H.M. and Forker, O.D. (1995) Analysis of Seasonal Milk Price Incentive Plan. Review of Agricultural Economics, 17, 383-393.

http://dx.doi.org/10.2307/1349581

[36] DairyNZ (2015) New Zealand Dairy Statistics.

http://www.dairynz.co.nz/TagListing/Index?tag=dairy%20industry%20statistics

[37] Bookstein, F.L. (Ed.) (1978) The Measurement of Biological Shapes and Shape Changes. Springer-Verlag, Berlin.

[38] Mitteroecker, P., Gunz, P., Windhager, S. and Schaefer, K. (2013) A Brift Review of Shape, Form, and Allometry in Geometric Morphometrics, with Application to Human Facil Morphology. Hystrix: The Italian Journal of Mammology, 24, 59-66.

[39] Klingenberg, C.P. (2011) MorphoJ: An Integrated Software Package for Geometric Morphometrics. Molecular Ecology Resources, 11, 353-357.

http://dx.doi.org/10.1111/j.1755-0998.2010.02924.x

[40] O’Higgins, P. and Jones, N. (1998) Facial Growth in Cercocebus torquatus: An Application of Three Dimensional Geometric Morphometric Techniues to the Study of Morphological Variation. Journal of Anatomy, 193, 251-272.

http://dx.doi.org/10.1046/j.1469-7580.1998.19320251.x

[41] Hammer, O. (2015) PAST.

http://folk.uio.no/ohammer/past

[42] Adams, D.C., Otarola-Castillo, E. and Sherratt, E. (2015) Geomorph: Software for Geometric Morphometric Analyses of 2D/3D Landmarks Data.

http://cran.r-project.org/web/packages/geomorph/

[43] Zelditch, M.L., Swiderski, D.L. and Sheets, H.D. (Eds.) (2004) Geometric Morphometrics of Biologists: A Primer. 2nd Edition, Elsevier Academic Press, San Diego, 443.

[44] Goodall, C. (1991) Procrustest Methods in the Statistical Analysis of Shape. Journal of the Royal Statistical Society, 53, 285-339.

[45] Zaiontz, C. (2015) The Data Analysis for This Paper Was Generated Using the Real Statistic Resource Pack Software.

www.real-statistics.com

[46] Klingenberg, C.P., Barluenga, M. and Meyer, A. (2002) Shape Analysis of Symmetric Structure: Quantifying Variation among Individuals and Asymmetry. Evolution, 56, 1909-1920. http://dx.doi.org/10.1111/j.0014-3820.2002.tb00117.x