Roseburia abundance associates with severity, evolution and outcome of acute ischemic stroke

MATERIALS AND METHODS
Study participants
This is a prospective observational cohort study. Patients with acute ischemic stroke were consecutively recruited from May 2018 to June 2019 with the following inclusion criteria: 1) aged 50 years or older; 2) local residents for over 6 months; 3) Magnetic Resonance Imaging (MRI)-confirmed ischemic stroke in the anterior circulation within 3 days of symptom onset; and 4) signed written informed consents. Exclusion criteria included: 1) cerebral hemorrhagic stroke; 2) a history of chronic inflammatory or immune diseases (e.g., rheumatoid arthritis, systemic lupus erythematosus, or inflammatory bowel disease); 3) a history of severe liver or kidney dysfunction, hematological diseases, and malignancies; 4) administration of probiotics, antibiotics, corticosteroids or immunosuppressants within the past 1 months; and 5) insufficient collection of fecal or blood samples.

Baseline characteristics and sample collection
We collected demographic information and medical histories from all participants by face-to-face interview. The etiology of ischemic stroke was classified by the Trial of Org 10172 in Acute Stroke Treatment (TOAST) criteria. Biochemical parameters including serum levels of total cholesterol (TC), high-density lipoprotein cholesterol (HDL), low-density lipoprotein cholesterol (LDL), fasting glucose, glycated hemoglobin, blood urea nitrogen (BUN), serum creatinine (Scr) and uric acid (UA) were collected after overnight fasting within 24 hours of admission and measured at the hospital central laboratory with laboratory staff blinded to clinical data. Stress hyperglycemia (SHG) was also included as a better biomarker of critical illness than absolute hyperglycemia. It was calculated using the following formula: fasting glucose/glycated hemoglobin ratio. Stroke severity was assessed by experienced neurologists on admission using the National Institute of Health Stroke Scale (NIHSS) score and retested at 24 hours, 3 days and 7 days. Patients were divided into two groups: minor stroke, who had admission NIHSS score ≤ 3, and non-minor stroke with admission NIHSS score > 3.

Sterile fecal containers and instructions were distributed to each study participant on admission. Approximately 2 g of fresh fecal samples were collected from each participant within 24 hours after admission and immediately (within 1 hour) stored at -80℃ until analysis.

Functional outcomes
Functional outcomes were quantified using the modified Rankin scale (mRS) score at 30 days and 1 year through routine telephone interview. Poor functional outcome was defined as mRS score > 2.

DNA extraction and high throughput sequencing
DNA extraction and sequencing were supported by the Shanghai Genesky Biotechnology Company (Shanghai, China) not knowing group assignment. According to the instructions, fecal genomic DNA was extracted from the fecal samples using the QIAamp® DNA Stool Mini Kit (Qiagen, Hilden, Germany). The V3-V4 hypervariable regions of the bacterial 16S rRNA gene were amplified by polymerase chain reaction (PCR) with the forward primer (5-CCTACGGGNGGCWGCAG-3) and the reverse primer (5-GACTACHVGGGTATCTAATCC-3). High throughput sequencing was performed on the Illumina Miseq platform using the 2×250 bp paired-end read protocol.

Bioinformatics and statistical analysis
The unique reads were clustered into operational taxonomic units (OTUs) by UPARSE with a 97% similarity cutoff. All OTUs were classified based on Ribosomal Database Project (RDP) Release 9 by Mothur. Within-individual (α) diversity (including observed species, Chao 1, ACE, Shannon, Simpson, and Coverage index) was used to measure the richness or evenness of taxa within each sample, and was analyzed by Mothur. Between-individual (β) diversity was provided for comparison of the taxonomic profiles between microbial communities. Unweighted and weighted UniFrac principal coordinate analysis (PCoA) based on OTUs were performed by R version 3.4.3 (Vegan package). Permutational multivariate analysis of variance (PERMANOVA; Adonis function) was carried out to examine whether there were statistical differences in bacterial community composition (β-diversity) between groups. Metastats analysis and linear discriminant analysis (LDA) effect size (LEfSe) were used to determine the significantly discriminative taxa between groups. Bacteria with significant differences (absolute value of logarithmic LDA score > 2) between the two groups were plotted on taxonomic bar plots. We also used BugBase to predict potential microbiome phenotypes, including aerobic, anaerobic, containing mobile elements, facultatively anaerobic, biofilm forming, gram-negative, gram-positive, potentially pathogenic, and stress tolerant.

