A total of 4,584 IS cases were recruited from five different hospitals in Jiangsu province, China, which included Jurong People’s Hospital, Nanjing Jiangning Hospital, Yixing People’s Hospital, the Affiliated Hospital of Xuzhou Medical University, and Huai’an First People’s Hospital from 2009 to 2021. Included patients were diagnosed according to the computed tomography (CT), magnetic resonance imaging (MRI), or angiography with symptoms lasting more than 24 h. Excluded patients were those who had non-atherosclerotic ischemic stroke, acute coronary syndrome, tumor, autoimmune disease, or severe kidney failure according to the patient’s history of disease, admission diagnosis and discharge diagnosis. Besides, we use the TOAST (Trial of ORG 10172 in Acute Stroke Treatment) criteria to group the IS subjects. The TOAST classification denotes five subtypes of ischemic stroke: large-artery atherosclerosis (LAA), cardioembolism (CE), small-vessel occlusion (SVO), stroke of other determined etiology (SOE), and stroke of undetermined etiology (SUE). 4,663 controls free of stroke were randomly sampled from local communities and matched to the cases on age (±2 years) and gender. Subjects with acute coronary syndrome, tumor, autoimmune disease, and severe kidney failure were also excluded from the study. Additionally, we chose 314 pairs of IS cases and controls to detect THBS1 mRNA expression from 2019 to 2021. IS cases had a definite TOAST type and qualified retention samples. The controls were selected from the cohorts in the same or adjacent areas of IS cases using age- and gender- matching method. For above 314 IS patients, we collected the National Institutes of Health Stroke Scale (NIHSS) score and modified Rankin scale (mRS) score of IS patients when discharge, and followed-up their NIHSS scores and mRS scores at 1, 3, and 6 months after discharge. Demographic and clinical characteristics of the study subjects in above two parts were listed in Table 1.

This research incorporated two prospective cohort studies. The community-based cohort study was conducted from 2009 to 2022, which recruited 4128 participants from Guanlin Town and Xushe Town, Yixing city (Jiangsu, China). 4098 baseline subjects without stroke were followed up until May 25, 2022 for stroke onset. The detailed information about this cohort has been described previously (Dong et al., 2021).

The hospital-based cohort study enrolled 3158 IS patients aged between 35 and 80 years from Yixing People’s Hospital and followed an average of 5.86-year to record their long-term death outcome based on the annual death data from Yixing Center for Disease Control and Prevention. The follow-up period ended on May 25, 2022. The four endpoints were all-cause death, stroke death, IS death, and HS death. Finally, we observed a total of 488 deaths, among which 245 died from stroke and over half of the stroke-induced deaths were attributed to IS (n = 161). Meanwhile, we collected outcome of the long-term death from local Centers for Disease Control and Prevention (CDC) to assessment the prognosis of IS patients comprehensively, including all-cause death, stroke death, and IS death. Followed up to May 13, 2022, we observed a total of 28 deaths. Of the 20 deaths from stroke, 13 deaths were attributed to IS.

Demographic and clinical characteristics of the study population in the two cohort studies were summarized in Supplementary Table 2.

