A Four-Factor Immunoscore System That Predicts Clinical Outcome for Stage II/III Gastric Cancer

Gastric cancer is one of the most common malignancies worldwide (1). Conventionally, tumor/lymph node/metastasis (TNM) staging has been used to assess the prognosis of gastric cancer patients; however, the clinical prognosis may vary significantly among patients with the same TNM stage due to its highly heterogeneous nature (2). Thus, researchers are intensely searching for complementary factors such as molecular signatures or genetic features (3). In addition to tumor-intrinsic factors, ongoing work is revealing extrinsic immune factors that affect tumor prognosis (4–7). However, due to the complicated characteristics of tumor immune response, the prognostic value of single immune factors are unreliable and have limitations when applied to large-scale populations. Thus, a comprehensive immunoscore system is needed to make more accurate prognoses for gastric cancer patients. A designated TNM-Immune (TNM-I) system has been defined in colorectal cancer by assessment of intratumor T cells, which is a prognostic factor superior to the American Joint Committee on Cancer (AJCC)/Union for International Cancer Control (UICC) TNM classification (8). However, an immunoscore system with prognostic significance for gastric cancer patients in clinical settings has not been identified so far.

The TNM-I system for colorectal cancer is based on density and location of CD3/CD8-positive cells, rather than on the function of T cells. However, the latter has an important impact on the prognosis of gastric cancer. A hallmark of tumor progression is the evasion of immune destruction, and PD-L1/PD-1 is one of the most important signaling pathways mediating tumor immune escape (9, 10). Pembrolizumab and nivolumab, monoclonal antibodies that block PD-L1/PD-1 interactions, are now used therapeutically for advanced melanoma and non–small cell lung cancer (NSCLC; refs. 11, 12). A completed phase Ib trial (Keynote-012) also showed an overall response rate (ORR) to pembrolizumab treatment was 33.3% among patients with gastric cancer (13). Because blockade of PD-L1/PD-1 can be effective in gastric cancer therapy, expression of PD-L1 and PD-1 was thought to have prognostic significance. However, the prognostic value of PD-L1/PD-1 remains controversial (14–23). One reason is that PD-L1 can be expressed by many cell types, such as tumor cells or immune cells; whether PD-L1 expressed by different cells has different prognostic values is unclear (14, 15, 18). Additionally, PD-L1 is regulated by intrinsic and extrinsic mechanisms; whether this has a different impact on the prognostic value of PD-L1 is not known (24). Gastric cancer is closely associated with inflammation, large numbers of preexisting CD8⁺ T cells are found at the tumor sites of gastric cancer patients. However, the prognostic effect of increased CD8⁺ T-cell infiltration in gastric cancer is still under debate (6, 7, 14, 25). These controversial viewpoints may be attributed to the complicated characteristics of tumor immune response. Therefore, a more reliable model incorporating multiple immune effectors is necessary and may provide better prognosis compared with single-factor predictors.

In the current study, expression of PD-L1 in tumor cells and immune cells, expression of PD-1 in immune cells, and infiltration of CD8⁺ T cells were evaluated in two independent patient cohorts with stage II–III gastric cancer and analyzed for clinical characteristics. We developed an immunoscore system based on the hazard ratios (HR) derived from the multivariable analyses. This approach can further divide patients with the same TNM staging into subgroups who are at low, moderate, or high risk, thereby providing more accurate prognoses and clues for the selection of immunotherapy options for gastric cancer patients.

This study was conducted on two independent cohorts of gastric cancer patients with stage II–III, who were followed up at the First Hospital of China Medical University. Patients were enrolled into two cohorts if they met the following criteria: Gastric cancer with stage 2–3; D2 lymphadenectomy; no neoadjuvant therapy before surgery; received 5-flurouracil-based standardized combination chemotherapy if relapsed. Newly collected specimen with more than 2 years of follow-up (from February 2012 to September 2012) were assigned to the discovery cohort (n = 51) and patients with 5 years of follow-up data (between April 2006 and July 2011) were selected according to the random number generated by the R/sample function to form validation cohort (n = 115). The median follow-up period was 31 months (range, 3–100 months), and 73 patients (43.98%) died during this period, with a 5-year overall survival (OS) probability of 50.1%. The study was approved by the Ethics Committee of China Medical University, and all procedures were conducted in accordance with ethical principles. Clinical information of all patients was retrieved from the Hospital Information System (characteristics are listed in Supplementary Table S1).

