Compositional Graphical Lasso Resolves the Impact of Parasitic Infection on Gut Microbial Interaction Networks in a Zebrafish Model

Microorganisms are ubiquitous in nature and manage critical ecosystem services, ranging from global nutrient cycling to human health. Microbes do not exist in a vacuum in nature, but instead are members of diverse communities known as microbiomes, wherein microbes may interact with other microbes of different species. These microbial interactions may be favorable or antagonistic and are crucial for the successful establishment and maintenance of a microbial community, and frequently result in important pathogenic or beneficial effect to the host or environment (Braga et al., 2016). A thorough understanding of how microbe interact with one another is critical to uncovering the underlying role microorganisms play in the host or environment. However, it is a highly challenging biological task as it is estimated that only 1% of bacteria are cultivatable. The inability to culture the majority of microbial species has motivated the use of culture-independent methods for microbiome studies in different environments (Faust and Raes, 2012; Tang, 2019).

Fortunately, recent innovations in in situ DNA sequencing provide opportunities to infer how microbes interact with one another or their environments. Modern studies on microbial interactions frequently rely on DNA sequencing techniques through the bioinformatic analysis of taxonomically diagnostic genetic markers (e.g., 16S rRNA) sequenced directly from a sample. The counts of the taxonomically diagnostic genetic markers can be used to represent the abundance of microbial species, e.g., Operational Taxonomic Units (OTUs) or phylotypes (e.g., genera), in a sample. Here, the frequency with which a taxon’s marker is observed in a sequence library represents its relative abundance in the community. When such abundance data are available from many communities, interactions among microbiota can be inferred through statistical correlation analysis (Faust and Raes, 2012). For example, if the relative abundances of two microbial taxa are statistically correlated, then it is inferred that they interact on some level. This approach has been used to document interaction networks in the healthy human microbiome (Faust and Raes, 2012), as well as free-living microbial communities (Freilich et al., 2010), and has been useful for generating hypotheses of host-microbiome interaction (Morgan et al., 2015).

Despite the great potential of microbiome interaction networks as a tool to advance microbiome research, the power of this approach has been limited by the availability of effective statistical methods and computationally efficient estimation techniques. Microbial abundance data possess a few important features that pose tremendous challenges to standard statistical tools. First, the data are represented as compositional counts of the 16S rRNA sequences because the total count of sequences per sample is predetermined by how deeply the sequencing is conducted, a concept named sequencing depth. The counts only carry information about the relative abundances of the taxa instead of their absolute abundances. Second, the sequencing depth is always finite and often varies considerably across samples in a microbiome dataset. Thus, the observed relative abundance of a taxon in a sample is only an estimator of its true relative abundance with the variance depending on sample-specific sequencing depth, causing the “heteroscedasticity” issue (McMurdie and Holmes, 2014). Third, the data are high-dimensional in nature. It is likely that the number of taxa is far more than the number of samples in any biological experiment.

When such abundance data are available, one common strategy to resolve the interactions among microbial taxa is to use correlation-type analyses. For example, after the sample correlations are calculated between the relative abundances of each pair of microbial taxa, a threshold is then applied such that an interaction is deemed present if the sample correlation exceeds the threshold. More recently developed methods have started to account for the compositional feature and aim to construct sparse networks for the absolute abundances instead of relative abundances, including SparCC (Friedman and Alm, 2012), CCLasso (Fang et al., 2015), and REBACCA (Ban et al., 2015), among others.

The above-mentioned network methods based on marginal correlations could lead to spurious correlations that are caused by confounding factors such as other species in the same community. To eliminate the detection of spurious correlations, interactions among taxa can be modeled through their conditional dependencies given the other taxa. The Gaussian graphical model is a classical approach to modeling the conditional dependency, in which the conditional dependency is determined by the nonzero entries of the inverse covariance matrix of a multivariate normal distribution used to model the data. Graphical lasso (Yuan and Lin, 2007; Banerjee et al., 2008; Friedman et al., 2008) and neighborhood selection (Meinshausen et al., 2006) are two commonly used methods to estimate sparse inverse covariance matrix under the high-dimensional Gaussian graphical model. However, when applied to microbiome abundance data, its multivariate normality assumption was violated by either the count or the compositional features of the data.

Several methods have been proposed to infer microbial conditional dependence networks based on Gaussian graphical models, such as SPIEC-EASI (Kurtz et al., 2015), gCoda (Fang et al., 2017), CD-trace (Yuan et al., 2019), and SPRING (Yoon et al., 2019). In order to transform the discrete counts to continuous variables and to remove the compositionality constraint, all these methods take the centered log-ratio transformation (Aitchison, 1986) on the observed counts as their first step. However, the centered log-ratio transformation suffers from an undefined inverse covariance matrix of the transformed data. To partially address this issue, these methods impose a sparsity assumption on the inverse covariance matrix and adding an $L_{1}$-norm penalty to their objective functions. Additionally, these methods ignore the heteroscedasticity issue as they simply treat the observed relative abundances as the truth. Ignoring the heteroscedasticity issue could impact downstream analysis including constructing microbial interaction networks (McMurdie and Holmes, 2014).

In this article, we provide a new statistical tool to help unleash the full potential of microbiome interaction networks as a research tool in the microbiome field. We adopt the logistic normal multinomial distribution to model the compositional count data (Aitchison, 1986; Billheimer et al., 2001; Xia et al., 2013). Compared to previous methods, this model accounts for the heteroscedasticity issue as the sequencing depth is treated as the number of trials in the multinomial distribution. Additionally, the additive log-ratio transformation applied to the multinomial probabilities results in a well-defined inverse covariance matrix in contrast to the centered log-ratio transformation. Based on this model, we develop an efficient algorithm that iterates between Newton-Raphson and graphical lasso for estimating a sparse inverse covariance matrix. We call this new approach “compositional graphical lasso”. We establish the theoretical convergence of the algorithm and illustrate the advantage of compositional graphical lasso in comparison to current methods under a variety of simulation scenarios. We further apply the developed method to the data from the Zebrafish Parasite Infection study (Gaulke et al., 2019) (see Section 1.3) to investigate how microbial interactions associate with parasite infection.

Helminth parasites represent a significant threat to the health of human and animal populations, and there is a growing need for tools to treat, diagnose, and prevent these infections. A growing body of evidence points to the gut microbiome as an agent that interacts with parasites to influence their success in the gut. To clarify how the gut microbiome varies in accordance with parasitic infection dynamics, the Zebrafish Parasite Infection Study was a recent effort (Gaulke et al., 2019) conducted at Oregon State University, which assessed the association of an intestinal helminth of zebrafish, Pseudocapillaria tomentosa, and the gut microbiome of 210 4-month-old 5D line zebrafish. Among these fish, 105 were exposed to P. tomentosa and the remaining 105 were unexposed controls. At each of the seven time points after exposure, a randomly selected group of 30 fish (15 exposed and 15 unexposed) were euthanized and fecal samples were collected. The parasite burden and tissue damage in P. tomentosa-infected fish were also monitored over 12 weeks of infection.