The case-control and cohort studies were all approved by the Research Ethics Committee of Nanjing Medical University (#2018571) and all participants signed the informed consent voluntarily.

Demographic characteristics, smoking and drinking habits, diseases history, and medication history were collected by questionnaire survey. All investigators were trained and qualified uniformly. Subjects who have ever smoked one pack of cigarettes a day in the past were considered smokers. Drinker was defined as individuals drinking alcohol at least three times daily. Physical examination including systolic blood pressure (SBP, mmHg) and diastolic blood pressure (DBP, mmHg) were measured at least three times.

Hypertension was defined as a self-reported history of hypertension, or elevated levels of blood pressure (SBP ≥ 140 mmHg, DBP ≥ 90 mmHg), or taking antihypertensive drugs recently. Those with fasting blood glucose level (GLU) ≥ 7.0 mmol/L, or self-reported history of diabetes, or taking hypoglycemic drugs were considered as diabetes. Dyslipidemia was defined as abnormal changes in lipid levels [total cholesterol (TC) ≥ 6.2 mmol/L, triglyceride (TG) ≥ 2.3 mmol/L, low-density lipoprotein-cholesterol (LDL-C) ≥ 4.1 mmol/L, high-density lipoprotein-cholesterol (HDL-C) < 1.04 mmol/L], or self-reported diagnosis of dyslipidemia, or currently taking lipid-lowering drugs.

Venous peripheral blood was collected from each subject after 8 h from the last meal into vacuum anticoagulation tubes with ethylene diamine tetraacetic acid dipotassium salt (EDTA-K2) and stored at –20°C. Leukocytes were obtained by gradient centrifugation from a subgroup of the study population. GLU was measured using the glucose oxidase method. TC was detected by the cholesterol oxidase-peroxidase method. TG was detected by the glycerophosphate oxidase-peroxidase method. HDL-C was detected by the catalase removal method and LDL-C was measured by sulfuric acid precipitation method.

The THBS1 gene (Gene ID: 7057, Locus NC_000015.10) locates on chromosome 15q14 and spans 18,388 bp, and consists of 22 exons. We searched the SNPs from the upstream 2 kb to the downstream 1 kb and selected tagging SNPs (tagSNPs) through the database of the Chinese Han population in Beijing (CHB) and China of the International Hap MAP Project. Three tagSNPs (tagSNPs rs2292305, etc. tagSNP rs2236741, and tagSNP rs2292304) would be available for candidate SNP selection. Included tagSNPs met the criteria of minor allele frequency (MAF) ≥ 0.05 and linkage disequilibrium (LD) r² ≥ 0.8. A functional candidate strategy was also applied to select potential functional SNPs on the bioinformatics effect prediction website SNPinfo Web Server.¹ The SNP function prediction results of rs2292305 showed no informed predictive biological function, therefore, we selected the closely linked tagSNP rs3743125 (r² = 0.932) with a prior predictive biological function as the substitute. We did not include rs2292304 because its probes and primers were not designed successfully. Finally, two tagSNPs of THBS1 gene, rs2236471 (C > T) and rs3743125 (G > A) were selected and genotyped in this study. Detailed biological information and function prediction were summarized in Supplementary Table 1.

DNA was extracted from the unfrozen venous peripheral blood using the protein precipitation method (Eaglink, EGEN2024, China) and then preserved at –20°C. Each DNA sample was quantified using the NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA). The polymerase chain reaction (PCR) TaqMan MGB probe assay was performed to amplify the two SNPs in the GeneAmp^® PCR system 9700 thermal cycle (Applied Biosystems, Foster City, CA). The results were post-read on the 7900HT real-time PCR system (Applied Biosystems, Foster City, CA) with the Sequence Detection System (SDS) 2.4 software. The successful call rates of two SNPs genotyping were both 100%.

White blood cells were separated from peripheral blood by gradient centrifugation and stored at –20°C. Peripheral leukocyte was reserved in RNA protective additive (Eaglink, EGEN2026, China) at –20°C. Total RNA was extracted using the Whole Blood RNA Extraction Kit (Yuan, Yu-BR02-1, China) and quantified using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA). Isolated RNA was reversely transcribed into cDNA using the PrimeScript™ RT Reagent Kit (Takara, RR047A, Japan). THBS1 gene mRNA expression of peripheral leukocyte was measured using SYBR Green quantitative real-time polymerase chain reaction (RT-qPCR), and the housekeeping gene Glyceraldehyde-3-Phosphate Dehydrogenase (GAPDH) was tested as an endogenous control. The qPCR reaction system included Platinum^® SYBR^® Green qPCR SuperMix-UDG (Invitrogen, 11733-046, USA), RNase free dH₂O, forward primer, reverse primer, and cDNA. Samples were incubated at 95°C for 5 min, followed by 40 cycles of 95°C for 10 s, 57°C for 20 s, and 72°C for 20 s on the QuantStudio™ 7 Flex Real-Time PCR System platform (Applied Biosystems, 4485700, CA). All samples were analyzed in three parallels, and cycle threshold values were recorded. The 2^–ΔΔCT method was used to calculate relative expression levels of THBS1 normalized by GAPDH. THBS1 mRNA’s forward primer sequence (5′–3′) was AGACTCCGCATCGCAAAGG, and the reverse primer sequence (5′–3′) was TCACCACGTTGTTGTCAAGGG. GAPDH mRNA’s forward primer sequence (5′–3′) was GGAGCGAGATCCCTCCAAAAT, and the reverse primer sequence (5′–3′) was GGCTGTTGTCATACTTCTCATGG.

