A framework for paired-sample hypothesis testing for high-dimensional data

Ioannis Bargiotas Argyris Kalogeratos Nicolas Vayatis

Centre Borelli, ENS Paris-Saclay Gif-sur-Yvette, France

Abstract-The standard paired-sample testing approach in the multidimensional setting applies multiple univariate tests on the individual features, followed by p-value adjustments. Such an approach suffers when the data carry numerous features. A number of studies have shown that classification accuracy can be seen as a proxy for two-sample testing. However, neither theoretical foundations nor practical recipes have been proposed so far on how this strategy could be extended to multidimensional paired-sample testing. In this work, we put forward the idea that scoring functions can be produced by the decision rules defined by the perpendicular bisecting hyperplanes of the line segments connecting each pair of instances. Then, the optimal scoring function can be obtained by the pseudomedian of those rules, which we estimate by extending naturally the Hodges-Lehmann estimator. We accordingly propose a framework of a two-step testing procedure. First, we estimate the bisecting hyperplanes for each pair of instances and an aggregated rule derived through the Hodges-Lehmann estimator. The paired samples are scored by this aggregated rule to produce a unidimensional representation. Second, we perform a Wilcoxon signed-rank test on the obtained representation. Our experiments indicate that our approach has substantial performance gains in testing accuracy compared to the traditional multivariate and multiple testing, while at the same time estimates each feature's contribution to the final result.

Index Terms-statistical hypothesis testing, paired-sample testing, p-value correction, Hodges-Lehmann estimator, pseudomedian, multidimensional data.

I. INTRODUCTION

In many situations, it occurs that the data instances of two samples are paired, as they correspond to two measurements of the same population of subjects. The problem of testing the homogeneity of two paired samples arises in various applications, ranging from medicine to finance. Especially in medicine, comparing a cohort before and after a specific treatment is very common. Scientists are interested in knowing if there is any statistical difference between the paired samples, and they investigate that either by multivariate hypothesis testing (e.g. using the multivariate Hotelling T 2 -test, HT2 for short) or by relying on a multiple testing (MT) strategy that performs a series of univariate tests. For instance, one can apply the MT strategy using the paired t-test, which comes with a normality assumption about the data. Another possibility is to use the well-known one-sample Wilcoxon signed-rank test (WSR), which is based on linear rank statistics. WSR is one of the most widely-used nonparametric approach in pairedsample hypothesis testing (a resource for several univariate tests can be found in [START_REF] Hollander | Nonparametric statistical methods/Ch. 2 The One-Sample Location Problem[END_REF]). Rank statistics have been introduced to avoid making assumptions about the distribution of a data sample under the null hypothesis. However, the MT strategy has raised a well-known scientific debate concerning mainly the significant increase of Type I error (false positive rate). Therefore, MT is usually followed by p-value adjustments that aim to control this increase [START_REF] Feise | Do multiple outcome measures require p-value adjustment?[END_REF]. Among the various proposed p-value adjustments, Bonferroni correction is the one used in the vast majority of studies, as it is simple to comprehend and easy to calculate, although it is slightly less powerful than others [START_REF] Gou | A class of improved hybrid hochberg-hommel type step-up multiple test procedures[END_REF].

Besides the natural interest in Statistics, the Machine Learning community has shown an increasing interest in hypothesis testing. For instance, in the case of two-sample testing, they have proposed meaningful solutions of how to learn an optimal decision function (i.e. a classifier) with which the dataset can be optimally scored and -in a second step-tested for any significant change between the samples [START_REF] Hediger | On the use of random forest for two-sample testing[END_REF][START_REF] Kim | Classification accuracy as a proxy for two-sample testing[END_REF][START_REF] Lopez-Paz | Revisiting classifier two-sample tests[END_REF][START_REF] Rosenblatt | Better-than-chance classification for signal detection[END_REF]. Many studies have stressed the importance of choosing a suitable classification criterion for this purpose. For example, the fact that the empirical AUC can be viewed as the extension of Mann-Whitney U -statistic in the multivariate setup [START_REF] Clémençcon | Empirical performance maximization for linear rank statistics[END_REF][START_REF] Cortes | AUC optimization vs. error rate minimization[END_REF], brought the conclusion about the superiority of maximizing the area under the ROC curve (AUC) criterion in non-parametric two-sample testing. This result led to a number of algorithms for two-sample testing for multidimensional data, based on the AUC maximization [START_REF] Bargiotas | Revealing posturographic profile of patients with parkinsonian syndromes through a novel hypothesis testing framework based on machine learning[END_REF][START_REF] Clémençon | AUC optimization and the two-sample problem[END_REF].

Despite the excessive need for paired-sample analysis, neither theoretical foundations nor practical recipes have been proposed for how such an optimal decision function can be defined for paired data samples. The general goal of this work is to rely on pairs of instances for defining an optimal scoring function that will allow the scoring of the dataset, the application of pairedsample testing, and the identification of data features that have significantly changed between the two samples. To this end, we present a framework for multidimensional paired sample testing, which is inspired by the testing procedures introduced in prior works for the two-sample setting [START_REF] Bargiotas | Revealing posturographic profile of patients with parkinsonian syndromes through a novel hypothesis testing framework based on machine learning[END_REF][START_REF] Clémençon | AUC optimization and the two-sample problem[END_REF][START_REF] Hediger | On the use of random forest for two-sample testing[END_REF]. At the core of our approach lies the definition we introduce for the pseudomedian hyperplane, formulated by the Hodges-Lehmann estimator, which essentially allows us to see the paired-sample testing problem from the classification viewpoint.

Motivational example

When researchers want to compare a cohort before and after a specific medical treatment using multiple features, they can employ already proposed multivariate hypothesis tests. Despite their usefulness, such tests often have limitations in several dimensions [START_REF] Bai | Effect of high dimension: by an example of a two sample problem[END_REF], or they might exhibit inadequate power when dealing with slightly difficult cases, such as heavy-tailed multivariate data [START_REF] Purdom | Error distribution for gene expression data[END_REF]. For this reason, the use of multivariate tests, e.g. the HT2 test and its variants, is less common. Researchers usually rely on an MT procedure to find those features that change significantly after the treatment. Nevertheless, the plethora of acquisition modalities and computed features challenge the hypothesis testing when following this approach. Many works either fail to correctly use MT by adjusting the p-values [START_REF] Williams | Low power, type II errors, and other statistical problems in recent cardiovascular research[END_REF], or they do not proceed to such adjustment because it is criticized for a significant increase in Type II error (false negative rate) [START_REF] Feise | Do multiple outcome measures require p-value adjustment?[END_REF][START_REF] Perneger | What's wrong with Bonferroni adjustments[END_REF]. For that reason, several biostatisticians recommend disclosing all the results (significant or non-significant) of the conducted analysis. The violation of this recommendation, the habitual misuse of those tests [START_REF] Thiese | The misuse and abuse of statistics in biomedical research[END_REF], the little agreement on when or how to adjust the p-values, combined with the relatively small available cohorts, may lead to false conclusions, lack of reproducibility, and, as a consequence, a significant delay to reach scientific consensus. Such problems led certain scientific communities to redefine statistical significance and to embrace debatable guidelines [START_REF] Benjamin | Redefine statistical significance[END_REF].

Contribution

This work proposes a framework for multivariate pairedsample hypothesis testing by aggregating information extracted independently from each pair of instances. For each pair, we compute a perpendicular bisecting hyperplane, which is perpendicular to the line segment connecting the two instances and passes from the segment's midpoint. By definition, any point on such a hyperplane is equidistant to the two associated instances, and also separates the data space in two parts, similar to what the midpoint does in the unidimensional setting. We propose an aggregation of those hyperplanes in the vain of a pseudomedian estimated by the Hodges-Lehmann estimator, which is directly related to the Wilcoxon signed rank (WSR) test, where it has been proven to be a consistent estimator for the unidimensional case. We call our test Multivariate WSR (MWSR), since its principle can be seen as a multivariate extension of the WSR test, where the Hodges-Lehmann pseudomedian hyperplane provides a scoring of the data. Although the perpendicular bisecting hyperplanes provide 'naive' linear decision rules that separate the two samples, yet their pseudomedian can be an adequate decision rule that can be used for scoring and ranking the instances. Our approach avoids the p-value adjustment, and the contribution of every feature to the final result can be measured using the linear coefficients of the pseudomedian hyperplane.

Finally, we provide empirical evidence that our test constitutes an important alternative to multivariate test or MT (especially for marginally separated cases) as it leads to clearly improved performance.

II. BACKGROUND A. Wilcoxon Signed Rank test (WSR)

Suppose we are in the univariate setting, the two paired samples are denoted by X, Y ∈ R, and (X i , Y i ) is the i-th out of the N pairs of data instances that corresponds to two measurements of the same subject, the i-th of the population. The WSR test has been designed to test whether a population is symmetric around a given median. In the context of the twosample testing, the WSR test is used to study the population