Previous analyses (Gaulke et al., 2019) of the Zebrafish Parasite Infection Study data have revealed that parasite exposure, burden, and intestinal lesions were correlated with gut microbial diversity. They also identified individual taxa whose abundance associated with parasite burden, suggesting that gut microbiota may influence P. tomentosa success. Numerous associations between taxon abundance, burden, and gut pathologic changes were also observed, indicating that the magnitude of microbiome disruption during infection varies with infection severity. However, it remains unclear how parasite success may disrupt or be modulated by the microbial interactions in the gut. Understanding how the microbial interactions associate with parasitic infection can help resolve potential drug or diagnostic test for parasitic infection.

To evaluate the performance of our method and to benchmark it against previously available tools, we take advantage of a unique data resource provided by the Tara Oceans Project on ocean planktons. This data set is particularly suitable for method comparison because it enjoys an experimentally validated sub-network of plankton interactome that has served as a gold standard for method benchmarking. Compared to other methods, compositional graphical lasso performs better in reconstructing the microbial interactions that are validated by the literature. In addition, it performs the best in picking up the keystone taxa, which are those with an excessive number of interactions with other taxa in the literature, such as Amoebophyra (Chambouvet et al., 2008), Blastodinium (Skovgaard et al., 2012), Phaeocystis (Verity et al., 2007), and Syndinium (Skovgaard et al., 2005). All these genera are well-described keystone organisms in marine ecosystems. Finally, compositional graphical lasso affords an opportunity to resolve novel modulators of community composition. For example, as one of a few described genera within the syndinean dinoflagellates—an enigmatic lineage with abundant diversity in marine environmental clone libraries, Euduboscquella (Bachvaroff et al., 2012) ranked high in the degree distribution uniquely by compositional graphical lasso.

To investigate the role of the gut microbiome plays in the parasite infections, we apply compositional graphical lasso to the Zebrafish Parasite Infection Study data. Interestingly, compositional graphical lasso identifies changes in interaction degree between infected and uninfected individuals for three taxa, Photobacterium, Gemmobacter, and Paucibacter, which are inversely predicted by other methods. Further investigation of these method-specific taxa interaction changes reveals their biological plausibility, and provides insight into their relevance in the context of parasite-linked changes in the zebrafish gut microbiome. In particular, based on our observations, we speculate on the potential pathobiotic roles of Photobacterium and Gemmobacter in the zebrafish gut, and the potential probiotic role of Paucibacter. Future studies should seek to experimentally validate the ecological roles of Photobacterium, Gemmobacter, and Paucibacter in the zebrafish gut, including their impacts on the rest of the microbial community and their roles in infection induced tissue damage.

Consider a microbiome abundance dataset with $n$ independent samples, each of which is composed of the observed counts of $K+1$ taxa, denoted by $\mathbf{x}_{i}=(x_{i,1},\ldots,x_{i,K+1})^{\prime}$ for the $i$-th sample, $i=1,\ldots,n$. Due to the compositional property of the data, the total count of all taxa for each sample $i$ is a fixed number, denoted by $M_{i}$. Naturally, a multinomial distribution is imposed on the observed counts:

where $\mathbf{p}_{i}=(p_{i,1},\ldots,p_{i,K+1})^{\prime}$ are the multinomial probabilities for all taxa and $\sum_{k=1}^{K+1}p_{i,k}=1$.

To apply the additive log-ratio transformation (Aitchison, 1986) on the multinomial probabilities, we choose one taxon, without loss of generality the $(K+1)$-th taxon, as a reference to which all the other tax are compared. The transformed multinomial probabilities are given by

Let $\mathbf{z}_{i}=(z_{i,1},\ldots,z_{i,K})^{\prime}$ for $i=1,\ldots,n$, and further assume that they follow an i.i.d. multivariate normal distribution

where $\boldsymbol{\mu}$ is the mean and $\boldsymbol{\Sigma}$ is the covariance matrix. Let $\boldsymbol{\Omega}=\boldsymbol{\Sigma}^{-1}$ be the inverse covariance matrix or the precision matrix.

The above model given in (1)–(3) is often referred to as the logistic normal multinomial model. In this model, a multinomial distribution is imposed on the compositional counts, which is the distribution of the observed data given the multinomial probabilities. In addition, to capture the variation of the multinomial probabilities across samples, we impose a logistic normal distribution on the multinomial probabilities as a prior distribution. We thereby obtain as our final model the logistic normal multinomial model, which is a hierarchical model with two levels.

The logistic normal multinomial model has a long history in modeling compositional count data and it has also been applied to analyze microbiome abundance data. For example, Xia et al. (2013) proposed a penalized regression under this model to identify a subset of covariates that are associated with the taxon composition. Our objective is different from Xia et al. (2013) as we aim to reveal the microbial interaction network by finding a sparse estimator of the inverse covariance matrix $\boldsymbol{\Omega}$ in (3). It is also noteworthy that Jiang et al. (2020) has the same objective as ours. However, Jiang et al. (2020) did not make full use of the logistic normal multinomial model as it focused on correcting the bias of a naive estimator of the $\boldsymbol{\Sigma}$ that does not require the logistic normal part of the model. By contrast, we aim to find an estimator of $\boldsymbol{\Omega}$ directly based on the logistic normal multinomial model.

From the logistic normal multinomial model in (1)–(3), we aim to derive an objective function of $\boldsymbol{\Omega}$ to estimate the microbial interaction network. To this end, we take a two-step procedure similar to SPIEC-EASI (Kurtz et al., 2015). In the first step, we find estimated values of $\mathbf{z}_{1},\ldots,\mathbf{z}_{n}$ based on the logistic normal multinomial model given $\mu$ and $\boldsymbol{\Sigma}$; in the second step, we find an estimate of $\boldsymbol{\Omega}$ based on the estimated values of $\mathbf{z}_{1},\ldots,\mathbf{z}_{n}$.

In the first step, we consider the posterior distribution of $\mathbf{z}_{1},\ldots,\mathbf{z}_{n}$ given the data $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$ and find the maximum a posteriori (MAP) estimates of $\mathbf{z}_{1},\ldots,\mathbf{z}_{n}$. For $i=1,\ldots,n$, the logarithm of the posterior density function of $\mathbf{z}_{i}$ given $\mathbf{x}_{i}$ is

