extfieldsofstudy: Biology Mathematics

text:
Narrowing the gap on heritability of common disease by direct estimation in case-control GWAS

One of the major developments in recent years in the search for missing heritability of human phenotypes is the adoption of linear mixed-effects models (LMMs) to estimate heritability due to genetic variants which are not significantly associated with the phenotype. A variant of the LMM approach has been adapted to case-control studies and applied to many major diseases by Lee et al. (2011), successfully accounting for a considerable portion of the missing heritability. For example, for Crohn's disease their estimated heritability was 22% compared to 50-60% from family studies. In this letter we propose to estimate heritability of disease directly by regression of phenotype similarities on genotype correlations, corrected to account for ascertainment. We refer to this method as genetic correlation regression (GCR). Using GCR we estimate the heritability of Crohn's disease at 34% using the same data. We demonstrate through extensive simulation that our method yields unbiased heritability estimates, which are consistently higher than LMM estimates. Moreover, we develop a heuristic correction to LMM estimates, which can be applied to published LMM results. Applying our heuristic correction increases the estimated heritability of multiple sclerosis from 30% to 52.6%.

Main text
The mystery of the "missing heritability" is a term commonly used to denote the gap between the expected heritability of many common diseases, as estimated by family and twin studies, and the overall additive (narrow-sense) heritability obtained by accumulating the effects of all singlenucleotide polymorphisms (SNPs) that have been found to be significantly associated with these conditions in genome-wide association studies (GWASs) [6][7][8][9] . Many diseases which comprise a considerable portion of the health-care burden display such a gap, including type-1 and type-2 diabetes, bipolar disorder, schizophrenia, Alzheimer's disease, multiple sclerosis and Parkinson's disease.
Researchers have proposed several hypothetical solutions to this mystery. These theories include rare causative variants, which are undetected by the current GWAS methodology, common variants with small effects, which do not pass the significance threshold and are therefore unaccounted for, gene-gene and gene-environment interactions which are overlooked by the additive model assumed by the GWAS scheme, epigenetic effects and more 7,9,10 .
Clearly, different theories have major implications for our understanding of human disease, and also dictate different strategies for discovery of the underlying genetic causes of disease. For example, identifying rare variants requires a focus on deep sequencing 11 , while detecting small effects of common variants requires increasing the sample sizes dramatically, or conducting massive meta-analyses 12 . Effective planning of genetic research can therefore only be guided by a satisfactory and well founded allocation of missing heritability to its various possible sources and causes. Additionally, accurate estimates of heritability can facilitate better personalized genetic risk predictions 13 .
Yang et al. 1 pioneered the use of LMMs to estimate heritability of continuous traits (e.g. height) from GWAS data, while accounting both for significantly and insignificantly associated SNPs, thus providing an estimate of the total heritability explained by common SNPs. Their method was applied to numerous phenotypes including human height, body mass index, von Willebrand factor 14 , gene expression 15 and intelligence 16 . The LMM method was later adapted to dichotomous disease phenotypes by Lee et al. 2 , assuming the well-known liability threshold model 17 .
They showed that for the case of a random (unascertained) population sample the LMM method can be used to estimate heritability by applying it directly to the observed phenotype as if it were continuous, and correcting the resulting "observed scale" estimate as suggested by Dempster and Lerner 17 .
The ascertained case-control scenario, which is relevant for most GWAS studies, proved more challenging, as the enrichment of cases violates the assumption of normality of the genetic component critical for LMM estimation. Lee et al. proposed a complex mathematical solution for this case, and demonstrated that it worked in limited simulations. They applied their method to three phenotypes from the Wellcome Trust Case-Control Consortium (WTCCC) 18 : Crohn's disease (CD), bi-polar disorder (BD) and type-1 diabetes (T1D), and it has since been applied to many major diseases including schizophrenia 4 , Alzheimer's disease, multiple sclerosis 5 and Parkinson's disease 3 . In all these cases, accounting for insignificant SNPs resulted in much higher estimates of heritability than the estimates obtained by accumulating SNPs which were found significant in GWAS.
The basic idea of heritability estimation methods in GWAS is that individuals who are more correlated genetically are more likely to have similar phenotypes. The strength of this connection depends on the heritability -higher heritability implies stronger connection. The genetic correlation of every pair of individuals can be estimated from their genotypes. Since the individuals in the study are unrelated, the correlations between their genotypes are typically small, but non-zero. By accumulating the information across all pairs of individuals in the study, one can leverage these minor differences to separate the phenotypic variance into its genetic and environmental components, resulting in an estimate of heritability.
Ascertained case-control studies pose a considerably harder challenge to deal with theoretically than non-ascertained (prospective or observational) studies, as the fact that cases are overrepresented in the study creates a wide range of artifacts. The normality of the genetic component, assumed by the liability threshold model, is violated by over-sampling of cases, as is the assumption that the genetic and environmental components are independent. Methods that critically rely on these assumptions, like LMM, can not be applied anymore, without requiring complex and potentially questionable adjustments.
We therefore developed our GCR method for estimating heritability of polygenic phenotypes in ascertained case-control studies. The key novelty of our method is that we model the selection process directly and account for ascertainment by conditioning on the selection of observed in-dividuals. This is contrary to the existing methods of Lee et al. 2 and Zhou et al. 19 , which apply methodologies suited for prospective studies of a continuous normal phenotype, and then attempt to correct the estimates so that they account for ascertainment.
Technically, we adopt the commonly used liability threshold model 2,17 , and derive analytically the relationship between the genetic correlation of any two individuals, and the similarity between their phenotypes, conditional on the fact that both individuals were selected for the study.
While this relationship has no closed-form expression, taking the first-order approximation yields a linear relationship between the expected phenotypic similarity and the genetic correlation which depends on the heritability. Our approach then entails performing a regression of phenotype similarities between pairs of individuals on their genetic correlation, as estimated from GWAS data.
The slope of the regression is then transformed to an estimate of the heritability (Online Methods).
We applied GCR to the same three phenotypes from WTCCC: CD, BD and T1D, and four additional phenotypes: type-2 diabetes (T2D), coronary artery disease (CAD), rheumatoid arthritis (RA) and hypertension (HT). Our method yields considerably higher estimates than the LMM method 2 for CD and BD (Table 1).
To explore the nature of this considerable gap between our estimates and previously published estimates of the additive heritability, we conducted a wide range of realistic simulations.
We simulated an end-to-end generative model, starting from minor allele frequencies (MAFs) and individual SNP effects, through genotypes and phenotypes, and finally the entire selection process (Online Methods).
Our simulations differ from all previously reported simulations of case-control studies in the context of heritability estimation in that they yield realistic, rather than degenerate, correlation structure among the individuals in the study. Since we simulated genotypes, the correlation between individuals is not restricted to a small subset of possible values as in the simulations of Lee et al. 2 . Additionally, when both heritability and ascertainment are high, cases tend to be more genetically similar than expected by chance, which was not accounted for in Lee et al.'s simulations.
In particular, their assumed correlation structure is highly degenerate, with as many as 99.98% of correlations being exactly 0. For further discussion of the different simulation setups see Online Methods and Supp. Material.
Our simulations showed that using the LMM method for heritability estimation in casecontrol studies yielded considerably negatively-biased estimates, while our method consistently generated unbiased estimates ( Fig. 1 and Supp. Figs. 1-7).
The simulations also demonstrated a strong relation between the bias of the LMM approach and the increased variance of genetic effects due to ascertainment. This allowed us to derive a heuristic correction for the LMM heritability estimates, which corrected its bias successfully in all our simulation scenarios (Online Methods, Fig. 2 and Supp. Figs. [12][13][14][15]. To validate our heuristic correction, we applied it to LMM heritability estimates for all seven WTCCC phenotypes and compared the corrected estimates to the ones obtained using GCR. The correlation between the corrected LMM estimates and GCR estimates was 0.98, indicating that our heuristic correction generalizes beyond our particular simulations (See Online methods and Supp. material for more details).
We then applied our correction to the results of several published studies which used the LMM approach. As expected, for studies with low ascertainment or phenotypes with low heritability, the corrected estimates were not substantially different from published estimates (Alzheimer's disease 5 , endometriosis 5 , schizophrenia 4 and Parkinson's disease 3 ). However, the corrected heritability estimate of multiple sclerosis -the most ascertained study we inspected -is 52.6%, compared to the uncorrected estimate of 30%. See Supp. material for more details.
An important aspect of heritability estimation is the inclusion of known (fixed) effects that should be accounted for, like known associated SNPs, known environmental effects or sex. However, estimation of fixed effects from ascertained data under the normality assumption might produce biased estimators of both the fixed effects and the heritability [20][21][22] . Conversely, our method can be rigorously extended to allow for fixed effects (Online Methods). Our simulations suggest that GCR produces accurate estimates of heritability in the presence of fixed effects (see Supp. material for more details). Specifically for the WTCCC phenotypes, previous analysis 18 suggested that there's little population structure (after removing individuals of non British descent). We validated this conclusion using the statistical test described in Patterson et al. (2006) 23 (see Supp. material for more details). We therefore included only sex as a fixed effect in our heritability estimation in Table 1. The estimates of both LMM and GCR for most phenotypes without inclusion of a sex fixed effect are not substantially different from the estimates in Table 1, with the exception of the heart-disease related phenotypes (results without fixed effect not shown). For CAD in particular, it is well known that sex is a major risk factor. Accordingly, the GCR estimate without fixed effect is 71.5%, dropping to about 61% once it is added, as seen in Table 1. GCR is computationally very efficient compared to the LMM approach (running time scales quadratically rather than cubically in the number of individuals), and so the running time is seconds rather than hours on WTCCC cohorts. As the typical size of GWAS continues to increase, this efficiency can prove critical in allowing GCR to remain computationally practical compared to other approaches. Fast running time is also useful in allowing us to use resampling approaches like the jackknife 24 for estimating confidence intervals and standard errors, rather than relying on complex and potentially inaccurate parametric approximations. All standard errors and confidence intervals we report are based on the jackknife.
We also experimented with a more involved version of our approach that included the second order term of the Taylor series expansion. This did not yield more accurate results than our first order GCR in our simulations (Supp. material).
Recently, Zuk et al. 10 have shown an intriguing general result on the connection between the derivative of the dependence of phenotypic similarity on the proportion of identity by descent (IBD) at the average population IBD, and the narrow-sense heritability, regardless of population structure and phenotype-genotype architecture. While Zuk et al.'s result uses IBD and our method relies on estimated genetic correlations, the latter is, in fact, an unbiased estimator of the former, which is unknown. Hence, our method can produce unbiased estimators of narrow-sense heritability, even in the presence of population structure, by multiplying the estimate by a factor of (1-average population IBD).
In conclusion, our new proposed GCR method for estimation of heritability in case-control GWAS, which is based on a regression of phenotype similarities on genotype correlations, is shown to be efficient and accurate. It improves on existing methodology and generates substantially higher estimates of heritability for two major diseases inspected: Crohn's disease and bipolar disorder. Moreover, we provide a heuristic correction for published LMM heritability estimates, which suggests that the heritability of multiple sclerosis is also considerably larger than previously thought.

Methods
Note: An implementation of our code for GCR estimation, simulations, and heuristic correction of LMM can be accessed at: https://sites.google.com/site/davidgolanshomepage/gcr Heritability estimation using genetic correlation regression Liability threshold model -notations. Denote K the prevalence of a condition in the population and P the prevalence in the study.
Under the liability threshold model, we assume that each individual i has an unknown liability l i = g i + e i where g i is a genetic random effect, which can be correlated across individuals, and e i is the environmental random effect, which is assumed to be independent of each other and of the genetic effects. Both effects are assumed to follow a Gaussian distribution with variances σ 2 g and 1 − σ 2 g respectively. A person is then assumed to be a case if her liability exceeds a threshold t = Φ −1 (1 − K), i.e. the phenotype y i is given by y i = I{l i > t}. This definition guarantees that the prevalence in the population is indeed K.
Selection probabilities. When the study is observational, the probability of being included in the study is independent of the phenotype. However, in a case-control study, the proportion of cases is usually greatly ascertained. To model this fact, we define a random indicator variable s i indicating whether individual i was selected to the study.
The commonly used "full" or "complete" ascertainment assumption 25 is P(s i = 1 | y i = 1) = 1. While this assumption can be relaxed, as discussed later, it simplifies subsequent analysis.
Suppose the population is of size n and that the expected size of the study is n s < n. The expected number of cases in the study is nK. Additionally, the proportion of cases in the study is P , so: yielding: denote p the probability of a control being included in the study (i.e. p = P(s i = 1 | y i = 0)). The expected number of controls in the study is n(1 − K)p . Additionally, the proportion of controls in the study is (1 − P ) so: Solving for p yields: From here it follows that the probability of being included in the study for a given individual (with unknown phenotype) is: Our results do not depend strictly on the full ascertainment assumption (that is, P(s i = 1 | y i = 1) = 1). The latter assumption can be relaxed such that for any p * , as long as the probability of being selected as a control is multiplied by the same probability. This can model any step prior to the selection procedure, for example the probability that an individual is approached by the health administration to begin with. For example, a nonascertained study (K = P ) might involve only a proportion p * of the population. In this case deriving the selection probability of a given individual yields p * as expected.
Heritability estimation. Next, consider a pair of individuals in the study, whose genetic effects are correlated and denote by ρ the correlation between their genetic effects.
Denote by Z ij the product of the standardized phenotypes: The variable Z ij can take three values: We write down the expected value of Z ij , conditional on the fact that s i = s j = 1 (the individuals are part of the study) and given ρ: We apply Bayes rule to the first of the three expressions on the right: P(y i = y j = 1 | s i = s j = 1; ρ) = P(s i = s j = 1 | y i = y j = 1; ρ)P(y i = y j = 1; ρ) P(s i = s j = 1; ρ) .
Under the full ascertainment assumption P(s i = s j = 1 | y i = y j = 1; ρ) = 1, and so P(y i = y j = 1 | s i = s j = 1; ρ) = P(y i = y j = 1; ρ) P(s i = s j = 1; ρ) . Similarly: and since a control is selected to the study with probability K(1−P ) P (1−K) , this boils down to: For the case of y i = y j , one individual is a case, and is automatically selected, while the other is a control and is selected with probability K(1−P ) P (1−K) . Hence: Using these results we get: Denote the numerator by A(ρ) and the denominator by B(ρ). We wish to approximate the last equation using a Taylor series around ρ = 0. Since the individuals are unrelated, the correlation is expected to be close to 0, and therefore such an approximation is expected to be good. Such an approximation would take the form: See Supp. Mat. for a discussion of a second order approximation.
Note that with ρ = 0, the phenotypes of the two individuals are i.i.d. and so A(0) = 0.
Therefore, the Taylor approximation can be simplified: Similarly, with ρ = 0 the events of being included in the study are i.i.d. for both individuals, so All that remains is to find A (0). Towards that end, we are interested in computing the probabilities of the three possible combinations of phenotypes: −∞ˆ∞ t f ρ,σ 2 g (l 1 , l 2 )dl 1 dl 2 , and P(y i = y j = 0; ρ, σ 2 where f ρ,σ 2 g is the multivariate Gaussian density, namely: with Σ denoting the covariance matrix of the liabilities, given explicitly by: The determinant of Σ is |Σ| = 1 − ρ 2 σ 4 g and its inverse is so the density function f ρ,σ 2 g can be written as: . Differentiating A(ρ) requires differentiating each of the three double integrals w.r.t. ρ: Setting ρ = 0 in the last expression yields: Explanation: we differentiate and set ρ = 0. By the chain rule, the derivative of any expression with ρ 2 is 0 at ρ = 0, and obviously the derivative of any expression which does not depend on ρ is 0. The only expression whose derivative is therefore not 0 at ρ = 0 is −2l 1 l 2 ρσ 2 g in the numerator of the exponent. The denominator of the exponent is 2 at ρ = 0, and so the derivative at ρ = 0 is Using these results we can write down A (0): and so: Hence, when the error of the approximation is small, the slope obtained by regressing Z ij on G ij is an unbiased estimator of P (1−P ) yields an unbiased estimator of σ 2 g -the liability scale heritability.
Extending the liability threshold model to include fixed effects. It is often desired to include fixed effects in the analysis of a complex phenotype. Such fixed effects might include external information such as sex, diet and exposure to environmental risks, but can also be genetic variants with known effects or estimates of population structure such as projections of several top principal components.
Since the liability threshold model is in fact a probit model, fixed effects can be included in the usual manner: where x i is a vector of the values of the relevant covariates and β is a vector of their respective effect sizes.
An individual is a case if l i > t, as before. However, an equivalent formulation would be to subtract the fixed effect from the threshold, rather than adding it to the liability: thus keeping the previous formulation of the liability as a sum of genetic and environmental effects.
Heritability estimation with known fixed effects. Assume first that the fixed effects are known, and so the t i 's are known. The probability of being a case, and of being included in the study, are no longer equal for all the observed individuals. We denote: the probability that the i'th individual is a case, and: the probability that the i'th individual is a case conditional on being selected for the study, where t i is computed using the fixed effects as described above. We now wish to derive the same first order approximation, while accounting for the newly introduced heterogeneity. We redefine: and follow the same steps as before (see supplementary materials for a full derivation). Using the first order Taylor approximation we conclude that the slope obtained by regressing Z ij on is an estimator of heritability on the liability scale.
Estimating heritability with unknown fixed effects. More often than not, the effects of relevant fixed effects are unknown and must be estimated from the data. However, estimating effect sizes under ascertainment in case-control studies is notoriously problematic. Specifically, under the threshold (probit) model, ignoring the ascertainment yields biased estimators.
A special exception is the case of logistic regression. In their seminal paper, Prentice and Pyke (1979) 26 proved that using a logistic regression to estimate fixed effects from ascertained data yields consistent estimators of these effects in the (unascertained) population, and that the ascertainment only biases the intercept.
We therefore suggest a two-step procedure for estimating heritability. First, we estimate the fixed effects using a logistic regression model. We then correct the effect of the ascertainment, and obtain the individual-specific thresholds. Lastly, we plug the thresholds into the estimation scheme described above.
More elaborately, by Bayes' formula: by the complete ascertainment assumption P(s i = 1 | y i = 1, x i ) = 1, and according to the selection scheme: We can thus solve for K i and express it is a function of P i : We then use logistic regression to obtainP i -a consistent estimator of P i , and use this estimate to obtain an estimate of K i , which is in turn used to estimate the threshold: Estimating the added variance due to fixed effects Lastly, the presence of fixed effects increases the variance of the liability, so σ 2 g no longer equals h 2 . The appropriate definition of heritability is now: where σ 2 t is the variance of the thresholds in the population, and so the estimate of σ 2 g can be transformed to an estimate of the heritability simply by dividing it by 1 + σ 2 t . We discuss how σ 2 t can be estimated from the data in the Supplementary material.

Simulations using a generative model
Our goal was to create a realistic simulation, covering a wide range of combinations of disease prevalence K, case sampling probability P and heritability σ 2 g , and for each one creating "natural" genotypes that realistically recreate the complex correlation structure induced by case-control sampling.
Given these parameters, our simulations proceeded as follows: 1. The MAFs of 10,000 SNPs were randomly sampled from U [0.05, 0.5].
2. SNP effect sizes were randomly sampled from N (0, 3. For each individual, we: (a) Randomly generated a genotype using the MAFs, and normalized it (according to Yang et al.'s model 27 ).
(b) Used the genotype and the effect sizes, to compute the genetic effect.
(d) Computed liability and phenotype.
(e) If the phenotype was a case -the individual was automatically included in the study.
Otherwise the individual was included in the study with probability K(1−P ) P (1−K) .

The genotypes of all included individuals were used to compute
This matrix G was then used to estimate heritability for the LMM using the GCT A software 27 and the correction of Lee et al. 2 ; and using GCR as described above.
We note that our choice to work with SNPs that are in linkage equilibrium was motivated by the analysis of Patterson et al. 23 . They show that for the purpose of generating correlation matrices, using SNPs in linkage disequilibrium (LD) is equivalent to using a smaller number of SNPs in linkage equilibrium. They also suggest a method for estimating the effective number of SNPs (i.e. the number of SNPs in linkage equilibrium leading to the same distribution of correlation matrices as a given set of SNPs in LD). Applying their method to WTCCC data suggests that the effective number of SNPs in linkage equilibrium is roughly one tenth of the actual number of SNPs. Hence, using 10,000 SNPs in our simulations is equivalent to simulating roughly 100,000 SNPs with realistic LD structure.

Heuristic correction for the LMM approach
Denote the estimate of the heritability on the liability scale, obtained through the method of Lee et al. 2 byσ 2 g , and the estimates obtained by our method byσ 2 g . Denote the true variance of the genetic effect under ascertainment by: We note a different analytical expression of σ 2 gcc is given in 2 . We validated correctness of our derivation of this expression numerically.
As detailed in the Supplementary Materials, our simulations demonstrated the following properties of these estimators: 1. They are both unbiased when there is no ascertainment.
2. Our estimateσ 2 g remains unbiased in all situations. However, in presence of ascertainment, σ 2 g is biased, and is not linear in the true σ 2 g for fixed K, P .
4. The bias ofσ 2 g worsens as the ascertainment factor K P decreases.
We therefore performed an extensive analysis of the relationship between K P and the bias of the "linearized" estimateσ 2 g σ 2 gcc σ 2 g in our simulations. Our analysis, as detailed in Supplementary Materials, led us to the following relationship betweenσ 2 g and the true underlying heritability σ 2 g : Correcting published results We have derived our heuristic correction using simulations wherein the true underlying heritability is known. However, contrary to previously used corrections, our correction is not a linear transformation of the estimate. When attempting to correct published estimates, we only know K, P andσ 2 g . We define our corrected estimateσ * 2 g to be the value of σ 2 g for which where the expectation is computed using the approximate relationship we derived previously. In other words, the corrected estimate of heritability is the value of heritability for which the expectation of the estimator is the observed estimate, where the expectation is calculated using our heuristic correction.    1 Heritability estimation using genetic correlation regression

Liability threshold model -notations
Denote K the prevalence of a condition in the population and P the prevalence in the study.
Under the liability threshold model, we assume that each individual i has an unknown liability l i = g i + e i where g i is a genetic random effect, which can be correlated across individuals, and e i is the environmental random effect, which is assumed to be independent of each other and of the genetic effects. Both effects are assumed to follow a Gaussian distribution with variances σ 2 g and 1 − σ 2 g respectively. A person is then assumed to be a case if her liability exceeds a threshold t = Φ −1 (1 − K), i.e. the phenotype y i is given by y i = I{l i > t}. This definition guarantees that the prevalence in the population is indeed K.

Selection probabilities
When the study is observational, the probability of being included in the study is independent of the phenotype. However, in a case-control study, the proportion of cases is usually greatly ascertained. To model this fact, we define a random indicator variable s i indicating whether individual i was selected to the study. The commonly used "full ascertainment" [11] assumption is P (s i = 1 | y i = 1) = 1. While this assumption can be relaxed, as discussed later, it simplifies subsequent analysis.
Suppose the population is of size n and that the expected size of the study is n s < n. The expected number of cases in the study is nK. Additionally, the proportion of cases in the study is P , so: nK n s = P, yielding: denote p the probability of a control being included in the study (i.e. p = P (s i = 1 | y i = 0)). The expected number of controls in the study is n(1 − K)p . Additionally, the proportion of controls in the study is (1 − P ) so: Solving for p yields: From here it follows that the probability of being included in the study for a given individual (with unknown phenotype) is:

Relaxing the "full ascertainment" assumption
In fact, the latter results do not depend strictly on the full ascertainment assumption (that is, P (s i = 1 | y i = 1) = 1). The latter assumption can be relaxed such that for any p * , as long as the probability of being selected as a control is multiplied by the same probability. This can model any step prior to the selection procedure, for example the probability that an individual is approached by the health administration to begin with. For example, a non-ascertained study might involve only a proportion p * of the population. In this case deriving the selection probability of a given individual yields p * as expected.

Heritability estimation
Next, consider a pair of individuals in the study, whose genetic effects are correlated and denote by ρ the correlation.
Denote by Z ij the product of the standardized phenotypes: The variable Z ij can obtain three values: We write down the expected value f Z ij , conditional on the fact that s i = s j = 1 (the individuals are part of the study) and given ρ: We apply Bayes rule to the first of the three summands on the right: Under the full ascertainment assumption P(s i = s j = 1 | y i = y j = 1; ρ) = 1, and so P(y i = y j = 1 | s i = s j = 1; ρ) = P(y i = y j = 1; ρ) P(s i = s j = 1; ρ) . Similarly: and since a control is selected to the study with probability K(1−P ) P (1−K) , this boils down to: For the case of y i = y j , one individual is a case, and is automatically selected, while the other is a control and is selected with probability K(1−P ) P (1−K) . Hence: Using these results we get: Denote the numerator by A(ρ) and the denominator by B(ρ). We wish to approximate the latter equation using a Taylor series around ρ = 0. Such an approximation would take the form: We later discuss a second order approximation as well.
Note that with ρ = 0, the phenotypes of the two individuals are i.i.d. and so A(0) = 0. Therefore, the Taylor approximation can be simplified: Similarly, with ρ = 0 the events of being included in the study are i.i.d. for both individuals, so B(0) = K 2 P 2 . All that remains is to find A (0). We are interested in computing the probabilities of the three possible combinations of phenotypes: −∞ˆ∞ t f ρ,σ 2 g (l 1 , l 2 )dl 1 dl 2 , and g is the multivariate Gaussian density, namely: with Σ denoting the covariance matrix of the liabilities, given explicitly by: The determinant of Σ is |Σ| = 1−ρ 2 σ 4 g and its inverse is and so the density function f ρ,σ 2 g can be written as: .
Deriving A(ρ) requires deriving each of the three double integrals w.r.t. ρ: Setting ρ = 0 in the last expression yields: Explanation: we differentiate and set ρ = 0. By the chain rule, the derivative of any expression with ρ 2 is 0 with ρ = 0, and obviously the derivative of any expression which does not depend on ρ is 0. The only expression whose derivative is therefore not 0 at ρ = 0 is −2l 1 l 2 ρσ 2 g in the numerator of the exponent. The denominator of the exponent is 2 at ρ = 0, and so the derivative is l 1 l 2 σ 2 Similarly: Using these results we can write down A (0): and so: Hence, when the error of the approximation is small, the slope obtained by regressing Z ij on G ij is an unbiased estimator of P (1−P ) yields an unbiased estimator of σ 2 g -the liability scale heritability.

Second order approximation
While the first order approximation yields very satisfactory results in our simulations, one can obtain a better estimator using a better approximation. The second term of the Taylor series takes the form: We have already derived: Now, B(ρ) is the probability that two individuals with genetic correlation ρ are included in the study: P (y i = y j = 0; ρ, σ 2 g ).

7
Using the previous derivatives of the double integrals yields: Computing A (0) requires computing the second derivative of the two-dimensional density at ρ = 0 : Hence: and the second order approximation can be written as: Since K and P are assumed to be known, the estimation problem boils down to a single-variable non-linear regression in σ 2 g .
2 Dealing with fixed effects

Extending the liability threshold model
It is often desired to include fixed effects in the analysis of a complex phenotype. Such fixed effects might include external information such as sex, diet and exposure to environmental risks, but can also be genetic variants with known effects or estimates of population structure such as projections of several top principal components.
Since the liability threshold model is in fact a probit model, these effects can be included in the usual manner: where x i is a vector of the values of the relevant covariates and β is a vector of their respective effect sizes.
An individual is a case if l i > t, as before. However, an equivalent formulation would be to substract the fixed effects from the threshold, rather than adding them to the liability: thus keeping the previous formulation of the liability as a sum of genetic and environmental effects.

Heritability estimation with known fixed effects
Assume the fixed effects are known, and so the t i 's are known. We define: and: to be the probability of the i'th individual being a case when accounting for fixed effects through the adjustment of the threshold t i .
We redefine: so now Z ij can obtain four possible values: We now wish to derive the same first order approximation before, while conditioning on the fixed effects.
Repeating the same steps as before while contitioning on t i , t j the expression for A(ρ) is now: Each phenotype was standardized to have mean 0, conditional on the relevant fixed effects. Additionaly, with ρ = 0 the phenotypes are independet and so the expected value of Z at ρ = 0 is 0. An immediate result is that A(0) = 0. To see this consider the Taylor expansion of but on the other hand: Therefore, A(0;ti,tj ) B(0;ti,tj ) = 0 and so A(0; t i , t j ) = 0. We conclude that the first order Taylor approximation is again of the form: To derive an explicit expression for A (0; t i , t j ) we need to differentiate the double integrals as before. When both i and j are cases, with thresholds t i , t j respectively, the double integral takes the form: differentiating it w.r.t. ρ and setting ρ = 0 yields: and differentiating the other double integrals yields similar results: and: We differentiate A(ρ; t i , t j ) w.r.t. ρ, and using the previous results and some algebra we get: As a sanity check, set P i = P j = P and t i = t j = t to get: which is the expression derived for the no fixed effects case.
Moreover, the the probability of inclusion of individual i in the study, conditional on the fixed effects is now: and so: Plugging the derived expressions for A (0; t i , t j ), B(0; t i , t j ), we conclude that regressing Z ij on yields an estimator of heritability on the liability scale.

Dealing with unknown fixed effects
More often than not, the effects of relevant fixed effects are unknown and must be estimated from the data. However, estimating effect sizes under ascertainment in case-control studies is notoriously problematic. Specifically, under the threshold (probit) model, ignoring the ascertainment yields biased estimators.
A special exception is the case of logistic regression. In their seminal paper, Prentice and Pyke (1979) [10] proved that using a logistic regression to estimate fixed effects from ascertained data yields consistent estimators of these effects in the (unascertained) population, and that the ascertainment only biases the intercept.
We therefore suggest a two-step procedure for estimating heritability. First, we estimate the fixed effects using a logistic regression model. We then correct the effect of the ascertainment, and obtain the individual-specific thresholds. Lastly, we plug the thresholds into the estimation scheme described above.
More elaborately, by Bayes' formula: by the complete ascertainment assumption P(s i = 1 | y i = 1, x i ) = 1, and according to the selection scheme: We can thus solve for K i and express it is a function of P i : We then use logistic regression to obtainP i -a consistent estimator of P i , and use this estimate to obtain an estimate of K i , which is in turn used to estimate the threshold:t i = Φ −1 1 −K i , and the estimates of the individual-wise thresholds are used for estimating the liability-scale heritability.

Estimating the added variance due to fixed effects
Lastly, the presence of fixed effects increases the variance of the liability, so σ 2 g no longer equals h 2 . The appropriate definition of heritability is now: where σ 2 t is the variance of the thresholds in the population, and so the estimate of σ 2 g can be transformed to an estimate of the heritability simply by dividing it by 1 + σ 2 t . Therefore, obtaining an estimate of the heritability in the presence of fixed-effects requires an estimate of σ 2 t in the population. To estimate σ 2 t we use the law of total variance: | y)) where y is the phenotype of the individual. Furthermore: and: Once we condition on a specific phenotype, both the expected and variance of T can be estimated from the data, as they are no longer affected by enrichment of cases in the data. Specifically, we use the estimated thresholdst i to estimate the expected value and the variance of the threshold for cases and controls, and plug the estimates into the equations above to obtainσ 2 t . We discuss how σ 2 t can be estimated from the data in the Supplementary material.

Description of the Lee et al. simulation scheme
Lee et al. perform simulations of ascertained case-control studies in the following manner: Phenotypes are simulated in blocks of 100 individuals Given the true genetic variance σ 2 g , liabilities for each are sampled from a multivariate Gaussian distribution with mean 0 and covariance matrix Σ = Gσ 2 g + I(1 − σ 2 g ), where G is given by: that is -individuals within each block have a genetic correlation of 0.05 while individuals in different blocks are perfectly unrelated. Individuals with liabilities higher than the threshold t are considered cases.
Cases are always included in the study. Controls are included in the study with probability K 1−K . Blocks are generated until 100 individuals are accumulated, and this process is repeated 100 times, so in total 10, 000 cases and controls are accumulated.

Discussion of problems with the simulations in Lee et
al.
There are several key aspects in which the simulations of Lee et al. differ from the true generative process of the data.
First, the resulting correlation matrix is highly degenerate, with most of the correlations being 0 -as much as 99.99\%. Since the underlying idea behind heritability estimation from unrelated individuals is leveraging the minor -but non-zero -correlations among unrelated individuals, restricting the simulations to largely degenerate correlation structures is highly unrealistic and counterproductive.
Second, in reality, correlations between individuals span a wide range of values. while in Lee et al.'s simulations the correlations are either 0 or 0.05. Hence negative correlations are impossible, and the expected correlation is larger than 0, contrary to the unrelated individuals assumption.
Thirdly, the selection procedure, in which cases are up-sampled compared to their prevalence in the population, results in cases being more similar genetically than controls. This is a well known phenomenon which has recently attracted considerable attention in the context of GWAS [14]. This is hardly captured by the simulation process of Lee et al. For example, cases from different blocks always have 0 genetic correlation between them in Lee et al.'s simulations, while in a realistic scenario they should have positive genetic correlation, whose magnitude depends on the heritability of the phenotype and the number of causative loci.
Lastly, for small enough prevalence, all blocks in a typical simulation would contain only a single individual. In this case the genetic effects would be completely interchangeable with the environmental effects since G = I and so heritability would be impossible to estimate. In reality, however, a very small prevalence should yield very closely related individuals, thus generating an intuitively easier estimation problem.
We also note that Lee et al.'s simulations use an unrealistic number of individuals -10,000.

Simulations using a generative model
To see if these problems have a major impact, we ran simulations using the full generative model. This was done as follows: 1. The MAFs of 10,000 MAFs were randomly sampled from U [0.05, 0.5].
2. SNP effect sizes were randomly sampled from N (0, 3. For each individual, we: (a) Randomly generated a genotype using the MAFs, and normalized it (according to Yang's model).
(b) Used the genotype and the effect sizes, to compute the genetic effect.
(e) If the phenotype was a case -the individual was automatically included in the study. Otherwise the individual was included in the study with probability K(1−P ) P (1−K) .
5. The genotypes of all included individuals were used to compute G = ZZ m where Z is the matrix of normalized genotypes.
6. We used GCTA [13] to estimate heritability with G as the genetic correlation matrix, and a 0 − 1 vector of the phenotypes as input.
We ran ten repetitions of this simulation for all combinations of σ    We then applied GCR to the same data. The results are portrayed in figure 3. As can clearly be seen from the figure, our estimators are linear in the underlying true heritability. Indeed, as can be seen in ????, applying the correction derived earlier yielded accurate estimators.

Distribution of estimates around the mean
Apart from being unbiased, GCR estimates display an equal or lower variance to the LMM estimates. This can be seen in figure 5.

Weak ascertainment
To test the robustness of our method to the different values of P and K, we ran a wide range of simulations, with less extreme ascertainment. The results are given in figure 6, showing how the downward-bias of the LMM estimates is smaller for less ascertained studies, and disappears completely when the studies are non-ascertained. in both cases the GCR yields unbiased estimates. True h^2 estimated h^2 (corrected) using matrix regression Figure 6: Comparing corrected LMM estimates (top) to GCR estimates (bottom) for studies with low-to-intermediate ascertainment.

non-ascertained simulations
We also applied our method to non-ascertained simulations, where both the LMM and GCR methods produce unbiased estimators, as can be seen in figure 7.

Simulations with population structure
To study the robustness of our method to the presence of subtle population structure, we modified our simulations to introduce population structure similar to the simulations in Lee et al. For every batch of 100 individuals, we randomly selected 5% of the SNPs and fix them for the entire batch. By doing so, the expected correlation within each batch is 0.05, as in Lee  The resulting GCR estimates were unbiased for all these scenarios, as can be seen in figure 8.

Simulations with different numbers of SNPs
The number of SNPs used in the simulations determines the variance of the realized genetic correlations between individuals. Since GCR relies on a first-order apporximation around zero correlation, high variance of genetic correlations might decrease the accuracy of GCR. To study the robustness of GCR to such changes, we reran our initial simulations using 1,000 or 100 SNPs. GCR still produced unbiased estimates. As expected, when the number of SNPs was smallest (100) and the heritability was highest (0.9), the results displayed cosiderably

Second order simulations
To study the accuracy of GCR, we implemented a version of GCR which utilizes the second order approximation as detailed before. We then compared the first-and second-order estimates of heritability on our initial set of simulated data. The resulting second order estimates were almost identical to the results obtained by the first order approximation (figure 10 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q

Simulations with fixed effects
To test that our method still generates accurate estimates of heritability even in the presence of fixed effects, we used the following simulation scenario: In addition to generating the genotypes as described earlier, for each simulated individual we generated a 0-1 "sex" variable with probability 0.5 for each value. We fixed the heritability to 0.5 and the proportion of cases in the study to 0.3, and defined the risk of individuals with sex = 0 to 0.005. The risk of individuals with sex = 1 was either 0.005, 0.01 or 0.02, i.e. the relative risk (RR) was 1, 2 or 4. The thresholds were adjusted to depend on the sex accordingly. Phenotypes were generated after accounting for the fixed effects by way of changing the threshold as described earlier, and individuals were selected for the study based only on their phenotypes (i.e. independently of their sex). We simulated 100 sets of genotypes/phenotypes/sex for each of the three possible RR values. We then estimated the heritability using LMM and GCR, each time with and without including the fixed effect in the analysis. The results are displayed in figure 11. Our simulations demonstrate that GCR correctly accounts for fixed effects while LMM estimates are biased with and without accounting for fixed effects.  Figure 11: Estimation of heritability with and without accounting for fixed effects using LMM and GCR. The true underlying heritability is 0.5 and the only method which is consistently accurate is GCR with fixed effects.

Inference
Standard inference for linear regression is not applicable in our case due to the breaking of several of the basic assumptions used in linear regression inference. For example, the errors are non-normal and not independent. We therefore use the jackknife [5] to estimate the variance of our estimates. Since our method is three to four orders of magnitude faster than LMM estimation, applying such a procedure is not considerably time-demanding. More specifically, we use the jackknife to estimate the standard deviation of the estimate, and use this standard deviation to construct a confidence interval (CI). At each jackknife iteration we remove one individual from the data, and remove all relevant entries to the genetic correlation matrix. We constructed 95% confidence intervals for our simulations and, these CIs covered the true value of the heritability in 98% of the simulations, indicating the jackknife is conservative in this case, and the true standard error of the estimates is usually smaller than the jackknife estimates.

Heuristic correction of LMM estimates
While we have demonstrated that our method can be used to yield accurate estimates of heritability, it might also be of interest to correct existing estimates of heritability which were obtained using LMM and so might be biased.
In their seminal work Dempster and Lerner [3] derive the relationship between additive heritability on the observed scale and heritability on the liability scale under the liability threshold model. However, D&L's work is done under a no-ascertainment assumption. Lee et al. [8] extend their result to ascertained case-control studies, and obtain the following relationship: Note that the relationship depends on the ratio of variances, so that the relationship is not linear. The variance of the genetic effects under ascertainment is given by: We note a different analytical expression of σ 2 gcc is given in Lee et al. We validated our derivation of this expression numerically. A script demonstrating the correctness of our derivation can be found on our webpage.
Despite the fact that this is a relationship between parameters, and we are interested in correcting estimates, we applied this correction to the estimates obtained in our simulations. Our reasoning was twofold -first, we hoped that the biases introduced by the estimation and the biases introduced by using the observed scale would be interchangeable, and thus correcting for the latter would expose the former, making its analysis easier. Second, from ??????????, it seems that LMM works well when heritability is low, regardless of K and P , and also that LMM works well when ascertainment is low. We speculated that this is due to σ 2 gcc being close to σ 2 g in such scenarios. The results of applying this correction to the simulations is shown in figure  12. As can be seen from the figure, the resulting estimates are still biased, except for the case of non-ascertained studies, but the dependency of the corrected estimates on the true underlying heritability seems linear. we speculated that the transformation indeed corrects for the effect of estimating under the observed scale (rather than under the liability scale), but do not correct for the additional effects introduced by the LMM estimation.
The resulting relationships between the corrected estimates and the true underlying heritabilities seemed linear, but with a slope which differs from 1, we speculated that the difference was due to the normality assumptions made by LMM, which are less accurate when ascertainment is strong. We therefore attempted to model the slope as a function of K and P . Since for the noascertained case K = P and the slope is 1, we speculated that the slope depends on K P , and that this size captures the degree of ascertainment. We therefore estimated the slope for each (K, P ) tuple and explored the dependence of the slope on the ratio K P . Using a Box-Cox plot [1] (figure 13) indicated that K P is the best transformation of the K P variable.  Figure 13: A Box-Cox plot using K P as the dependent variable and the slope as the covariate indicates that a square-root transformation is a good fit.
We then used linear regression to estimate the relationship: slope(K, P ) = 1.3 − 0.3 K P We concluded by applying this correction on top of the previous correction (i.e. dividing the estimates by the slope). The results of this correction are given in figure 14. While the results are not as well behaved as when using GCR, it seems that applying these two corrections one on top of the other reduces the bias considerably for all simulated ascertainment scenarios. While after correcting, the results are unbiased, in terms of variance, GCR still outperforms LMM

Applying the correction to published data
So far, we have derived our heuristic correction using simulations wherein the true underlying heritability is known. However, when attempting to correct published estimates, we only know K, P andσ 2 g LM M . We define our corrected estimateσ 2 g * to be the value for which E[σ 2 g LM M ; σ 2 g =σ 2 g * ] =σ 2 g LM M , i.e. the value of underlying heritability for which the expected value of the estimator is the observed estimate, where the expectation is calculated using our heuristic correction.
Confidence intervals are derived by applying the same procedure to the top and bottom limits of a 95% confidence interval based on the published standard deviation of the estimate.

Testing the correction on WTCCC data
To validate our suggested heuristic correction, and to make sure it generalizes beyond our simulations, we computed LMM estimates of heritability for all 7 WTCCC case-control studies, and applied the correction to the resulting estimate. We also applied our method to the data and compared the result to the corrected estimates. The results are highly correlated (ρ = 0.97) and can be seen in figure 16.

Applying our correction to other published results
We applied our correction to several published case-control studies [8,6,7,4]. The corrected heritability estimates are given in  [8,6,7,4] 6 Accounting for imperfect linkage-disequilibrium Yang et al. [12] suggest that the genetic correlation matrix, estimated based on the genotyped data, is in fact not the true underlying correlation matrix, as the latter should be computed only using causal SNPs. Since causal SNPs are unknown, and often not genotyped, the estimated genetic correlation is a noisy and biased estimate of the true genetic correlation. To correct for this effect, Yang et al. suggest using a different correlation matrix where G is the estimated genetic correlation matrix, and β is a constant correcting the bias between the true and estimated genetic correlations. Using extensive simulations using real data and realistic assumptions they estimate β at roughly 0.875. It is easy to show that using G * in LMM estimation yields an estimate of heritability which is 1 βĥ 2 , whereĥ 2 is the heritability estimate. In the case of GCR, we regress the product of normalized phenotypes Z ij on G ij . Hence, plugging in G * ij = βG ij as the covariate has the exact same effect of increasing the estimated heritability by a factor of 1 β .

Quality Control
Following Lee et al. (2011) we applied a stringent quality control (QC) process on the WTCCC data to avoid detecting spurious heritability due to genotyping differences between cases and controls or between the different control groups. We removed SNPs with MAF>5%, with missing rate >1% and SNPs which displayed a significantly different missing rate between cases and controls (p-value<0.05). We also removed SNPs which deviated from Hardy-Weinberg (HW) equilibrium in the control groups (p-value <0.05). Additionally we removed SNPs which displayed a significant difference in frequency between the two control groups. Sex chromosomes were excluded from the analysis. We removed all the individuals appearing in the WTCCC exclusion lists. These include duplicate samples, first or second degree relatives, individuals which are not of European descent and other reasons. In addition we removed individuals with missing rate >1% and all individual pairs with an estimated relationship >0.05 based on the correlation matrix. The last step is done to ensure individuals in the study are not closely related.

No evidence of population structure
Lee et al. [8] include several top principal component vectors as fixed effects in an attempt make sure heritability estimates are not inflated due to population structure in the WTCCC samples. We suggest that this procedure is not required for two reasons. First, we note that [2] perform a thorough investigation of population structure in the WTCCC samples, and conclude that once individuals of different ancestry are removed from the samples, the effect of population structure is very minute. Second, [9] present a straightforward statistical test to decide how many principal components should be included as covariates in an association study. However, their method assumes that the individuals are randomly sampled from the population. In our case, cases are over-sampled in the study, and so LD regions hosting a causative variant might be mistaken for population structure. One simple way to overcome this problem is to dilute the number of SNPs, making sure each LD block is represented by a small number of SNPs. It is then unlikely that the differences between the case and control groups due to caustive variants would be captured by a top ranking principal component.
We applied the method of [9] do the GWAS data using every k'th SNP for increasing values of k. If there's a real population structure, we expect the number of significant principal components to be non-zero asymptotically. However, this was not the case, in agreement with [2] (data not shown).
We note that including principal components which tag such regions hosting causative variants would result in an underestimation of the heritability, as such a vector would be highly correlated with the phenotype. It is therefore highly undesired to include principal components which do not tag actual population structure.


extfieldsofstudy: Medicine Biology

text:
Enterococcus faecium Biofilm Formation: Identification of Major Autolysin AtlAEfm, Associated Acm Surface Localization, and AtlAEfm-Independent Extracellular DNA Release

ABSTRACT Enterococcus faecium is an important multidrug-resistant nosocomial pathogen causing biofilm-mediated infections in patients with medical devices. Insight into E. faecium biofilm pathogenesis is pivotal for the development of new strategies to prevent and treat these infections. In several bacteria, a major autolysin is essential for extracellular DNA (eDNA) release in the biofilm matrix, contributing to biofilm attachment and stability. In this study, we identified and functionally characterized the major autolysin of E. faecium E1162 by a bioinformatic genome screen followed by insertional gene disruption of six putative autolysin genes. Insertional inactivation of locus tag EfmE1162_2692 resulted in resistance to lysis, reduced eDNA release, deficient cell attachment, decreased biofilm, decreased cell wall hydrolysis, and significant chaining compared to that of the wild type. Therefore, locus tag EfmE1162_2692 was considered the major autolysin in E. faecium and renamed atlAEfm. In addition, AtlAEfm was implicated in cell surface exposure of Acm, a virulence factor in E. faecium, and thereby facilitates binding to collagen types I and IV. This is a novel feature of enterococcal autolysins not described previously. Furthermore, we identified (and localized) autolysin-independent DNA release in E. faecium that contributes to cell-cell interactions in the atlAEfm mutant and is important for cell separation. In conclusion, AtlAEfm is the major autolysin in E. faecium and contributes to biofilm stability and Acm localization, making AtlAEfm a promising target for treatment of E. faecium biofilm-mediated infections.

calization of Acm, a unique adhesin described in E. faecium to bind to collagen types I and IV. Finally, we show the presence of eDNA at the cell septum of the wild-type strain and between the chains of the AtlA Efm -deficient mutant. DNase treatment resulted in cell separation, a finding which suggests that eDNA contributes to the chaining phenotype in the AtlA Efm -deficient mutant.

RESULTS
Identification of putative E. faecium autolysins. Six putative autolysin-encoding genes were selected by bioinformatic analysis, and single-crossover mutants were constructed. Locus tag EfmE1162_2692 has the highest similarity and sequence coverage compared to AtlA in E. faecalis (51% amino acid similarity and 80% sequence coverage) (see Table S1 in the supplemental material). Growth curves of these six mutants in brain heart infusion (BHI) broth showed no difference compared to the wild type; however, in tryptic soy broth (TSB) medium with 1% glucose (TSBg), all strains grew to a lower optical density at 600 nm (OD 600 ) than with growth in BHI medium ( Fig. S1A and B). Mutant 0587::pWS5 showed earlier log phase in TSBg than did other strains.
The putative autolysin locus tag EfmE1162_2692 is resistant to Triton X-100, exerts peptidoglycan hydrolysis activity, and is involved in cell separation and cell adhesion. Several functional assays were performed to characterize the autolysis activity of the six predicted autolysins. In early stationary phase, bacterial cells are known to be more vulnerable to cell lysis when exposed to Triton X-100 while major autolysin mutants show less lysis (18). Triton X-100 lysis assays of E1162 wild type and the six putative autolysin mutants revealed that with the exception of mutant 2692::pWS9, the OD 600 of the wild type and the mutants was decreased by 60%. A decrease of only 20% in the OD 600 was detected in mutant 2692::pWS9, confirming a resistance to lysis (Fig. 1). Therefore, we renamed locus tag EfmE1162_2692 to atlA Efm .
Although the observed phenotype of 2692::pWS9 is similar to that of other major autolysins, we attempted to construct a markerless mutant of atlA Efm (see Materials and Methods), but unfortunately, this was unsuccessful. As an alternative, we selected from a mariner transposon library in E1162 a mutant (atlA Efm ::genta) that contained an insertion of a gentamicin resistance gene at position 278 relative to the start of the gene. This insertion of a single resistance gene is considerably smaller than pWS9 and even less likely to result in polar effects on adjacent genes that encode a type 2 phosphatidic acid phosphatase family protein (upstream) and a leucyl-tRNA synthetase (downstream). However, to confirm that all observed phenotypes were attributed only to insertional inactivation of atlA Efm , we complemented atlA Efm ::genta in trans, resulting in atlA Efm ::gentaϩatlA Efm .
Growth curves of atlA Efm ::genta in BHI medium showed no difference compared to the wild type, while in TSBg, atlA Efm :: genta reached a higher OD in stationary phase (see Fig. S1A and B in the supplemental material). As shown in atlA Efm ::pWS9, at-lA Efm ::genta is also resistant to lysis when analyzed in the Triton X-100 assay (Fig. 1).
AtlA Efm contains six LysM domains which are conserved in other major autolysins and have been shown to be able to bind and hydrolyze peptidoglycan from the cell wall, causing cell death (19). To confirm hydrolysis of peptidoglycan by AtlA Efm , supernatants of the wild type, all predicted autolysin mutants, and the complemented strain atlA Efm ::gentaϩatlA Efm were subjected to zymography. In Fig. 2, three autolytic bands of 75, 50, and 37 kDa were observed. Mutants 2692::pWS9 and atlA Efm ::genta lacked the 75-kDa band compared to the wild type and the complemented strain, a band which is the expected size of the major autolysin AtlA Efm . The second band of approximately 50 kDa was missing from the mutant 1800::pWS6. Fifty kilodaltons is the expected size of locus tag EfmE1162_1800 (46.8 kDa), indicating that this gene possibly encodes a secondary autolysin. None of the screened mutants lacked the third autolytic band (37 kDa).
By hydrolyzing peptidoglycan, autolysins can play a role in cell separation and growth. Chaining of bacterial cells, a result of impaired cell partitioning, has been observed in other species when the major autolysin gene is inactivated. In overnight cultures, long chains of atlA Efm ::genta were observed, implying the cell separation defect (see Fig. S2A to H in the supplemental material). No chaining was observed in the other predicted autolysin mutants or in the complemented atlA Efm mutant.
Since in other bacterial species, autolysins are also implicated in early attachment of cells, we subsequently studied early adhesion on polystyrene of the wild type and the six putative autolysin mutants. Only atlA Efm ::genta had substantially fewer cells attached predicted autolysin mutants, and the complemented AtlA Efm -deficient strain were incubated in the Bioscreen C system at 37°C with continuous shaking, and absorbance at 600 nm (OD 600 ) was recorded every 15 min for 21 h. On the x axis, time is indicated in hours (h); on the y axis is the percentage of surviving cells in suspension. (33%) than the wild type ( Fig. 3), while attachment was restored in the complemented strain atlA Efm ::gentaϩatlA Efm .
In conclusion, atlA Efm ::genta not only was less vulnerable to Triton X-100 and less capable of initial adherence but also reduced peptidoglycan hydrolysis activity and exhibited long chaining, supporting that this gene serves as the actual major autolysin in E. faecium.
AtlA Efm is a prerequisite for normal biofilm development in E. faecium. We studied the effect of autolysin mutants on biofilm formation in three different models: polystyrene, semistatic, and flow cell models. In the polystyrene model, all predicted autolysin mutants were tested, and only atlA Efm ::genta formed 80% less biofilm than the wild type while biofilm formation was restored in the complemented mutant (see Fig. S3 in the supplemental material).
Due to this result, only atlA Efm ::genta was tested in the other assays. In the semistatic model, biofilm formation of atlA Efm :: genta was considerably reduced compared to that of the wild type. Average thickness was reduced by 80% while biomass was reduced by 90% relative to that of the wild type. This phenotype was complemented in atlA Efm ::gentaϩatlA Efm (Fig. 4A to C).
In a flow cell model, biofilm formation of the AtlA Efm mutant was also attenuated. In this model, wild-type cells firmly attached to the slide, formed small clumps, and developed into thick biofilm covering the whole slide within a few hours. In contrast, atlA Efm ::genta was impaired in initial attachment and subsequently formed only small and thin clumps, which were not stable enough to develop to mature biofilm at 17 h of growth (see Movie S1A in the supplemental material). Again, the wild-type  atlA Efm ::gentaϩatlA Efm (C) was grown for 24 h on poly-L-lysine-coated glass in TSBg at 120 rpm at 37°C. Cells were stained with propidium iodide (red). The average thickness and biomass of biofilms were measured at five random positions and analyzed with Comstat/Matlab software. The average thickness was 1.86 m for the wild type, 0.14 m for atlA Efm ::genta, and 0.84 m for atlA Efm ::gentaϩatlA Efm . The mean biomass was 0.73 m 3 /m 2 for the wild type, 0.03 m 3 /m 2 for atlA Efm ::genta, and 0.45 m 3 /m 2 for atlA Efm ::gentaϩatlA Efm . The difference between the wild type and the mutant was significant (P Ͻ 0.005). Pictures were taken at ϫ63 magnification with a 2.5 optical zoom. In the flow cell model, biofilm was grown for 17 h in a Stovall flow cell system in 1:10 TSB with 1% glucose (0.13 ml/min) at 37°C. Pictures were taken at ϫ40 magnification with a 2.5 optical zoom for the wild type (D), a 3.5 optical zoom for atlA Efm ::genta (E), and a 2.5 optical zoom for atlA Efm ::gentaϩatlA Efm (F). Cells were stained with Syto 9 (green) and propidium iodide (red). phenotype was restored when the mutant was complemented ( Fig. 4D to F).
AtlA Efm is cell surface localized at the septum and is secreted in biofilm. In other bacterial species, the major autolysin has been identified both in the supernatant of planktonic cells and localized to the bacterial surface (11,20). We used AtlA Efm -specific antibodies to detect cell surface expression of AtlA Efm in planktonic cells. Flow cytometry analysis demonstrated that AtlA Efm was detected at the wild-type cell surface, with higher expression in log phase than stationary phase (data not shown). AtlA Efm antibodies were specific, since in all of the following experiments, incubations with secondary antibodies alone yielded no signal.
In addition, confocal microscopy showed that AtlA Efm was also localized at the septum of cells in log phase, supporting a role in cell separation, as already concluded from the chaining phenotype in the mutant (Fig. 5A). As expected, no AtlA Efm signal was observed from the mutant (Fig. 5B). In the complemented strain ( Fig. 5C), the same septum localization as that identified in the wild type was observed. In 24-h biofilms, AtlA Efm was not specifically localized at the septum of the cells but was identified between cells in the matrix and randomly around a subset of cells. In the mutant, again, no AtlA Efm was detected, while the complemented strain exhibited the same localization of AtlA Efm as the wild type ( Fig. 5D to F).
AtlA Efm contributes to eDNA release and eDNA-dependent biofilm formation. By exerting lysis activity, autolysins can contribute to eDNA release in planktonic cells and also in biofilm, thereby contributing to biofilm formation and stability. The presence of eDNA in E. faecium wild type, atlA Efm ::genta, and atlA Efm :: gentaϩatlA Efm was analyzed by isolation of DNA present in the supernatant of overnight cultures. As shown in Fig. S4 in the supplemental material, supernatant of atlA Efm ::genta contained less DNA than the wild type and the complemented strain. To further investigate the role of eDNA in biofilm, we studied the effect of DNase in a semistatic biofilm model. Incubation with DNase at a t of 0 reduced biofilm biomass and average biofilm thickness by 94% and 91%, respectively, in the wild type and atlA Efm :: gentaϩatlA Efm after 24 h of biofilm growth (Fig. 6A, B, E, and F). These results show that eDNA is important in biofilm formation. Unexpectedly, incubation of the mutant with DNase leads to separation of bacterial chains into single cells and increased cell attachment but not clump and biofilm formation ( Fig. 6C and D). Similar to the biofilm semistatic model, chains were separated within 5 min during planktonic growth in an eight-well chamber slide when DNase was added to atlA Efm ::genta (Movie S1B).
To further investigate how DNase treatment results in chain separation, we incubated 4-h and 24-h biofilms and overnightcultured planktonic cells with double-stranded DNA (dsDNA) antibodies. Using confocal microscopy, we were able to detect DNA at the cell septum in wild-type 4-h early biofilm ( Fig. 7A and B) and between the chains of overnight planktonic growth of at-lA Efm ::genta ( Fig. 7C and D), results which concur with our observations of the effect of DNase on cell separation in the chained atlA Efm ::genta and suggest that eDNA contributes to the chaining phenotype in the AtlA Efm -deficient mutant. Immunoelectron microscopy clearly showed the chaining phenotype in the AtlA Efmdeficient mutant, confirmed by ␣-dsDNA gold incubation labeling of eDNA at the cell septum of atlA Efm ::genta chains from overnight cultures ( Fig. 7E and F). In both methods, incubation with only secondary antibodies of the same isotype did not give specific signals (see Fig. S5E and F in the supplemental material). In 4-h and 24-h atlA Efm ::genta biofilms, DNA was also observed between the cells ( Fig. S5B and C).
The effect of exogenous AtlA Efm on phenotypes of wild-type and mutant strains was studied using confocal microscopy. For this, the wild type, atlA Efm ::genta, and atlA Efm ::gentaϩatlA Efm were incubated with or without AtlA Efm in BHI in an eight-well chamber. When AtlA Efm was added, as with DNase, chains started to be released into separate cells in the AtlA Efm -deficient mutant, but unlike DNase, clumping of cells occurred in all strains. It is important to mention that this phenomenon was observed in BHI, a medium in which E. faecium forms less-established biofilms (see Movie S1C in the supplemental material).
These results indicate that AtlA Efm plays different roles in biofilm, where it contributes to biofilm stability via eDNA release, and during planktonic growth, where it may induce cell clumping.
AtlA Efm influences surface protein localization and binding to collagen types I and IV. Some autolysins are known to act as adhesins and may bind extracellular matrix (ECM) proteins (microbial surface components recognizing adhesive matrix molecules [MSCRAMMs]). Therefore, we verified the role of AtlA Efm in binding to ECM proteins using a whole-cell enzyme-linked immunosorbent assay (ELISA). When the wild type, atlA Efm :: genta, and atlA Efm ::gentaϩatlA Efm were incubated with ECM proteins, the AtlA Efm -deficient mutant, unlike the wild type and the complemented strain, specifically did not bind to collagen types I and IV (Fig. 8). To confirm whether AtlA Efm is responsible for this phenotype, an ELISA with the purified protein was conducted; however, no binding of purified AtlA Efm to collagen types I and IV was observed (data not shown).
Since AtlA Efm seemed not directly involved in collagen binding, we hypothesized that AtlA Efm contributed to cell surface expression of another collagen-binding surface protein. Nallapareddy et al. (21) described Acm, a collagen-binding adhesin, biofilm, cultured without (A) and with (B) DNase, atlA Efm ::genta biofilm cultured without (C) and with (D) DNase, and atlA Efm ::gentaϩatlA Efm biofilm cultured without (E) and with (F) DNase. Twenty-four-hour biofilm was grown in a semistatic model in TSBg at 120 rpm at 37°C. Cells were stained with propidium iodide. When DNase was added, the average thickness was 0.068 m for the wild type, 0.099 m for atlA Efm ::genta, and 0.10 m for atlA Efm ::gentaϩatlA Efm . The mean biomass was 0.019 m 3 /m 2 for the wild type, 0.037 m 3 /m 2 for atlA Efm ::genta, and 0.035 m 3 /m 2 for atlA Efm :: gentaϩatlA Efm . The difference between biofilms with and without DNase was significant for the wild type and the complemented AtlA Efm -deficient strain (P Ͻ 0.05). Pictures were taken at ϫ63 magnification with a 2.5 optical zoom. . Cells were incubated with ␣-dsDNA and goat ␣-mouse Alexa 488 (green). The bacterial membrane was stained with FM 95-5 (red). Transmission electron microscopy of atlA Efm :: genta overnight cultures (E and F). Cells were incubated with ␣-dsDNA and rabbit ␣-mouse of the same isotype. Immunogold particles are seen linked to wire-like structures between the cells as well as to polymeric structures outside cells, which likely represents released DNA.
which, so far, is the only E. faecium surface protein that specifically binds to collagen types I and IV.
Confocal microscopy using ␣-Acm antibodies revealed that Acm was localized at the septum and at the poles in the wild type and the complemented strain, while in atlA Efm ::genta, Acm was localized only at the septum ( Fig. 9A to C). Control antibodies from the same isotype gave no labeling (see Fig. S6B in the supplemental material). Absence of exposed Acm at the poles of the AtlA Efm -deficient mutant most probably resulted from the defect in cell separation in this strain and may explain the reduced binding to collagen types I and IV. Specificity of ␣-Acm antibodies was confirmed, since flow cytometry revealed, unlike with the E1162 wild type, no Acm expression in strain E980, a commensal isolated with a transposon insertion in acm (Fig. S6A).
To confirm the later hypothesis, the AtlA Efm -deficient mutant was tested for collagen type I and IV binding in a whole-cell ELISA in the presence of DNase, which, as demonstrated above, resolves the long chains in this mutant. As expected, incubation of the atlA Efm ::genta mutant with DNase restored binding to collagen types I and IV (Fig. 8).

DISCUSSION
In this study, we identified and characterized the major autolysin in E. faecium, AtlA Efm , and demonstrated that AtlA Efm , as in other species, is directly involved in biofilm formation, cell separation, cell wall hydrolysis, and eDNA release (22). Furthermore, the AtlA Efm -deficient mutant displayed attenuated cell surface localization of Acm.
Based on sequence similarity, AtlA Efm is composed of three domains, a central catalytic domain, responsible for glucosaminidase activity (23), a C-terminal domain, composed of LysM modules that afford peptidoglycan affinity and possibly target autolysins to the division septum and poles (24), and a T/E-rich N-terminal domain, the function of which is still unknown (23). Like the major autolysin of E. faecalis (14), AtlA Efm contains a lysozyme-like superfamily domain and six LysM domains which are not present in any of the other putative autolysins analyzed. Lysozyme-like superfamily domains are known to be involved in the hydrolysis of beta-1,4-linked polysaccharides (25). Moreover, atlA Efm is present in all sequenced E. faecium strains, suggesting a highly conserved and essential function.
Five types of experiments indicated that AtlA Efm is the major autolysin of E. faecium. First, the AtlA Efm -deficient mutant dem-onstrated considerably less cell lysis in the presence of Triton X-100 than did the wild type. Second, in the zymography, a lytic band of 75 kDa, which is identical to the predicted molecular mass of AtlA Efm , was identified as AtlA Efm . Third, the AtlA Efm -deficient mutant showed decreased early adhesion compared to that of the wild type. Fourth, devious cell separation by chain formation was observed only in the AtlA Efm -deficient mutant, and finally, the AtlA Efm -deficient mutant had 80% less biofilm production than the wild type in a polystyrene assay. All phenotypes were able to be restored by complementation in trans. These five lines of evidence demonstrate that AtlA Efm plays a major role in E. faecium E1162 autolysis and should be considered the major autolysin. In addition to AtlA Efm , we identified a second autolysin (locus tag EfmE1162_1800) by zymography. In E. faecalis, a second autolysin, designated AtlB, has lytic properties and can act as a surrogate of AtlA but is less efficient in cell separation, probably due to the presence of fewer LysM domains (26). Also in E. faecium, locus tag EfmE1162_1800 has only two LysM domains, compared to the six domains in AtlA Efm . The exact function of this second autolysin remains to be investigated.
Confocal microscopic observations revealed that the AtlA Efmdeficient mutant formed long chains of cells compared to the wild type in planktonic growth and biofilm. The formation of long chains of cells upon inactivation of the major autolysin has also been observed previously (27,28) and implies that autolysin plays a role in cell separation. Our findings that AtlA Efm was localized at the septum of cells in log phase and that addition of exogenous AtlA Efm resolved chained cells into single cells support a role of AtlA Efm in cell separation. Schlag et al. (19) elucidated the mechanism underpinning localized peptidoglycan hydrolysis in staphylococci. Although wall teichoic acid prevents S. aureus Atl from binding to the cell wall, binding of Atl through its amidase domain at the cell septum is increased, presumably due to a lower wall teichoic acid concentration (19). This mechanism explains why autolysin-derived peptidoglycan hydrolases are active primarily at the site of cell separation (29,30). In E. faecalis, cell chaining in an AtlA mutant disappeared after 24 h of culture, indicating that eventually AtlB can cleave the cell septum in the absence of AtlA (31). However, this does not apply to E. faecium, since chaining was also observed in overnight-cultured atlA Efm ::genta and in all other conditions tested.
An unexpected finding was that addition of DNase to the AtlA Efm -deficient mutant also resulted in chains resolving into single cells. This means that for the first time we have demonstrated that DNA plays a role in cell-cell interactions in E. faecium. In E. faecalis, DNA can be secreted by metabolically active cells in the absence of lysis, and this actively secreted eDNA is present at the cell septum and important in early biofilm formation (32). As in E. faecalis, we were able to demonstrate that in wild-type E. faecium and in the AtlA Efm -deficient mutant, eDNA is present at the cell septum. Specifically, the fact that this was also seen in the AtlA Efm -deficient mutant supports the hypothesis that also in E. faecium, eDNA was able to be released in the absence of a major autolysin and is involved in cell separation. Together, these results indicate that both peptidoglycan and eDNA are involved in cellcell interactions upon cell division that allow both AtlA Efm and DNase to separate the bacterial chains. The exact organization of peptidoglycan and eDNA at the cell septum and the possible interaction between these two components are currently not known and remain to be elucidated.
We also demonstrated the presence of AtlA Efm in the biofilm matrix. As in E. faecalis, AtlA Efm might be involved in autolysis of a subpopulation in the biofilm by a fratricidal mechanism (14). In E. faecalis, this is mediated by the protease GelE and coordinated by a two-component system. Such a regulatory mechanism of AtlA Efm expression is currently unknown in E. faecium.
Incubation of a single layer of AtlA Efm -deficient mutant cells with purified AtlA Efm resulted in clump formation, suggesting that secretion of AtlA Efm contributes to this process. This indicates that AtlA Efm may have different functions during planktonic and biofilm lifestyles of E. faecium.
Biofilm formation is a multistep process that requires bacterial adhesion to a substrate followed by cell-cell interactions, resulting in multiple layers of biofilm. Apart from impaired initial adherence of the AtlA Efm -deficient mutant, which may be attributed to the chaining phenotype (33)(34)(35), the biofilm deficiency also seemed to result from reduced stability of the matrix in biofilm clumps. This is most probably due to reduced amounts of eDNA but might also be due to a lack of cohesive activity of AtlA Efm when secreted in biofilm.
Apart from biofilm formation, cell attachment through binding to ECM proteins is an important step in colonization. We observed reduced binding of the AtlA Efm -deficient mutant to collagen types I and IV, most likely due to altered localization of Acm. Peptidoglycan hydrolases, in general, have been described to be involved in surface protein assembly, specifically proteins related to secretion systems in Gram-negative bacteria (36). However, this is the first documentation of the role of a major autolysin in protein surface localization in Gram-positive bacteria. This observation implies an essential contribution of AtlA Efm to E. faecium colonization, since Acm seems to be the major collagen-binding surface protein. Collagen type I is the most common collagen in the human body and is found in all tissues, while collagen type IV is present exclusively in the basal lamina (37), for instance in the gut, which is an important niche of E. faecium. The results presented here provide insights into how AtlA Efm , a newly identified major autolysin in E. faecium, plays an important role in cell lysis and separation, biofilm formation, cell adhesion, and Acm localization and thus may contribute to the pathogenesis of E. faecium infections. Expression, distribution, and activation of the autolysin are probably tightly controlled in E. faecium, and these mechanisms need to be elucidated.
Nosocomial E. faecium infections are caused by multiresistant strains, which tend to even become panresistant. In addition, current drugs fail to eradicate biofilms and new drugs are barely needed. Although animal models did not show a reduction in virulence of AtlA-deficient mutants in E. faecalis (38), we expect AtlA Efm to behave differently. First, genetic inactivation of the major autolysin AtlA Efm affects not only biofilm formation but also binding to the ECM collagen types I and IV. Second, after prolonged growth, no apparent replacement of AtlA Efm function by the secondary autolysin, as in E. faecalis, occurs. Therefore, AtlA Efm may serve as a promising candidate for development of new targeted strategies in the treatment of E. faecium biofilmrelated infections (39).

MATERIALS AND METHODS
Identification of putative E. faecium autolysins. A bioinformatic screen was performed to identify putative autolysins in the E. faecium E1162 genome (40). Using BLASTp (http://www.ncbi.nlm.nih.gov/), six genes were selected for further characterization based on amino acid similarity to atlA from E. faecalis V583 (see Table S1 in the supplemental material).
Bacterial strains, plasmids, growth conditions, and determination of growth curves. E. faecium and Escherichia coli strains used in this study are listed in Table S2 in the supplemental material. The ampicillinresistant E. faecium strain E1162, recently sequenced in our lab, was used throughout this study (40). Unless otherwise mentioned, E. faecium was grown in BHI at 37°C. For biofilm assays, TSBg was used. E. coli DH5␣, BL21 (Invitrogen), and EC1000 (41) were grown in Luria-Bertani (LB) medium. Where necessary, antibiotics (Sigma-Aldrich, Saint Louis, MO) were used at the following concentrations: gentamicin, 300 g/ml (E. faecium) and 25 g/ml (E. coli); spectinomycin, 300 g/ml (E. faecium) and 100 g/ml (E. coli); erythromycin, 50 g/ml with added lincomycin at 50 g/ml (E. faecium) and 150 g/ml (E. coli); and ampicillin, 150 g/ml (E. coli). Determination of growth curves in triplicate was described by Zhang et al. (42) using a Bioscreen C instrument (Oy Growth Curves AB).
Gene disruption and mutant complementation. Predicted autolysin genes were disrupted and complemented as described before (43). Genespecific primers and modified pWS3 vectors are mentioned in Table S3 in the supplemental material, resulting in single-crossover mutant strains in Table S2. To confirm the phenotype of 2692::pWS9, we tried to construct a markerless mutant using the approach previously described by us (42); however, complete or partial markerless deletion of atlA Efm was unsuccessful. Therefore, we screened 400 clones from a mariner transposon mutant library for reduced biofilm formation in polystyrene plates and selected a mutant that carried the gentamicin resistance gene aac(6=)aph(2Љ) as an insertion in atlA Efm , atlA Efm ::genta (42). Finally, atlA Efm :: genta was complemented using a modified pMSP3535 vector (44) (Table S2). For that, the entire atlA Efm gene was amplified using primers atlA Efm -comp-BamHI-F and atlA Efm -comp-XhoI-R (Table S3) and cloned in pMSP3535, resulting in pMSP3536. pMSP3536 was introduced in atlA Efm ::genta by electroporation, resulting in atlA Efm ::gentaϩatlA Efm .
Autolysis assay. To monitor bacterial lysis, overnight cultures of wildtype E. faecium, all single-crossover mutants, atlA Efm ::genta, and the in trans-complemented strain atlA Efm ::gentaϩatlA Efm were inoculated at an initial OD 660 of 0.05 into 5 ml BHI and incubated at 37°C. Cells at an OD 660 of 0.5 were washed with phosphate-buffered saline (PBS) and resuspended in 0.02% Triton X-100. Cell suspensions were incubated in the Bioscreen C system at 37°C with continuous shaking, and absorbance at 600 nm was recorded every 15 min for 21 h. Each experiment was performed in triplicate.
Zymography. A zymographic analysis of cell wall-associated hydrolases was conducted as described by Liu et al. (45) and Ju et al. (46) with some modifications. Briefly, extracellular murein hydrolases of the wild type, all mutants, and the complemented strain atlA Efm ::gentaϩatlA Efm were collected from 50 ml BHI medium at an OD 660 of 0.3. To prepare substrates for the zymogram, 300 ml of wild-type overnight BHI culture was used to extract the cell wall from the bacteria. Lytic activity as a result of cell wall hydrolysis appeared as a white band on a blue background.
Initial polystyrene adherence assay. The initial adherence assay was performed in triplicate as described previously with minor modifications (47). E. faecium strains cultured overnight in BHI were diluted up to an OD 660 of 0.05 and grown until an OD 660 of 0.3. Two-hundred-microliter bacterial suspensions were added in triplicate to the polystyrene plate. Bacteria were stained with crystal violet and eluted with 96% ethanol.
Biofilm polystyrene assay. The biofilm polystyrene assay was performed in triplicate as described previously with some modifications (47). In brief, plate-grown bacteria were resuspended in TSBg and the OD 660 was measured. Cell suspensions were diluted to an OD 660 of 0.01 and incubated for 24 h. The plates were washed and stained as described above.
Biofilm semistatic model and CLSM. For the biofilm semistatic model, overnight bacterial cultures were diluted to an OD 660 of 0.01 in 6 ml TSBg and added to a coverslip coated with poly-L-lysine (0.45 m; diameter, 12 mm; Becton Dickinson) inside a well from a six-well polystyrene plate (Corning Inc.). When applied, 1.5 g/l of DNase (Roche) was added at the start of the biofilm (t ϭ 0). Biofilms were grown at 37°C for 24 h at 120 rpm. After 24 h, the coverslips were washed once with 0.85% NaCl, and the biofilms were chemically fixed in 8% glutaraldehyde (Merck) for 20 min and washed again with 0.85% NaCl. The biofilms were stained with 15 g/ml propidium iodide (PI) in 0.85% NaCl for 15 min. After incubation, the stain was removed and coverslips were transferred to glass microscope slides and analyzed by a confocal laser scanning microscope (CLSM) (Leica SP5), equipped with an oil plan-Neofluar ϫ63/1.4 objective. PI was excited at 633 nm. Z stacks were taken with an interval of 0.42 m. Pictures were analyzed with LAS AF software (Leica), and biofilm was quantified using Comstat (48)/Matlab R2010b software (the MathWorks). The average thickness and biomass of the biofilms were measured at five randomly chosen positions. Flow cell biofilm model. Dynamic biofilms were studied in a Stovall flow cell system (Life Science, Inc., Greensboro, NC). Twenty-four-hour TSB cultures were diluted until an OD 660 of 0.01 and inoculated in the chambers. Biofilms were grown in TSB diluted in PBS (1:10 [vol/vol]) with 1% glucose under a flow of 0.13 ml/min during 17 h. Biofilm development was scanned at regular intervals of 7 min (ϫ40 objective) with a DFC360 FX digital camera kit SP5 (Leica). After 17 h, the flow was increased to 0.5 ml/min to wash away loose cells. Final biofilms were stained with live/dead stain (BAC light kit; Invitrogen), and images were acquired using a confocal microscope (Leica SP5). eDNA isolation. Total eDNA from the E1162 wild type, atlA Efm :: genta, and atlA Efm ::gentaϩatlA Efm was isolated from BHI overnight cultures with isopropanol and ethanol precipitation as described previously (14).
Cloning, expression, and purification of recombinant AtlA Efm . The E1162 atlA Efm gene was amplified with atlA Efm -exp-BamHI-F and atlA Efm -exp-NotI-R primers (see Table S3 in the supplemental material) and cloned in frame into the expression plasmid (Table S2), expressed and purified as previously described by Bardoel et al. (49).
Flow cytometry. AtlA Efm and Acm were detected at the surface of E. faecium E1162, atlA Efm ::genta, atlA Efm ::gentaϩatlA Efm , and E. faecium E980 grown in BHI until an OD 660 of 0.3 or stationary phase by flow cytometry as described previously (47). ␣-AtlA Efm and ␣-Acm rabbit immune serum (1:100) was used to detect AtlA Efm and Acm.
AtlA Efm , Acm, and eDNA localization. AtlA Efm , Acm, and eDNA localization in E. faecium was performed on single cells and in the semistatic biofilm model. For single cells, 5 l of bacteria grown to an OD 660 of 0.3 or in overnight culture was placed on a coverslip coated with poly-L-lysine and dried at 55°C for 20 min. For both conditions, cells were washed with 3 ml PBS, fixed with 3% paraformaldehyde for 15 min, and washed again with PBS. To identify AtlA Efm , Acm, and eDNA, cells were incubated with ␣-AtlA Efm or ␣-Acm IgG antibody (1:500 in PBS with 1% bovine serum albumin [BSA]) or monoclonal ␣-dsDNA (Abcam, Cambridge, MA) (1: 1,000 in PBS with 1% BSA) for 1 h on ice and subsequently washed with PBS. A secondary IgG Alexa Fluor 488 goat ␣-rabbit or ␣-mouse antibody (Life Technology) (1:500 in PBS with 1% BSA) was added and incubated for an additional 1 h on ice to detect binding to first antibodies. Cells were washed once more and incubated with 5 g/ml FM 5-95 dye (Invitrogen) for 2 min on ice. Fluorescence was analyzed in the confocal microscope (Leica SP5). Alexa 488 and FM 5-95 were excited at 488 nm. As a control, bacteria were treated as described above but omitting the first antibodies (␣-AtlA Efm , ␣-Acm, and ␣-dsDNA).
eDNA localization in the AtlA Efm -deficient mutant by transmission immunoelectron microscopy. Transmission immunoelectron microscopy was performed as described previously with some modifications (50). eDNA in atlA Efm ::genta overnight culture was labeled with 1:50diluted specific monoclonal ␣-dsDNA mouse immune serum (Abcam, Cambridge, MA) in PBS containing 1% bovine serum albumin. For the negative control, this primary antibody step was omitted. Grids were examined using a FEI Tecnai T20 transmission electron microscope at a magnification of ϫ30,000 to ϫ65,000, and images were recorded using a digital 4k-by-4k charge-coupled device (CCD) camera.
Effect of AtlA Efm and DNase on E. faecium strains. The effect of At-lA Efm and DNase on the phenotype of the E1162 wild type, atlA Efm ::genta, and atlA Efm ::gentaϩatlA Efm was analyzed during bacterial growth under a confocal microscope. The wild type, atlA Efm ::genta, and atlA Efm :: gentaϩatlA Efm were grown in BHI at 37°C for 24 h. After overnight growth, cells were diluted to an OD 660 of 0.005 in BHI and inoculated with and without purified AtlA Efm (5 g) or DNase (1.5 g/l) in a Lab-Tek chambered eight-well coverglass (Thermo, Fisher Scientific). Cell growth was recorded at regular intervals of 3 min (ϫ40 objective) with the DFC360 FX digital camera kit SP5 (Leica).
ELISA. Binding of the E1162 wild type, atlA Efm ::genta, and atlA Efm :: gentaϩatlA Efm to immobilized fibrinogen, fibronectin, collagen type I, collagen type III, and collagen type IV was assayed by a whole-cell ELISA as described previously (51).