Z = {Z 1 , ..., Z N } comprising the paired differences Z i = Y i - X i , for i = 1, ..., N . Assumptions 1)
The differences Z 1 , ..., Z N are independent and identically distributed (i.i.d.), i.e. Z i ∼ F , where F represents their continuous probability density function. 2) F should be a symmetric distribution around its median θ, which is usually referred to as effect size (e.g. the effect size of a certain medical procedure) [START_REF] Hollander | Nonparametric statistical methods/Ch. 2 The One-Sample Location Problem[END_REF]. Formally, this suggests that:

θ: F (θ + x) = F (θ -x), ∀x ∈ R.

Null and alternative hypotheses

The type of change the WSR detects is a distribution shift from X to Y by a quantity θ, which is defined as the median of the differences contained in Z. Thinking in regard to the pairwise differences Z i , essentially converts the paired-sample testing problem to one-sample testing over the population Z. Consequently, the null and the alternative hypotheses of the WSR paired-sample test can be derived from those of the one-sample test. The two-sided WSR test [START_REF] Wilcoxon | Individual comparisons by ranking methods[END_REF] considers: H 0 : θ = 0, the Z i 's are symmetric around θ = 0;

H 1 : θ ̸ = 0, the Z i 's are symmetric around θ ̸ = 0. (1)
Under the null hypothesis H 0 , there is no difference between the two samples, and therefore the considered treatment or intervention caused a zero shift (θ = 0) in location. This simply asserts that the Z i 's are symmetrically distributed around 0.

Procedure

The WSR statistic T + Z is calculated from the absolute values |Z 1 |, ..., |Z N | of the differences and their signs. The absolute values are sorted in ascending order, and let R i denote the rank of |Z i |, i = 1, ..., N , in that latter order. Let the positive sign indicator variables ψ i , i = 1, ..., N, where

ψ i = 1 if Z i > 0, and 0 if Z i < 0. The product ψ i R i is the positive signed rank of Z i . The WSR statistic T +
Z is then the sum of the positive signed ranks:

T + Z = N i=1 ψ i R i . (2)
According to the two-sided version of the test, and at a given significance level α, if the X and Y samples differ only by a location shift θ, then the (Z i -θ)'s have a distribution that is symmetric around 0. The default procedure is for testing with θ = 0. To test with θ = θ 0 ̸ = 0, one needs to subtract a given θ 0 of interest from each Z i difference to form a modified sample Z (θ0) = {Z 1 -θ 0 , ..., Z N -θ 0 }. Then, T + Z (θ 0 ) is computed as the sum of the positive signed ranks for Z (θ0) .

B. Hodges-Lehmann estimator as effect size

The effect size θ associated with the WSR statistic T + Z , is estimated by the pseudomedian of the differences, which is in turn estimated by the Hodges-Lehmann estimator θ [START_REF] Jr | Estimates of Location Based on Rank Tests[END_REF]:

θ = median 1 2 (Z i + Z j ); ∀i ≤ j = 1, ..., N . (3)
In words, this is the median of all the M = 1 2 N (N +1) pairs of midpoints of the differences. The vector

W Z = [ 1 2 (Z i + Z j )] ij , ∀i ≤ j = 1, .
.., N is known as the vector of Walsh averages [START_REF] Walsh | Some significance tests for the median which are valid under very general conditions[END_REF]. Supposing there are no ties and no zeros among the Z i 's, the number of positive Walsh averages (W + Z ) is equal to the WSR statistic T + Z [START_REF] Høyland | Robustness of the Hodges-Lehmann estimates for shift[END_REF]. Without loss of generality, suppose that F satisfies the fact that the median of differences (θ) and the pseudomedian are unique. Then, the Hodges-Lehmann estimator θ is a consistent estimator of the pseudomedian, which can generally be different from the median θ [START_REF] Jr | Estimates of Location Based on Rank Tests[END_REF]. However, when F is symmetric, the median and the pseudomedian coincide.

As we mentioned above, θ is associated with the WSR test. When θ = 0, then the distribution of the statistic T + Z is symmetric around its mean rank M

2 . An empirical estimator of θ is the amount θ that should be subtracted from each Z i so that the value of T + Z (θ 0 ) is as close to M 2 as possible. Briefly, one needs to calculate the amount that Z should be shifted in order to become a sample with median 0. Recall that, under Assumption l and 2, each of the (Z i -θ)'s comes from a population median that is equal to 0.

C. Approaches for multivariate paired-sample testing

Here, we present the basic ideas behind two of the most widely-used approaches for multivariate paired-sample testing, where X, Y ∈ R d , d > 1. One option is to use MT with p-value correction. For instance, to test with a target significance level α, the Bonferroni correction considers d independent tests, one per dimension, and uses in each of them the 'corrected' significance α d [START_REF] Abdi | Bonferroni and šidák corrections for multiple comparisons[END_REF]. Another alternative is the multivariate Paired Hotelling T 2 -test (HT2) [START_REF] Hotelling | The generalization of student's ratio[END_REF], which is a generalization of paired t-test. Its null hypothesis asserts that the mean vectors µ X , µ Y , of the population are equal against the alternative hypothesis that these mean vectors are not equal, hence H 0 :

µ X = µ Y against H 1 : µ X ̸ = µ Y .
Same as in the univariate case, the differences between the paired instances are Z i = Y i -X i , for i = 1, ..., N and the hypothesis testing becomes

H 0 : µ Z = 0 against H 1 : µ Z ̸ = 0 where µ Z = µ Y -µ X .
Under the assumption of normality for Z, the HT2 test statistic is given by:

T 2 = N ZT Σ -1 Z Z, (4)
which is proportional to the distance between m Z and 0, where Z is the sample mean vector and Σ -1 Z is the inverse covariance matrix (or precision matrix).

III. PROPOSED PAIRED-SAMPLE TESTING FRAMEWORK

Testing in multiple dimensions

Let us start by establishing that, in the unidimensional case,

X Y distance 1 rule 1 distance 2 rule 2 distance 3 rule 3 X Y distance 1 rule 1 distance 2 rule 2 distance 3 rule 3
Pseudomedian Pseudomedian (a) Example in 1d. This is the typical paired-sample setting. We draw as dotted lines the intervals of the differences Z 1 , Z 2 , Z 3 between the instances of three pairs. The midpoints and the pseudomedian rules are shown as continuous lines vertical to the axis.

X Y distance 1 rule 1 distance 2 rule 2 distance 3 rule 3 X Y distance 1 rule 1 distance 2 rule 2 distance 3 rule 3 Pseudomedian Pseudomedian (b) Example in 2d.
Here, each distance associated to a pair of instances is shown as a dotted line, while the 1d midpoint rule is now a hyperplane (continuous lines) that is a perpendicular bisector to the line segment connecting the pair, and hence pass from the segment's midpoint. by adding θ 2 to the median of X, we can get what we term as a pseudomedian decision rule separating the space in two parts. This is a notion that we generalize to propose our paired-sample framework for multiple dimensions, and the Multivariate WSR test (MWSR). For d > 1, we deal with the point-cloud of X, Y ∈ R d , thus the Euclidean distance Z i = ||Y i -X i || 2 between pairs of instances (X i , Y i ) can be seen as the analogue to the difference Z i = Y i -X i we saw earlier in 1d. Further, each separating rule associated with a midpoint of the 1d case, now becomes a (d-1)-dimensional perpendicular bisecting hyperplane: it is perpendicular to the line segment connecting the pair and also passes from its midpoint. Each such hyperplane is computed by taking into account only one specific pair of instances, yet it splits the space in two parts, and therefore it can be seen as a decision rule that could hopefully classify the data in two parts, the X and the Y part. Fig. 1 shows an insightful visual example: Fig. 1a presents the 1d setting, as discussed in the earlier sections: the paired samples lie on a line, each difference Z i is a distance on that line, and midpoints can be computed for each (X i , Y i ) pair. Fig. 1b shows how our approach operates in multiple dimensions. The decision rule for a pair (X i , Y i ) is its perpendicular bisecting hyperplane C i :

w T i x + b i = 0, x ∈ R d
, where w i are the coefficients of the hyperplane, and b i is its intercept. With the known elements, it is easy to compute the exact w i and b i ; then, sign(C i (x)) classifies an input datapoint x, while the signed distance from x to C i can be seen as a score for x based on this hyperplane.

In simple terms, we propose to see the scoring function derived from the hyperplane of the i-th pair as an analogue to the half of the empirical effect size calculated by this specific pair of instances. Furthermore, the aggregation of all the N rules can provide a more robust scoring function. This aggregation can be computed by directly defining the pseudomedian decision rule in multiple dimensions: The linear coefficients w i and the intercept b i of each C i rule are used for estimating the Hodges-Lehmann pseudomedian hyperplane through the ŵ, using the same procedure as θ (see Sec. II-B). The two-step MWSR testing procedure We compute a linear decision rule for each pair (X i , Y i ), hence we gather N such models. In order to extend the procedure presented in Sec. II-B to the general d-dimensional setting, we need to define first the Walsh model-wise average