where $\propto$ denotes that two quantities are equal up to a term not depending on $\mathbf{z}_{i}$. In the above derivation, we ignored the marginal density function of $\mathbf{x}_{i}$ that does not depend on $\mathbf{z}_{i}$, which does not affect the estimation of $\mathbf{z}_{i}$.

By independence between all the samples, the logarithm of the posterior density function of $(\mathbf{z}_{1},\ldots,\mathbf{z}_{n})$ given the data $(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})$ can be written as (again, ignoring a term independent of $\mathbf{z}_{1},\ldots,\mathbf{z}_{n}$):

Given the values of the multivariate normal parameters $\boldsymbol{\mu}$ and $\boldsymbol{\Omega}$, one can maximize (4) with respect to $(\mathbf{z}_{1},\ldots,\mathbf{z}_{n})$. This leads to the MAP estimator $(\hat{\mathbf{z}}_{1},\ldots,\hat{\mathbf{z}}_{n})$.

In the second step, we find a sparse inverse covariance estimator of $\boldsymbol{\Omega}$ based on the estimate values $(\hat{\mathbf{z}}_{1},\ldots,\hat{\mathbf{z}}_{n})$. Hereby, we use the graphical lasso estimator, which minimizes the $L_{1}$ penalized negative log-likelihood function as follows:

It turns out that the above two-step procedure is equivalent to minimizing an overall objective function with respect to both $\mathbf{z}_{1},\ldots,\mathbf{z}_{n}$ and $(\boldsymbol{\mu},\boldsymbol{\Omega})$:

In other words, $\mathbf{z}_{1},\ldots,\mathbf{z}_{n}$, $\boldsymbol{\mu}$, and $\boldsymbol{\Omega}$ are all treated as unknown parameters in the minimization of the objective function (6).

It is noteworthy that similar ideas to the above two-step procedure have been used in existing microbial interaction network estimation methods. For example, SPIEC-EASI is also a two-step procedure. In its first step, the abundance counts are converted into their centered log-ratio transformed data; in its second step, either graphical lasso or neighborhood selection is applied to estimate a sparse inverse covariance matrix based on the centered log-ratio transformed data in the first step.

Although both our method and SPIEC-EASI can be regarded as two-step procedures, we underline two important distinctions between them. First, our method is based on the additive log-ratio transformation and SPIEC-EASI uses the centered log-ratio transformation. As mentioned in the Introduction, the centered log-ratio transformation suffers from an undefined inverse covariance matrix of the transformed data while the additive log-ratio transformation does not. Second, in the first step, our method accounts for the sequencing depths [$M_{i}$’s in (1)] and further the uncertainty of the observed relative abundances, and thus addresses the heteroscedasticity issue. However, the heteroscedasticity issue is ignored in SPIEC-EASI.

The objective function (6) naturally includes three sets of parameters $(\mathbf{z}_{1},\ldots,\mathbf{z}_{n})$, $\boldsymbol{\mu}$, and $\boldsymbol{\Omega}$, which motivates us to apply a block coordinate descent algorithm. A block coordinate descent algorithm minimizes the objective function iteratively for each set of parameters given current values of the other sets. Given the initial values $(\mathbf{z}_{1}^{(0)},\ldots,\mathbf{z}_{n}^{(0)})$, $\boldsymbol{\mu}^{(0)}$, and $\boldsymbol{\Omega}^{(0)}$, a block coordinate algorithm repeats the following steps cyclically for iteration $t=0,1,2,\ldots$ until the algorithm converges.

1. 1.

   Given $\boldsymbol{\mu}^{(t)}$ and $\boldsymbol{\Omega}^{(t)}$, find $(\mathbf{z}_{1}^{(t+1)},\ldots,\mathbf{z}_{n}^{(t+1)})$ that maximizes (6).

2. 2.

   Given $(\mathbf{z}_{1}^{(t+1)},\ldots,\mathbf{z}_{n}^{(t+1)})$ and $\boldsymbol{\Omega}^{(t)}$, find $\boldsymbol{\mu}^{(t+1)}$ that maximizes (6).

3. 3.

   Given $(\mathbf{z}_{1}^{(t+1)},\ldots,\mathbf{z}_{n}^{(t+1)})$ and $\boldsymbol{\mu}^{(t+1)}$, find $\boldsymbol{\Omega}^{(t+1)}$ that maximizes (6).

Below, we will present the details of this algorithm in each iteration. For the initial values $(\mathbf{z}_{1}^{(0)},\ldots,\mathbf{z}_{n}^{(0)})$, we use their maximum likelihood estimators from the multinomial distribution, i.e.,

If $x_{i,k}=0$ for some $i$ and $k$, we add a small constant to it to evaluate the log ratio. For the initial value $\boldsymbol{\mu}^{(0)}$, we have a closed form minimizer of $\boldsymbol{\mu}$ for (6) given the values of $(\mathbf{z}_{1},\ldots,\mathbf{z}_{n})$, which is $\boldsymbol{\mu}=\bar{\mathbf{z}}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{z}_{i}$. Therefore, we set the initial value as $\boldsymbol{\mu}^{(0)}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{z}_{i}^{(0)}$. Finally, for the initial value $\boldsymbol{\Omega}^{(0)}$, we use the estimate of the graphical lasso algorithm taking the sample covariance matrix computed from $\mathbf{z}_{1}^{(0)},\ldots,\mathbf{z}_{n}^{(0)}$ as input.

In step 1, given $\boldsymbol{\mu}^{(t)}$ and $\boldsymbol{\Omega}^{(t)}$, minimizing the objective function (6) with respect to $(\mathbf{z}_{1},\ldots,\mathbf{z}_{n})$ is equivalent to minimizing the following objective function with respect to each $\mathbf{z}_{i}$ separately, for $i=1,\ldots,n$:

The above objective function is a smooth and convex function in $\mathbf{z}_{i}$ with the following positive definite Hessian matrix

where $e^{\mathbf{z}_{i}}=(e^{z_{i,1}},\ldots,e^{z_{i,K}})^{\prime}$ and $\text{diag}(e^{\mathbf{z}_{i}})$ is the diagonal matrix with the diagonal elements $e^{z_{i,1}},\ldots,e^{z_{i,K}}$. Therefore, we apply the Newton-Raphson algorithm to find the minimizer numerically. In addition, we implement a line search procedure in each Newton-Raphson iteration following the Armijo rule (Armijo, 1966). This procedure ensures sufficient decrease in the objective function at each iteration to prevent possible divergence of the algorithm.

Step 2 is similar to the initialization step, in which $\boldsymbol{\mu}$ has a closed-form solution and is updated as $\bar{\mathbf{z}}^{(t+1)}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{z}_{i}^{(t+1)}$ from the current numerical values of $(\mathbf{z}_{1}^{(t+1)},\ldots,\mathbf{z}_{n}^{(t+1)})$ that are computed from the Newton-Raphson algorithm in step 1.