The missing values of TC (1.5%), HDL (1.5%), LDL (1.5%), fasting glucose (3.7%), glycated hemoglobin ratio (3.7%), and UA (5.9%) were interpolated with the median. Propensity score-matched (PSM) analysis was used to obtain matched pairs of samples from the minor stroke group and the non-minor stroke group. In the PSM algorithm, the corresponding propensity score of the grouping variable (minor or non-minor) was calculated for each patient with a 1:1 nearest-neighbor matching algorithm with a caliper width of 0.2 of the propensity score, with age, sex, and coronary heart disease (CHD) as covariates. Spearman’s rank correlation coefficient was used to explore the correlation of different genera with biochemical parameters, NIHSS scores obtained at different timepoints and functional outcomes. We used linear mixed-effects models with random intercepts and slopes to test whether the relative abundance of discriminative taxa (e.g., genus Roseburia) or Firmicutes to Bacteroidetes ratio (F/B ratio) or gram-negative/gram-positive ratio account for the evolution of NIHSS scores through the first 7 days of hospitalization. Since the NIHSS score was highly skewed, the natural logarithm transformation [ln (NIHSS + 1)] was applied. Grand-mean centering for continuous covariates with meaningless 0 values (such as age) was performed. Multivariable logistic regression analyses were also used to evaluate the associations between the relative abundance of discriminative taxa and functional outcomes at 30 days and 1 year. The resulting p values were adjusted using the Benjamini-Hochberg false discovery rate (FDR) correction. Two-sided p value < 0.05 was considered significant.

材料与方法
研究对象
本研究为前瞻性观察队列研究。于2018年5月至2019年6月连续招募急性缺血性脑卒中患者，纳入标准如下：1）年龄≥50岁；2）本地居住时长超过6个月；3）症状发作3天内经磁共振成像（Magnetic Resonance Imaging, MRI）证实为前循环缺血性脑卒中；4）签署书面知情同意书。排除标准包括：1）出血性脑卒中；2）慢性炎症性或免疫性疾病病史（如类风湿关节炎、系统性红斑狼疮、炎症性肠病）；3）严重肝肾功能不全、血液系统疾病及恶性肿瘤病史；4）过去1个月内使用过益生菌、抗生素、糖皮质激素或免疫抑制剂；5）粪便或血液样本采集量不足。

基线特征与样本采集
通过面对面访谈收集所有参与者的人口统计学信息与病史。缺血性脑卒中的病因采用ORG 10172急性脑卒中治疗试验（Trial of Org 10172 in Acute Stroke Treatment, TOAST）标准进行分类。于入院后24小时内隔夜空腹后采集血清总胆固醇（total cholesterol, TC）、高密度脂蛋白胆固醇（high-density lipoprotein cholesterol, HDL）、低密度脂蛋白胆固醇（low-density lipoprotein cholesterol, LDL）、空腹血糖、糖化血红蛋白、血尿素氮（blood urea nitrogen, BUN）、血清肌酐（serum creatinine, Scr）及尿酸（uric acid, UA）等生化指标，由医院中心实验室中对临床数据设盲的实验人员完成检测。应激性高血糖（stress hyperglycemia, SHG）作为优于绝对高血糖的危重症生物标志物，通过空腹血糖/糖化血红蛋白比值计算得出。脑卒中严重程度由经验丰富的神经科医师于入院时采用美国国立卫生研究院卒中量表（National Institute of Health Stroke Scale, NIHSS）评分进行评估，并于入院后24小时、3天及7天复测。将患者分为两组：入院NIHSS评分≤3分的小卒中组，以及入院NIHSS评分>3分的非小卒中组。

