Quantifying (non)parallelism of gut microbial community change using multivariate vector analysis

Abstract Parallel evolution of phenotypic traits is regarded as strong evidence for natural selection and has been studied extensively in a variety of taxa. However, we have limited knowledge of whether parallel evolution of host organisms is accompanied by parallel changes of their associated microbial communities (i.e., microbiotas), which are crucial for their hosts' ecology and evolution. Determining the extent of microbiota parallelism in nature can improve our ability to identify the factors that are associated with (putatively adaptive) shifts in microbial communities. While it has been emphasized that (non)parallel evolution is better considered as a quantitative continuum rather than a binary phenomenon, quantitative approaches have rarely been used to study microbiota parallelism. We advocate using multivariate vector analysis (i.e., phenotypic change vector analysis) to quantify direction and magnitude of microbiota changes and discuss the applicability of this approach for studying parallelism, and we compiled an R package for multivariate vector analysis of microbial communities (‘multivarvector’). We exemplify its use by reanalyzing gut microbiota data from multiple fish species that exhibit parallel shifts in trophic ecology. We found that multivariate vector analysis results were largely consistent with other statistical methods, parallelism estimates were not affected by the taxonomic level at which the microbiota is studied, and parallelism might be stronger for gut microbiota function compared to taxonomic composition. This approach provides an analytical framework for quantitative comparisons across host lineages, thereby providing the potential to advance our capacity to predict microbiota changes. Hence, we emphasize that the development and application of quantitative measures, such as multivariate vector analysis, should be further explored in microbiota research in order to better understand the role of microbiota dynamics during their hosts' adaptive evolution, particularly in settings of parallel evolution.


| INTRODUC TI ON
Parallel evolution, the repeated evolution of similar traits in independent lineages in response to similar selective pressures, is a widespread phenomenon and provides strong evidence for natural selection (Colosimo et al., 2005;Elmer et al., 2010;Losos et al., 1998;Rosenblum et al., 2017;Steiner et al., 2009). Yet, the extent of parallelism varies considerably across levels of biological organization (e.g., genotype vs phenotype) and across taxa (Bolnick et al., 2018). Substantial variation in parallelism can even be found among closely related populations adapting to seemingly similar habitats (Stuart et al., 2017). Traditionally, (non)parallel evolution has been regarded, and classified, as a binary phenomenon (evolution is parallel or not).
However, it has recently been argued that by considering parallel evolution as a quantitative continuum, we will be better able to identify and understand the genetic and ecological factors that affect the extent of parallelism (Bolnick et al., 2018).
The study of parallelism has recently been extended to hostassociated microbial communities and in particular the gut microbiota, the microbial community inhabiting a host's gut . To investigate microbiota parallelism, it can be useful to adopt both theoretical and methodological approaches developed for studying parallel evolution (e.g., Rennison et al., 2019). Gut microbial communities are highly diverse (Brooks et al., 2016;Human Microbiome Project Consortium, 2012;Youngblut et al., 2019) and affect host physiology in many ways (e.g., nutrient metabolism; Turnbaugh et al., 2006). The gut microbiota of an increasing number of host species is being characterized, and we are obtaining a more comprehensive picture of the extensive diversity of host-associated microbes (Song et al., 2020;Tarnecki et al., 2017). The gut microbiota is shaped by host genetics and ecological factors (Benson et al., 2010;Goodrich et al., 2014;Li et al., 2017;Spor et al., 2011;Sullam et al., 2012), and can impact the ecology and evolution of their hosts Zepeda Mendoza et al., 2018).
Study systems in which closely related populations or species have independently adapted to similar ecological niches are particularly useful for studying the evolutionary ecology of host-associated microbial communities (e.g., Härer et al., 2020). In these systems, one can ask whether phenotypic or ecological changes that have occurred repeatedly in multiple host populations (i.e., parallel evolution), are associated with parallel changes in microbial communities (i.e., microbiota parallelism). We would like to emphasize that microbiota parallelism solely describes repeatability in the direction and magnitude of change of microbial communities, but not necessarily their parallel evolution. There is now growing interest in determining whether parallel adaptation of hosts is associated with parallel microbiota changes, as parallelism among independent gut microbial communities suggests changes could be predictable and adaptive (Delsuc et al., 2014;Härer et al., 2020;Song et al., 2020).
Integration of microbiota data from a range of host populations that have repeatedly and independently adapted to similar ecological niches (e.g., Song et al., 2020) provides a powerful opportunity to investigate the ecological and evolutionary dynamics of host-microbe interactions.
Gut microbiota parallelism is predicted if hosts are adapting to similar trophic niches, since diet is known to be a major factor shaping gut microbial communities Smits et al., 2017;Turnbaugh et al., 2009). However, several ecological and genetic factors could further promote or hinder parallelism; these include similarity of host ecology, physiology and genetics, as well as, differential environmental exposure, and mode of microbial transmission (see Discussion for more details). This begs the question: Does the gut microbiota change in a predictable manner during their hosts' parallel adaptation to similar trophic niches, and if so, what factors affect the likelihood of observing parallelism? To address questions of gut microbiota predictability and parallelism, we need hypothesis-driven tests leveraged in systems with well-characterized host ecology and repeated patterns of niche shifts. It is also imperative to employ quantitative statistical metrics. Here we suggest that multivariate vector analysis, a quantitative method, can be used to estimate the degree of microbiota parallelism, which might allow identifying the key ecological and evolutionary processes shaping variation in microbial communities. This method, originally termed 'phenotypic change vector analysis' was developed by Collyer and Adams for studying magnitude and direction of multivariate phenotypic change (Adams & Collyer, 2009;Collyer & Adams, 2007), and it has previously been applied to study variation in phenotypic (Stuart et al., 2017) and gut microbiota (Rennison et al., 2019) parallelism in threespine stickleback fish. In this method, vectors connecting the multivariate means (centroids) of microbial communities are estimated for pairs of host populations.
The resulting vectors are then compared in a pairwise fashion for all host populations (Figure 1; see Section 2 for more details). This approach provides information not only on the direction (angle between vectors), but also the magnitude (vector length) of microbiota divergence. However, it is the angle that quantifies parallelism, the smaller the angle between two population pairs, the more parallel the pattern of divergence ( Figure 1) (Bolnick et al., 2018;Stuart et al., 2017). Crucially, when integrated with additional ecological or genetic data, this quantitative approach allows direct tests of factors that affect direction and magnitude of microbiota changes.
Parallel divergence in trophic ecology is well-documented in several teleost fishes; e.g., threespine stickleback (Bell & Foster, 1994;Taylor & Mcphail, 1999), African and Neotropical cichlids (Elmer et al., 2014;Muschick et al., 2012), lake whitefish (Bernatchez et al., 1999) and Trinidadian guppies (Reznick et al., 1996) (Table A1). To exemplify the utility of multivariate vector analysis for quantifying parallelism in compositional and functional changes of gut microbial communities, we reanalyzed published 16 S rRNA gene sequencing data sets from these model systems. We discuss how estimates of magnitude of divergence and parallelism can be used to identify factors that affect microbial communities associated with many host lineages and give recommendations on the use of this approach. We further acknowledge current limitations of using multivariate vector analysis for studying microbiota parallelism and emphasize the need for further development of this approach in microbiota research.
When applied to a wide range of host organisms, we argue that multivariate vector analysis has the potential to give a unique insight into the microbiota dynamics associated with their hosts' adaptation to different ecological niches.