In step 3, given $(\mathbf{z}_{1}^{(t+1)},\ldots,\mathbf{z}_{n}^{(t+1)})$ and $\boldsymbol{\mu}^{(t+1)}=\bar{\mathbf{z}}^{(t+1)}$, the objective function for $\boldsymbol{\Omega}$ can be simplified as

where $\mathbf{S}^{(t+1)}=\frac{1}{n}\sum_{i=1}^{n}(\mathbf{z}_{i}^{(t+1)}-\bar{\mathbf{z}}^{(t+1)})(\mathbf{z}_{i}^{(t+1)}-\bar{\mathbf{z}}^{(t+1)})^{\prime}$. It is obvious that minimizing the objective function (8) becomes a graphical lasso problem (Yuan and Lin, 2007; Banerjee et al., 2008; Friedman et al., 2008). It is well known that the graphical lasso objective function is a convex function in $\boldsymbol{\Omega}$ (Banerjee et al., 2008) and efficient algorithms have been developed for its optimization (Friedman et al., 2008). In this paper, we implement this step using the graphical lasso algorithm included in the huge (Zhao et al., 2012) package in R.

The above block coordinate descent algorithm iterates between Newton-Raphson and graphical lasso and is designed specifically to optimize the objective function (6) for compositional count data. Therefore, we name this algorithm the compositional graphical lasso algorithm, and the entire approach the compositional graphical lasso method including both the model and the algorithm for the analysis of microbiome abundance data.

Unfortunately, the objective function (6) is not necessarily a convex function jointly in $(\mathbf{z}_{1},\ldots,\mathbf{z}_{n})$, $\boldsymbol{\mu}$, and $\boldsymbol{\Omega}$. However, we have shown that it is convex in each of the three subsets of its parameters. The convergence property of such an optimization problem has been studied in the literature. For example, Tseng (2001) studied the convergence property of a block coordinate descent method applied to minimize a nonconvex function with certain separability and regularity properties. We will establish the convergence property of the compositional graphical lasso algorithm following Tseng (2001).

Recall that our algorithm treats the three sets of parameters $(\mathbf{z}_{1},\ldots,\mathbf{z}_{n})$, $\boldsymbol{\mu}$, and $\boldsymbol{\Omega}$ as three blocks and optimizes for each block iteratively. In addition, as in Tseng (2001), the objective function (6) can be regarded as the sum of two parts, the first of which is an inseparable but differentiable function given by

and the second of which is a sum of separable and differentiable functions given by $\ell_{1}(\mathbf{z}_{1},\ldots,\mathbf{z}_{n})+\ell_{2}(\boldsymbol{\Omega})$, where

Tseng (2001) established the convergence property of a block coordinate descent algorithm under regularity conditions on $\ell_{0}$, $\ell_{1}$, and $\ell_{2}$.

To present the main convergence property of the compositional graphical lasso algorithm, let’s review the definition of a cluster point in real analysis. A cluster point of a set $\mathcal{A}\subset\mathbb{R}^{n}$ is a real vector $\mathbf{a}\in\mathbb{R}^{n}$ such that for every $\delta>0$, there exists a point $\mathbf{x}$ in $\mathcal{A}\setminus\{\mathbf{a}\}$ satisfying that $\|\mathbf{x}-\mathbf{a}\|_{2}<\delta$. Obviously, any limit point of the set $\mathcal{A}$ is a cluster point. Furthermore, define a cluster point of the compositional graphical lasso algorithm to be a cluster point of the set $\{(\mathbf{z}_{1}^{(t)},\ldots,\mathbf{z}_{n}^{(t)},\boldsymbol{\mu}^{(t)},\boldsymbol{\Omega}^{(t)}):t=0,1,2,\ldots\}$, which are minimizers found at each iteration $t$. Then, the following theorem presents a theoretical property for every cluster point of our algorithm as follows.

The proof of Theorem 1 can be found in the supplementary materials. This theorem guarantees that a cluster point, usually a limit point, of the compositional graphical lasso algorithm, is at least a stationary point. It is noteworthy that there exists a global minimizer for the objective function in (6) because we have proved the coerciveness of the objective function in the proof (Calafiore and El Ghaoui, 2014, Lemma 8.4). Therefore, to achieve global optimization in practice, one can run the algorithm multiple times starting with different initial values and choose the one solution that yields the smallest objective function among the multiple ones.

In addition, the values of the objective function at each iteration, i.e.,
$\{\ell(\mathbf{z}_{1}^{(t)},\ldots,\mathbf{z}_{n}^{(t)},\boldsymbol{\mu}^{(t)},\boldsymbol{\Omega}^{(t)}):t=0,1,2,\ldots\}$, will always converge. This is because that the objective function is bounded below (as $\ell_{0}+\ell_{2}$ is bounded below as shown in the proof and $\ell_{1}>0$ by definition) and our algorithm results in non-increasing objective function values between two iterations. Therefore, the values of objective function will always converge to a limit point. In practice, we have always observed the numerical convergence of both the minimizers and the values of the objective function after a certain number of iterations.

There is a large body of literature on the selection of a tuning parameter in the variable selection framework. Common approaches can be broadly categorized into three types: criteirion-based methods, prediction-based methods, and stability-based methods. Criterion-based methods such as Akaike information criterion (AIC) (Akaike, 1974) and Bayesian information criteria (BIC) (Schwarz et al., 1978) balance the model complexity and the goodness of fit, prediction-based methods such as cross validation (Stone, 1974; Geisser, 1975) and generalized cross validation (Golub et al., 1979) aim to minimize the expected prediction error of the selected model on independent datasets. Stability-based methods such as stability selection (Meinshausen and Bühlmann, 2010) and the Stability Approach to Regularization Selection (StARS) (Liu et al., 2010) select a model with high stability under subsampling or bootstrapping the original data.

In this work, we apply StARS to select the tuning parameter $\lambda$ in our objective function (6). In StARS, we draw $N$ subsamples without replacement from the original dataset with $n$ observations, each subsample of size $b$. For each value of the tuning parameter $\lambda$, we obtain an estimate of $\boldsymbol{\Omega}$, i.e., a network for each subsample. Then, we measure the total instability of these resultant networks across the $N$ subsamples. The total instability of these networks is defined by averaging the instabilities of each edge across the $N$ subsamples over all possible edges, where the instability of each edge is estimated as the twice the sample variance of the Bernoulli indicator of whether this edge is selected or not in each of the $N$ subsamples.

Starting from a large penalty which corresponds to the empty network, the instability of networks increases as $\lambda$ decreases. StARS stops and selects the tuning parameter to be the minimum value of $\lambda$’s with which the instability of the resultant networks is less than a threshold $\beta>0$. In principle, StARS selects a tuning parameter so that the resultant network is the densest among networks with a total instability less than a threshold $\beta$ without violating some sparsity assumption. The selected network is the “densest on the sparse side,” as it starts with the empty network and stops when the instability first across the threshold.