Standard indirect immunoperoxidase protocols were used for the immunohistochemistry assay. Briefly, embedded tumor tissues were sectioned to 4-μm thickness, dewaxed, and rehydrated. Heat-induced antigen retrieval in citrate buffer was conducted followed by 3% hydrogen peroxide. After blocked with bovine serum albumin (BSA), sections were incubated with indicated primary antibodies PD-L1 (#13684, clone E1L3N, CST), PD-1 (ZM-0381, clone MRQ-22, OriGene), CD8 (MAB-0021, clone c8/144B, Maixin Biotech) at 4°C overnight. Sections were further incubated with the UltraSensitive SP (Mouse/Rabbit) IHC Kit (9710, Maixin BiotechChina) and visualized by staining with 3,30-diamino-benzidine tetrahydrochloride (Maixin Biotech) and hematoxylin. For double-color IHC, sections were stained with DouMaxVision (BCIP/AEC) IHC Kit (Kit-9998, Maixin Biotech). Briefly, specimens were incubated with primary mouse antibodies PD-1 (ZM-0381, clone MRQ-22, OriGene) at 25°C for 60 minutes, and visualized by the AP (alkaline phosphatase)-BCIP/NBT substrate system. Subsequently, specimens were blocked with double-stain enhancer and incubated with primary rabbit antibodies CD8 (RMA-0514, clone sp16, Maixin Biotech) at 25°C for 60 minutes, and visualized by the horseradish peroxidase (HRP)–3-amino-9-ethylcarbazole (AEC) detection system.

All specimens were examined by two independent pathologists in a blinded manner based on staining percentage and intensity of positive cells. The staining pattern of PD-L1 was defined as negative or positive if <1% or ≥1% in TCs (13), respectively, or if <5% or ≥5% in ICs (26), respectively, according to reference. Staining of CD8 cells was defined as less or more if infiltration <20% or ≥20% in the tumor site; PD-1 was quantified by staining intensity in a 2-grading system as followed: PD-1^low if intensity = 0 or 1 (Fig. 1D and E), PD-1^high if intensity = 2 (Fig. 1F). The staining intensity of PD-1 was defined as 0 = negative, 1 = weak (with canary yellow staining), and 2 = intermediate–strong (with brown staining being similar to germinal center T cells in reactive tonsil tissues) according to reference (27, 28). For analysis of heterogeneity of these immune antigens, a series of optimal experimental processes have been carried out to reduce the deviation. First of all, H&E stains of a number of wax blocks from the same patient's sample were examined carefully by pathologist. The most representative block, which covered multiple heterogeneous regions, was selected to make a whole tissue section for IHC staining. Second, expression of immune antigens was reviewed in total tumor area and results were estimated as relative percentage staining and intensity staining. To minimize the impact of this spatial heterogeneity, tumor samples were collected from gastric cancer patients with stages 2 to 3 who underwent resection without receiving radiochemotherapy. Thus, the interference from late stage of disease and treatment of radiochemotherapy were excluded.

Fresh gastrectomy tissues were kept in precooled Hank's solution immediately, minced, and digested with medium containing collagenase IV (Gibco, 17104-019, 1 mg/mL) and DNase I (sigma, D4527, 200 IU/mL) at room temperature for 1 hour. After digestion, single-cell suspension was obtained by filtered and centrifuged. For flow cytometric analysis, cells were stained with the following Abs: anti-human CD3 (BD, 555332), anti-human CD4 (BD, 555349), anti-human CD8 (BD, 555369 or 300908), anti-human PD-1 (biolegend, 329906 or 329910), anti-human PD-L1 (BD, 557924) in BD Accuri C6. Isotype control Ab staining was used as a negative control.

Frozen specimens were sectioned to 4-μm thickness, fixed for 20 minutes with 4% PFA (paraformaldehyde). After being blocked with 5% BSA (bovine serum albumin), sections were incubated with indicated primary antibodies: rabbit anti-human CD8 (ab4055, Abcam), rat anti-human CD8 (ab60076, abcam), mouse anti-human PD-1 (ab52587, abcam), and rabbit anti-human PD-L1 (#13684, CST) at 4°C overnight, followed by sections-paired secondary antibody: Alexa Fluor Plus 488 goat anti-rabbit (A-11034, Thermo Fisher), Alexa Fluor Plus 555 goat anti-mouse (A-32727, Thermo Fisher), Alexa Fluor Plus 647 goat anti-rat (A-21247, Thermo Fisher) at 25°C for 2 hours. Isotype abs were used as negative control. Slides were mounted with DAPI mounting medium.

Characteristics of the patients were compared using a t test for continuous variables or a Spearman test for categorical variables. OS was calculated from the time of surgery till death or the last follow-up visit (August 11, 2014). Survival analysis was performed using a Kaplan–Meier method, and the differences were assessed by a two-tailed log-rank test. To evaluate the effect of the expression of PD-1, PD-L1, and CD8 on the OS, univariate and multivariate analyses using a Cox proportional hazard regression classification were carried out, and HRs were estimated with 95% confidence interval (95% CI) limits. The prognostic accuracy of our predict model compared with the TNM stage model was conducted by receiver operating characteristic (ROC) analysis.

Factors found to be associated with the OS by multivariate analysis were enrolled for the development of the prognostic system: TC PD-L1, IC PD-L1, IC PD-1, CD8, Lauren Classification, and AJCC stage (as listed in Table 1). Hazard ratios (HR) of each factor from multivariate analyses were used to derive coefficient scores calculated as follows: HRs of each prognostic factor were divided by the smallest one (1.72) and rounded to the nearest integer value (5). Coefficient scores of four immune effectors were then added up to assign each patient with an immunoscore index (from 0 to 5), on which the final system was based.

The clinical characteristics of two independent cohorts of gastric cancer patients with stage II–III are summarized in Supplementary Table S1. For the discovery cohort, 41 males and 10 females were enrolled with the median age at the diagnosis being 58 years (range, 41–76 years). The validation cohort has included 90 males and 25 females, with the median age at the diagnosis being 60 years (range, 25–79 years). The majority of tumors were located in lower or middle part of the stomach, only 10.24% were located in upper one-third part of the stomach. Together, 77.7% patients had lymphatic invasion, and 80.7% patients belonged to stage III according to the criteria released by the Seventh AJCC. The characteristics of patients in current study were consistent with those of our previous study (5) and those of Chinese gastric cancer patients (29).

We detected the expression of PD-L1 and PD-1 in ICs infiltrating the tumor site. As shown by IHC staining, most PD-L1⁺ ICs or PD-1^hi ICs were detected in IC aggregates infiltrating into the tumor tissues (Fig. 1B and F), whereas the remainder were scattered within the tumor region (Fig. 1C and G). Typical staining of PD-L1 in ICs appeared to be membranous with variable cytoplasmic staining (Fig. 1B and C), whereas staining of PD-1 in ICs was membranous-like (Fig. 1F and G). PD-L1⁺ ICs were observed in 33.33% (17 out of 51) and 17.4% (20 out of 115) of samples from the discovery and validation cohorts, respectively. Immune cells categorized as PD-1^hi were detected in 33.33% (17 out of 51) and 30.43% (35 out of 115) patients in the discovery and validation cohorts, respectively. Frequencies of immune-factor expression in tumor-infiltrating cells (Fig. 2F–H) were confirmed by flow cytometry (Fig. 2A–E) and in situ multicolor immunofluorescence staining by confocal analysis (Fig. 2I–M).

Patients with PD-L1⁺ ICs had shorter OS compared with those with PD-L1^– ICs (P = 0.0798 in the discovery cohort, P = 0.0191 in the validation cohort; Fig. 3A and D). In the validation cohort, the 5-year survival rate was 30% in the PD-L1⁺ IC group versus 51.18% in the PD-L1^– IC group, and the median survival was 22 months in the PD-L1⁺ IC group versus 80 months in the PD-L1^– IC group after 5 years of observation. Patients with PD-1^hi ICs had poorer OS compared with those with PD-1^low IC (P = 0.0012 in the discovery cohort, P = 0.0246 in the validation cohort), and the 5-year survival rate was 35.96% in the PD-1^hi IC group versus52.52% in the PD-1^low IC group in the validation cohort (Fig. 3B and E). In the validation cohort, median survival was not reached in the PD-1^low group, but was 33 months in the PD-1^hi group after 5 years of observation. No correlation was found between PD-L1 or PD-1 expressed by ICs and the clinical characteristics of the patients.

CD8⁺ T cells play a pivotal role in antitumor immunity; however, their prognostic value is controversial (6, 7, 14, 25). We detected infiltrations of CD8⁺ T cells in tumor tissues by IHC staining. The density of CD8⁺ T cells in patients varied widely (Fig. 1H and I). Infiltration of >20% of CD8⁺ T cells (referred to as CD8^More) was observed in 41.57% patients (25 out of 51) in the discovery cohort and 38.26% patients (44 out of 115) in the validation cohort. Patients with CD8^More tended to have better OS compared to those with CD8^Less, but the differences were not significant (P = 0.377 in the discovery cohort, P = 0.187 in the validation cohort; Fig. 3C and F), and the 5-year survival rate was 58.45% in the CD8^More group versus 40.71% in the CD8^Less group in the validation cohort. In the validation cohort, median survival was not reached in the CD8^More group, whereas the CD8^Less group had a median survival of 49 months after 5 years of observation. We found that patients with CD8^More/PD-1^low had significantly better survival than those with CD8^Less/PD-1^hi (Fig. 4A and B), which might explain the previously contradictory value of intratumor CD8⁺ T cells in the prognosis of gastric cancer. In addition, the positive staining areas of CD8 and PD-1 overlapped, which was confirmed by flow cytometry analysis on isolated tumor-infiltrating cells (Fig. 4C and D), IHC staining in tumor serial sections (Fig. 4E and F), double IHC staining (Fig. 4G and H), and in situ immunofluorescence staining of coexpression of CD8 and PD-1 (Fig. 4I–L). Together, these clinical data demonstrated that PD-1 expression could distinguish subsets of intratumor CD8⁺ T cells with possibly different function in gastric cancer, and the infiltration of CD8⁺ T cells with PD-1^hi was associated with worse prognoses. No correlation was observed between expression of CD8 protein and other clinical characteristics of patients.

Expression of PD-L1 in TCs of gastric cancer patients was on cell membranes, with variable cytoplasmic staining (Fig. 5A). Tumor cells were 31.37% and 41.74% PD-L1⁺ in the two cohorts, using the 1% standard threshold of PD-L1 in TCs, as was used to detect 40% positive PD-L1 expression in the Keynote-012 clinical trial. We found the typical distribution of PD-L1⁺ TCs (Fig. 5A, panel 2) to be restricted primarily to the lymphocyte infiltration region in gastric cancer (Fig. 5A, panel 3). Tumor cell expression of PD-L1⁺ was correlated with CD8^More, as confirmed by the Spearman test (Supplementary Table S2).

Patients with PD-L1⁺ TCs had significantly better OS compared with those with PD-L1^– TC (P = 0.0375 in the discovery cohort, P = 0.0029 in the validation cohort; Fig. 5B and C). In the validation cohort, the 5-year survival rate was 66.11% in the PD-L1⁺ TC group versus 33.88% in the PD-L1^– TC group. Median survival had not been reached after 5 years in the PD-L1⁺ TC group, whereas the PD-L1^– TC group had a median survival of 37 months. Also, patients with PD-L1⁺TC and CD8^More had the best OS (Fig. 5D), which represented an adaptive immune status in gastric cancer. In addition, PD-L1⁺ TC was correlated to better Borrmann classification and AJCC stage as shown by Supplementary Table S2.

Clinical characteristics and immune effectors of total 166 patients were included and subjected to univariable analysis. Lauren classification (diffuse type: HR = 2.099, P = 0.009; mix: HR = 2.838, P=0.004), differentiation (moderate: HR= 1.199, P = 0.751; poorly: HR = 2.215, P = 0.125), AJCC (stage III HR = 4.778, P = 0.001), TC PD-L1 (negative HR = 2.474, P = 0.001), IC PD-L1 (positive HR = 1.944, P = 0.011), IC PD-1 (positive HR = 2.236, P = 0.001), and CD8 (CD8^Less HR =1.513, P = 0.097) were selected for the multivariable analyses. Two tumor characteristics (Lauren classification, AJCC stage) and four immune variables (TC PD-L1, IC PD-L1, IC PD-1, and CD8) had independent prognostic value for the OS of gastric cancer patients, according to the results of the multivariable analyses (Table 1). However, differentiation was not found to be significant.

Although four factors—expression of PD-L1 on TCs or ICs, expression of PD-1 on ICs, and numbers of CD8⁺ ICs—had independent prognostic significance, their complicated interaction during antitumor immune responses prevented any one of these factors from providing accurate predictions of patients' survival. Therefore, we develop a comprehensive immunoscore system for the prediction of patients' survival. The coefficient of each immune factor was first calculated according to their HR values gained from the multivariable Cox hazard analysis (Materials and Methods; Table 1). Then each patient was assigned an immunoscore index by summing the HR values of four immune variables. Finally, all the patients were included in a Kaplan–Meier survival analysis, according to their immunoscores index. Patients could be separated into different risk groups (Fig. 6A). Among them, 87% with score 0 for the observation period were still alive. In contrast, 75% patients with a high score (3–5) died during the period. The immunoscore system could split patients into different risk groups who were in the same stage II or stage III (Fig. 6B and C), or even substage (Supplementary Fig. S1), which may improve the prognostic accuracy of TNM staging. To verify the accuracy of our four-factor prognostic model compared with the TNM staging prognostic model, ROC analysis was performed. Our four-factor model yielded an area under the ROC curve (AUC) of 0.856 for prediction of mortality at 5 years, which was superior to TNM staging with an AUC of 0.676 (Fig. 6D).

Gastric cancer is a highly heterogeneous cancer. Although patients can receive gastrectomy, the clinical outcome may vary significantly among patients with same stage. Thus, AJCC stage alone does not provide accurate outcome prediction for operable gastric cancer patients. Accumulating evidence has suggested that tumor progression is affected not only by its intrinsic characteristics, but also by extrinsic immune effectors (30). In addition to the infiltration of T cells, immune checkpoint molecules such as PD-L1 and PD-1 play a pivotal role in cancer (31). However, the impact of immune checkpoint molecules on prognoses for gastric cancer patients remains to be elucidated (14–17), due to not only the lack of uniform criteria, but also the complex characteristics of the antitumor immune response. In the current study, we analyzed the expression of PD-L1 and PD-1, the infiltration of CD8⁺ T cells, and the association between expression levels of these immune effectors and patients' OS, in order to build a comprehensive immunoscore system based on the HRs of immune variables with prognostic significance. Our prognostic model yielded an AUC of 0.856 for prediction of mortality at 5 years, which was superior to TNM stage with an AUC of 0.676.

An immunoscore system developed by Galon J and colleagues has demonstrated that immune cell infiltration, compared with intrinsic tumor characteristics, has a better prognostic value for colorectal cancer (8, 30). However, only the amount and location of CD3/CD8 lymphocyte infiltration, rather than the function of these immune cells, was included in this system. Another five-feature–based immunoscore is available to predict the recurrence and survival of gastric cancer. Although CD45RO and CD66b were added to the system, critical factors that could reflect the interaction between tumor cells and immune cells were not included (32). In addition, the tumor microenvironment (TME) was divided into four stratifications according to the expression of PD-L1 in tumor cells and presence of TILs, but the prognostic value of this system was unclear for gastric cancer patients (33). In contrast, our immunoscore system contained not only the amount but also function of the intratumor immune cells that could reflect interaction between tumor cells and immune cells, and it could fill in several gaps of the TME classification. First, we noted that the expression of PD-L1 in TC was derived from the infiltration of CD8⁺ T cells, and patients with PD-L1⁺CD8⁺, an adaptive subgroup, had the best prognosis in gastric cancer. In contrast, PD-L1⁺ ICs were found to be associated with worse survival, rather than PD-L1⁺ TCs, which was therefore more critical for predicting response to immune checkpoint inhibitors. Second, PD-1 should be combined with CD8 when analyzing prognostic significance. Lastly, although PD-L1, PD-1, and CD8 had independent prognostic significance in gastric cancer, single factors were not sufficient for accurate prognoses of cancer patients—the comprehensively immunoscore system was necessary.

Expression of PD-L1 in TCs was regulated by two major mechanisms (34): extrinsic immune-induced PD-L1 expression (24) and intrinsic oncogenic activation mechanism (35–38). As the majority of gastric cancer are associated with inflammation (3), TC PD-L1 expression in gastric cancer could be upregulated by extrinsic inflammatory factors. Consistently, our IHC results showed that PD-L1 expression in TC is mainly detected in areas surrounded by gatherings of immune cells in gastric cancer tissues, with correlations between TC PD-L1 and more infiltration of CD8⁺ T cell confirmed by the Spearman test. Patients with PD-L1⁺ tumor cells and CD8^More T cells have the best prognosis compared to other groups, indicating an adaptive immune status existing in gastric cancer patients, consistent with reports in melanoma and hepatocellular carcinoma (39, 40), although PD-L1 has been reported to correlate with the worse survival of gastric cancer patients (14). A possible explanation for the contrary results is the different patient populations enrolled. We enrolled 166 gastric cancer patients with stages II–III who underwent resection without receiving preoperative neoadjuvant therapy, thus PD-L1 staining were not influenced by the chemotherapy. In Thompson's study, however, only 34 patients with gastric cancer or gastro-oesophageal junction (GEJ) cancer were enrolled, and 44% patients had received neoadjuvant chemotherapy and/or radiation before the detection of PD-L1. Moreover, PD-L1 expression in TC was reported to correlate with Epstein–Barr virus or microsatellite instability in previous studies (41), both of which are associated with better survival in gastric cancer (42, 43). Additionally, we also detected the PD-L2 level in 166 gastric cancer patients; however, its prognostic value was not significant.

PD-L1 is more commonly expressed in TILs than that in TCs, and PD-L1 expressed in TILs, but not in TCs, was relevant to response to the PD-L1 mAb (MPDL3280A) treatment in a range of cancers (44). In our study, PD-L1⁺ ICs, rather than PD-L1⁺ TCs, correlated with poorer OS. Immune cells positive for PD-L1 could be dendritic cells (DC) or macrophages in the tumor microenvironment, which exert immunosuppressive function (45, 46). Consistently, in the Keynote-012 clinical trial, overall responses are higher in patients with PD-L1⁺ ICs, whereas the PD-L1⁺ TCs were not useful in predicting the response, which indicated that PD-L1⁺ ICs were more critical to the suppression of antitumor immunity. As for the prognostic impact of PD-1, it was correlated with worse prognosis in Japanese patients with stage II/III gastric cancer, but associated with better outcomes in Caucasians with stage I–IV gastric cancer (15, 17). Here, we demonstrated for the first time that patients with high intensity of PD-1 correlated with poorer survival in Chinese gastric cancer patients, whereas there was no correlation between PD-1 and patient's survival if patients were divided by positive and negative expression of PD-1. Of note, our result showed that patients with CD8^MorePD-1^low had significantly better OS compared with those with CD8^lessPD-1^hi, which may explain the controversial prognostic prediction of CD8 in gastric cancer (6, 7, 14, 25). CD8⁺ T cells with PD-1^hi/CTLA4^high are partially exhausted T cells, and strongly correlate with response to PD1 mAb therapy (47). This provided molecular evidence that the intensity of PD-1 expression could impact the function of intratumor CD8⁺ T cells, which supported our conclusion that CD8 and PD-1 should be combined to evaluate their prognostic significance for gastric cancer patients.

Taken together, we built an immunoscore system based on hazard ratios of four immune variables. This comprehensive immunoscore system can improve the accuracy of TNM staging for predicting survival and is an essential complement for the AJCC staging system for operable gastric cancer patients with stages II–III. Future studies should measure other critical immune markers for inclusion in the system, such as Bim and colleagues (48). Further development of immune signatures will facilitate choosing appropriate immunotherapy options for gastric cancer patients.

No potential conflicts of interest were disclosed.

Conception and design: T. Wen, Z. Wang, Y. Li, X. Qu, Y. Liu

Development of methodology: T. Wen, X. Yang

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): T. Wen, Z. Wang, Y. Li, Z. Li, S. Wang, J. Qu, K. Hou, L. Xu, C. Li, J. Wang, J. Liu, L. Chen

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): T. Wen, Y. Li, Z. Li, X. Che, Y. Fan, J. Qu, K. Hou, W. Zhou, X. Qu, Y. Liu

Writing, review, and/or revision of the manuscript: T. Wen, Z. Wang, Y. Li, X. Che, X. Qu, Y. Liu

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): T. Wen, Z. Wang, J. Zhang, Y. Liu

Study supervision: T. Wen, Y. Li, X. Che, X. Qu, Y. Liu

This work was supported by the Key National Science and Technology Major Project of China (No. 2013ZX09303002 to Y.P. Liu), the National Natural Science Foundation of China (Nos. 31300743 and 31770963 to T. Wen, 81673025 to L. Xu, 81602098 to Y.B. Fan), Science and Technology Plan Project of Liaoning Province (No. 2015020734 to T. Wen, 2014225013 and 2014226033 to Y.P. Liu), the Key Laboratory Project of the Education Department of Liaoning Province (No. LZ2014037 to Y.P. Liu), The National Key Research and Development Program of China (No. 2017YFC1308900 to Y.P. Liu), National Science and Technology Major Project of the Ministry of Science and Technology of China (No. 2017ZX09304025 to X. Qu), and the Health and Family Planning Commission of Liaoning Province (No. LNCCC-D07-2015 to T. Wen).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.