Two HQ Y. pestis strains were collected in April, 2017 during daily surveillance, and eight more ones were isolated from rodents and their body fleas that were captured in the following large-scale field investigation (Guo et al., 2018). The mouse traps were set near the burrows of host animal, and in every day the traps were checked for 3 days, as well as to search for dead or diseased animals around the investigation area. After fumigating by diethyl ether, the body hairs of the captured animals were combed and body fleas were uncovered. Fleas were collected in normal saline, and their gastric contents were inoculated on plates. Blood samples from dogs and rodents were centrifuged to separate the serum for F1-Antibody detection by indirect hemagglutination test (IHA). Samples of liver and spleen from the captured live rodents, and samples of gland, heart, lung, liver, spleen, bone marrow from dead rodents were used to be cultured first and then be detected the F1 antigen by reverse indirect hemagglutination test (RIHA).

Ten Y. pestis isolates were collected from nine rodents and one parasitic flea in four villages in LJ during 2006–2018. Among the nine strains detected, seven strains were isolated from rodent corpses, and two strains were isolated from live rodents during the plague outbreak (Table 1).

We selected 10 HQ isolates (Y. pestis strains were isolated from seven live rodents and three fleas) and 10 high virulence Y. pestis strains that were isolated from LJ (the adjacent city) for DNA extraction. The extracted DNA was prepared for sequencing after recovery (the preserved Y. pestis strains were inoculated in added blood (rabbit) LB medium and cultured at 28°C for 24 h), activation, and proliferation (the activated strains were inoculated in the new medium and the culture conditions were same as before). A QIAGEN DNeasy Blood & Tissue (No. 69506) kit (QIAGEN Shanghai, China) and a Promega Wizard Genomic DNA Purification Kit (A1120) were used to extract the DNA for next-generation sequencing and third-generation sequencing, respectively.

Whole-genome sequencing was performed on an Illumina MiSeq platform with 150-bp paired-end libraries, and the raw sequencing reads were trimmed using Trimmomatic (Bolger et al., 2014). After filtering, clean reads with target mean coverage >80 × and quality score >20 remained. The clean reads were assembled de novo using SPAdes (Bankevich et al., 2012). Isolates LJ935 and HQ16 were also sequenced on a Pacific Biosciences (PacBio) RS II platform at the Beijing Genomics Institute (BGI, Shenzhen, China) to obtain the complete genome sequence. About 100,000 PacBio long reads were generated, which corresponds to 158 × coverage of the estimated Y. pestis genome. The average read length was 6,774 and 9,142 for LJ935 and HQ16, respectively. Sequence reads were assembled using Unicycler¹.

We sequenced the genomes of 10 HQ isolates and 10 LJ isolates (Yulong wild rodent plague focus), and compared them with 368 previously published genomes that represented the global diversity of Y. pestis (Table 1 and Supplementary Table 1). All assemblies were aligned against a reference genome of strain CO92 (accession no. NC_003143.1) using MUMmer (v3.0) (Kurtz et al., 2004) to generate whole-genome alignments and to identify single nucleotide polymorphisms (SNPs) in the core genome. The raw sequencing reads were mapped to the assemblies to evaluate the SNP accuracy using SOAPaligner (Li R. et al., 2009) as described previously (Yang et al., 2019). This process excludes unreliable SNPs located in repeated regions with low-quality scores (<20) or supported by few reads (<10 paired-end reads). Maximum likelihood phylogeny was inferred using PhyML (Guindon et al., 2010) from an alignment of core genome SNPs using the genome of strain CO92, which was the phylogenetically most related genome to our isolates at the time of analysis. The phylogeny was rooted using phylogroup 0.PE7 genomes as outgroups. Confidence in individual branching relationships was assessed using 100 bootstrap replicates and the tree was visualized using FigTree software². Because the core genomes in different datasets were variable, which could affect SNP calling, we recalled SNPs for the newly sequenced genomes using the same pipeline to obtain higher resolution. We identified 52 high quality SNPs from the datasets of the 20 new isolates. These SNPs were used as the character set in the Applied Maths Bionumerics 6.6 software³, which was used to construct a minimum spanning tree (MSTree) with hypothetical intermediate nodes. The MSTree was fully parsimonious for all 20 Y. pestis isolates and was confirmed using Y. pestis Z176003 (1.IN2 phylogroup strain, accession no. NC_014029.1) as the outgroup.