To evaluate the performance of compositional graphical lasso, we conduct simulation studies and compare it with other network estimation methods.

Given that our goal is to estimate the true network, i.e., $\boldsymbol{\Omega}$ in (3), we consider the following three types of true precision matrices $\boldsymbol{\Omega}=(\omega_{kl})_{1\leq k,l\leq K}$, which are different in the pattern of edge distributions as well as the degree of connectedness.

1. 1.

   Chain network: $\omega_{kk}=1.5$, $\omega_{kl}=0.5$ if $|k-l|=1$, and $\omega_{kl}=0$ if $|k-l|>1$. A node is designed to be connected to its adjacent nodes, and the connectedness of nodes is balanced.

2. 2.

   Random network: $\omega_{kl}=1$ with probability $3/K$ for $k\neq l$. A node is connected to all other nodes randomly with a fixed probability, set to be $3/K$. Similar to the chain structure, the connectedness of nodes is balanced.

3. 3.

   Hub network: All nodes are randomly split into $\lceil K/20\rceil$ disjoint groups, and a hub node $k$ is selected from each group. For any other node $l$ in the same group, $\omega_{kl}=1$. All the remaining entries of $\boldsymbol{\Omega}$ are zero. Here, nodes are partitioned into the same group at random, but is then designated to be connected to the hub node at certain. The degree of connectedness among nodes is extremely unbalanced in this case: the hub nodes are connected to all the other nodes in its group (around 20 nodes) and all the other nodes are only connected to the hub node in its group, i.e., just one node.

In addition to the true network, we also consider two other factors that are expected to influence the result. The first factor is the sequencing depth, $M_{i}$, in the multinomial distribution (1). We simulate $M_{i}$ from uniform distributions, the detail of which will be discussed in the following subsections. The second factor is the degree of variation in the logistic normal distribution (3). For each of the three types of precision matrices, we consider an additional factor by multiplying a positive constant $c$ to $\boldsymbol{\Omega}$ so that the true precision matrix is $c\boldsymbol{\Omega}$. We choose $c=1$ and $c=1/5$ separately and call the two settings low and high compositional variation, respectively.

The data are simulated following the logistic normal multinomial model in (1)–(3). We first simulate $\mathbf{z}_{i}\sim N(\boldsymbol{\mu},\boldsymbol{\Sigma})$ independently for $i=1,\ldots,n$; then, we perform the inverse log-ratio transformation (also know as the softmax transformation, the inverse transformation of (2)) to obtain the multinomial probabilities $\mathbf{p}_{i}$ for $i=1,\ldots,n$; last, we simulate multinomial counts $\mathbf{x}_{i}$ from a multinomial distribution with sequencing depth $M_{i}$ and probabilities $\mathbf{p}_{i}$. Throughout this simulation study, we fix $n=100$ and $K=200$.

The simulation results are based on 100 replicates. For each replicate, we apply compositional graphical lasso, neighborhood selection, and graphical lasso separately to obtain a sparse estimator of $\boldsymbol{\Omega}$. For neighborhood selection and graphical lasso, we first obtain an estimate of $\mathbf{z}_{1},\ldots,\mathbf{z}_{n}$ from the multinomial distribution via the additive log-ratio transformation

and then apply neighborhood selection and graphical lasso directly on the estimates $\tilde{\mathbf{z}}_{1},\ldots,\tilde{\mathbf{z}}_{n}$ by treating them as surrogates for their true counterparts, i.e., $\mathbf{z}_{1},\ldots,\mathbf{z}_{n}$. These methods are almost identical to SPIEC-EASI although we replaced the centered log-ratio transformation in SPIEC-EASI with additive log-ratio transformation for a fair comparison.

To compare the performance of the three methods in terms of network recovery, all three methods are applied with a sequence of tuning parameter values, and their true positive rates (TPR) and false positive rates (FPR) in terms of edge selection are recorded for each value of $\lambda$. An ROC curve is plotted from the average TPR and the average FPR over the 100 replicates for each of the tuning parameters.

In addition, we apply StARS to select an optimal tuning parameter $\lambda$. Following the recommendation in Liu et al. (2010), we set the threshold for the total instability to be $\beta=0.05$, the size of each subsample $b=7\sqrt{n}$, and the number of subsamples $N=50$. Once the optimal tuning parameter is determined by StARS, we fit the whole dataset with the selected tuning parameter and evaluate the resultant network using three criteria: precision, recall, and F1 score, which are defined as

where TP, FP, and FN are numbers of true positives, false positives, and false negatives, respectively.

In this subsection, we evaluate the performance of compositional graphical lasso on dense data, in which most of the simulated counts are nonzero. This corresponds to a simulation setting where the sequencing depths $M_{i}$’s are relatively larger than the number of taxa $K+1$. Still, to evaluate the effect of the sequencing depth on the performance of network estimation methods, we simulate $M_{i}$ from two uniform distributions, Uniform$(20K,40K)$ and Uniform$(100K,200K)$, and call the two settings low and high sequencing depth in this subsection, respectively.

Figure 1 presents the ROC curves for compositional graphical lasso (Comp-gLASSO), neighborhood selection (MB), and graphical lasso (gLASSO), from which we can see that compositional graphical lasso dominates its competitors in terms of edge selection in all settings. In particular, the advantage of the compositional graphical lasso over neighborhood selection and graphical lasso is the most obvious when the compositional variation is high and the sequencing depth is low, no matter which type of network structure is considered. On the contrary, the three methods perform very similarly for all types of network structures when the compositional variation is low and the sequencing depth is high. The difference between compositional graphical lasso and the rest is intermediate for the other two settings when both compositional variation and sequencing depth are high or when both are low. Comparing graphical lasso and neighborhood selection, they tend to perform more similarly although graphical lasso seems to outperform neighborhood selection in some settings with a small margin.

The above observations agree with our expectation about how the two factors, compositional variation and sequencing depth, affect the comparison between the methods. Recall that neighborhood selection and graphical lasso replace the true values of $\mathbf{z}_{1},\ldots,\mathbf{z}_{n}$ by their estimates/surrogates $\tilde{\mathbf{z}}_{1},\ldots,\tilde{\mathbf{z}}_{n}$ as in (12) without taking into account the estimation accuracy or uncertainty of these surrogates. First, a higher sequencing depth leads to more accurate surrogates $\tilde{\mathbf{z}}_{1},\ldots,\tilde{\mathbf{z}}_{n}$; therefore, it is not surprising to see that the three methods perform more similarly when the sequencing depth is high. Second, a higher compositional variation results in a higher variation in $\mathbf{z}_{i}$’s and further in $\mathbf{p}_{i}$’s that are multinomial probabilities. Since neighborhood selection and graphical lasso ignore the multinomial component in the model, it is also not surprising to see that their performance is deteriorated by a high compositional variation.