W C = 1 2 (C i + C j ) ij , ∀i ≤ j = 1, .
.., N , as the average of the coefficients and intercepts for pairs of linear models. For the linear decision rules that we employ, this is: Having specified how we can pass from averaging instances to averaging model parameters, we can now express the Hodges-Lehmann estimator in the form as before:

Ĉ * = median 1 2 (C i + C j ); ∀i ≤ j = 1, ..., N , (5)
which computes the median coefficient-wise. The estimated Ĉ * is the empirical optimal scoring function where the final classification scores ( Ŝ1 , Ŝ2 ) are calculated for X and Y , respectively. Note that the operation of finding Ĉ * from a set of N decision rules can be seen as model aggregation, which in this case is not a generic one (e.g. like the voting used in Ensemble Learning), but rather a special aggregation, directly linked to the WSR test that we employ next. The last step amounts to performing a WSR test on Ŝ1 and Ŝ2 . In addition to the given p * -value and the effect size θ * , the user also gets informed about the contribution of each dimension to the final result (through the coefficients of the linear model Ĉ * ).

The feature importance can be easily obtained considering that the more orthogonal is a decision rule (hyperplane) to a certain feature axis, the higher the associated coefficient will be. Furthermore, the sign indicates the positive or negative relationship with the response. The procedure of the proposed framework is detailed in Alg. 1.

IV. EXPERIMENTS

The experiments in this section compare the performance of our testing method, the MWSR test, to the classical multiple testing (MT) procedure with Bonferroni p-value adjustment, and the HT2 test. We use several synthetic scenarios and one real clinical dataset. The synthetic scenarios have increasing difficulty and allow us to check the limits of the proposed test. The comparison is two-fold: i) in terms of sensitivity in detecting a difference between the paired samples (p-value, effect size), and ii) in terms of estimating accurately the contribution of each feature to the test's outcome, which is helpful for the user to interpret the results.

Algorithm 1

The MWSR paired-sample testing framework

Input: X, Y ∈ R N ×d are the 2 • N paired samples; Output: Ĉ * , ( Ŝ * 1 , Ŝ * 2 ), p * -value, θ * , I * First step: Compute a scoring for i = 1, .., N do Ci ← perpendicular_bisector(Xi, Yi) end for k ← 1; M ← 0N×N for i = 1, .., N do for j = i, .., N do W C,k ← 1 2 (Ci + Cj) ▷ the Walsh average of hyperplanes k ← k + 1 end for end for Ĉ * ← median(WC )
▷ the pseudomedian aggregate, see Eq. 5

Ŝ * 1 , Ŝ * 2 ← get_scores( Ĉ * (X, Y )) ▷ classification-based scoring Second step: Paired-sample test over the computed scores p * -value, θ * ← WSR( Ŝ * 1 , Ŝ * 2 )
▷ p-value and effect size

I * ← w( Ĉ * ) ▷ feature importance index return Ĉ * , ( Ŝ * 1 , Ŝ * 2 ), p * -value, θ * , I *

Synthetic datasets

Synthetic data are simulated by pairing two samples coming from two Gaussian distributions, with feature-wise correlation:

R X (k) ,Y (k) = cov(X (k) , Y (k) ) σ X (k) σ Y (k) = 0.5, (6)
where X (k) , Y (k) are vectors representing the paired samples with only the k-th feature, cov(•, •) is their covariance, and σ X (k) , σ Y (k) are the respective standard deviations. To produce the dataset for each scenario, we set:

• a fixed population size, N = 30 pairs of data instances;

• the dimensions d = {10, 20, 30, 60}, mimicking the number of features that a usual study may have;

• a standard deviation std = {1, 2};

• the first 90% of the dimensions to have no statistical difference between the two samples, and hence to be randomly drawn (separately) from N (0, std);

• the last 10% of the dimensions to present the same difference in mean value, hence producing a linear shift. We allow this shift to also increase progressively from 0 to 1 to investigate the detection sensitivity of the methods. Given a scenario, we generate 20 cases and apply all statistical tests. We report the average performance, namely the percentage of cases with a significant shift obtained by each statistical test, as a function of the size of the shift (referred to as average difference in the figures). Moreover, we provide results regarding the inferred feature importance. Both elements should be examined jointly to validate the acquired results further.

Real dataset

We extend our empirical validation by employing a typical real clinical dataset, with a relatively low population and multiple features. It concerns posturographic assessment for subjects with Parkinson's syndrome (PS). This dataset, initially used in [START_REF] Bargiotas | Predicting future falls of parkinsonians using posturography and random forest[END_REF] 2d). The performance is presented as a function of the separation distance between the two distributions (average difference on x-axis).

Neurology department of the HIA, Percy hospital in Clamart, France, who were diagnosed with PS. The subjects underwent a posturography assessment using a force platform (here a Wii Balance Board (Nintendo, Kyoto, Japan)) that captures the trajectory of the center of pressure (CoP) exerted by the entire body over time when an individual stands on them. The assessment comprises two examinations with a 6-month difference in time, which are the paired samples we use in our experiment. Each time their postural stability was recorded for 25 seconds while maintaining an upright position on a force platform with eyes open.

To characterize subjects' postural control, the dataset provides 16 features that had been previously proposed as indicators of postural stability [START_REF] Quijoux | A review of center of pressure (COP) variables to quantify standing balance in elderly people: Algorithms and open-access code[END_REF]. In detail: Percentiles (95% and 5%) (cm), Range (cm), Variance (cm 2 ), Mean Instant Velocity (cm/s), Acceleration (cm/s 2 ) and Frequency (Hz) below which 95% of the signal energy is found, for both X-medio-lateral (ML) and Y-antero-posterior (AP) axes, confidence ellipse area (cm 2 ) that covers the 95% of the points of the trajectory and the angular deviation (in degrees • ).

Results and discussion

In all the synthetic cases in Fig. 2, the proposed MWSR test is always by far superior to the MT and the multivariate HT2 test in detecting the difference between the paired samples. The MWSR is more sensitive than the compared tests, as it manages to detect much smaller shifts. Moreover, it reports no false positive results in the case of no shift in location. Worth noting that when d > N , the classical HT2 test fails completely due to the non-existence of the inverse of the sample covariance matrix (see in Fig. 2 the four cases with d ≥ 30). The difference in performance between the MWSR and the MT increases in higher dimensions as well as when the variance is higher (see Fig. 2, e-h), which indicates a robust behavior to that kind noise. Note that the MWSR has no particular computational overhead; it runs in less than 0.05s for the small real dataset, and this increases slightly with N .

Concerning the feature importance indexes obtained by the MWSR, they are effectively recognizing the last 10% of the features that contribute to the distribution shift (see in Fig. 3 the features x55-x60). This is made clear by the higher weights (i.e. larger stacks) assigned to those features.

In the real dataset, the MWSR test found that the Parkinsonians changed significantly after 6 months, while the MT did not detect any significant change.

Despite the above notable advantages, the proposed methodology is designed and limited for linear shifts in location between the paired samples, same as the univariate WSR test. Moreover, a point that needs attention is to normalize the data so that the feature importance indexes are better estimated and are easier to interpret. The weight assigned to a particular feature is influenced by the other features, particularly when there are correlations among them.

V. CONCLUSIONS

We presented a sound strategy for paired-sample testing for multidimensional data, which outperforms the classical multivariate approaches (Hotelling T 2 test, the multiple testing with p-value adjustment) in our experiments. Our method has two additional advantages: i) it is generally simple and understandable, and ii) it offers a tool that allows the user to interpret the actual contribution of every feature to the final result. Beyond the illustrative preliminary experiments in real and synthetic datasets that we presented in this paper, we intend to investigate the relative efficiency of this test to other setups. We also plan to investigate theoretically the regime in which this method is qualitatively superior.

Figure 1 :

 1 Figure1: An illustrative example of how we extend the WSR paired-sample hypothesis testing from one dimension (a) to multiple dimensions with the MWSR test (b). We see each pair's midpoint as a decision rule separating the data space, and we generalize this in multiple dimensions by means of perpendicular bisecting hyperplanes.

Figure 3 :

 3 Figure3: The relative feature importance in the paired-sample testing result of different approaches, and their corresponding average performance for N = 30 pairs coming from two Gaussian distributions with std = 1, in d = 60 dimensions (corresponds to the case of Fig.2d). The performance is presented as a function of the separation distance between the two distributions (average difference on x-axis).

Features

  A p-value per feature obtained by MT