The 25 variable number of tandem repeat (VNTR) loci initially reported by Pourcel et al. (2004) were used to type the 20 newly isolates. Based on the VNTR locus information, we extracted the VNTR profiles of the HQ and LJ isolates and compared them with the VNTR profiles of three Xinjiang strains of 1.IN4 (XJ2651, XJ2650, XJ2656) reported previously (Zhang et al., 2018).

We aligned the 20 newly isolates to the selected reference genome, LJ935, using BWA software (v0.7.17) (Li and Durbin, 2010) to identify indels <50-bp long, and visualized them using SAMtools (Li H. et al., 2009). The indel loci and the 300-bp sequences flanking the indels were extracted and compared with the nucleotide sequences in the NCBI database⁴ to obtain the alleles. The annotation information was obtained according to the position of the alleles (Cui et al., 2020). All the assemblies were used in Prokka (Seemann, 2014) for gene annotation, and the annotation results (gff3 files) were used in Roary (Page et al., 2015) to identify the pan-genome and generate a matrix of the presence/absence of the genes in each sample as described previously (Yang et al., 2019). Gene rearrangements in the two complete maps (LJ935 and HQ16) and assemblies were identified by performing a multiple genome alignment with rearrangements using Mauve⁵ (Darling et al., 2004). Finally, a copy number variation (CNV) analysis was performed using the pipeline described in see section “Genomic Variation Detection.” Briefly, the paired-end reads were mapped to the LJ935 genome using BWA (Li and Durbin, 2009). The CNVs were predicted and corrected by CNGpro (v1.25.0) in each 100-bp window (Brynildsrud et al., 2015; Brynildsrud, 2018). In the data statistical analysis, all the results were classified by region (Lijiang or Heqing) and host state (dead or live rodent).

Total RNA was isolated from the four strains isolated from HQ (HQ16 and HQ32) and LJ (LJ236 and LJMS) using the Qiagen RNeasy^® Mini kit (No. 74104) and RNAprotect Bacteria Reagent (No. 76506), and purification was performed according to the manufacturer’s instructions. Sequencing libraries were generated using NEBNext^® UltraTM RNA Library Prep Kit for Illumina^® (NEB, United States) following manufacturer’s recommendations and index codes were added to attribute sequences to each sample. The prepared libraries were sequenced on an Illumina Novaseq platform and 150-bp paired-end reads were obtained. Read quality statistics were evaluated using FastQC⁶.

Paired-end raw reads were trimmed using Trimmomatic (Bolger et al., 2014). Aligned reads were annotated and counted using HTSeq (v0.6.1) (Anders et al., 2015). Differentially expressed genes (DEGs) between the two groups were detected using DESeq2 with the following cutoffs: Benjamini-Hochberg adjusted P-value (p.adjust) < 0.05 and absolute log₂(fold change) >1 (Deng et al., 2019). The Pearson correlation coefficient (R) of RNA-seq replicates was computed based on the DESeq2 regularized logarithm (rlog) normalized read counts (Love et al., 2014). The DEGs were annotated based on reference genome of Yersinia pestis YN1683 (accession no. NZ_ADTD00000000.1). Rockhopper (McClure et al., 2013) was used to identify new intergenic region transcripts, and Blastx was compared with the Non-Redundant Protein Sequence Database to annotate the newly predicted transgenic regions. The DEGs were functionally annotated by sequence homology searches against the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (Kanehisa et al., 2010) and the Swiss-Prot database⁷. In this study, we focused on quantifying only known mRNA sequences, so unknown genes or isoforms were not considered (Kim et al., 2020).

The health status of the captured rodents in HQ was judged by observing their manner, reaction, activity, and hair as described previously (Guo et al., 2018). During the dissection, the rodent lymph nodes, heart, lung, liver, and spleen tissues were examined pathologically, and subcutaneous hemorrhage, and pleural and peritoneal effusion were noted. The spleen and liver samples were sectioned and stained with hematoxylin and eosin and observed under a microscope.