Figure 2 presents the recall, precision, and F1 score from 50 replicates of the estimated network resulted from the tuning parameter selected by StARS. The first observation would be that the precisions of both compositional graphical lasso and graphical lasso are much worse than their recalls, whereas the precisions and recalls are more comparable for neighborhood selection. Interestingly, StARS results in a much more sparse network for neighborhood selection than the other methods under the same stability threshold, suggesting that fewer edges selected by neighborhood selection are stable enough (within the instability threshold). When it comes to method comparison, compositional graphical lasso has much higher recall than neighborhood selection in most settings, but have comparable or lower precision in most of the settings with high sequencing depth. The network from compositional graphical lasso has higher F1 score than the ones from neighborhood selection in most settings, except when sequencing depth is high and compositional variation is low for chain and hub networks. In addition, the network from compositional graphical lasso has higher or comparable precision, recall, and F1 score than the ones from graphical lasso in all settings. Similar to the observations from the ROC curves, the advantage of compositional graphical lasso is more obvious with a low sequencing depth or a high compositional variation.

In this subsection, we evaluate the performance of compositional graphical lasso on sparse data, in which a range of proportions of the compositional counts are simulated as zero. We only present the simulation results for the chain network, given that the simulation results for the other network types are similar. To examine how our method performs for different sparsity levels, we simulate $M_{i}$ from four uniform distributions, Uniform$(8K,16K)$, Uniform$(4K,8K)$, Uniform$(2K,4K)$, and Uniform$(K,2K)$, respectively. Note that the sparsity level of the simulated data also depends on the compositional variation (see Section 3.1). In a typically simulated dataset with high compositional variation, the sparsity level is around 40%, 50%, 60%, and 70% with the four uniform distributions for $M_{i}$; in a typically simulated dataset with low compositional variation, the sparsity level is around 5%, 10%, 25%, and 40% with the four uniform distributions for $M_{i}$. In both cases, we refer to these four sparsity levels as 1, 2, 3, and 4, with 1 the least sparse and 4 the most sparse. Due to the limited space, we place the figures resulted from this simulation in the supplementary materials.

Figure S1 presents the ROC curves for compositional graphical lasso (Comp-gLASSO), neighborhood selection (MB), and graphical lasso (gLASSO). Similar to Figure 1, we can see that compositional graphical lasso dominates its competitors in terms of edge selection across all sparsity levels. As the sequencing depths for sparse data are relatively small compared to the ones used for dense data, the advantage of the compositional graphical lasso over neighborhood selection and graphical lasso is obviously and consistently observed in all simulation settings. In addition, graphical lasso and neighborhood selection tend to perform more similarly although graphical lasso seems to outperform neighborhood selection with a small margin. This result is consistent with that observed for dense data in Section 3.2.

Figure S2 presents the recall, precision, and F1 score from 50 replicates of the estimated network resulted from the tuning parameter selected by StARS. The first observation is that all methods perform worse for sparse data than for dense data from the comparison between Figure S2 and Figure 2; all methods perform worse as the sparsity level increases. Similar to Figure 2, Figure S2 implies that compositional graphical lasso outperforms graphical lasso and neighborhood selection, with a higher recall, precision, and F1 score consistently observed in most settings. The only exception is that, when data are extremely sparse with a low compositional variation, graphical lasso seems to outperform compositional graphical lasso by a very small margin as both methods are not working well. For sparse data, StARS results in an almost empty network for neighborhood selection, resulting in zero recall, precision, and F1 score for most settings. This agrees with our observation from the simulation results for dense data that neighborhood selection tends to produce a much sparser network than compositional graphical lasso and graphical lasso when its tuning parameter is selected by StARS.

To better understand the ocean, the largest ecosystem on the earth, the Tara Oceans Project aims to build the global ocean interactome that can be used to predict the dynamics and structure of ocean ecosystems. To achieve this, the Tara Oceans Consortium sampled both plankton and environmental data at 210 sites from the world’s oceans using the 110-foot research schooner Tara during the Tara Oceans Expedition (2009-2013). The data collected was later processed using sequencing and imaging techniques. One unique advantage of the Tara Oceans Project is that it has generated a list of 91 genus-level marine planktonic interactions that have been validated in the literature (Lima-Mendez et al., 2015). Though this list only comprises interactions between microbes that represent a small fraction of the total marine eukaryotic diversity and is therefore far from complete, it could serve as partial ground truth for us to evaluate the interactions identified by different methods. Thus, our major goal is to use the Tara Oceans Project as a benchmark study in order to compare the performance of different methods in constructing the planktonic interactions.