Figure 4 :

 4 Figure 4: The relative importance of each feature as indicated by our approach and the MT on the posturographic dataset (d = 16 features).

  , includes 30 subjects (mean age: 79.6±4.4 years) from the

			Dimensions: 60, N:30, Variables std: 1			Dimensions: 60, N:30, Variables std: 1
	Stacked Average w of every Feature	0.00 0.02 0.04 0.06 0.08 0.10	0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0 Average Difference	Stacked Significant Features Percentage(%)	0 100 200 300 400 500 600	0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 Average Difference 0.95 1.0	. . . . . .	x30 x54 x55 x56 x57 x58 x59 x60 x1
	x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 (a) Average hyperplane coefficients ( ŵ) for each feature (d = 60) x13 x14 x15 x16 x17 x18 x19 x20 x21 x22 x23 x24 x25 x26 x27 x28 x29 x30 x31 x32 x33 x34 x35 x36 x37 x38 x39 x40 x41 x42 x43 x44 x45 x46 x47 x48 x49 x50 x51 x52 x53 x54 x55 x56 x57 x58 x59 x60	x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 (b) Percentage of times (%) that features are found significant after MT x13 x14 x15 x16 x17 x18 x19 x20 x21 x22 x23 x24 x25 x26 x27 x28 x29 x30 x31 x32 x33 x34 x35 x36 x37 x38 x39 x40 x41 x42 x43 x44 x45 x46 x47 x48 x49 x50 x51 x52 x60 x59 x58 x57 x56 x55 x54 x53 with Bonferroni correction (d = 60)

The average performance of paired-sample testing approaches in synthetic datasets with N = 30 pairs of instances coming from two Gaussian distributions, both with either std = 1 or 2 (first and second row, respectively). The performance is presented as a function of the separation distance between the two distributions (average difference on x-axis). On the y-axis it appears the % of significant results detected by each method, and on the x-axis the progressive difference in the mean value for the 10% of the dimensions (d).


Stucki David

Brites Daniela

Coscolla Mireia

Trauner Andrej

Stefan Niemann

Rutaihwa Liliana

Borrell Sonia

Marie Ballif

Sebastien Gagneux
email: sebastien.gagneux@unibas.ch

Jeljeli Leïla

Prof. Sebastien Gagneux

M Bo

F.N S Ak

Mycobacterium tuberculosis lineage 4 comprises globally distributed and geographically restricted sublineages

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Ecologists distinguish between generalists and specialists depending on the width of an organism's ecological niche1. In infectious diseases, the niche of a given pathogen is determined by host range and the agent's capacity to survive in the environment2. Some microbes are obligate pathogens restricted to one or several host species3,4, others are mainly free-living and only occasionally pathogenic5. Little is known on the niche width of obligate human pathogens3. The causative agent of tuberculosis, known as the Mycobacterium tuberculosis complex (MTBC), is an obligate pathogen that comprises seven phylogenetic lineages adapted to humans and two lineages adapted to various wild and domestic animal species6. Some human-adapted MTBC lineages have received particular attention. For example, Lineage 2, which includes the Beijing family of strains, has repeatedly been associated with drug resistance7. Lineage 2 likely originated in East Asia8,9, and has recently been expanding in some parts of the world10. By contrast, Lineages 5 and 6 (also known as Mycobacterium africanum West Africa I and II), and Lineage 7 are largely restricted to West-and East Africa, respectively11,12. The observation that the human-adapted MTBC population is phylogeographically structured has led to the hypothesis that the different lineages might be adapted to particular human populations13. Support for this notion comes from the observation that sympatric host-pathogen associations in human tuberculosis remain stable over time, even in metropolitan settings where host and pathogen populations intermix14-17. Moreover, sympatric host-pathogen associations are perturbed in HIV coinfected patients14, indicating that in the context of reduced host immune-competence, the different lineages can successfully infect and cause disease irrespective of the host genetic background.

Contrary to the other main human-adapted MTBC lineages, Lineage 4 occurs at significant frequencies on all inhabited continents18. It is hence geographically the most widespread cause of human tuberculosis19. Yet, the reasons for this global success are unknown. Lineage 4 has been shown to exhibit enhanced virulence in macrophage and animal models of infection, albeit with much variation between different Lineage 4 strains19,20. Moreover, molecular epidemiological studies have reported considerable variation in the transmission success of different Lineage 4 strains in clinical settings19. These observations suggest that Lineage 4 is genetically and phenotypically diverse, and this diversity might determine the epidemiology of different Lineage 4 subtypes in different parts of the world. The purpose of this study is to get a better understanding of the global population structure of Lineage 4 and the evolutionary forces that have contributed to the success of Lineage 4 across the world. For this we combined large-scale single nucleotide polymorphism (SNP)-typing with targeted whole-genome sequencing of a global collection of Lineage 4 clinical isolates.

Results

MTBC Lineage 4 comprises 10 separate sublineages

We first analyzed 72 published genome sequences of Lineage 4 clinical strains from global sources21,22. These strains harbored 9,455 variable single nucleotide positions which divided Lineage 4 into 10 sublineages (L4.1.1 to L4.10 in Fig. 1a and Supplementary Fig. 1). We used four complementary approaches to validate these sublineages. First, we performed a principal component analysis, which showed a clear separation of seven sublineages (L4.1.1, L4.1.3, L4.1.2, L4.2, L4.3, L4.4, L4.10; Supplementary Fig. 2). Sublineages L4.5, L4.6.1/Uganda and L4.6.2/Cameroon were less clearly separated. Second, we found that the mean pairwise genetic distance between pairs within the sublineages was significantly lower than between sublineages (276 SNPs versus 602 SNPs, Wilcoxon rank sum test, p < 0.0001, Supplementary Fig. 3). Overall, the mean pairwise SNP distance between any two strain pairs was 565 SNPs. Third, we calculated pairwise fixation indexes (F ST ) to evaluate the degree of population differentiation. All F ST values between the sublineages were larger than 0.33 (Supplementary Table 1), indicating that these populations are separated. Fourth, we mapped previously reported phylogenetic markers onto our genome-based phylogenetic tree15,23-28. Most of these markers were congruent with our sublineage definition (Supplementary Fig. 1).

Sublineages differ in their phylogeographic distribution

Because the MTBC exhibits limited sequence variation and no signficant ongoing horizontal gene exchange, SNP homoplasies are extremely rare, making SNPs ideal phylogenetic markers29. We further scrutinized the 9,455 variable positions among the 72 MTBC Lineage 4 genomes, and found 51 to 277 specific for one of each of the 10 sublineages. All of these variable positions were mutually exclusive, i.e. they showed no homoplasy. We selected a subset of these sublineage-specific SNPs and used these to screen a global collection of 3,366 Lineage 4 clinical isolates from 100 countries using various genotyping platforms30-35. First, we developed a novel sublineage-specific multiplexed SNP-typing assay using the Luminex platform as previously reported36, and used that method to screen 2,001 isolates (Supplementary Table 2). In addition, we screened 741 isolates using the Sequenom MassARRAY platform (Supplementary Table 3)37, and 624 isolates by PCR and Sanger sequencing (Supplementary Table 4). Overall, 3,181/3,366 (94.5%) Lineage 4 isolates were successfully assigned to a sublineage (Supplementary Table 5). An additional 92/3,366 (2.7%) isolates harbored the reference allele for all sublineages, indicating they belonged to one or several additional and unknown sublineages. For the remaining 93/3,366 (2.8%) isolates, no classification could be obtained for various technical reasons. Among the 3,181 Lineage 4 isolates assigned to one of the 10 sublineages, L4.3/LAM was the most frequent, accounting for 20.3%, followed by L4.6.1/Uganda (14.2%), L4.10/PGG3 (11.9%), L4.4 (10.1%), and L4.1.2/Haarlem (9.9%) (Fig. 1b).

Mapping the proportion of each sublineage by country showed that the sublineages differed in their geographical distribution (Fig. 2). Specifically, L4.1.2/Haarlem, L4.3/LAM and L4.10/PGG3 occurred globally (Fig. 3a, Supplementary Fig. 4). By contrast, L4.1.3/Ghana, L4.5, L4.6.1/Uganda and L4.6.2/Cameroon occurred at high frequencies in specific regions of Africa or Asia, and were almost completely absent from Europe and the Americas (Fig. 3b). The geographical spread of the three remaining sublineages was intermediate (Fig. 2, Supplementary Figs. 4 and5). L4.1.1/X mainly occurred in the Americas and in lower proportions in few countries of Southern Africa, Asia and Europe. L4.2 and L4.4 occurred in high proportions among isolates from particular countries in Asia and Africa, but were largely absent from the Americas (Fig. 2, Supplementary Figs. 4 and5). A similar pattern of sublineage distribution was observed when normalizing by TB prevalence38 and country surface area (Supplementary Fig. 6).