Three isolates of HQ were selected randomly and compared with their phylogenetically nearest neighbors (three isolates from LJ), one strain from each of three other distinct Y. pestis lineages (Yunnan domestic rodent strain, Jianchuan wild rodent strain, and EV76 strain), and one Yersinia pseudotuberculosis strain. The following biochemical tests were performed: sugar–alcohol assays (maltose, glycerol, melibiose, arabinose, and rhamnose); nitrate reduction test; V virulence factor test, pigmentation (pgm) virulence factor test, and nutritional requirement detection. For the sugar–alcohol assays, Y. pestis strains were diluted and dripped into fermentation tubes, which were cultured at 37°C for 7 days. For the nitrate reduction test, the Y. pestis strains were cultured in nutrient liquid containing nitrate at 28°C for 48 h, and the biochemical reactions were observed after adding two Gliss reagents. For the V virulence factor test, diluted bacterial suspensions (10⁴ CFU per ml) were cultured in calcium-deficient medium or LB medium at 37°C for 72 h, then transferred to 28°C for 24 h. The determination of V was according to the number of bacterial colonies on the two mediums. For the pgm virulence factor test, the bacterial suspension was diluted to the desired bacterial titer (10³ CFU per ml) and evenly coated on Congo red plates, then incubated at 28°C for 3 days. For the nutritional requirement detection, the bacteria were incubated in glutamate-deficient medium or phenylalanine-deficient medium and cultured at 28°C for 4 days, then colony growth was observed.

The field investigation of the HQ plague epidemic was conducted from 17 April to 28 May 2017, which surrounded Damachang Village, HQ County within a radius of 5 km. Six villages were included and 7,090 individuals were checked through household inspection and quarantine, and no one with high fever or who was seriously ill was found (Guo et al., 2018). In total, 354 live rodents were captured and 242 body fleas were collected for Y. pestis detection, and no dead or diseased animals were found around the investigation area. Ten Y. pestis strains were isolated in HQ; six isolates were from Eothenomys miletus (E. miletus), one was from Apodemus chevrieri (A. chevrieri), one was from Ctenophthalmus quadratus, and two were from Neopsylla specialis specialis (N. specialis specialis). The surveillance found that the rodent plague had appeared only in the cultivated land, with a radius of 300 meters around the earliest isolation site (Figure 1). The surveillance of plague in Yunnan was performed according to standard protocol in National Scheme of Plague Surveillance of China. More details about the plague surveillance programmes are provided in http://www.gov.cn/yjgl/2005-08/30/content_28245.htm.

All the rodents captured in HQ were live rodents and they remained healthy with no obvious clinical signs of disease; that is, they exhibited normal daily activity, had shiny hair, bright eyes, and agile movements, which indicated a healthy state. Anatomical observation showed that the rodents had slight subcutaneous hyperemia, normal sized glands, and a small amount of red pleural and peritoneal effusion. The colors of the heart, lung, liver, and spleen tissues were normal, and there was no swelling, nodules, or necrosis in these organs. Three rodents positive for bacterial culture were selected randomly and their livers and spleens were sectioned. Microscopic examination of these sections revealed lesions in different levels of the liver tissue. Some lesions were lymphocyte infiltrations in portal areas and hepatocyte cords, and some were both lymphocyte infiltration and patchy necrosis of hepatocytes, and dilation and congestion of hepatic sinusoids (Figure 2). No obvious pathological changes were found in spleen sections. Pathological–anatomic observation of the liver and spleen tissue sections suggested that the Y. pestis isolated from HQ might have low virulence for E. miletus and A. chevrieri.

To define the phylogenetic position of HQ live rodent isolates in Y. pestis genealogy, we sequenced the genomes of 10 HQ isolates and 10 LJ isolates (Yulong wild rodent plague focus), and compared the sequences with 368 published genome sequences that represented the global diversity of Y. pestis (Table 1 and Supplementary Table 1). We built a maximum likelihood tree of 388 Y. pestis genomes on the basis of 5,458 SNPs. The HQ and LJ genomes grouped together in phylogenetic tree and formed a monophyletic lineage that was located between the known phylogroups, 1.ORI and 1.IN3 (Figure 3A). In a previously reported phylogenetic tree that was based on the diversity of 25 VNTR loci, three Y. pestis strains that were isolated from Xinjiang Autonomous Region, China also were located between 1.ORI and 1.IN3 and were attributed to 1.IN4 group (Zhang et al., 2018). However, by comparing the VNTR profiles of the HQ and LJ isolates with the VNTR profile of the 1.IN4 strains, we found at least 14 different VNTR loci among them (Supplementary Table 2). Therefore, the genetic evidence indicated that the HQ and LJ strains formed a new phylogenetic group, which was designated as 1.IN5. Notably, the 10 HQ isolates (all from live rodents) were clustered together (Figure 3B), composing a unique genotype within the 1.IN5 group. Although all 10 LJ genomes were in the 1.IN5 group, they were isolated mainly from rodent corpses, suggesting the LJ strains had higher virulence than the HQ strains.