As the partial ground truth is a list of genus-level interactions, we choose to analyze the genus-level abundance data, which are aggregated from the original OTU abundance data downloadable from the Tara Oceans Project data repository (http://taraoceans.sb-roscoff.fr/EukDiv/). As a benchmark study, we only include the 81 genera that are involved in the list of gold-standard interactions in our analysis. In addition, we discard the samples with too few sequence reads (less than 100), resulting in 324 samples left in our analysis.

Similar to the simulation study, we apply compositional graphical lasso, graphical lasso, and neighborhood selection to estimate the interaction network among the 81 genera. We first pick the genus Acrosphaera, which has the largest average relative abundance among those genera not involved in the gold-standard list, and use this genus as the reference taxon for all three methods. Then, we apply each method with a sequence of 70 decreasing tuning parameter values, resulting in a sequence of interaction networks starting from an empty network. Finally, we apply StARS to find the optimal tuning parameter, in which the parameters $\beta$, $b$, and $N$ are set the same as in the simulation study.

First, we compare the three methods in terms of how highly they rank the literature validated interactions amongst their top reported edges. Specifically, we start with a large value of the tuning parameter that results in an empty network, then decrease the tuning parameter so that the network becomes denser, and stop until the network has about 200 edges (out of a total of 3240 possible edges). At each value of the tuning parameter, we plot the number of literature validated interactions included in the network versus the total number of edges of the network, resulting in a step function for each method (Figure 3(a)). From Figure 3(a), we observe that compositional graphical lasso identifies slightly more literature validated interactions than graphical lasso until the total number of edges arrives 175 and graphical lasso identifies one more literature validated interaction afterwards. Neighborhood selection selects much fewer literature validated interactions than either compositional graphical lasso or graphical lasso. These observations imply that compositional graphical lasso slightly outperforms graphical lasso in reconstructing the literature validated interactions, while its advantage over neighborhood selection is much more obvious.

Second, we compare the overall topologies of the three interaction networks with tuning parameters selected by StARS. We find that compositional graphical lasso, graphical lasso, and neighborhood selection identify 749, 921, and 190 edges, respectively, with the same instability threshold used in StARS. This agrees with our observation in the simulation study that the network from neighborhood selection is much sparser than those from compositional graphical lasso and graphical lasso. The degree distributions from the networks estimated by the three methods are shown in Figure 3(b). The center of the three degree distributions are ranked as neighborhood selection, compositional graphical lasso, and graphical lasso in the ascending order, which is also reflected in the densities of the three interaction networks.

Third, we investigate the high-degree nodes, i.e., hub genera, in the interaction networks. High-degree nodes are often thought to represent taxa that elicit an effect on a large number of other members of their community, such as keystone taxa and generalistic parasites. It is observed that there are a few hub genera that have an excessive number of interactions with other genera in the literature, such as Amoebophrya, Blastodinium, and Parvilucifera, which are referred to as benchmark hub genera. Although the literature validated interactions are rather incomplete, it is still of interest to evaluate how well the three methods pick up those benchmark hub genera. Since the density of networks from the three methods are rather different, it is hard to compare the degrees of the hub genera from the three networks directly, but it is reasonable to compare the ranks of those degrees within each degree distribution. The method that generates lower ranks (degree of genera ranked in descending order) for those hubs in their degree distribution are believed to pick up the benchmark hub genera better.

A list of 7 benchmark hubs (which has degree $\geq 5$) along with their degrees from the incomplete network constructed from the literature is shown in Table 1, followed by the corresponding ranks of those genera in the degree distributions from each of the three methods and their corresponding degrees in the parentheses. We can see that compositional graphical lasso generates lower ranks than graphical lasso for all 7 genera, while neighborhood selection generates lower ranks than compositional graphical lasso for 3 genera, and the opposite for the other 4 genera. Overall, compositional graphical lasso performs the best in picking up the benchmark hub genera among the three methods.

We find that compositional graphical lasso predicts a high degree for genera that are known to act as keystone species or parasitize other taxa, many of which are also high-degree in the literature validated network. For example, Phaeocystis elicits a high degree in both literature validated and compositional graphical lasso networks. This genus is a well-described keystone organism in some marine ecosystems (Verity et al., 2007), where it causes large phytoplankton blooms and plays an important role in the global cycling of carbon and sulfur (Verity and Smetacek, 1996). Similarly, these two networks predict a high degree for taxa that are known parasites, such as Blastodinium (Skovgaard et al., 2012), Amoebophyra (Chambouvet et al., 2008), and Syndinium (Skovgaard et al., 2005). In some instances, compositional graphical lasso uniquely reveals high-degree genera even when compared to the literature validated network, such as in the case of the parasite Euduboscquella (Bachvaroff et al., 2012), which is ranked the 5th by compositional graphical lasso but only has 1 interaction in the literature. Given the fact that the literature validated network only captures a subset of interactions that exist in nature (i.e., those interactions that have been explicitly tested), we posit that compositional graphical lasso affords an opportunity to resolve novel modulators of community composition, such as keystone taxa and generalistic parasites, that follow-up experiments can validate.

The Zebrafish Parasite Infection Study, recently conducted at Oregon State University, used a zebrafish helminth intestinal infection model to resolve how the gut microbiome and parasite burden covary during infection (Gaulke et al., 2019). This study quantified the infection burden of an intestinal helminth of zebrafish, Pseudocapillaria tomentosa, and the gut microbiome in 210 4-month-old zebrafish. Half of these fish were exposed to the parasite, P. tomentosa, and the other half were unexposed. Given that not all exposed fish would be ultimately infected, parasite burden was measured more accurately as the total number of worms in their fecal samples. In addition, the same fecal samples were used to measure the abundance of the gut microbiome species via 16S rRNA sequencing, which resulted in gut microbiome data for 207 fish. Among these fish, 81 were infected after being exposed to P. tomentosa and the other 126 fish were not infected, as indicated by the total number of worms. Our major goal is to evaluate how the microbial interaction network associates with successful parasite infection in the gut.

Similar to the Tara Oceans Project, we choose to analyze the genus-level abundance data. We discard those genera that have a nonzero abundance in less than 5% of the samples, resulting in 42 genera left in our analysis. In addition, we follow the same strategy to define a reference genus as in Jiang et al. (2020), i.e., combining those OTUs that do not have a genus-level taxonomic classification into a pseudo-genus and using it as the reference genus. The analysis is conducted in the same way as in the Tara Oceans Project, except that the methods (compositional graphical lasso, graphical lasso, and neighborhood selection) are applied separately for the uninfected and the infected fish, resulting in three interaction networks for each group of fish. Analyzing the uninfected and infected fish separately allows us to compare the microbial interaction networks between the two groups.

First, we compare the overall topologies of the interaction networks with tuning parameters selected by StARS. We find that compositional graphical lasso, graphical lasso, and neighborhood selection identify 262, 312, and 78 edges, respectively, for uninfected fish, and 241, 259, and 60 edges, respectively, for infected fish. Comparing the two groups of fish, the interaction networks for the uninfected fish are slightly denser than the infected fish. The comparison among the three methods agrees with our observation in the simulation study and in the Tara Oceans Project that the network from neighborhood selection tends to be much sparser than those from compositional graphical lasso and graphical lasso. The degree distributions based on the three methods are shown in Figure 4, which again suggests high similarity between the two groups of fish. Similar to the Tara Oceans Project, neighborhood selection results in the lowest median degree, followed by composition graphical lasso and graphical lasso, regardless of whether the fish are infected or not.

Second, we further investigate the degree of each node in the interaction networks. This analysis helps identify those high-degree nodes, i.e., hub genera, which are indicative of keystone taxa in the microbial community. Due to the limited space, we refer to Table S2 in the supplementary materials that presents all 42 genera with their degrees in each network in descending order. Although the networks are similar in density between the two groups of fish, we have observed a substantial difference between the degrees of individual nodes. To see this we sort the nodes in each network based on their degrees and compare the top ten nodes between networks. Taking compositional graphical lasso as an example, only three out of the top ten nodes overlap between uninfected and infected fish: Paicibacter, Photobacterium, Fusobacterium. This suggests that the network structures between the uninfected and infected fish are more distinctive than what their overall topologies might appear. The hub genera identified by different methods are quite different from each other too. For the group of uninfected fish, only three out of the top ten nodes overlap between the networks from the three methods: Aeromonas, Photobacterium, and Rheinheimera; similarly, for the group of infected fish, only three out of the top ten nodes overlap between the networks from the three methods: Fusobacterium, Paucibacter, and Yersinia.

Third, we compare the networks generated from uninfected and infected fish to clarify the relationship between infection and the interactions that exist in the gut microbiome, motivated by the hypothesis that intestinal dysbiosis is defined not only by changes in community composition, but also by how members of the community ecologically interact. In order to understand this relationship, we identify microbes whose interaction set changes between infected and uninfected hosts then compare these results across methods. One of the bacterial taxa that compositional graphical lasso identifies as increasing in interaction degree among infected individuals, while other methods estimate a decrease in detected interaction degree, is the genus Gemmobacter. Our earlier work (Gaulke et al., 2019) showed that the abundance of this taxon positively links to parasite exposure, indicating that it is more ecologically successful in infected versus uninfected intestines. The analysis from compositional graphical lasso indicates that this increase in Gemmobacter relative abundance in infected individuals is coincident with an increase in the ecological interactions between Gemmobacter and other members of the gut microbiome, as it interacts with a greater number of other taxa in infected individuals. While Gemmobacter is not a particularly well studied genus, prior work links it to dysbiotic diseases in a variety of hosts (Bates et al., 2022; Ni et al., 2020). Our observation of an increase in its interaction with other taxa in infected individuals strengthens the hypothesis that Gemmobacter may act as an agent of dysbiosis, as its increase in its relative abundance not only links to a disease context, but its more integrated role in the microbial community (i.e., it potentially impacts a larger number of taxa in the gut).

The results from compositional graphical lasso also uniquely suggest that the Paucibacter genus in uninfected individuals displays a higher number of identified interactions relative to infected individuals. The other methods we evaluated find that Paucibacter’s interactions are relatively consistent between infected and uninfected individuals. The results from compositional graphical lasso are notable because the genus Paucibacter includes clades that have been identified as being conserved in the healthy zebrafish gut microbiome (Sharpton et al., 2021), presumably because these taxa are critical to the commensal microbial community. Based on our observations, namely that disruption to the broader microbiome by parasite infection appears to disrupt the set of interactions between Paucibacter and other taxa in the community, we hypothesize that Paucibacter is a keystone member of the healthy zebrafish gut microbiome and that parasite exposure induced dysbiosis occurs in part by disrupting Paucibacter in the gut. Collectively, our observations reveal patterns of interaction between infected and uninfected individuals that clarify the potential ecological role of Gemmobacter and Paucibacter, which future research should seek to empirically test.

In addition to these observations about Gemmobacter and Paucibacter, compositional graphical lasso also finds insightful patterns about Photobacterium, a genus which contains pathogens of marine fish (Romalde, 2002; Osorio et al., 2018), though its role as a pathogen in zebrafish is less clear. However, as prior work has cultured Photobacterium from the zebrafish gut (Cantas et al., 2012) and other studies have demonstrated that Photobacterium elicits toxicity to zebrafish cells in culture (Osorio et al., 2018), observations collectively point to Photobacterium as a potential pathogen or pathobiont of zebrafish. Consistent with this speculation, prior work found that Photobacterium statistically explained variation in P. tomentosa worm burden as well as worm infection induced intestinal inflammation (Gaulke et al., 2019). We sought to determine if our interaction network analysis could clarify the role of Photobacterium in the zebrafish gut microbiome, especially as it relates to parasite infection.

Interestingly, Photobacterium is identified as an interaction hub (i.e., a relatively high degree node) in all networks assembled from all methods we evaluated. However, compositional graphical lasso uniquely reveals that the number of interactions linked to Photobacterium increases in infected fish (19 edges) as compared to uninfected fish (16 edges), whereas the other approaches observe the opposite pattern (graphical lasso: 11 versus 17 edges; neighborhood selection: 3 versus 5 edges). These differences in the Photobacterium subgraph across approaches also includes variation in the specific taxa that Photobacterium is inferred to interact with. To better understand changes in the set of Photobacterium interactions, we extracted the first-degree sub-graph connected to Photobacterium that each of the three methods produced for both the infected and uninfected host interaction networks. Compositional graphical lasso alone estimates eight interactions in the Photobacterium-specific infected host network, and six of these eight genera have previously been linked to parasite exposure, disease burden, or histopathology score. Two of these parasite etiology-linked taxa are Aeromonas and Gemmobacter, genera which are notable because they contain microbes that elicit pathogenic or pathobiotic phenotypes. For example, the Aeromonas genus includes well-characterized and abundant opportunistic pathogens, such as A. hydrophila (Saraceni et al., 2016), and the Gemmobacter genus, as noted above, has been shown to opportunistically increase abundance (Huang et al., 2020) in dysbiotic fish gut communities which display impeded overall gut function. Given that Photobacterium proliferates in the infected host gut, these observations suggest that Photobacterium may work alongside and even promote the growth of other pathobiotic taxa to disrupt the composition of the gut microbiome and induce dysbiosis upon infection.

Based on these collective observations, we hypothesize that Photobacterium is an intestinal pathobiont of the zebrafish gut and contributes to dysbiosis in infected fish. Photobacterium is a relatively prevalent genus in the zebrafish intestines, even in uninfected fish. However, P. tomentosa infection links to an increase in the relative abundance of Photobacterium, which also positively associates with intestinal hyperplasia (Gaulke et al., 2019), possibly due to the cytotoxic effect of Photobacterium. Given that Photobacterium is also an interaction hub whose influence on the community increases in infected fish, at least according to compositional graphic lasso, infection induced changes in Photobacterium relative abundance may drive additional changes in the success of other taxa in the microbiome, including opportunistic pathogens and pathobionts like members of Aeromonas and Gemmobacter. Going forward, future studies should seek to experimentally validate our novel hypothesis about the pathobiotic role of Photobacterium in the zebrafish gut, including its impact on the rest of the microbial community and its role in infection induced tissue damage. If our hypothesis is accurate, Photobacterium may serve as an important model taxon for discerning how the gut microbiome and helminths interact to impact infection outcomes.

Collectively, these results provide evidence for our hypothesis that in the case of intestinal helminth infection, dysbiosis may be defined by not only a change in the composition of the gut microbiome, but a restructuring in how key members of the microbial community interact with the rest of the microbiota. Notably, these patterns were only observed by compositional graphical lasso, as the other methods did not observe unique evidence of such infection-associated inversions for Gemmobacter, Paucibacter, or Photobacterium.

Finally, we visualize the six genus interaction networks, three for infected fish and three for uninfected fish (Figure 5). For a better visualization, we only keep the top $100$ edges for the networks resulted from compositional graphical lasso and graphical lasso as they are much denser than the ones from neighborhood selection. Their edges are ranked in the exactly same way as in the Tara Oceans Study, first by selection probability and then by edge weight (see Figure S3 in the supplementary materials for details). For all networks, darker blue implies higher magnitude in the absolute value of partial correlation.