Populations that occupy a broader variety of environments may exhibit a wider geographic distribution. Humans differ in their susceptibility to TB39, and human genetic diversity may thus determine the width of the ecological niche accessible to different MTBC genotypes40,41. The geographical restriction of particular MTBC genotypes might reflect local adaption of these pathogen variants to the corresponding human host populations13,15. Such a sympatric host-pathogen association in human TB is compatible with the "local" sublineages observed here, and supports the notion that these sublineages represent ecological specialists. By contrast, the three "global" sublineages could represent generalists capable of infecting and causing disease in many different human populations. This notion was supported by the fact that the three generalist sublineages L4.1.2/Haarlem, L4.3/LAM and L4.10/PGG3 were observed in 49, 47 and 47 countries, respectively, whereas the specialist sublineages L4.1.3/Ghana, L4.5, L4.6.1/Uganda and L4.6.2/Cameroon were only found in few countries each (3, 7, 9 and 10 countries, respectively). The country frequencies of the remaining three sublineages L4.1/X, L4.2/Ural and L4.4 were intermediate [START_REF] Coll | A robust SNP barcode for typing Mycobacterium tuberculosis complex strains[END_REF][START_REF] Fenner | HIV infection disrupts the sympatric host-pathogen relationship in human tuberculosis[END_REF] and 26 countries, respectively) (Supplementary Fig. 4).

The different geographical distribution of generalist and specialist sublineages could be due to intrinsic biological factors, extrinsic factors such as human migration, or both. Hence we next performed various population genomic analyses to explore the genomic characteristics of these Lineage 4 generalists and specialists, as well as the role of human migration in the global spread of the most successful generalist sublineage.

Genomic features of generalist and specialist sublineages

The geographic and niche distribution of populations can be correlated with their genetic variability or with that of their ancestors. One possible reason for the restricted host range of the specialist sublineages might be historical, i.e. the ancestor populations of the extant specialist populations may have harbored more deleterious mutations, restricting their host range. To assess this possibility, we characterized the mutations which contributed to the divergence of the different sublineages; these mutations are variants that have become fixed during the evolution of these sublineages. We focused on the substitutions that occurred in all isolates of any of the generalist sublineages (L4.1.2/Haarlem, L4.3/LAM and L4.10/ PGG3) and compared them to the substitutions that occurred in all isolates of any of the three specialist sublineages (L4.6.1/Uganda, L4.5 and L4.6.2/Cameroon). We identified nonsynonymous SNPs predicted to have a functional effect using SIFT42. We found that overall, the specialist and generalist sublineages showed a similar proportion of fixed substitutions (among all substitutions) predicted to impact gene function (23.0% versus 20.6%, χ 2 test p=0.62; Supplementary Table 6), suggesting that the mutational load of the ancestor populations did not differ significantly between generalists and specialists.

Small populations with restricted geographic ranges are expected to have reduced levels of genetic diversity43. Thus, one possible restriction to niche expansion by specialist sublineages could be that these sublineages have low genetic diversity precluding adaptation to new hosts. We characterized the genetic diversity associated with the process of diversification in Lineage 4 generalists and specialists, focusing on L4.3/LAM and L4.6.1/ Uganda, globally the most frequent generalist and specialist sublineages of Lineage 4 in our dataset, respectively (Fig. 1b). We analyzed the whole-genome sequences of 293 L4.3/LAM clinical strains representing the global diversity of this sublineage. These were selected from a global collection of 2,132 L4.3/LAM isolates based on standard genotyping data (Supplementary Table 7, Supplementary Figs. 789). For L4.6.1/Uganda, we analyzed wholegenome sequences of 203 clinical strains from Uganda and several neighboring countries (Supplementary Table 7, Supplementary Figs. 10 and11). This sample included 28 L4.6.1/ Uganda strains identified through screening of 13,067 publically available MTBC whole genome sequences (see Online Methods) 44-56. Comparing the genetic diversity between these two bacterial populations showed that L4.3/LAM was significantly more diverse than L4.6.1/Uganda (mean number of 395 SNPs between pairs compared to 215 SNPs, respectively; Wilcoxon rank sum test p<0.0001), consistent with the expected difference between generalists and specialists43.

Antigenic diversity in Lineage 4 sublineages

We previously reported that in the human-adapted MTBC, experimentally confirmed human T cell epitopes were conserved57,58. This is unlike many other pathogens where genomic regions encoding antigens tend to be diverse as a result of antigenic variation linked to immune escape59. When we assessed the evolutionary conservation of 1,226 experimentally confirmed human T cell epitopes60 in L4.6.1/Uganda by calculating their dN/dS, we found that these epitopes were significantly more conserved than the non-epitope regions of the corresponding T cell antigens (Wilcoxon rank sum test, p<0.0001, Fig. 4). This result was consistent with our previous findings for the MTBC overall57,58. However for L4.3/LAM, we saw the opposite, i.e. the T cell epitopes showed a significantly higher dN/dS than the non-epitope regions (Wilcoxon rank sum test, p<0.0001, Fig. 4). To test whether a high dN/dS in T cell epitopes is characteristic of the generalist sublineages, we analyzed the genomes of 228 L4.2/Haarlem strains and 301 L4.10/PGG3 strains identified by screening of 13,067 publically available genomes (Supplementary Table 7, Supplementary Figs. 12 and13). We found that in contrast to L4.3/LAM, the epitope regions in these generalist sublineages were more conserved than the corresponding non-epitope regions, i.e. similar to L4.6.1/Uganda and the MTBC overall57,58 (Fig. 4). Consistent with previous reports57,58,61, essential genes62 were significantly more conserved than nonessential genes in all sublineages, except L4.3/LAM in which the dN/dS of essential and nonessential genes were not significantly different (Fig. 4) One of the limitations of our dN/dS analyses was that despite a large number of genomes analyzed, within individual sublineages, the mean number of pair-wise differences in regions encoding T cell epitopes was very small (Supplementary Table 8), limiting the accuracy of dN/dS inferences for epitopes. Hence, we assessed T cell epitope diversity by comparing the number of epitopes affected by nonsynonymous variants in the different sublineages (Fig. 5). We found that in all four sublineages, the majority of epitopes were completely conserved, consistent with our previous findings for the MTBC overall57,58. However, each of the three generalist sublineages showed significantly more epitopes harboring at least one amino acid change when compared to the specialist sublineages L4.6.1/Uganda (Fig. 5, χ 2 tests p<0.0001 for all comparisons). It is possible that this comparably higher epitope diversity in generalists might reflect interactions with broader host populations.

The epitopes interrogated in our analysis were encoded by a total of 304 antigens. The number of antigens containing nonsynonymous variation in epitopes was 60, 26, 46 and 48 antigens in L4.3/LAM, L4.6.1/Uganda, L4.2/Haarlem and L4.10/PGG3, respectively. When excluding nonsynonymous mutations present only in one strain in each sublineage (which likely represent transient mutations), the number of antigens dropped to 20, 11, 12, 24 in L4.3/LAM, L4.6.1/Uganda, L4.2/Haarlem and L4.10/PGG3, respectively (Supplementary Table 9). Interestingly, 10 of those antigens exhibited independent parallel, nonsynonymous variation in epitope regions in the different sublineages (Supplementary Table 9). Of those antigens, Pst1, an adhesin promoting phagocytosis63, and FbpB, the precursor of the secreted antigen 85-B64, had already been pointed out as encoding diverse epitopes by a previous study, in which several MTBC lineages were compared57. (Supplementary Table 9). Other antigens exhibiting parallel nonsynonymous changes by different sublineages include known immunodominant, secreted antigens such as Mpb6464, MPT3265 and MPT7066 and three latency-associated antigens (Rv1733c67, Rv3034c, Rv262868, Supplementary Table 9).

Origin and global spread of the L4.3/LAM sublineage

Irrespective of the putative biological differences between the Lineage 4 sublineages, human migration could also have led to variation in the global distribution of MTBC lineages. Because the most successful sublineage of Lineage 4 was also frequently found in Europe, we hypothesized that the global success of L4.3/LAM was driven by European migration and colonization. To test this hypothesis, we first determined the most likely geographical origin of the most recent common ancestor of L4.3/LAM using two methods for reconstruction of ancestral states69. By both methods, Europe was predicted as the most likely place of origin of L4.3/LAM (100% and 99.6%, respectively) (Fig. 6a, Supplementary Fig. 14). Moreover, the ancestral geographical regions reconstructed for subsequent nodes in the phylogeny were consistent with the spread of L4.3/LAM from Europe to other parts of the world (Fig. 6a). Finally, we found that L4.3/LAM strains from Europe were genetically more diverse than L4.3/LAM strains from other continents, which further supports a European origin for this sublineage (Fig. 6b, Kruskall-Wallis test p<0.0001; Fig. 6c).

Discussion

Our findings show that the global success of Lineage 4 is a consequence of both biological and social phenomena. Specifically, we found that Lineage 4 is genetically diverse, and that this diversity is phylogeographically structured. The phylogeography of Lineage 4 supports an ecological distinction between globally represented generalists and geographically restricted specialists. Our in-depth population genomic analyses of one specialist and three generalist sublineages showed that even though the majority of human T cell epitopes were completely conserved in all four sublinages, the proportion of epitopes with amino acid substitutions was significantly higher in generalists. Finally, we demonstrate a likely European origin for L4.1/LAM, the most frequent and globally widespread generalist sublineage of Lineage 4.