入院时向每位研究对象发放无菌粪便容器及采集说明。每位参与者于入院后24小时内采集约2g新鲜粪便样本，并立即（1小时内）储存于-80℃直至后续分析。

功能结局
通过改良Rankin量表（modified Rankin scale, mRS）于30天和1年通过常规电话访谈量化功能结局。不良功能结局定义为mRS评分>2分。

DNA提取与高通量测序
DNA提取与测序工作由上海天昊生物科技有限公司（中国上海）完成，该公司对分组信息不知情。按照试剂盒说明书，使用QIAamp® DNA粪便微量试剂盒（Qiagen, 德国希尔登）从粪便样本中提取粪便基因组DNA。采用正向引物（5'-CCTACGGGNGGCWGCAG-3'）与反向引物（5'-GACTACHVGGGTATCTAATCC-3'），通过聚合酶链式反应（polymerase chain reaction, PCR）扩增细菌16S rRNA基因的V3-V4高变区。采用Illumina Miseq平台，以2×250 bp双端读长方案进行高通量测序。

生物信息学与统计分析
将唯一读段通过UPARSE以97%相似性阈值聚类为操作分类单元（operational taxonomic units, OTUs）。所有OTUs基于核糖体数据库项目（Ribosomal Database Project, RDP）第9版，通过Mothur进行分类。个体内（α）多样性（包括观测物种数、Chao1、ACE、Shannon、Simpson及Coverage指数）用于衡量每个样本内类群的丰富度或均匀度，通过Mothur完成分析。个体间（β）多样性用于比较不同微生物群落的分类学特征。基于OTUs的未加权与加权UniFrac主坐标分析（principal coordinate analysis, PCoA）通过R版本3.4.3（Vegan包）实现。采用置换多元方差分析（permutational multivariate analysis of variance, PERMANOVA；Adonis函数）检验各组间细菌群落组成（β多样性）是否存在统计学差异。采用Metastats分析与线性判别分析效应量（linear discriminant analysis effect size, LEfSe）确定组间具有显著鉴别性的类群。将两组间差异显著（对数LDA评分绝对值>2）的细菌绘制于分类学柱状图中。本研究还使用BugBase预测潜在的微生物组表型，包括需氧、厌氧、携带移动元件、兼性厌氧、生物膜形成、革兰氏阴性、革兰氏阳性、潜在致病性及耐应激。

TC（1.5%）、HDL（1.5%）、LDL（1.5%）、空腹血糖（3.7%）、糖化血红蛋白比值（3.7%）及UA（5.9%）的缺失值采用中位数插补法处理。采用倾向得分匹配（propensity score-matched, PSM）分析从非小卒中组与小卒中组中获取匹配样本对。在PSM算法中，以年龄、性别及冠心病（coronary heart disease, CHD）为协变量，采用1:1最近邻匹配算法，卡尺宽度为倾向得分的0.2倍，为每位患者计算分组变量（小卒中或非小卒中）对应的倾向得分。采用Spearman秩相关系数探索不同菌属与生化指标、不同时间点采集的NIHSS评分及功能结局的相关性。本研究使用带有随机截距和斜率的线性混合效应模型，检验鉴别性类群（如罗斯氏菌属（Roseburia））的相对丰度、厚壁菌门/拟杆菌门比值（F/B比值）或革兰氏阴性/革兰氏阳性比值是否可解释住院前7天内NIHSS评分的变化。由于NIHSS评分呈高度偏态分布，采用自然对数变换[ln(NIHSS + 1)]进行处理。对0值无实际意义的连续协变量（如年龄）进行总体均值中心化。还采用多变量logistic回归分析评估鉴别性类群的相对丰度与30天及1年功能结局的关联。所得P值采用Benjamini-Hochberg错误发现率（false discovery rate, FDR）校正进行调整。双侧P值<0.05被认为具有统计学显著性。