We used EpiData 3.1 software (The EpiData Association, Odense, Denmark) for duplicate entry and consistency check of the collected data. Continuous variables were presented as median [inter-quartile range (IQR)] for non-parametric data. Categorical variables were presented as frequencies and percentages. For group-wise comparisons, the Kruskal–Wallis test or Mann–Whitney test was used for continuous variables with abnormal distribution. The chi-square test (χ²) was used to compare the differences in categorical variables between case and control groups. The Fisher’s exact test was used to estimate whether the genotype frequencies in controls and IS group met the Hardy–Weinberg equilibrium (HWE) law. Binary logistic regression was applied to calculate the odds ratios (ORs) and corresponding 95% confidence intervals (CIs) for the association of THBS1 variants and IS with adjustment for covariates (age, gender, smoking, drinking, hypertension, diabetes, and dyslipidemia). We used Cox proportional hazard regression to estimate the association with hazard ratios (HRs) and 95% CIs as well as after adjustment for covariates in the cohort study. Kruskal-Wallis H test was conducted to test the trend in mRNA levels among different groups of SNP genotypes. We used Spearman’s rank correlation to evaluate the correlations between THBS1 mRNA expression and NIHSS scores and MRS Scores after discharge in IS cases. We also used restricted cubic splines (RCS) with four knots at the 20th, 40th, 60th, and 80th centiles to flexibly model the association of THBS1 mRNA expression with the risk of long-term deaths in IS patients.

Haplotype association analyses as outlined by Schaid et al. (2002) were performed to test the associations of statistically inferred Haplotype with IS weighted with their estimated probability. All data analyses were carried out using SAS software 9.4 (SAS Inc., Cary, N.C., USA) and R 4.1.1 version.² A two−tailed P-value < 0.05 was considered statistically significant.

Table 1 shows the detailed demographic and clinical characteristics of 4,584 IS cases and 4,663 matched controls in the genetic case-control study. The median age of IS group (66 years) was comparable with that of the control group (66 years, P = 0.355). IS group had higher proportions of male (59.4%) and smokers (22.1%) while a lower proportion of drinkers (13.4%) than control group (41.6, 20.1, and 19.6%; P < 0.05). There were significant differences in SBP, DBP, GLU, TC, HDL-C, and LDL-C levels (P < 0.001) between the IS cases and controls. IS group had higher prevalence of hypertension (84.0%), diabetes (27.4%), and dyslipidemia (41.3%) than controls (50.1, 14.2 and 38.0%; P ≤ 0.001).

In the case-control study for transcriptome level analysis, IS group presented higher level of SBP but lower levels of GLU, TC, and HDL-C than control group (Table 1). Age and gender were both matched for case and control. IS group had lower proportions of smokers (15.3%) and drinkers (9.9%) than control group (26.1% and 29.9%; P ≤ 0.001). IS group had higher prevalence of hypertension (86.6%) while a lower prevalence of dyslipidemia (23.2%) than control group (72.9 and 51.6%; P < 0.001). Table 1 also listed the characteristics of participants in the cohort study for the risk of IS incidence. The clinical characteristic of the follow-up study for the long-term death after IS were presented in Supplementary Table 2.

As shown in Table 2, the genotype and allele distributions of SNP rs2236741 followed HWE both in the control and IS groups (both P-values > 0.05). Even though we have double-checked the genotyping results and controlled the quality by comparing the genotype frequencies of cases and controls of each plate, the allele frequencies of rs3743125 did not accord with the HWE in controls (P = 0.024) but accord with HWE in IS cases (P = 0.098).