Our observation that Lineage 4 is phylogenetically diverse is in line with previous findings27,70, and highlights the importance of large and globally representative samples when studying the population structure of human pathogens. We found that Lineage 4 comprises at least 10 sublineages, which differ in their geographical distribution. The phylogeography of these sublineages is consistent with an ecological separation into specialists and generalists, with some sublineages showing an intermediate geographical distribution. Our phylogenetic analyses also showed that the three generalist sublineages identified within Lineage 4 were not monophyletic (Fig. 1a), indicating that generalism was acquired multiple times independently during the evolution of Lineage 4. Specialist sublineages also emerged multiple times, which is consistent with local adaptation to separate human populations13.

One could argue that the reason for specialist sublineages being geographically restricted is they diverged later than the generalist sublineages during the evolution of Lineage 4, and thus had insufficient time to spread globally. However, based on recent findings by Comas et al.71, the African specialist sublineages already existed at least several centuries ago, perhaps even several millennia ago, depending on the age of the most recent common ancestor of the MTBC that has been estimated between 70'000 years9,21 and 6'000 years72,73. Thus, this timespan should have offered ample opportunity for the specialist sublineages to become more geographically widespread.

The genetic diversity of the specialist sublineage L4.6.1/Uganda was significantly lower than that of the generalist L4.3/LAM, as expected from populations with restricted geographical ranges43. Concomitantly, the diversity of T cell epitopes in the specialist sublineage L4.6.1/ Uganda was also significantly lower than in any other of the three generalist sublineages analyzed. Whether the low genetic diversity of the specialist sublineage has hindered the adaptation of these strains to other human populations or reflects a restricted niche due to the lack of opportunity for spreading will need to be explored in future studies.

In all sublineages analyzed, the large majority of T cell epitopes were completely conserved, which is in agreement with previous reports for the MTBC overall57,58. This suggests that both these generalists and specialists do not use antigenic variation as a main mechanism of immune evasion. Despite this general trend, we found that some antigens have acquired nonsynonymous mutations in parallel in the different sublineages, suggesting that variation in these particular antigens might be beneficial. For example, acquiring nonsynonymous variation may allow particular antigens to be recognized by T cell receptors of different human populations, which might be beneficial in the presence of different human HLA alleles58. This could also provide an explanation for the differences in the degree of variation in T cell epitopes of the generalist and specialist sublineages, as generalist sublineages are expected to interact with a broader range of HLA alleles. Alternatively, some nonsynonymous mutations in epitopes might reflect escape from human T cell recognition58. More work is needed to determine if and how the limited diversity in T cell epitopes in the MTBC is linked to adaption to different host populations and/or immune escape.

Two independent phylogeographic analyses predicted Europe as the most likely geographical origin for the most recent common ancestor of L4.3/LAM. A European origin for L4.3/LAM was further supported by our finding that strains belonging to this sublineage were more genetically diverse in Europe compared to Africa, Asia and America. Taken together, these results suggest a role for Europeans for the spread of L4.3/LAM across the globe. Given the high frequency of L4.3/LAM in Europe (Fig. 2, Fig. 3a), particularly in TB patients from the Iberian Peninsula and in Latin America74,75, Portuguese and Spanish exploration, trade and conquest over the last centuries may have contributed to the global dissemination of this sublineage76.

Of note, the Americas lack specialist sublineages, including the three African specialist sublineages, despite centuries of slave trade. Importantly, this also applies to MTBC Lineage 5 and 6 (i.e. M. africanum) which today are largely limited to parts of West Africa11, the source of most of African slaves shipped to the Americas. Even if these lineages did reach the Americas at the time, they later might have been replaced by generalist sublineages from Europe including L4.3/LAM, following the massive influx of Europeans to the Americas during the 19 th and early 20 th centuries77, a time when the European TB epidemic was at its peak73,78. Importantly, these human migrations can be viewed as natural experiments, in which diverse human populations came into contact with different MTBC genotypes. As mentioned above, there is evidence that the African specialist sublineages of Lineage 4 already existed in sub-Saharan Africa centuries ago71. Following European contact, generalist sublineages were introduced to Africa and today, a significant proportion of human tuberculosis in Africa is caused by L4.1/LAM and other generalists (Figs. 2 and3a). By contrast, no significant spill-over of African specialist sublineages has occurred into Europe or American populations of European ancestry. Three of the 10 sublineages showed an intermediate pattern of geographical distribution.

Independent of the open question as to whether they might represent generalists or specialists, it is interesting to note that none of these three sublineages were found at significant frequency and proportion in Europe. They might therefore represent generalist of a non-European origin. Deeper analyses are needed to shed more light on these sublineages.

Our study is limited in that many of the MTBC isolates analyzed come from convenience samples and might therefore not be representative of a particular country. However, for the analysis of sublinage distributions by SNP-genotyping, we included more than 3,000 clinical isolates from 100 countries, which should reduce any potential selection bias. For the deep genomic analyses, we selected strains basesd on a large and diverse collection of classical genotyping patterns, and in addition, screened <13,000 MTBC whole genome sequences available in public repositories. As a further limitation, some isolates in our collection were obtained from patients who recently emigrated from a high tuberculosis incidence region into a low-incidence country. However, we excluded cases from ongoing transmission and focused on immigrants with reactivation disease, i.e. they were most likely infected in their country of origin before moving abroad. Moreover, we used country of birth for all analyses as opposed to country of tuberculosis diagnosis.

In conclusion, our findings indicate that the global success of Lineage 4 partly results from the different evolutionary strategies adopted by different sublineages. These strategies reflect an ecological distinction between specialists and generalists. The specialist sublineages are adapted to their sympatric host populations and geographically restricted. The generalist sublineages exhibit a broader ecological niche and are geographically widespread. Moreover, Europeans contributed to the global spread of the most successful generalist sublineage of Lineage 4. Our results highlight the ecological and epidemiological relevance of the deep phylogenetic diversity within the MTBC79. More generally, exploring potential differences between specialists and generalists in other pathogens will improve our understanding of the biology and epidemiology of infectious diseases.

available WGS data from several published studies8,45-56,75 and other whole genome data available in the public domain. We used KvarQ80 to screen for the L4.6.1/Uganda-specific SNPs described below. Starting from 13,067 genome sequences and excluding all clustered isolates except for one representative of each cluster, we identified 28 additional L4.6.1/ Uganda genome sequences which we included in our analysis of a total of 203 genomes (Supplementary Table 7). For genomic analysis of L4.1.2/Haarlem and L4.10/PGG3 strains, we screened the same 13,067 isolates (plus our own collection) for clade-specific SNPs of these two sublineages. We identified 505 genome sequences of strains of L4.1.2/Haarlem and 748 sequences of L4.10/PGG3. After excluding problematic sequences and strains in transmission clusters (criteria see above), we used 228 strains of L4.1.2/Haarlem and 301 strains of L4.10/PGG3. H37Rv was used as outgroup for all sublineage phylogenies except L4.10/PGG3, for which an isolate of L4.1.2/Haarlem was used (H37Rv belong to L4.10/ PGG3).

Whole genome sequencing, variant calling and filtering

WGS of new MTBC isolates was performed using Illumina chemistry (MiSeq, HiSeq2000/2500, NextSeq; paired end or single end). Illumina MiSeq-generated sequencing reads were clipped for adapters with Trimmomatic81 before mapping. We used a previously described pipeline for the mapping of short sequencing reads to the reference genome (a reconstructed hypothetical MTBC ancestor) with BWA 0.6.221. SNPs were called with SAMtools 0.1.19, and excluded if the coverage was less than 10% or more than 200% of the average coverage of the genome, if not supported by at least two reads on each strand, or if the quality was less than 30. All SNPs were then annotated using H37Rv reference annotation (AL123456.2) with Annovar82 and customized scripts. SNPs in regions annotated as "PE/PPE/PGRS", "maturase", "phage", "insertion sequence" were excluded. Additionally, we excluded SNPs in genes with previously identified repetitive regions58. Small insertions and deletions called by BWA/SAMtools as "INDEL" were not considered for the analyses. The presence of large genomic deletions reported previously15,28,74 was assessed by manually inspecting BAM alignment files from BWA mappings in Artemis for the presence of reads at the genomic regions with described deletions. Alternatively, we used a new testsuite in KvarQ80 to check for reads aligning to 25 bp query sequences of the corresponding deletion.

Phylogenetic and population genetic analyses for the definition of sublineages