| Data acquisition
We obtained 16 S rRNA gene sequencing data from six published studies of parallel evolution in teleost fishes: herbivorous and carnivorous African cichlids (Baldo et al., 2017), benthic and limnetic Neotropical Midas cichlids , benthic and limnetic threespine stickleback from British Columbia (Rennison et al., 2019), freshwater and estuarine threespine stickleback from Oregon (Steury et al., 2019), benthic and limnetic lake whitefish (Sevellec et al., 2018), and low-predation and highpredation Trinidadian guppies (Sullam et al., 2015). Sample sizes for each population are indicated in Table A1 and information on sequencing platform, amplicon region and NCBI archiving are listed in Table A2. For each data set, all samples were included in a single sequencing run. To improve readability, we will refer to different host lineages as populations, whether they are populations of the same species or distinct species. We tested for gut microbiota parallelism (i) among population pairs, and (ii) between an outgroup and several focal populations (Figure 1). Outgroups were selected based on phylogenetic information (see Table A1). Marine threespine stickleback colonized freshwater environments around 10,000-12,000 years ago (Bell & Foster, 1994). Hence, a marine population was selected as the outgroup since it represents the ancestral state. For Midas cichlids, the species A. citrinellus from Lake Nicaragua represents the outgroup to all crater lake species investigated, since crater lakes were colonized from the two great lakes of Nicaragua (L. Managua and L. Nicaragua) within the last 5000 years (Kautt et al., 2020). In African cichlids, we selected a species from Barombi Mbo as the outgroup to several species from Lake Tanganyika based on a recent phylogeny by Irisarri et al. (2018). In this study system, the focal carnivorous populations F I G U R E 1 Illustration of vector analysis for determination of gut microbiota parallelism. Vectors connect the population means (centroids) between population pairs (a) or between an outgroup and several focal populations (b). The phylogenetic trees in (a, b) represent schematics to provide general information on the phylogenetic relationships of the studied populations. For the population pair comparisons, the two populations used to calculate a vector (γ) were most closely related to each other providing repeated and independent cases of ecological divergence. The populations or species included for each analysis are listed in Table A1. Angles between vectors provide a quantitative measure of parallelism (c) and range from anti-parallel, to orthogonal to parallel (adopted from Bolnick et al., 2018). Angles between multivariate vectors were measured based on PCoA scores and represent a quantitative measure of gut microbiota parallelism. Vector lengths (L) provide information on the magnitude of gut microbiota changes. Centroids are shown as bold symbols (squares and circles) whereas individual data points are shown as faint symbols to illustrate differences in data distribution. Note that the direction of vectors is important and should be consistent across comparisons to obtain biologically meaningful results (e.g., always from ecotype A to ecotype B within a study system).
were chosen to avoid phylogenetic clustering; the closest relatives of all focal populations differ in trophic ecology (e.g., herbivores or omnivores) (Baldo et al., 2017). Further information on the populations used can be found in Table A1. Sequence data were downloaded from the NCBI Sequence Read Archive (SRA); information on sequencing platforms, sample sizes and accession numbers are provided in Tables A1 and A2. Data was converted from SRA to FASTQ format using the fastq-dump function of the SRA Toolkit v2.9.6-1 (https://ncbi.github.io/sra-tools/).

| Gut microbiota analysis
Forward reads had higher sequence quality than reverse reads, and read lengths varied across studies due to differences in sequencing technology, which led to non-overlap of reads for some studies.
Hence, we only used forward reads to achieve higher consistency in analysis across the different data sets. Forward reads were imported into the open-source bioinformatics pipeline QIIME2 (Bolyen et al., 2019). Sequence quality control was done with the plugin DADA2 (Callahan et al., 2016) and a phylogenetic tree was produced with FastTree 2.1.3 (Price et al., 2010); read numbers before and after DADA2 filtering are provided in the Dryad database (see Data Accessibility Statement). Taxonomy was assigned against the 16 S rRNA gene Silva database version 132 (Quast et al., 2013) using the feature-classifier classify-sklearn plug-in in QIIME2 (Pedregosa et al., 2011). Taxonomic assignment was not done for Trinidadian guppies (Sullam et al., 2015) as this study used a different region of the 16 S rRNA gene (V1-V3). Rarefaction depths and sequence lengths varied across data sets (Table A1), ASV sequences can be obtained from the Dryad database (see Data Accessibility Statement).

Encyclopedia of Genes and Genomes (KEGG) orthologs and Enzyme
Commission numbers were predicted with the PICRUSt2 plugin in QIIME2 (Douglas et al., 2020;Kanehisa et al., 2012) with a maximum nearest-sequenced taxon index (NSTI) cutoff of 2. Across the study systems, more than 90% of ASVs (92.7-99.5%) were below this cutoff, except for the study on whitefish where the proportion was slightly lower (85.3%; Table A3). Mean and median NSTI scores ranged from 0.45-2.689 and 0.072-0.366, respectively. Based on distance matrices for all these different metrics, principal coordinate analyses (PCoA) were performed and PCoA scores were used as input for multivariate vector analyses.

| Multivariate vector analysis
We quantified compositional and functional gut microbiota parallelism using multivariate vector analysis and compiled our code into an R package that can be obtained from github (https://github.com/ andre as-haere r/multi varve ctor). We largely followed the methodology reported in Rennison et al. (2019), which to date is the only study that has used this approach for studying gut microbiota parallelism. We found different degrees of parallelism for gut microbiota function than reported in the original study. This is likely due to differences in data processing and analysis pipelines, emphasizing the need to standardize data analysis when making inferences across studies. Multivariate vectors were calculated by connecting the population means (centroids) of PCoA scores, either between population pairs or between an outgroup and focal populations (as depicted in Figure 1). The dimensionality of the data sets used to estimate angles for gut microbiota composition and function is listed in Tables A4 and A5, respectively. Angles were measured between these vectors and were calculated for all possible pairwise comparisons in each data set (e.g., between all three benthic-limnetic population pairs in threespine stickleback from British Columbia, Canada).
The direction of vectors was held consistent, e.g., from benthic to limnetic across all comparisons within a study system, representing a repeated measure of evolutionary divergence between populations.
Yet, we would like to mention that this does not necessarily reflect the direction of evolutionary change in all of our study systems (i.e., ancestral to derived). Angles for gut microbiota composition  (Tables A4 and A5). Hence, in the main text we only present Bray-Curtis dissimilarity for gut microbiota composition and MetaCyc pathway abundances based on Bray-Curtis dissimilarity for the inferred functional metagenome.
The angles provide us with a quantitative measure of parallelism within and across study systems. Smaller angles (below 90°) indicate parallelism, angles around 90° indicate orthogonal change and larger angles (above 90°) indicate anti-parallelism ( Figure 1c). A more detailed discussion on different interpretations of the distribution of angles, but also on the limitations of this method can be found in Bolnick et al. (2018) and in Watanabe (2022). To statistically test for gut microbiota parallelism, previous studies proposed to either regard changes as parallel when angles do not deviate from 0° (Bolnick et al., 2018) or when angles are significantly smaller than 90° (Rennison et al., 2019). To give us an initial simplistic indication of parallelism patterns, we used one-sample t-tests with an angle of 90° as the null expectation for non-parallelism since (almost) all data were normally distributed based on Shapiro-Wilk tests (Shapiro & Wilk, 1965). The only exception was gut microbiota composition of threespine stickleback population pairs from Rennison et al. (2019), for which we performed a one-sample Wilcoxon signed-rank test (Wilcoxon, 1945). However, due to the non-independence of pairwise angles (Watanabe, 2022), we also quantified parallelism by calculating distributions of random angles in multidimensional space (which is centered at 90°) and using Monte Carlo simulations (with 10 5 iterations) to test for significant parallelism, or by performing a Rayleigh test which is used to examine the unimodal concentration of directional vectors (Mardia et al., 1979;Watanabe, 2022). We compare the results obtained using these different methods and discuss their interpretation.
To quantify the magnitude of gut microbiota changes, we calculated means of vector lengths (meanL) for each population pair.
Correlation analyses were based on non-parametric Spearman rank correlation coefficients, as not all data were normally distributed. For the population pair comparisons, we also tested whether the taxonomic level at which microbial communities are studied affect the magnitude and direction of gut microbiota change based on the multivariate vector analysis. Taxonomic assignment was done in QIIME2 (more information provided above), biom tables were created for each study system at different taxonomic levels of the bacterial communities (phylum, class, order, family, genus, species). From each biom table, we produced a Bray-Curtis distance matrix and calculated PCoA scores. Estimates of angles and vector lengths were calculated based on these PCoA scores, similar to the other analyses mentioned above.

| Gut microbiota parallelism across population pairs
First, we tested for parallelism across population pairs where vectors connect two closely related populations, representing cases of repeated divergence. Levels of gut microbiota parallelism varied considerably within and across teleost fish model systems for parallel evolution ( Figure 2). Within study systems, we detected a wide range of angles (e.g., 38.5-88.8° for gut microbiota composition in African cichlids; Figure 2a), and parallelism estimates were also highly variable across study systems. When comparing mean angles against a null expectation of 90°, statistically significant parallelism was only found for gut microbial composition (mean: 68. 9°; one-sample t-test: p < .001, t = −6.566) and function (mean: 38.9°, p < .001, t = −12.513) among herbivorous and carnivorous African cichlids (Figure 2a, b). We detected suggestive evidence for parallelism of gut microbial composition in benthic and limnetic threespine stickleback from British Columbia (mean: 81°, one sample Wilcoxon signed rank test: p = .087), but the sample size of this study was very small (n = 3). Similar results were obtained when using Monte Carlo simulations to compare mean angles against a multidimensional null distribution, we detected significant parallelism for African cichlids' gut microbiota composition and function (p < 1 × 10 −5 for both tests), and suggestive evidence in gut microbiota composition of benthic and limnetic threespine stickleback from British Columbia (p = .087) as well as in gut microbiota function of freshwater and estuarine threespine stickleback from Oregon (p = .052; Figure A3). These findings were further supported by Rayleigh tests, in which significant concentrations of angles were detected for gut microbiota composition (S = 113.01, p < 1.26 × 10 −9 ) and function (S = 160.47, p = 6.29 × 10 −21 ) only in African cichlids. When performing multivariate vector analysis based on different taxonomic levels of gut microbial communities (species to phylum), we found that estimates of F I G U R E 2 Significant gut microbiota parallelism among population pairs was detected in African cichlids for (a) gut microbiota composition and (b) function. In the outgroup comparisons, all three study systems showed significant gut microbiota parallelism for composition (c) and function (except for African cichlids; d). For all three study systems, outgroups inhabit distinct water bodies from the focal populations (Table A1). Numbers of comparisons are indicated next to the name of each study system; the populations used for each analysis as well as samples sizes for each population are stated in Table A1. Here, we show the results of testing mean angles against 90°, and no statistical tests were conducted for threespine stickleback population pairs from lakes and estuaries in Oregon and benthic and limnetic Midas cichlids from Nicaraguan crater lakes as we only had one comparison for each data set (a, b).
direction and magnitude were generally consistent, and independent of taxonomic resolution (Figure 3).

| Gut microbiota parallelism across focal populations compared to outgroup
When comparing an outgroup to focal populations adapted to similar ecological niches (Figure 1b), we found strong evidence for compositional and functional gut microbiota parallelism across all three study systems; almost all angles were smaller than 90° (Figure 2c, d).
By comparing mean angles against a null expectation of 90°, gut microbiota changes associated with carnivory in African cichlids from Lake Tanganyika (compared to an herbivorous outgroup from Lake Barombi Mbo) were significantly parallel for composition Angles for gut microbiota composition and function were strongly correlated across study systems (Spearman's ρ = 0.726, p < .001; Figure 4c). Angles did not differ between gut microbiota composition and function (p = .21, t = −1.279; Figure 4d); even when only considering comparisons with angles below 90° for both measures (p = .723, t = 0.358). There was no correlation between gut microbiota parallelism (angles) and mean magnitude of gut microbial change (meanL), for composition (ρ = −0.08, p = .691) or function (ρ = 0.073, p = .716; Figure A5c, d).

| DISCUSS ION
Multivariate vector analysis offers the opportunity to quantify direction and magnitude of microbiota changes, thereby allowing identification of factors that shape microbial communities across host populations. The purpose of our study is to exemplify and advocate the use of this approach in microbiota research, compare different statistical approaches to test for parallelism, discuss current limitations, and give recommendations on methodological aspects.
To this end, we applied multivariate vector analysis to investigate the gut microbiota of several classic teleost fish model systems that exemplify parallel ecological shifts. While this manuscript is focused on the study of parallelism, multivariate vector analysis can also be used to study the extent of convergence or divergence of microbiota changes, when incorporating information on the direction of evolutionary change in host lineages (Bolnick et al., 2018). It should be noted that multivariate vector analysis has only been used in one microbiota study (Rennison et al., 2019), and more conceptual and methodological research is needed to determine how characteristics of microbiota data sets affect parallelism estimates in order to comprehensively interpret the biological implications of such results.
We would also like to emphasize that the small set of study systems included here are meant to exemplify the use of this quantitative approach. The goal was not to identify or disentangle the effects of specific factors that shape gut microbiota parallelism, but we encourage future studies to apply quantitative analyses to larger data sets in order to investigate general patterns and processes using the methods presented here.

| Multivariate vector analysis to quantify microbiota parallelism
When studying microbiota parallelism, one admittedly simplified prediction is that parallel adaptation of hosts to a novel diet trans- There has been no conclusive evidence presented to suggest gut microbiota parallelism in Nicaraguan Midas cichlids , lake whitefish (Sevellec et al., 2018) or Trinidadian guppies (Sullam et al., 2015).
The variation in gut microbiota parallelism described in previous studies might be biologically real or could be due to differences in the methodological approaches used, which limits our ability to make inferences across systems. Different statistical methods were used to analyze the gut microbiota data of these species; four of the studies used permutational analysis of variance (PERMANOVA) to determine the effects of ecotype, the environment, and their interaction to infer parallelism. One study on African cichlids further inferred parallelism based on PCoA scores of the first axes (Baldo et al., 2017). A drawback of such approaches is that gut microbiota changes are scored as parallel or non-parallel in a binary manner, based on a given significance threshold. Such methods also do not provide estimates of the magnitude of change. Thus, we lack quantitative information on the extent and variation of (non)parallelism ( Figure 1c). Multivariate vector analysis allows quantification of variation across independent population pairs. When applied to large data sets, the method can facilitate identification of factors that underlie variation in parallelism. This is important as we know that a multitude of genetic, ecological, and physiological factors affect the composition of microbial communities in different host lineages, including teleost fishes (Amato et al., 2019;Baldo et al., 2017;Bletz et al., 2016;Bolnick, Snowberg, Hirsch, Lauber, Org, et al., 2014;Burns et al., 2017;Zhang et al., 2016). When applied to our six case studies, we found similar results as other statistical methods; there was only strong evidence of parallelism among population pairs of African cichlids, and weaker evidence for benthic-limnetic threespine stickleback (Figure 2a, b). The consistency of these results highlights the suitability of multivariate vector analysis for studying microbiota parallelism, with the added value of quantitative data that can be used in direct hypothesis testing.

| Methodological considerations
We tested whether methodological aspects, such as the choice of metric for determining differences in bacterial community composition (e.g., Bray-Curtis dissimilarity, weighted and unweighted UniFrac) or the taxonomic level at which microbial communities are F I G U R E 3 Estimates of direction (left column) and magnitude (right column) of gut microbiota change did not differ substantially with the taxonomic level used for multivariate vector analysis.
studied, affect estimates of direction and magnitude of gut microbiota change. We detected strong correlations among angle estimates generated using the different metrics ( Figure A1), and statistical significance of parallelism was generally consistent across different metrics for both taxonomic composition (Table A4) and inferred metagenome function (Table A5). This suggests that the parallelism estimates from multivariate vector analysis are not strongly affected by the choice of diversity metric. It should be noted that multivariate vector analysis is not restricted to analyzing data from principal coordinates analysis (PCoA); data from other ordination methods, such as non-metric multidimensional scaling (NMDS), can also be used as input for the analysis. NMDS data have been used in a parallelism analysis of the threespine stickleback gut microbiota, and the results are qualitatively similar to the PCoA analysis reported here (Rennison et al., 2019). Further, results were robust across taxonomic levels (Figure 3), which was a bit surprising as the composition of microbial communities might be more stochastic at lower taxonomic levels. Thus, predictability of microbiota change (i.e., parallelism) could have been expected to be stronger at higher taxonomic levels. However, our results suggest that researchers should be able to tailor the taxonomic level to the needs of their particular study.
Previous studies utilized different approaches for significance testing of angles from multivariate vector analysis (Bolnick et al., 2018;Rennison et al., 2019). One study suggested that to infer parallelism, methods should be based on the distribution of random angles, which depends on data dimensionality (Watanabe, 2022), rather than testing against a certain angle (e.g., 90°). Fortunately, in highly multidimensional space, such as that of microbiota data, random angles approximately follow a normal distribution centered around 90° (Watanabe, 2022). Hence, we argue that one-sample t-tests comparing the mean angles of our empirical data against the null expectation of 90° is a useful first test of parallelism. However, given that the 90° value is only a rough approximation, we implemented two additional statistical analyses to test for significant parallelism. We used Monte Carlo simulations to compare the estimated mean empirical angles against a multidimensional null distribution.
We further implemented Rayleigh tests which test for the unimodal concentration of directional vectors (Mardia et al., 1979). Across these three statistical methods, we obtained largely consistent results. Since each statistical approach tests for a different property (see Watanabe, 2022 for a more detailed discussion), the use of multiple approaches helps obtain a more comprehensive and robust picture of whether changes of the gut microbiota are consistent with a parallel, orthogonal, or anti-parallel pattern of change ( Figure 1c).
Multivariate vector analysis was developed for studying a range of phenotypic traits (Adams & Collyer, 2009), and has been used to study morphological and behavioral parallelism (Stuart et al., 2017).
Only one study has applied this approach for studying gut microbiota F I G U R E 4 Levels of gut microbiota (non)parallelism were strongly correlated for gut microbiota composition and function for population pairs (a) and the outgroup comparisons (c); the dashed line has a slope of 1. There was only suggestive evidence for a difference in angles between gut microbiota composition and function for population pairs (b) but not for the outgroup comparison (d). Information on the populations used for each analysis are stated in Table A1. †p < .1.
to obtain a comprehensive picture of microbiota parallelism in a phylogenetically and geographically broad range of host lineages.

| Progress and prospects in implementing multivariate vector analysis
Multivariate vector analysis quantifies changes in microbial communities and enables testing the effects of ecological or genetic factors on parallelism. Across six model systems of parallel evolution, we observed extensive variation in gut microbiota parallelism (Figure 2), and we sought to conduct preliminary analyses to explore some of the factors that might affect parallelism. For example, we found that changes in gut microbiota function might be more parallel than in taxonomic composition, which is in line with findings of individual studies (Figure 4) Rennison et al., 2019). This pattern could be explained by the functional redundancy within microbial communities . Variation in taxonomic composition might be produced by historical contingency, priority effects or microbial dispersal ability (Costello et al., 2012;Martinez et al., 2018). Yet, taxonomically distinct bacterial taxa can provide similar metabolic functions, potentially causing stronger signatures of parallelism in gut microbiota function when hosts adapt to similar ecological niches. At the same time, parallelism estimates for gut microbiota composition and function were strongly correlated ( Figure 4), indicating that similar factors shape the extent of parallelism. We suggest that work integrating parallelism estimates for both composition and function from a variety of disparate host taxa is needed to explore this pattern further. But, it is important to be cautious when interpreting results on inferred metagenome function using tools such as PICRUSt2, particularly in non-model organisms (Douglas et al., 2020); the reliability of metagenome prediction highly depends on the database of available bacterial genomes, which can vary across host organisms (Sun et al., 2020).
Shifts in gut microbial communities were found to be most parallel among population pairs when the overall magnitude of divergence (meanL) was greatest ( Figure A5); this raises the question of whether substantial divergence in the gut microbiota might indicate shared adaptive changes. Strong selection pressures may lead to more determinism in microbial community assembly and, consequently stronger microbiota parallelism. In contrast, when there is little microbiota divergence, differences might be produced mainly by stochastic processes, e.g., priority effects and drift (Martinez et al., 2018). These processes are unlikely to generate similar microbiota shifts, and theoretical work suggests that parallelism may only be seen when the selection landscape is highly parallel . For morphological traits and genetic divergence, it has been shown that the degree of environmental variation among evolutionary replicates directly predicts the magnitude of parallelism based on multivariate vector analysis (Stuart et al., 2017). If more similar selection pressures translate to more parallel host phenotypic and ecological change, this could lead to similar shifts in microbial communities. Accordingly, host adaptation to similar diets, when accompanied by changes in gut morphology, is expected to promote gut microbiota parallelism (Baldo et al., 2017;Muegge et al., 2011), whereas adaptation to different diets could lead to gut microbiota anti-parallelism. Yet, strong parallelism could also be explained by variation in environmental factors, as seen in our outgroup comparisons (Figure 2c, d). In these settings, focal populations and the outgroups live in strongly differentiated environments (Barluenga & Meyer, 2010;Ormond et al., 2011;Torres-Dowdall et al., 2017), suggesting that similar patterns of divergence (parallelism) might be primarily driven by abiotic (physicochemical properties) and biotic (microbial communities) differences between environments, in addition to or instead of host trophic ecology. This was further supported by findings in Midas cichlids and threespine stickleback, where angles did not appear to be smaller for comparisons that only included the outgroup and a certain ecotype (i.e., benthics or limnetics) than for comparisons that included the outgroup and both ecotypes ( Figure A7). Yet, one should be careful when drawing conclusions from these results since sample sizes were small. Future work that estimates the similarity of abiotic and biotic factors among host populations will be key to determining whether parallel selective regimes generally translate to parallelism in changes of microbial communities.
Variation in microbial communities can be strongly associated with host phylogeny and the extent of genetic divergence (Benson et al., 2010;Brooks et al., 2016;Goodrich et al., 2014;Li et al., 2017;Youngblut et al., 2019); thus, increasing genetic (and phylogenetic) distance among hosts (Smith et al., 2015) could affect the likelihood of observing microbiota parallelism. Stronger parallelism might be predicted for host lineages that split earlier, with sufficient time to diverge ecologically. Results from Neotropical and African cichlids hint at such an association; the very recently diverged crater lake Midas cichlids from Nicaragua (<5000 years ago) (Kautt et al., 2020) showed no evidence for gut microbiota parallelism. In contrast, there was strong evidence in much older African cichlid lineages, where divergence times are in the range of millions of years ( Figure 2) (Baldo et al., 2017). Threespine stickleback and whitefish colonized freshwater lakes after the last ice age (<12,000 years ago) and formed ecologically distinct species pairs adapted to different niches (Bernatchez et al., 1999;Matthews et al., 2010;Schluter & Mcphail, 1992). Yet, we only detected some evidence for parallelism in benthic-limnetic stickleback (Figure 2), suggesting that an association with host divergence time is not necessarily expected in general. Parallelism was also stronger in the outgroup comparisons (Figure 2), which could be driven by stronger genetic divergence of the outgroup compared to focal populations (but also by environmental differences, see previous paragraph). However, for two of the study systems (Midas cichlids and threespine stickleback) the outgroup and the focal populations split very recently (<12,000 years), whereas population pairs of African cichlids split much earlier. Hence, the stronger parallelism in the outgroup comparisons cannot be explained solely by host divergence time. Again, a diverse sampling of host taxa, potentially including intra-as well as interspecific host comparisons, will be key to robustly test whether divergence time is a key factor determining patterns of gut microbiota parallelism.
There are numerous possible uses of multivariate vector analysis to investigate changes in host-associated microbial communities as many genetic and ecological factors may determine parallelism.
For example, it would be particularly interesting to compare patterns among trophic specialists and generalists. Specialists might have a more stable gut microbiota compared to the gut microbiota of generalists, which is expected to fluctuate more over time (Baniel et al., 2021;Smits et al., 2017). Thus, this variation in gut microbiota plasticity (Kolodny & Schulenburg, 2020) could affect observed parallelism.
The mode of gut microbiota transmission could be another interesting factor to consider using multivariate vector analysis. The acquisition of gut microbiota from environmental sources (i.e., horizontally) vs transfer from mother to child (i.e., vertically) could affect the likelihood of observing parallelism, depending on whether these sources are shared among hosts (Mulder et al., 2009;Smith et al., 2015).
Further, multivariate vector analysis could also be leveraged in experimental studies to investigate parallelism associated with short-term microbiota shifts in response to environmental changes, for example, after a host organism is exposed to a novel diet or a certain pathogen.
It would be interesting to compare the extent of microbiota parallelism over short (ecological) and longer (evolutionary) timescales. To obtain a comprehensive understanding of the evolutionary ecology of the gut microbiota, future studies should make an effort toward combining quantitative measures with a diverse range of genetic and phenotypic host data, as well as environmental data.

| CON CLUS IONS
Our study exemplifies the use of multivariate vector analysis for studying microbiota dynamics, and discusses potential advantages compared to more commonly used statistical approaches. By using a common analytical framework, multivariate vector analysis allows quantification of the direction and magnitude of microbiota change. When applied across a broad range of host taxa, these estimates can be leveraged to examine general patterns of repeatability. Combining quantitative estimates with host-associated and environmental data offers the possibility to improve our knowledge of the eco-evolutionary processes that shape microbial community dynamics. Identification of these factors will improve our ability to predict (putatively adaptive) shifts in microbial communities. Hence, we encourage further adoption of quantitative measures for studying microbiota dynamics during adaptive evolution of their hosts, particularly in settings of parallel evolution.

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
All sequencing data included in this study has been obtained from data sets published in the NCBI Sequence Read Archive How to cite this article: Härer, A., & Rennison, D. J. (2022).
Quantifying (non)parallelism of gut microbial community change using multivariate vector analysis. Ecology and Evolution, 12, e9674. https://doi.org/10.1002/ece3.9674 A PPE N D I X A F I G U R E A 1 Parallelism estimates among different gut microbiota composition metrics were strongly correlated across comparisons of (a) population pairs and (b) outgroups vs focal populations. Angles calculated from the multivariate vector analysis are indicated on the x-and y-axes.

F I G U R E A 2
Parallelism estimates among different gut microbiota function metrics were strongly correlated across comparisons of (a) population pairs and (b) outgroups vs focal populations. Angles calculated from the multivariate vector analysis are indicated on the x-and y-axes.

F I G U R E A 3
Histograms showing the distribution of angles among population pairs for gut microbiota composition (left column) and function (right column). The vertical dashed line represents the mean angle for each study system, and the multidimensional distribution of random angles is illustrated by the red curves, the shape of a curve is determined by the number of PCoA axes included in each analysis.

F I G U R E A 4
Across population pair comparisons, angles for gut microbiota function were significantly smaller than for composition (p < .001, t = 5.5881) when only considering comparisons that showed evidence of gut microbiota parallelism (angles <90° for composition and function).

F I G U R E A 5
The strength of parallelism (angles) and the mean magnitude of gut microbiota change (meanL) were correlated for composition and function across population pair comparisons (a, b); the more parallel gut microbiota changes were (i.e., smaller angles), the higher the magnitude of change. No correlation was observed for outgroup vs focal population comparisons (c, d).

F I G U R E A 6
Histograms showing the distribution of angles between an outgroup and focal populations for gut microbiota composition (left column) and function (right column). The vertical dashed line represents the mean angle for each study system, and the multidimensional distribution of random angles is illustrated by the red curves, the shape of a curve is determined by the number of PCoA axes included in each analysis.

F I G U R E A 7
For the outgroup comparisons, we tested whether parallelism is stronger among focal populations of the same ecotype compared to focal populations of different ecotypes for gut microbiota composition and function in threespine stickleback from British Columbia (a, c) and Midas cichlids (b, d). The three different categories include angles of "outgroup-benthic vs outgroup-benthic", "outgroup-benthic vs outgroup-limnetic" or "outgroup-limnetic vs outgrouplimnetic" comparisons. We did not detect evidence that angles were smaller when only comparing among the same ecotype, suggesting that parallelism is mainly driven by the shared adaptation to a common environment among the focal populations, more so than by the adaptation to benthic and limnetic niches specifically.
TA B L E A 1 Identity of populations/species selected from each dataset that were included in our study, and information on their trophic ecology or habitat types as well as sample sizes. Note: Please note that one species of African cichlids (C. coloratus) is planktivorous; it was included in our analyses to provide another comparison since it showed gut microbiota variation from the carnivorous P. straeleni in the same direction as other herbivorous species along the major axis of differentiation (Baldo et al., 2017). Different subsets were used in the population pair and outgroup vs focal species analyses. For African cichlids, the herbivorous species from Barombi Mbo (S. steinbachi) represents an outgroup to the carnivorous species from Tanganyika based on a recent phylogeny (Irisarri et al., 2018). For threespine stickleback and Midas cichlids, outgroups are inferred ancestral populations that repeatedly colonized novel environments (Bell & Foster, 1994;Kautt et al., 2020). Sequencing read lengths and rarefaction depths used for all analyses are stated in bold for each study system and type of comparison.

TA B L E A 1 (Continued)
TA B L E A 2 Information on sequencing platforms, amplified region and data archiving for all datasets included in this study.

TA B L E A 3
Inferred metagenome function was predicted with the PICRUSt2 plugin in QIIME2 with a maximum nearest-sequenced taxon index (NSTI) cutoff of 2, and reads above this value were discarded for these analyses. Note: Statistical significance was determined at the .05 level and results were highly consistent across metrics, only for threespine stickleback from British Columbia and lake whitefish did we detect some suggestive evidence (p < .1) for parallelism for one of the three metrics, but not in the other two. We performed one-sample t-tests when data was normally distributed; for data with non-normal distribution one-sample Wilcoxon signed-rank tests were used (indicated by asterisks). The dimensionality of each data set is indicated in brackets.
p values between .05 and 0.1 are indicated in bold and italics, p values below .05 are indicated in bold.

TA B L E A 4
Comparison of mean angles and statistical tests for parallelism (angles <90°) for taxonomic composition of the gut microbiota across three different metrics: Bray-Curtis dissimilarity, weighted and unweighted UniFrac. Note: Statistical significance was determined at the .05 level and results were consistent across metrics for all study systems. We performed one-sample t-tests when data was normally distributed; for data with non-normal distribution one-sample Wilcoxon signed-rank tests were used (indicated by asterisks). The dimensionality of each data set is indicated in brackets.
p values between .05 and 0.1 are indicated in bold and italics, p values below .05 are indicated in bold.

TA B L E A 5
Comparison of mean angles and statistical tests for parallelism (angles <90°) for inferred metagenome function across three different metrics based on bray-Curtis dissimilarity: MetaCyc pathway abundances (PWAB), enzyme commission numbers (EC) and Kyoto encyclopedia of genes and genomes orthologs (KO).