Stage II or III breast cancer patients who received NAT were enrolled in this study. The clinical and pathologic characteristics were reviewed retrospectively from medical records. The presence of estrogen receptors (ER), progesterone receptors (PR), and Her2 were determined by immunohistochemical staining. The ER or PR status was considered negative when less than 1% of the tumor cells showed positive staining. For Her2 staining, a score of 0 or 1+ was considered negative; specimens with a score of 2+ were further tested with fluorescence in situ hybridization analysis. The tumor histological grade was defined using the Nottingham combined histological grading system. This study was approved by the institutional review board (IRB number: 201704009RINC).

At the initial diagnosis (defined as before NAT), a 10-mL sample of blood was collected and stored in an EDTA-containing tube. Then, all patients were treated with NAT and received breast surgery. After NAT and breast surgery (defined as after NAT), another 10 mL of blood was sampled. Within three hours of blood sampling, the plasma was extracted after centrifugation at 1000× G for 10 minutes then stored at -80°C (21). Cell-free DNA was extracted using a QIAamp Circulating Nucleic Acid Kit (Qiagen, Germantown, MD, USA) according to the manufacturer’s protocol.

The library was constructed using a QIAseq Targeted DNA Panel with a customized gene list. The customized panel was designed to amplify the coding regions of the following genes: TP53, PIK3CA, Her2, GATA3, CDH1, PTEN, AKT1, ESR1, S100A7-9, ZNF703, B2M, CCND1, GATA3 and c-MYC. According to the manufacturer’s protocol, 10 ng of DNA was digested briefly into small fragments by a fragmentation enzyme at 32°C and 72°C. The DNA fragments were added to the QIAseq IL-N7 adapters, followed by target enrichment polymerase chain reaction (PCR) using the QIAGEN IL-Forward primer and the targeted DNA Panel primers. Finally, the library was amplified with universal PCR. The DNA library was then checked by using an Agilent Chip High Sensitivity DNA kit. KAPA library quantification kits were used to quantify the final concentration. The final DNA library was sequenced with the following Illumina platforms: Illumina MiSeq Reagent Kit v2, 2 x 150 bp reads or Illumina NextSeq 550 system Mid-Output Kit, 2 x 150 bp reads.

Previously, we have constructed an analytic pipeline of post‐NGS bioinformatics (22). First, BWA software (version 0.5.9) was used to align the raw sequencing data to the reference human genome [Feb. 2009, GRCh37/hg19; SAMtools (version 0.1.18)]. Picard (version 1.54) was used to perform the necessary data conversion, sorting, and indexing. GATK was used for variant calling with the Mutect2 and VariantFiltration parameters. Finally, ANNOVAR was used to annotate the genetic variants. Pathogenic and likely pathogenic variants were defined according to the American College of Medical Genomics and Genetics (ACMG) guidelines (23). The presence of ctDNA was determined by the presence of pathogenic and likely pathogenic variants, which are also considered tumor mutations. For variants of uncertain significance, if the prevalence of the variants in the normal population was less than 0.01 in a genomic database (1000 Genomics, ESP6500 and ExAC) and predicted to be deleterious by computer software (SIFT, PolyPhen2, and CADD), then they were classified as “highly suspected deleterious”. The above filtering analyses removes germline variants as much as possible (24); these variants are highly suspected to originate from tumors, so the detection of these variants could be considered indicative of ctDNA.