A phylogenetic tree was generated with all Lineage 4 genomes, plus several reference genomes from all other MTBC lineages. Pairwise SNP distances were calculated using MEGA583 and the ape-package in R84. Fixation indices (F ST ; estimation of population separation) were calculated using Arlequin 3.5.85. Statistical significance was obtained by permutating haplotypes between sublineages. Principal Component Analysis (PCA) was done with Jalview86. Naming of Lineage 4 sublineages was adapted to Coll et al.27 whenever possible. However, in that publication, no criteria for the definition of sublineages were given, and not all sublineages were identified as such. We therefore added continuous numbers for the clades which were not defined by Coll et al. The new sublineages defined in this study are L4.1.3/Ghana and L4.10/PGG3 (the latter including L4.5, L4.8 and L4.9 according to nomenclature by Coll et al.). The full phylogenetic tree, including previous markers and spoligotyping family names is shown in Supplementary Fig. 1.

Identification of sublineage-specific SNPs

The alignment of all SNPs from the initial 72 MTBC Lineage 4 strains was imported into Mesquite87, in parallel with the phylogenetic tree generated from the same data in MEGA5. We used the "Trace Character History" module of Mesquite to map polymorphisms to clades. The full dataset of reconstructed positions was exported, and sublineage-specific SNPs were extracted as nucleotide differences between internal nodes of the phylogeny.

SNP-typing to screen for MTBC Lineage 4 sublineages

We developed a new SNP-typing assay to screen clinical isolates for the defined Lineage 4 sublineages. For this, we selected one "diagnostic" SNP per sublineage using previously defined methods and criteria36. Oligonucleotides were designed for a 10-plex MOL-PCR assay based on the Luminex xTag platform (Luminex, Austin, USA) (Supplementary Table 2)36. DNA extracts from clinical MTBC isolates were then screened with either i) the new MOL-PCR assay, ii) standard PCR amplification and subsequent Sanger sequencing of the region up-and downstream of the sublineage-specific SNP (see Supplementary Table 4 for PCR and sequencing primers), iii) a real-time PCR melting curve assay using the same SNPs (Supplementary Table 2), or iv) the MassARRAY platform (Sequenom, San Diego CA, USA) using phylogenetically redundant SNPs (Supplementary Table 3). The set of SNPs used in the MassARRAY typing scheme covered only six of 10 sublineages. Hence, all Lineage 4 isolates without any of the six SNPs with the mutant allele defined by MassARRAY typing (n=49) were subjected to the Luminex-assay described above. For all isolates, patient place of birth was used as country information. We obtained sublineageclassification data for 3,273 (97.2%) of a total of 3,366 isolates (Supplementary Table 5).

Spatial analysis and data presentation

For each country with Lineage 4 sublineage data available, sublineage proportions (compared to all Lineage 4 isolates from the same country) were calculated and mapped to a world map with ArcGIS ArcMap 10.0 (Esri, Redland, USA). A shapefile with country boundaries was used from DIVA-GIS, which is freely available. Categories for number of countries were defined as 0, 1-3, 4-10 and >10 countries. For individual sublineage "heat maps", countries with less than 3 isolates were not included. For the additional maps shown in Supplementary Fig. 6, sublineage proportions were normalized by the TB prevalence in the country as estimated by WHO38, and the area of the country. Other figures were generated with the ggplot2 library in R and GraphPad Prism 6.02 (GraphPad Software, San Diego, USA). Statistical analyses were performed with R or GraphPad Prism.

SIFT-analyses of functional effects of fixed sublineage SNPs

Analysis of SNPs fixed in each of the 10 sublineages were assessed for predicted functional consequence with the "Sorting Intolerant From Tolerant" (SIFT) in the software SIFT4G (v2.1)42 and the pre-compiled Mycobacterium tuberculosis database "GCA_000195955.2.22". Conservation levels of SNPs in the pre-compiled database had been obtained by comparing Mycobacterium tuberculosis H37Rv proteins to all proteins in the UniRef90 database. We pooled SNPs fixed in the generalist sublineages and the specialist sublineages, respectively, and excluded the L4.1.2/Ghana sublineage, as whole genome sequences of only two, very closely related isolates were available. Gene categories were analyzed based on the classification by Tuberculist88.

Phylogenetic reconstruction and population genetic analyses

The final alignment of polymorphic positions in all strains was used to estimate phylogenies with Bayesian methods using MrBayes 3.2.589 for L4.3/LAM and L4.6.1/Uganda sublineages (Fig. 6, Supplementary Fig. 10). For the Bayesian analysis we used a gamma rate distribution estimated from our dataset and a burn-in equal to 1/10 the number of generations; after the burn-in phase every 100 th tree was saved. Two parallel Markov chains were run in each of two runs. Tree length, log-likelihood score and alpha value of the gamma distribution were inspected for stationarity before termination of MrBayes. Trees were generated with standard parameters. A consensus tree was used for further analyses. Additionally, we used MEGA583 to generate Maximum Likelihood phylogenetic trees (Supplementary Figs. 9, 11, 12 and 13). We used the general time reversible (GTR) model of evolution, and 500 pseudoreplicates for bootstrapping confidence levels. Positions with gaps in more than 50% of taxa were ignored. Tree figures were generated using FigTree version 1.4.2. Pairwise SNP distances were calculated with the ape-package and the dna.dist function in R version 3.2.2, using raw counts of mutations and pairwise deletions for sites with gaps. For the comparison of pairwise number of SNP distributions overall (L4.3/LAM and L4.6.1/Uganda) and between continents for L4.3/LAM, a mean SNP distance to all isolates of the same population was calculated for each isolate, and a distribution of the mean pairwise distance plotted. Wilcoxon rank sum and Kruskall-Wallis tests were used to test for differences between continents as data were assumed to not be normally distributed. Average pairwise nucleotide diversities per site (π) were calculated as the average number of pairwise mismatches among a set of sequences divided by the total length of the interrogated sequences in base pairs (equation 4.21 in Ref. [START_REF] Hartl | Principles of population genetics[END_REF]). Confidence intervals for π were obtained by bootstrapping (1000 replicates) by re-sampling with replacement the nucleotide sites of the original alignments of polymorphic positions using the function sample in R. Lower and upper levels of confidence were obtained by calculating the 2.5th and the 97.5th quantiles of the π distribution obtained by bootstrapping. Code details are available upon request.

Antigenic diversity in human T cell epitopes

Experimentally confirmed human MTBC T cell epitope sequences were retrieved from the Immune Epitope Database on the 24 th of April 2015. Only linear epitopes from the MTBC (ID: 77643) tested in human T cell assays, with no MHC restrictions were selected (1,730 epitopes). The sequence of each epitope was blasted using blastP91 against the reference strain (H37Rv) to obtain genomic coordinates. Epitopes with no coordinates in H37Rv or for which no accurate coordinates could be determined (due to multiple hits) and epitopes in repetitive regions such as PE/PPE genes, phages-related genes and transposases were excluded, rendering a final set of 1,226 epitopes. Those epitopes are distributed across 304 antigens and have some overlapping sequences. In order to proceed with the sequence analysis, alignments were obtained by concatenating all epitope sequences after excluding sequence redundancy. Alignments of non-epitope containing antigens were obtained by excluding the regions described as epitopes from each respective antigen. To assess how other regions of the genome are evolving, alignments for essential and non-essential genes were also obtained62.

Alignments of epitopes and non-epitope containing antigens, essential and nonessential genes, were used to calculate pairwise dN/dS ratios for L4.3/LAM, L4.6.1/Uganda, L4.10/ PGG3 and L4.1.2/Haarlem sublineages. The dN/dS measures were calculated using all polymorphic sites within each sublineage and reflect therefore both within-sublineage substitutions and transient polymorphisms. Pairwise dN and dS values within each sublineage were calculated using the R package seqinr using the kaks function. To avoid having undetermined pairwise dN/dS values due to dN or dS being zero, a mean dN/dS was then calculated per sequenced isolate by dividing its mean pairwise dN by its mean pairwise dS with respect to all other sequenced isolates within each sublineage. The statistical differences between epitopes and non-epitope regions of antigens within each sublineage were accessed by using Wilcoxon rank sum tests with continuity correction implemented in R version 3.2.2.

Reconstruction of geographical origin of L4.3/LAM

The software RASP69 was used to reconstruct the hypothetical geographic origin of the MTBC L4.3/LAM ancestor genotype. The Bayesian phylogeny of 294 isolates (including H37Rv as outgroup) and the corresponding continent of birth of the patient were loaded as distribution. We used the S-DIVA (a parsimony based method) as well as the Bayesian Binary Method (BBM) implementation in RASP. A set of trees from MrBayes89 was used to correct for phylogenetic uncertainty in the S-DIVA analysis. Populations were defined according to country of birth of the patients and according to the United Nations definition. The isolates from Turkey, Libya, Algeria and Morocco were in the category "Europe and Mediterranean". RASP reconstruction was done without the outgroup (H37Rv). As we observed a single strain (from Ukraine) with a distinct, basal position in the phylogeny, we also performed a sensitivity analysis by excluding that isolate for the RASP analysis. With both methods, BBM as well as S-DIVA, the changes in proportions of continents were minor. With BBM, the proportion of "Europe/Mediterranean" for the "L4.3/LAM ancestor" decreased to 98.8%, and with S-DIVA, the proportion of "Europe/Mediterranean" decreased to 99.0% when excluding this basal isolate.