No significant association of the two tagSNPs at THBS1 with IS was observed in the case-control study (Supplementary Table 3). The adjusted ORs (95%CIs) for the addictive model of rs2236741 and rs3743125 were 0.955 (0.872–1.047) and 0.955 (0.890–1.025) after adjustment for covariates. No significant association was observed for the dominant and recessive models of the two SNPs with IS (Table 2). Further stratification analyses by age, gender, smoking, drinking, hypertension, diabetes, and dyslipidemia did not reveal any significant association (Supplementary Tables 4, 5). There was no significant association between THBS1 variants and TOAST subtypes of IS in the case-control study (Supplementary Table 6).

Compared with C-G haplotype of rs2236741 and rs3743125, the haplotype C-A, T-G, and T-A were identified to have no significant association with IS (Supplementary Table 7), even after adjustment for covariates [Adjusted ORs (95%CIs): 0.974 (0.894–1.062), 1.103 (0.829–1.469), and 0.933 (0.845–1.030), respectively] (Table 3). In addition, the haplotypes C-A, T-G, and T-A were not associated with LAA and SVO (Supplementary Table 8).

After an average of 10.15-year follow-up for 4098 subjects free of stroke at baseline, 319 incident IS cases (10.30%) were observed and the incidence density was 63.08 (per 10⁴ person-year). The variants of rs2236741 and rs3743125 were not associated with the incident risk of IS (Supplementary Table 9), even after adjustment for covariates [Adjusted HRs (95%CIs) for the addictive model: 1.025 (0.821–1.281) and 0.984 (0.829–1.168)] (Table 4).

The variants of rs2236741 and rs3743125 were not associated with the long-term death after IS (Supplementary Table 10). Adjusted HRs (95%CIs) of the addictive model of rs2236741 and rs3743125 for all-cause death, stroke death, IS death, and HS death were 1.014 (0.848–1.212), 1.052 (0.820–1.349), 0.990 (0.723–1.356), 0.948 (0.553–1.626), and adjusted HRs (95%CIs) of the addictive model of rs3743125 for all-cause death, stroke death, IS death, and HS death were 1.125 (0.982–1.289), 0.985 (0.809–1.199), 1.030 (0.810–1.310), 0.887 (0.583–1.351), respectively. The dominant and recessive models of the two SNPs had no association with the long-term death either (Table 5).

The mRNA expression level of THBS1 in IS cases was approximately equal to that in controls (1.01 vs. 0.99, P = 0.833). Further subgroup analyses by age, gender, smoking, drinking, hypertension, diabetes, and dyslipidemia did not observed significant differences in THBS1 mRNA expression was detected between IS cases and controls (Supplementary Table 11). Additionally, no significant difference in THBS1 mRNA expression was detected between the three TOAST subgroup of IS and control group (Figure 1). As shown in Figure 2, the THBS1 mRNA expression levels significantly differed across rs3743125 GG, GA, and AA carriers in the control group (P = 0.016), but not in IS group (P = 0.454).

RCS regression analyses did not identify significant linear or non-linear correlation between THBS1 mRNA expression and the risk of all-cause death, stroke death, and IS death in IS patients (Figure 3). The adjusted HRs (95%CIs) for all-cause death, stroke death, and IS death were 0.962 (0.842–1.099), 0.986 (0.908–1.071), 0.975 (0.821–1.158), respectively (Supplementary Table 12). Furthermore, THBS1 mRNA expression has no significant relevance to NIHSS scores of discharges (ρ = –0.042, P = 0.462), 1 month after discharge (ρ = 0.064, P = 0.445), 3 months after discharge (ρ = 0.044, P = 0.601), and 6 months after discharge (ρ = 0.065, P = 0.441) in IS cases (Figure 4). Similarly, there is no correlation between THBS1 mRNA expression and mRS scores of discharges (ρ = –0.010, P = 0.864), 1 month after discharge (ρ = 0.030, P = 0.716), 3 months after discharge (ρ = 0.043, P = 0.608), and 6 months after discharge (ρ = 0.066, P = 0.434) in IS cases (Figure 5).