Since the Her2, c-Myc, CCND1 and S100A genes can be amplified in some breast cancer tumors, we decided to use copy number variations (CNV) to indicate the presence of ctDNA (25–27). Copy number variations were analyzed by OncoCNV (https://github.com/BoevaLab/ONCOCNV) according to the authors’ instructions. The baseline control consisted of the ctDNA BAM files of 14 healthy people. The ctDNA BAM files from the breast cancer patients were compared to the BAM files from the control population by using OncoCNV’s default cghseg segmentation algorithm (28). The sequencing region of each targeted gene was divided into several segments. When the mean of all segments of each gene was significantly different from the baseline, such as when the copy number predicted was greater than three copies or fewer than one copy from the baseline, we considered that to indicate a CNV alteration, which indicated the presence of ctDNA.

The chi-squared test and Fisher’s exact test were used to calculate the significance of the variance between each group. Survival was estimated by Kaplan-Meier analysis. Cox proportional hazards regression analysis was used to estimate the hazards ratios of RFS with a corresponding 95% confidence interval (CI) for various factors. All p values are two-sided, and p-values less than 0.05 were considered statistically significant.

First, to confirm the accuracy of the NSG-based deep sequencing, we checked whether this method could distinguish the existence of low-abundance mutants from background errors arising from the polymerase chain reaction (PCR) or sequencing process. We constructed a TP53 mutant (NM_000546.6: c.844C>A) as a reference sample; then we utilized this TP53 mutant with serial concentrations of 100%, 10%, 1%, and 0.1% to test whether the experimental method could detect these mutants at these concentrations ( Supplementary Methods ). The results demonstrated that the signal from the 0.1% mutant was significantly higher than background errors ( Supplementary Figure S1A ), suggesting that NGS testing accurately detected mutants present at 0.1%. In addition, the mutation level could be measured with a linear fashion (R² = 0.9997, Supplementary Figure S1B ).

Second, in deep cell-free analyses, another source of variants that makes it hard to distinguish cancer mutations is clonal hematopoiesis of indeterminate potential (CHIP) (29–31). The CHIP mutations mostly occur in the DNMT3A, TET2, PPM1D, ASXL1 and TP53 genes (29), whereas pathogenic variants of breast cancer were most prevalent in TP53, PIK3CA, MAP3KA1, CDH1, and PTEN (32). Variants most likely to be indistinguishable from CHIP were located in TP53. Twenty-two tumors from the pre-neoadjuvant core biopsy tumors were available for DNA extraction and sequencing ( Supplementary Table S1 ). Among them, 6 patients had TP53 variants, and their TP53 variants co-existed in the ctDNA and DNA from tumor biopsies ( Supplementary Table S1 and Supplementary Figure S2 ), suggesting the TP53 variants origin from breast cancer, not CHIP mutations.

A total of 95 patients were enrolled in this study. The median age was 50.0 years old. Forty-one patients had ER(+) Her2(-) breast cancer, 29 patients had Her2(+) breast cancer, and 25 patients had triple-negative breast cancer (TNBC). Before NAT, tumors with T1, T2 and T3-4 size classifications were found in three, 54 and 38 patients of each population, respectively. Eighty-two patients had positive axillary lymph nodes. According to standard clinical practice, ER(+) Her2(-) breast cancer patients with large tumors were treated with NAT. Out of the 95 patients, 77 patients received anthracycline while 80 patients received taxane in their NAT regimens. All Her2(+) patients received adjuvant anti-Her2 target therapy (27 patients receiving trastuzumab, one another receiving trastuzumab/pertuzumab and the other receiving trastuzumab-DM1). After NAT, 13 patients achieved a pCR of their primary breast tumors; 82 patients did not have pCR. Among the 13 pCR patients, there was one ER(+) Her2(-), six Her2(+) and six TNBC patients. The frequency of pCR was significantly higher in patients with Her2(+) breast cancer or TNBC than ER(+)Her2(-) patients (p = 0.002). CtDNA was detected in 60 patients before NAT and 31 patients after NAT. All of the clinical and pathologic characteristics are shown in Table 1 .

Among the 95 patients, 19 patients were found to have ctDNA before and after NAT; 41 patients had ctDNA only before NAT, 12 patients had ctDNA only after NAT, and 23 patients had ctDNA neither before nor after NAT ( Supplementary Table S2 ). The most common genetic variants were in the TP53 (n = 28), followed by PIK3CA (n = 16), CDH1 (n = 15), and Her2 (n = 7) genes. Eighteen patients had altered CNVs in their ctDNA, including of AKT1, CCND1, CDH1, c-MYC, Her2, PIK3CA, S100A, and ZNF703, either before or after NAT ( Supplementary Table S2 and Figure 1 ). Before NAT, Patient #73 ( Figure 1A ) and Patient #24 ( Figure 1B ) exhibited copy number gains of the S100A and Her2 genes in ctDNA, respectively; after NAT, the copy numbers of these genes in ctDNA returned to normal levels. Patient #3 ( Figure 1C ) had a new copy loss of the PTEN gene after NAT. We observed gains of Her2 and c-MYC in patient #27 ( Figure 1D ) before NAT that were only partially resolved after NAT.

Patients who had ctDNA before NAT tended to have a larger tumor size than those who did not have ctDNA before NAT (mean 5.0 cm vs. 4.3 cm, p = 0.104). However, the presence of ctDNA after NAT did not correlate with the tumor size or LN numbers after NAT. Although the difference was not statistically significant, patients with pCR had a lower detection of ctDNA after NAT than patients with no pCR (patients with pCR vs. absence of pCR: 15.4% vs. 35.4%, p = 0.132). Additionally, the presence of ctDNA was not correlated with the immunophenotype of breast cancer.

The median follow-up time of the entire cohort was 5.1 years, and the 5-year recurrence-free survival (RFS) was 58% (95% CI 48.0 – 68.0%). For clinical factors, Kaplan-Meier analysis showed that the residual tumor size after NAT and N classification after NAT were prognostic factors for RFS; patients who achieved pCR tended to have a better RFS than patients who did not achieve pCR ( Figures 2A–C and Table 2 ). On the other hand, patients with ctDNA after NAT had a significantly inferior RFS (p < 0.001, Figure 2D ). Other factors, such as age, ctDNA detection before NAT, immunophenotype, initial tumor size before NAT and N classification before NAT and adjuvant chemotherapy did not influence RFS. RFS was similar between patients with and without TP53, PIK3CA and CDH1 mutations ( Table 2 and Supplementary Table S3 ).

We then analyzed the clinical and pathologic characteristics of patients with and without ctDNA after NAT, and no difference was found between the two patient groups ( Supplementary Table S4 ). After incorporating the residual tumor size, N classification after NAT, pCR and ctDNA after NAT, multivariate analysis showed that an N3 classification and ctDNA positivity after NAT were independent risk factors that predicted tumor recurrence (N3, hazard ratio (HR) 3.352, 95% CI 1.267 – 8.870, p = 0.015; ctDNA, HR 4.135, 95% CI 2.014 – 8.491, p < 0.0001). Other factors did not significantly impact RFS ( Table 2 ).

Next, we analyzed the 72 patients with detected ctDNA, either before or after NAT. Patients with ctDNA positivity after NAT had a significantly inferior RFS compared to those without detectable ctDNA ( Supplementary Figure S3 , p<0.001). After adjusting for tumor size (after NAT), N classification (after NAT) and pCR, multivariate analysis with the Cox model revealed that ctDNA positivity after NAT was the most significant risk factor that predicted tumor recurrence (HR 8.02, 95% CI 3.24 – 19.86, p < 0.0001) ( Supplementary Table S5 ).

The median RFS of all the patients with ctDNA positivity after NAT was 1.19 years. When stratified by the immunophenotypes, ctDNA positivity after NAT was associated with a significantly inferior RFS for ER(+) breast cancer or TNBC patients and a trend of higher recurrence rates for patients with the Her2 subtype ( Figures 3A–C ). The median RFS of ER(+) breast cancer, Her2 (+) breast cancer and TNBC patients with ctDNA positivity after NAT were 0.90, 2.52 and 0.74 years, respectively.

For the entire cohort, the presence of ctDNA after NAT was a significant risk factor associated with recurrence in both patients who achieved and did not achieve pCR ( Figures 3D, E , all p < 0.001). Because pCR was previously reported as a surrogate marker for survival in patients with Her2(+) and TNBC (6), we analyzed these patient subgroups. Between the two patient populations, pCR was related to a trend of improved survival compared to absence of pCR (HR 3.328, 95% CI 0.777 – 14.243, p = 0.105, Supplementary Table S6 ). Multivariate analysis showed that advanced nodal status and ctDNA after NAT were independently correlated with high risk (N2-3, HR 3.753, 95% CI 1.146–12.297, p = 0.029; ctDNA, HR 3.123, 95% CI. 1.139 – 8.564, p = 0.027), and pCR status did show a not significant correlation with recurrence ( Table 3 ). A potential reason for this phenomenon is that pCR only represents the therapeutic efficacy of local breast tumor and the ctDNA may indicate that an occult lesion is present that is not effectively treated with NAT. In our study, 13 patients achieved pCR after NAT, and among those patients, two exhibited ctDNA positivity after NAT. One TNBC patient (case #50) received neoadjuvant docetaxel/epirubicin (four cycles) and achieved pCR for her primary breast and axillary tumors. However, she had hepatic metastases at 6 months after mastectomy ( Supplementary Figure S4 ). The other patient (case #5) had Her2-positive breast cancer and received neoadjuvant docetaxel/trastuzumab (four cycles) and epirubicin/cyclophosphamide (four cycles). The pathology showed no residual tumors. Trastuzumab was continuously maintained for one year. At the end of trastuzumab treatment (13 months after mastectomy), a cerebellar metastasis was found. The other 11 patients who achieved a pCR did not have ctDNA after NAT nor did they experience recurrence or metastasis.