Figure 1 .

 1 Figure 1. Definition and global frequency of Lineage 4 sublineages. (a) We defined 10 sublineages based on the analysis of 72 MTBC Lineage 4 genome sequences published previously21,22. Sublineages were labeled according to Coll et al.27 (whenever possible) and previous designations based on spoligotyping (see Supplementary Fig. 1). Black triangles indicate sublineages identified as specialists, black circles indicate generalists. Filled shapes indicate sublineages, for which we performed deep genomic analyses. (b) Global proportion of each sublineage. A total of 3,366 MTBC Lineage 4 isolates were screened for sublineage-specific SNPs. L4.3/LAM was the most frequent sublineage globally.

Figure 2 .

 2 Figure 2. Global distribution of Lineage 4 sublineages. Pie charts showing proportions of the 10 Lineage 4 sublineages among all MTBC Lineage 4 isolates in each country. Circle sizes correspond to the number of isolates analyzed per country. A total of 3,366 MTBC Lineage 4 isolates were included. Color codes are as in Fig. 1.

Figure 3 .

 3 Figure 3. Country-specific proportions of sublineages reveal generalists and specialists. (a) The generalist sublineages L4.1.2/Haarlem, L4.3/LAM and L4.10/PGG3 were found globally at high proportions. (b) The locally restricted specialist sublineages L4.1.3/Ghana, L4.5, L4.6.1/Uganda and L4.6.2/Cameroon occurred at high frequencies in only a few countries and were restricted to certain geographical regions. Intensity of red indicates proportion of the sublineage among all Lineage 4 isolates in each country. Countries with fewer than three isolates in total are shown as "no data" and are filled white. A total of 3,366 Lineage 4 isolates were included in this analysis. The color scale for all sublineages is as indicated in Panel a, except for sublineage L4.1.3/Ghana (separate scale shown).

Figure 4 .

 4 Figure 4. Pair-wise ratios of rates of nonsynonymous to synonymous substitutions (dN/dS) in generalist and specialist sublineages for different gene categories. Abbreviations: Epi -experimentally confirmed human T cell epitopes; nEpi -non-epitope regions of T-cell antigens, both obtained from the Immune Epitope Database60; Essessential genes62; nEss -non-essential genes62. Wilcoxon rank sum tests: L4.6.1/Uganda (N=203) Epi vs nEpi, W=4952, p<0.001; L4.6.1/Uganda (N=203) Ess vs nEss, W=1415, p<0.001; L4.3/LAM (N=293) Epi vs nEpi, W=74540, p<0.001, L4.3/LAM (n=293) Ess vs nEss W=45067, p-value=0.29; L4.1.2/Haarlem (N=228) Epi vs nEpi, W=6561, p<0.001, L4.1.2/Haarlem (N=228) Ess vs nEss W=13369, p<0.001; L4.10/PGG3 (N=301) Epi vs nEpi, W= 27335, p<0.001, L4.10/PGG3 (N=301) Ess vs nEss W= 3103, p<0.001.

Figure 5 .

 5 Figure 5. Frequency distribution of the number of epitopes with nonsynonymous variants in generalist and specialist sublineages. A total of 1,226 T cell epitopes were included in the analysis. The number above each bar corresponds to epitope counts. Generalist sublineages L4.3/LAM, L4.1.2/Haarlem and (L4.10/PGG3. Specialist sublineage L4.6.1/Uganda. Tests: L4.6.1/Uganda vs L4.3/LAM Χ 2 = 27.04, p<0.001; L4.6.1/Uganda vs L4.1.2/Haarlem Χ 2 =15.75, p<0.001; L4.6.1/Uganda vs L4.1.2/PGG3 Χ 2 = 68.24, p<0.001.

Figure 6 .

 6 Figure 6. Genome-based phylogeny and diversity by continent of 293 strains of the L4.3/LAM sublineage. (a) Bayesian phylogeny with label colors indicating continent of strain origin: blue, Europe/ Mediterranean; red, Sub-Saharan Africa; yellow, America; pink, Asia. Numbers on nodes indicate posterior probabilities. Pie charts indicate reconstructed ancestral geographical regions of the internal nodes. The hypothetical L4.3/LAM-ancestor is labeled and a European origin for this ancestor was supported using a Bayesian Method (shown) and a Maximum Parsimony method (Supplementary Fig. 14). The pie colors correspond to the colors of the taxa labels. (b) Boxplot of pairwise genetic distances (number of polymorphisms) of L4.3/LAM strains by continent (p-values from Wilcoxon rank sum test). (c) Nucleotide diversity per site (π), measured by continent. Error bars indicate 95% confidence intervals. MTBC isolates from countries of the continent group "Oceania" (UN category; including Australia and New Zealand, Melanesia, Micronesia and Polynesia) were excluded for the genetic diversity analysis in panels B and C due the low number of samples.

Nat Genet. Author manuscript; available in PMC 2017 April 30.

Acknowledgements

We thank Sophie Lecher, Sarah Li, and Julia Zallet for technical support. Calculations were performed at sciCORE (http://scicore.unibas.ch/) scientific computing core facility at the University of Basel. This work was supported by the Swiss National Science Foundation (grants 310030_166687 (S.G.) and 320030_153442 (M.E.), Swiss HIV Cohort Study grant 740 (L.F.)), the European Research Council (309540-EVODRTB, S.G.), TB-PAN-NET (FP7-223681, S.N.), PathoNgenTrace (FP7-278864-2, S.N.) projects, SystemsX.ch (S.G.), the German Center for Infection Research (DZIF, S.N.), the Novartis Foundation (S.G.), and the Natural Science Foundation of China (91631301, Q.G.).

Data Availability Statement

All data generated or analyzed during this study are included in this published article (and its supplementary information files). Sequencing reads have been submitted to the EMBL-EBI European Nucleotide Archive (ENA) Sequence Read Archive (SRA) under the study accession number PRJEB11460.

Online Methods

Mycobacterial isolates

For the definition of Lineage 4 sublineages, we used 72 whole genome sequences of MTBC Lineage 4 and reference sequences of the other MTBC lineages published previously21,22 (Supplementary Table 7). These represented the largest collection of Lineage 4 wholegenome sequences available at that time. For the SNP-screening of clinical isolates for sublineage-classification, we used a retrospective global collection of 3,366 MTBC Lineage 4 isolates from 100 countries (Supplementary Table 5)15,30-35. All isolates had previously been identified as MTBC Lineage 4 by SNP-typing, genomic deletion analysis or spoligotyping. Approximately one third of these isolates were from patients who migrated to another country (1,106; 32.9%), and we used country of birth of the patient as a proxy for the origin of the MTBC strains. Two thirds of the isolates (2,260; 67.1%) were from countries where both country of isolation and country of birth were identical. Isolates of L4.6.1/Uganda from Uganda were genotyped in our previous work34. For the in-depth population genomic analysis of L4.3/LAM, we included previously published genomes21,27,44, and generated whole genome sequences of additional strains selected from a large collection of 2,132 MIRU-VNTR-genotyped isolates representing the global diversity of L4.3/LAM (Supplementary Fig. 7). Starting from 500 whole genome sequences, we excluded sequences with bad quality (sequencing coverage < 15x, proportion of homozygous variant calls <85%), isolates in transmission clusters (defined as isolate pairs differing by ≤12 SNPs) and strains with unknown country of origin, resulting in whole genome sequencing (WGS) data for 293 L4.3/LAM strains, which were included in the final analysis (Supplementary Table 7). For the in-depth population genomic analyses of L4.6.1/ Uganda, we generated WGS data from 175 isolates of the L4.6.1/Uganda genotype, selected for maximal geographic diversity and from previous studies34. Moreover, to further increase geographic coverage and genetic diversity among L4.6.1/Uganda strains, we analyzed all

Supplementary Material

Refer to Web version on PubMed Central for supplementary material.

Authors

David Stucki #1,2 , Daniela Brites #1,2 , Leïla Jeljeli #3,4 , Mireia Coscolla 1,2 , Qingyun Liu 6 , Andrej Trauner 1,2 , Lukas Fenner 1,2,5 , Liliana Rutaihwa 1,2 , Sonia Borrell 1,2 , Tao Luo 7 , Qian Gao 6 , Midori Kato-Maeda 8 , Marie Ballif 1,2,5 , Matthias Egger 5 , Rita Macedo 9 , Helmi Mardassi 4 , Milagros Moreno 10 , Griselda Tudo Vilanova 11 , Janet Fyfe 12 , Maria Globan 12 , Jackson Thomas 13 , Frances Jamieson 14 , Jennifer L.