Detection of SARS-CoV-2 antibodies in serum and dried blood spot samples of vaccinated individuals using a sensitive homogeneous proximity extension assay

A homogeneous PCR-based assay for sensitive and specific detection of antibodies in serum or dried blood spots (DBS) is presented and the method is used to monitor individuals infected with or vaccinated against SARS-CoV-2. Detection probes were prepared by conjugating the recombinant spike protein subunit 1 (S1), containing the receptor binding domain (RBD) of SARS-CoV-2, to each of a pair of specific oligonucleotides. The same was done for the nucleocapsid protein (NP). Upon incubation with serum or DBS samples, the bi- or multivalency of the antibodies (IgG, IgA or IgM) brings pairs of viral proteins with their conjugated oligonucleotides in proximity, allowing the antibodies to be detected by a modified proximity extension assay (PEA). Anti-S1 and anti-NP antibodies could be detected simultaneously from one incubation reaction. This Antibody PEA (AbPEA) test uses only 1 µl of neat or up to 100,000-fold diluted serum or one ø1.2 mm disc cut from a DBS. All 100 investigated sera and 21 DBS collected prior to the COVID-19 outbreak were negative, demonstrating a 100% specificity. The area under the curve, as evaluated by Receiver Operating Characteristic (ROC) analysis reached 0.998 (95%CI: 0.993–1) for samples taken from 11 days after symptoms onset. The kinetics of antibody responses were monitored after a first and second vaccination using serially collected DBS from 14 individuals. AbPEA offers highly specific and sensitive solution-phase antibody detection without requirement for secondary antibodies, no elution step when using DBS sample in a simple procedure that lends itself to multiplex survey of antibody responses.

with rates of asymptomatic casesinformation that is vital for developing policies to stem transmission [1]. Since the outbreak of the COVID-19 pandemic in 2019, there has been a rapid development of assays to assess immunity to this severe acute respiratory syndrome, caused by SARS-CoV-2. Most assays involve a solid-phase, such as in lateral flow assays (LFA), enzyme linked immunosorbent assay (ELISA), and chemiluminescence immunoassay (CLIA), a form of ELISA. ELISAs are the most commonly used methods for antibody detection in clinical laboratories and generally rely on dual-binding of patient antibodies, using virus antigen for capture and labeled secondary antibodies for detection. Among their advantages are the ability to detect both titers and isotypes of the antibodies. Frequently encountered problems in ELISA are weak signals and limited dynamic ranges. CLIA assays have excellent sensitivity with high signal-to-noise ratios. The Abbott (IL, USA) Architect SARS-CoV-2 IgG test, targeting antibodies directed against NP, is one of the most commonly used CLIAs, using magnetic beads coated with recombinant viral antigen and a luminescent substrate for automated analysis. The sensitivity for detection of virus-specific IgG reached 100% at 17 days after symptom onset or 14 days after PCR positivity for the virus [2]. Meso scale discovery (MD, USA; MSD) offers one of the most sensitive methods with a specificity of 97.4% and a sensitivity of 96.2% for samples collected 14 or more days after onset of symptoms [3]. Detection of antibodies directed against spike (S) or subunit 1 (S1) proteins is made after 50,000-fold dilution of serum or plasma.
Homogenous proximity extension assays (PEA) are protein detection reactions that differ in several respects from more standard immune assays for protein detection [4]. Like sandwich ELISAs, PEA depends on target recognition by pairs of affinity reagents, but no solid phases are needed for capture and no washes are required. Proteins are accurately quantified from levels below pg/ml and over wide dynamic ranges. Typically sample volumes of only 1 µl or less are used for multiplex assays of a hundred or more target proteins [5]. The assays generally use pairs of oligonucleotide-conjugated antibodies that are brought into proximity by binding the same target protein molecules in a sample. This induced proximity allows oligonucleotides to produce amplifiable DNA strands via DNA polymerization. The amplicons then serve as surrogate markers for the targeted proteins in real-time PCR or DNA sequencing reactions. By conjugating barcoded oligonucleotides to pairs of antigens rather than to target-specific antibodies, the assay technique was adapted for detection of specific antibody responses against those antigens, as antibodies in the sample will have the effect of bringing pairs of antigen molecules into proximity in an approach termed AbPEA. A related solution phase assay, agglutination-PCR (ADAP), has been developed for detection of total levels of SARS-CoV-2 antibodies [6]. This assay consists of four reaction steps and a very good sensitivity of 98% and specificity of 99.55% was reported.
Serological testing for SARS-CoV-2 antibodies typically uses serum or plasma samples, collected by venipuncture. Collection of DBS by finger prick has proved to be a useful means of monitoring SARS-CoV-2 antibody levels, as samples can be collected at home and sent by post for analysis [8,9]. The Siemens (Munich, Germany) SARS-CoV-2 total antibody assay achieved a 100% and 94% agreement between results for plasma and DBS samples for negative and positive patient samples, respectively. The assay used 6 mm punches from the DBS [10]. A concordance of 98.4% was shown between DBS and venous plasma samples using DBS sample volumes of 150 µl blood by the Architect SARS-CoV-2 IgG assay [11]. Similarly, there was 96.8% and 81.3% agreement for positive and negative samples, respectively, comparing serum and DBS with 3 mm punches by the Euroimmun (Lübeck, Germany) ELISA [12]. The rapid semi-quantitative Roche (Basel, Switzerland) Elecsys Immunochemistry assay can reveal SARS-CoV-2 antibodies from DBS samples in 18 min, but the sensitivity for detection of infected individuals was only 71.4% [2,13]. In contrast, the ADAP method achieved 100% sensitivity and specificity using DBS samples, but antibodies had to be eluted from the DBS for 90 min and then concentrated [7]. The elution of antibodies from the DBS samples is a step that is in common for the published assays and takes more than 1 h.
The AbPEA approach, presented herein, consists of three simple reaction steps of incubation with liquid or dried sample, extension/ prePCR and readout via real-time PCR without washing steps. The reaction uses 1 µl of diluted or undiluted patient serum samples or one ø1.2 mm cut from a DBS to reveal the immune status of sampled individuals. The high specificity and sensitivity, as well as broad dynamic range, of this detection method are demonstrated using sets of prepandemic and patient serum samples, as well as DBS samples collected over time from vaccinated individuals.

Production of recombinant SARS-CoV-2 spike receptor-binding domain (RBD)
To express recombinant S1-RBD, a codon-optimized cDNA, encoding amino acids 319-541 of the SARS-CoV-2 spike protein fused to the Nterminal signal peptide of immunoglobulin kappa, was synthesized by Eurofins Genomics (Ebersberg, Germany). C-terminal sortase and hexahistidine tags were inserted by PCR with nested primers (see Supplementary Table S3 for full protein sequence). The cDNA was cloned downstream of the CAG promoter in an expression vector containing Epstein-Barr Virus-derived sequences for episomal replication using Gibson assembly [14] using the HiFi DNA Assembly kit (New England Biolabs,Ipswich, MA, USA).
Recombinant S1-RBD was produced in 293 F cells (Catalog number: R79007, Thermo Scientific,Waltham, MA, USA) cultured in Freestyle 293 medium and transfected at 2 × 10 6 cells/ml, using FectoPro transfection reagent (Polyplus transfection, Illkirch, France) in a volume of 200 ml in a 1000 ml E-flask at 130 rpm in a Multitron 4 incubator (Infors, Velp, the Netherlands). Transfected cells were kept at 37ºC and the supernatant was harvested after 72 h and filtered using Polydisc AS 0.45 µm (Whatman; Cytiva, Marlborough, MA, USA) and loaded onto a 5 ml HisTrap Excel column (Cytiva). After sample loading, the column was washed with 10 column volumes of 500 mM sodium chloride, 30 mM sodium phosphate buffer with 30 mM imidazole before elution of the protein using the same phosphate buffer with a 15 column volumes linear gradient between 30 and 300 mM imidazole. Pooled fractions were concentrated using 10 kDa Vivaspin concentrators (10 kDa MWCO; Sartorius, Göttingen, Germany), passed over a HiPrep 26/10 desalting column (Cytiva) in phosphate-buffered saline and finally concentrated again. Integrity and purity of the RBD preparation was checked by Coomassie-stained SDS-PAGE ( Supplementary Fig. 1).

Patient serum samples
Two sets of serum samples were collected. The first consisted of 20 positive and 20 negative sera, previously analyzed by the Abbott SARS-CoV-2 IgG test targeting antibodies directed against NP (Catalog number. G92798R02). This cohort was used to develop the assay. A second set of samples consisted of 100 negative samples collected in 2019 before the COVID-19 pandemic, and 94 serum samples collected between March 1 and November 30, 2020. The latter samples were from individuals with symptoms of COVID-19 and diagnosed with SARS-CoV-2 infection by reverse transcription PCR (RT-PCR) or by characteristic symptoms. In addition, patient samples had been analyzed using the Abbott and ELISA tests for SARS-CoV-2 IgG antibodies at the Department of Clinical Microbiology, Laboratory medicine, Norrlands University hospital-region Västerbotten, Umeå University, Sweden. These patient serum samples were collected between 3 days and 8 months after diagnosed or suspected COVID-19.

DBS from SARS-CoV2 vaccinated individuals
DBS samples were collected from 101 individuals by unsupervised self-sampling. For 16 of these, DBS were obtained before vaccination and then weekly for 8 weeks, thereafter at 3, 4, 5 and 6 months after the first and second vaccinations. For this collection, one drop of blood of about 10-50 µl from a finger prick was spotted on Whatman™ paper (903 Protein Saver Card, Catolog number: 10531018, GE Healthcare Life Sciences, Sweden). DBS samples were dried for more than 60 min at room temperature (RT), and then stored at 4 • C. A circular piece of 1.2 mm in diameter was cut using a Uni-core™ puncher (Catalog number: WHAWB100074, GE Healthcare Life Sciences, Sweden) and used for AbPEA. DBS stored at either RT, 4 • C or 20 • C for 1 week yielded consistent results, as did repeated testing of the same samples, (data not shown). An earlier study measuring protein levels in DBS demonstrated that closely similar levels were detected in different samples cut from the same DBS [15].
Collections of all samples were approved by the Swedish Ethical Review Authority (Dnr 2020-06575 and 2021-03952). All DBS donors signed informed consent documents.

Protein-oligonucleotide conjugation
The sequences of all oligonucleotides in this study have been published previously [16], and those ones used here are described in supplementary Table S1. Thirty micrograms of S1-RBD and NP was activated with the 2 mM crosslinker dibenzylcyclooctyne-NHS ester (DBCO-NHS ester) (Catalog number: 761524, Sigma Aldrich) at RT for 30 min. Unreacted crosslinker was removed using a 7 K MWCO Zeba spin desalting column (Thermo Scientific). Both activated proteins were then divided into two aliquots and each aliquot was incubated at 4 • C overnight with either one of a pair of azide-modified, barcoded oligonucleotides (Click-FWD1 and Click-Univ-REV for the S1-RBD protein, and Click-FWD2 and Click-Univ-REV for the NP). Oligonucleotides used for covalent coupling to antigen have an azide residue at the 5 ′ ends. One of the oligonucleotides in each pair was assay-specific and 63 nt long, while the other protein aliquot was conjugated to a common 37 nt oligonucleotide called Click-Univ-REV. An assay-specific 89-mer oligonucleotide (hyb-REV), comprising 32-and 9-nt segments complementary to Click-Univ-REV and to Click-FWD, respectively, was then allowed to hybridize to the probe with the oligonucleotide Univ-REV. In addition, goat IgG was conjugated with both Click-FWD3 and Click-REV3, and used as an extension and PCR control. The extension primers included uracil residues in place of some T residues, as indicated in supplementary Table S1. The molecular beacon included a fluorophore, FAM, and a quencher, DABSYL, at either end.

Proximity extension assay
One µl serum sample or dilutions thereof were combined with 3 µl incubation solution of PBS buffer containing 5 mM EDTA, 100 µg/ml salmon sperm DNA (Sigma Aldrich), 0.1% BSA, 1 mM biotin, 100 nM goat IgG, 0.05% tween 20 and 1.33 nM of 3-plex PEA reagents in the form of pairs of oligonucleotide-conjugated S1-RBD and NP, as well as a goat IgG conjugate, serving as extension control. The combined samples and reagents were incubated at 37 • C for 1 h or 4 • C overnight with similar results. After incubation, 96 µl extension solution containing 1X Hypernova buffer (BLIRT S.A. Gdańsk, Poland), 1.5 mM MgCl 2 , 0.2 mM of each dNTP, 1 µM of universe forward and reverse primers, 2 U/ml Platinum™ Taq DNA polymerase (Invitrogen, Carlsbad, CA, USA) and 5 U/ml Hypernova DNA polymerase (Blirt Cat No. RP235) was added, thereby diluting the reactions 25-fold. The extension/pre-PCR reactions were incubated at 50 • C for 20 min, followed by 5 min heat-activation at 95 • C and then 17 cycles of pre-PCR of 95 • C for 30 s, 54 • C for 1 min, and 60 • C for 1 min using a primer pair common for all extension products. For the subsequent qPCR detection, 2.5 µl of the extension/ pre-PCR products were transferred to a 96 or 384-well plate and combined with 7.5 µl qPCR mix to a final concentration of 1X PCR buffer (Invitrogen), 2.5 mM MgCl 2 , 0.2 mM of each dNTP, 1.6 µM ROX reference Dye (Invitrogen), 0.3 µM molecular beacon (Biomers, Ulm, Germany), 0.9 µM of the appropriate amplicon-specific primer pair, 0.1 U uracil N-glycosylase and 0.5 U recombinant Taq polymerase (Invitrogen). Quantitative real-time PCR was run with an initial incubation at 25 • C for 30 min, followed by incubation at 95 • C for 5 min, and then 30 cycles of 15 s at 95 • C, and 1 min at 60 • C. To control for variation between experiments and microtiter plates, a standard serial dilution of rabbit anti-S1-RBD IgG (from 10 nM to 1 pM), extension/PCR control [15], as well as positive patient and negative patient controls, were also included in each 96 well microtiter plate.
For DBS, disks 1.2 mm in diameter were cut from the filter papers used for collecting DBS samples and added directly to reactions or first eluted in a 100 µl volume from which 1 µl was incubated with the oligonucleotide-conjugated antigens. Elution of antibodies from DBS was performed by adding DBS paper disks into 100 µl elution buffer (PBS containing 0.5% tween20, 1% BSA and 1% Halt protease inhibitor cocktail (Thermo Scientific)) with shaking on a shaker set at 900 rpm for 60 min at RT [17]. The AbPEA assay for DBS samples was identical to that for serum samples as described above. Due to short supply of Hypernova DNA polymerase, Onetaq Hotstart DNA polymerase (Catalog number: M0481S, New England Biolabs) was used, which functions for both DNA extension and replication. The solution added for extension and pre-PCR contained 100 mM Tris-HCl pH 8.5, 50 mM KCl, 3 mM MgCl 2 , 0.75% TritonX-100, 0.2 mM dNTP, 5% DMSO, 5% glycerol, 1 µM forward and reverse primer and 1 unit Onetaq Hotstart DNA polymerase. The performance of Onetaq Hotstart DNA polymerase was equivalent to that of Hypernova DNA polymerase for both serum and DBS samples in corresponding buffer.
For DBS samples collected before the pandemic, pre-vaccination and after a first vaccination, AbPEA was performed directly without elution and dilution. For DBS samples collected after a second vaccination or after infection and first vaccination, AbPEA was carried out after elution and 100-fold dilution, and the Ct values were recalculated by subtracting 6.64 from detected Ct of 100-fold dilution to compensate for the sample dilution.

Antibody detection by Abbott, Meso scale discovery and an in-house ELISA
The Abbott SARS-CoV-2 IgG assay and MSD V-PLEX SARS-CoV2 panel 11 were analyzed according to instructions by the suppliers (Abbott 6R86-22 and MSD K15437U, respectively). An in-house ELISA was performed as previously described [18].

Comparison between antibody detection using MSD and AbPEA
Eight serum samples with different antibody levels were selected and subjected to an MSD assay using V-Plex IgG indirect serology. In order to compare the assay sensitive and dynamic ranges between MSD and Ab-PEA, serum samples were diluted 50-, 500-, 5000-, 50,000-or 500,000fold. 50 µl of the diluted serum samples were used for MSD, while 1 µl samples were used for AbPEA. A serial dilution of rabbit anti-S1-RBD was also included in both assays.

Data and statistical analyses
The inter-assay coefficient of variation (CV) of the anti-S1-RBD assay was determined by analyzing a serial dilution of rabbit anti-S1-RBD IgG in duplicate on 10 different occasions. From the standard curve a high and a low concentration of rabbit antibodies (1 nM and 10 pM, respectively) were selected and used to calculate the overall mean, standard deviation (SD) and %CV using the formula: Overall %CV = (SD of plate Ct means / mean of plate means) x 100. The average %CV from the low and high concentration was used as the inter-assay %CV. The intra-plat%CV was investigated using goat IgG conjugated with both forward and reverse oligonucleotides, serving as an internal control. This reagent produces detectable extension products irrespective of the presence of SARS-CoV-2 antibodies in the sample and was included in all reactions. All graphs and figures were plotted in Rstudio.
Receiver operating characteristics (ROC) analysis was used to estimate the performance of anti-S1-RBD antibody detection by AbPEA. ROC analysis represents the rates of false negatives and positives detected by the assay across variable cutoffs and provides a means to determine the optimal cutoff point for the assay. ROC curve graphics was generated and the area under the curve (AUC) was calculated in Rstudio using the pROC package [19]. The optimal ΔCt cutoff values were determined as the ΔCt cutoff yielding the highest Youden statistics, using the Youden Index formula: J = Se + Sp-1, where Se is sensitivity and Sp is specificity, both varying between 0 and 1.

Assay sensitivity, dynamic range and optimal serum dilution
The PEA technique is commonly used for multiplex protein detection in, for instance, serum or plasma samples [16]. Here the technique was adapted for detection of antibodies against specific antigens. A schematic illustration of the AbPEA is shown in Fig. 1. The analytical sensitivity of the assay was first evaluated by investigating dilutions of a commercial purified rabbit anti-S1-RBD antibody preparation as shown in Fig. 2 A. Next, serial dilutions of convalescence serum samples were investigated, as illustrated in Fig. 2B and C. At 1 pM of the rabbit anti-S1-RBD antibodies in the 1 µl samples (0.15 pg) the signal over background was greater than 2.2 Ct (this ΔCt is calculated as Ct value of negative control minus Ct value of positive samples), corresponding to 4.6-fold over background levels. Such serial dilution of rabbit anti-S1-RBD antibodies was included in each microtiter plate in the following experiments to control for variation between plates and experiments. The variation of the standard curves across 10 plates and experiments performed at different times was analyzed with very good reproducibility as demonstrated in Fig. 2A. The inter-assay coefficients of variations (CV) were 8.13%, 6.02%, 5.03%, 3.25%, and 5.5% at 10 nM, 1 nM, 100 pM, 10 pM and 1 pM of the rabbit anti-S1 antibody detection, respectively. An intra-assay CV variation of 1.3% was calculated using the internal control of goat IgG conjugated with both forward and reverse oligonucleotides.
A set of 40 sera was used to investigate the performance of the assay, 20 of the samples (P1-P20) were positive for antibodies against NP using the Abbott SARS-CoV-2 IgG kit and 20 were negative controls (H1-H20). Five randomly selected anti-S1-RBD positive sera (P1, P3, P5, P7, P10) and 3 negative controls (H1, H3, H5) were used in initial analyses. One µl of the undiluted samples as well as up to 10,000-fold diluted positive sera and up to 100-fold diluted negative sera were analyzed for both anti-S1-RBD and anti-NC antibodies. As shown in Fig. 2B, anti-S1-RBD antibody levels in 10,000-fold diluted patient samples P1 and P5 were quite high with ΔCt values of 3.9 and 5.6, respectively, whereas the limits of detection for P3 and P10 were reached already at a 1000-fold dilution, and for P7 at a 100-fold dilution. The signals for serum samples from control subjects with no known antibodies against the virus were similar to that for the no template controls (NTC). The detection signals for anti-NP antibodies were lower than those for anti-S1-RBD antibodies (Fig. 2C). The highest detectable dilutions of anti-NP antibodies of around 500-fold was seen in samples P1 and P5, and 100-fold for samples P3 and P10, while sample P7 could only be diluted 10-fold. Thus, the assay for anti-NP antibody detection was more than 20-fold less sensitive than that for anti-S1-RBD antibodies, and the dynamic range extended over only two orders of magnitude. It is not known whether the lower sensitivity reflects differences in antibody levels or in Fig. 1. Schematic illustration of AbPEA for antibody detection. A) Preparation of the probes. Two portions of the recombinant S1 or NP antigens were each conjugated to a pair of oligonucleotides. The pairs of oligonucleotides contain a mutually complementary region, and specific barcode each. One of them also has a target sequence for a molecular beacon. B) AbPEA tests for detecting antibodies to SARS-CoV-2. Upon incubation with serum samples or cuts of DBS from previously infected or vaccinated individuals, the di-of multimeric nature of antibodies in the sample bring pairs of oligonucleotide-conjugated protein probes in proximity. This allows their attached oligonucleotides to hybridize and be extended by a DNA polymerase. The extension products incorporate barcodes from both oligonucleotides conjugated to a particular antigen. The resulting DNA extension products are then amplified using universe primers to generate sufficient amount amplicons (pre-PCR). Finally, antigen-specific amplicons are quantified by real-time PCR using specific primers and molecular beacons to record at which PCR cycle (Ct) a signal exceeds a threshold value. ΔCt values are calculated as the Ct of the negative control minus that of the sample. the performance of the AbPEA. The background in assays for NP antibodies was quite high. The levels of antibodies directed against S1-RBD and NP correlated among these samples ( Fig. 2B and C). Titration tests showed that undiluted sera with high levels of antibodies against either the S1-RBD or NP proteins resulted in signals at plateau levels, while 10fold dilutions of serum samples generated signals below plateau levels for all samples. Therefore, 10-fold serum dilutions were used in all the following analyses unless otherwise stated.
All 20 serum samples negative for NP antibodies according to the Abbott test had AbPEA values for anti-S1-RBD similar to those of the NTC, while the 20 positive samples had ΔCt values between 9.65 and 13.6, except for samples P9 and P12, for which the signals were similar   3. Pilot tests to establish AbPEA assays for anti-S1-RBD and anti-NP specific antibodies. 20 serum samples were used that had tested positive by the Abbott SARS-CoV-2 IgG test (labeled as P) and 20 tested negative (labeled as H). A) AbPEA detection of anti-S1-RBD antibodies. B) AbPEA detection of anti-NP antibodies. The y axes display Ct (empty column) and ΔCt (solid column) values for the detection of anti-S1-RBD (blue) or anti-NP antibodies (green). C) and D) ROC curve analysis of the anti-S1-RBD and NP antibody tests. The table at the bottom displays the data of the Abbott SAR-CoV-2 IgG test. The asterisks mark OD values close to the cutoff of 1 of the Abbott test.
to those for the NTC (Fig. 3A). Detection of anti-NP antibodies showed results similar to those for anti-S1-RBD antibody reactivity, although signals were generally lower (Fig. 3B). Consistent with the AbPEA results for anti-S1-RBD antibodies, anti-NP antibodies were not detected in samples P9 and P12. A ROC analysis was performed for the above anti-S1-RBD antibody tests to identify an optimal cutoff between negative and positive samples ( Fig. 3C and D). The ROC analysis shows the relationship between rates of false positive and negatives at different cutoff values for the recorded antibody levels using AbPEA. In general, the closer the AUC is to the maximum, the better the assay distinguishes among cases and controls. The inserted table in Fig. 3C shows the sensitivity and specificity of the anti-S1-RBD assay depending on the choice of ΔCt cutoff value. Samples P11 and P18 were clearly positive by AbPEA for both anti-S1-RBD and anti-NP, although they were close to the cutoff value according to the Abbott test for anti NP antibodies (See table at the bottom of Fig. 3). An optimal cutoff point was identified by finding the maximal Youden Index (J) [20]. Using a cutoff of ΔCt 2, the S1-RBD antibody assays achieved 95% sensitivity with no false positives. This sensitivity was estimated on the assumption that all samples with positive Abbott test results were indeed true positives, but as mentioned the Abbott test measures antibodies directed against NP, while AbPEA is targeting anti-S1-RBD antibodies [21].

Test of AbPEA performance using a set of patient serum samples
The AbPEA assay for anti-S1-RBD protein antibody reactivity was applied to analyze sera collected before the COVID-19 pandemic from a cohort of 100 individuals and 94 serum samples collected from individuals previously scored as positive for COVID-19, either by RT-PCR for SARS-CoV-2 RNA (68 samples) or according to characteristic symptoms. In addition, data were also available from analyses by an inhouse ELISA for all samples and by Abbott test for most of the samples. The sample donors had an age range spanning from 18 to 97-years, and they were collected from as early as 3 days after diagnosis to after more than 8 months. On the basis of the pilot test results, 1 µl of 10-fold diluted sera was used for initial AbPEA analyses. However, for samples with a ΔCt of less than 4, as well as samples collected within one month of infection, undiluted sera were also used for the analysis. All 100 samples collected in 2019 before the pandemic had Ct values of 21 or more (23.7 on average), similar to those of the NTC, demonstrating that the assay had a specificity of 100% (Fig. 4). Among 50 samples collected one month after onset of symptoms, 49 were positive with ΔCt over 6. One sample only yielded positive results when undiluted serum was used. All these samples were also positive by the Abbott and inhouse ELISA tests. Among 15 samples collected between 11 and 31 days after symptom onset, only one was negative by both AbPEA and ELISA. Four out of 5 samples collected on the 10th day were positive by AbPEA, while 3 of them were positive by ELISA. Among 23 samples collected between days 4-9, 5 were positive by AbPEA, but none according to ELISA (results are included in supplementary Table S2). Furthermore, detection signals varied over a broad range, spanning from 2.4 to 13.6 ΔCt, corresponding to amplicon levels ranging from 5.3-to 12,400-fold over background on a linear scale.
The assay performance in these cohorts was investigated by ROC analysis. A cutoff value of ΔCt 2 maximized the diagnostic sensitivity and specificity and preserved a 100% specificity using the Youden Index. The ROC analyses revealed an AUCs of 0.94 (95%CI: 0.91-1) for all 94 collected patient samples from 3 days (Fig. 4), and 0.998 (95%CI: 0.993-1) for samples collected from 11 days.

Comparison of AbPEA to other methods
As shown in Fig. 5A the results of AbPEA analyses of patient samples correlated slightly better with those of the in-house ELISA, both assays detecting S1-RBD antibodies. In contrast, the Abbott assay targets NP antibodies. The Abbott data presented in Fig. 5B were from 66 sera. All samples taken after 31 days of symptom were positive by all three methods. One of the significant differences between the three assay techniques was that ΔCt values from AbPEA tests correlated better with the time of sample collection, with higher values one month after symptom onset than within one month. Among 29 sera taken 3-10 days after symptom onset, 9 were positive by AbPEA and 5 positive by ELISA. Two samples that were positive by Abbott, but negative by both AbPEA and ELISA were collected 4 or 10 days after symptoms (data in Supplementary materials).
The assay sensitivity of AbPEA was further compared to one of the most sensitive assays, Meso Scale Discovery (MSD), by titration experiments. A selection of 8 serum samples with distinct antibody levels according to earlier AbPEA analyses were used for comparison. Sera were serially diluted from 50-to 500,000-fold for MSD and 10-to 100,000-fold for AbPEA. One µl and 50 µl of each dilution, respectively, was used for AbPEA and MSD. In the case of AbPEA, SARS-CoV-2 S1- Fig. 4. AbPEA analyses of sera from 94 COVID-19 positive patients and 100 negative control samples. A) The cutoff of ΔCt 2 was used to distinguish samples positive or negative for S1 antibodies. Color codes represent the samples collection times after symptom onset. B) ROC analysis of the AbPEA analysis for S1-RBD-specific antibodies yielded an AUC of 94.4% for all samples collected from 3 days (blue) and 99.8% for samples collected from 11 days after symptom onset (red).
RBD antibodies could be detected using as little as 0.1 nl of high-titer serum samples, such as samples 1 and 2, while a 1 µl sample volume also gave positive results for samples with low antibody levels, taken 8 days after symptom debut (results are included in supplementary  Table S3). Accordingly, samples of 1 µl of undiluted serum should be used in AbPEA to maximize diagnostic sensitivity, although a 1:10 dilution is helpful to distinguish among samples with high antibody levels. In MSD, a 50,000-fold dilution or 0.1 nl of serum was optimal for most samples according to the manufacturer. However, at this dilution, positive samples with low antibody levels from 8 days after symptom debut were negative.

Monitoring anti-S1-RBD antibody levels pre-and post-vaccination by AbPEA analysis of DBS
Anti-S1-RBD antibody levels in self-sampled DBS collected prior or at different time points after a first and a second vaccination were examined. In total, 101 healthy individuals from the medical faculty of Uppsala University participated in this survey. Sixteen of the participants collected blood samples prior to or on the day of vaccination, as well as at weeks 1, 1.5, 2, 3, 4, 5 and 6 following their first and then also a second vaccination. 1.2 mm diameter discs (corresponding to approximately 0.5 µl fresh blood) were added to an incubation solution  either directly or after elution and using 1 µl of 100-, 400-and 1000-fold dilutions. DBS collected from 21 individuals in 2015 and 2019 were used as negative controls. In addition, repeated tests of DBS samples taken at the same time and stored either at RT, 4 • C or − 20 • C for 1 week, as well as collected 3-5 weeks after a second vaccination yielded consistent ΔCt value. Furthermore, the same results were obtained when analyzing anti-S1 antibodies from ø1.2 mm DBS and 0.5 µl fresh blood from the same individual, indicating good efficiency of antibody detection from DBS (data not shown).
Anti-S1-RBD antibody levels were measured by AbPEA in a series of DBS collected from two individuals, one vaccinated individual without previous infection (Fig. 6A) and another having received one vaccination after two infections (7 and 5 months before the vaccination) (Fig. 6B). In the first case, antibody levels increased rapidly during the first two weeks after a first vaccination and reached the highest level in the 3rd week, remaining at similar levels even after a second vaccination when undiluted samples were analyzed (blue line). These results were at plateau levels for the assay, and therefore 100-fold. 400-fold and 1000fold dilutions of the samples were also analyzed. A significant increase in antibody levels was observed at one week after a second vaccination, which was sustained for over three weeks and decreased slightly after one month. For the previously infected person, antibody levels were already high before vaccination. A significant increase of antibodies after a first vaccination was observed after sample dilution, but only a limited further increase was detected after the second vaccination. The results suggested that a 100-fold dilution was sufficient for accurate antibody detection after a second vaccination. However, at a 100-fold dilution, antibody levels were not detectable 1.5 week (11 days) after a first vaccination. To investigate antibody levels in DBS over a wide concentration range, both punched-out samples without elution, as well as eluted and 100-fold diluted samples were used in the following experiments.
Antibody levels were measured in consecutive samples from 14 vaccinated individuals using punched-out disc directly without elution and 100-fold diluted DBS samples after elution. Antibody levels varied among the individuals, but the changes across two vaccinations were similar, increasing to a maximum at two weeks after a first vaccination, and sustained before increasing further after a second vaccination (Fig. 7). Antibody levels increased already to a maximum at one week after the second vaccination, and remained at similar levels for about three months before slowly decreasing.
In addition, AbPEA analyses of S1-RBD antibody responses were performed on DBS collected from 21 individuals in 2015 or 2019 before the pandemic and 141 DBS collected from 101 individuals during the pandemic, including 27 DBS collected before vaccination, 37 DBS collected after a first vaccination and 78 DBS obtained more than one month after a second vaccination (Fig. 8). The average AbPEA results increased by more than 8.76 ΔCt values or 430-fold by 3 weeks after a first vaccination, compared to levels before vaccination. The second vaccination resulted on average in a further boost of 4.2 ΔCt values, or another 20-fold. Antibody levels of individuals who had been infected before or after vaccination (red circles in Fig. 8) were higher than those of individuals who had not had the disease.

Discussion
A homogenous PCR-based assay, AbPEA, is presented, serving to sensitively and conveniently detect antibody responses to the SARS-CoV-2 proteins RBD and NP after infection or vaccination. The assay is based on the proximity ligation assay technology for protein detection [4], which has been further developed and commercialized as proximity extension assays (Olink Proteomics, www.olink.com). In AbPEA the dior multimeric nature of antibodies serves to bring into proximity pairs of antigen proteins, conjugated with oligonucleotides. After incubation, the reactions are diluted with the addition of a DNA polymerase. This allows oligonucleotides on pairs of conjugates that remain in proximity by virtue of having been bound by the same antibodies to be enzymatically extended, generating amplifiable reporter DNA strands with no requirement for washes.
A closely related DNA-assisted antibody detection technique, agglutination-PCR (ADAP), has previously been applied for detection of anti-SARS-CoV-2 antibodies both in serum/plasma and DBS samples from patients [6,7]. The assay is similar to AbPEA in that pairs of oligonucleotides are conjugated to S1-RBD protein molecules, which are brought into proximity upon binding to antibody, although a proximity ligation rather than extension mechanism is used to create the amplicon that is then measured via PCR. Their assay consists of four steps and uses 4 µl sample, while AbPEA generates results in three reaction steps using 1 µl sample. For analysis of DBS by ADAP, antibodies are first eluted in 1000 µl buffer and then concentrated.
In this report, the assay sensitivity and dynamic range of AbPEA were investigated, and the assay was benchmarked to other methods for measuring specific antibody responses. Furthermore, the changes in antibody levels were monitored during the course of two consecutive vaccinations using serially self-collected DBS. The performance of AbPEA was evaluated using both liquid serum samples and DBS from Anti-S1-RBD antibodies could be also detected from DBS collected from 11 days after vaccination. Assay sensitivity of AbPEA was similar to that of commercial assays from Abbott and MSD using the same serially diluted samples. AbPEA could detect anti-S1-RBD antibody levels in high-titer serum samples even after 10,000-fold dilution corresponding to 0.1 nl aliquots of serum. To achieve sensitive detection early in immune responses undiluted samples were used, while distinction among further elevated antibody levels after boosting required that samples were diluted.
By using both undiluted and 10-fold diluted serum samples, antibody concentrations could be monitored over more than 3 logs. Since AbPEA measures how antibodies can bring oligonucleotide-conjugated antigen molecules into proximity, the assay can detect all classes of immunoglobulins (Ig). Since no anti-Ig reagents are required, the same assay can be applied to assess antibody responses to a given protein in any species, for example for veterinary purposes (data not shown).
Most data presented here concern detection of anti-S1-RBD antibodies. However, one of the defining features of the PEA assay mechanism is that it lends itself to multiplexing. The assays developed by Olink Proteomics measure 96 or even more proteins and controls in individual samples. The reason for this far greater potential for multiplexing PEA reactions compared to sandwich ELISA reactions is that the assays are designed so that only the correct pairs of reagents can give rise to detectable amplicons. The AbPEA therefore has excellent potential to allow measurements of antibody responses to large numbers of antigens in the same assay. Such measurement of antibody repertoires can prove of diagnostic value in many contexts. Besides immune responses to infectious agents or vaccines, as demonstrated here, it is also of interest to survey responses to autoantigens, and antibodies to interferon have for example been linked to poor outcome in COVID-19 [22]. In preliminary experiments a four-plex panel was used to simultaneously assess antibody levels directed against S1-RBD, NP, IFNα1 and IFNω1 (data not shown). Consistent results were obtained when measuring antibody responses to S1-RBD and NP in individuals who had been infected with the virus. Furthermore, significant levels of anti-NP antibodies could be detected from DBS collected from individuals that had been infected by SARS-CoV2 after 2 doses of vaccinations against the S1-RBD protein (data not shown). The anti-NP assay showed a higher background and a relatively low signal compared to that for S1-RBD, however, and further optimization is warranted. Applying AbPEA to self-sampled DBS, antibody responses were detected by 1.5 weeks after a first vaccination. In order to demonstrate the further increased antibody levels after a second vaccination, DBS portions were eluted and diluted 100-fold.
In summary, AbPEA allows convenient detection of antibody levels across wide dynamic ranges in a convenient format with no requirement for washes and with a potential for extensive multiplexing.  , and among vaccinated individuals sampled were obtained pre-vaccination (PreV), 4-6 weeks after a 1st vaccination (1stV before 2ndV) and 1-8 month after 2nd vaccination (After 2nd) as reflected in bee swarm plots. Red circles indicate individuals who also had been infected by the virus. For DBS samples collected before the pandemic, pre-vaccination and after 1st vaccination, AbPEA were performed directly on paper punches without elution and dilution. For DBS samples collected after a 2nd vaccination or after infection and 1st vaccination, AbPEA was performed after elution and 100-fold dilution, but the Ct values were calculated by subtracting 6.64 from detected Ct of 100-fold dilution. P-values were determined by a Wilcoxon test. B) ROC values were calculated using data after the first vaccination. manuscript.

Declaration of Competing Interest
UL is a founder and shareholder of Olink Proteomics, having rights to the PEA technology. All other authors declare that they have no competing financial interests or personal relationships that could influence the work reported in this paper.

Data availability
Data will be made available on request.


Paleogenomes Reveal a Complex Evolutionary History of Late Pleistocene Bison in Northeastern China

Steppe bison are a typical representative of the Mid-Late Pleistocene steppes of the northern hemisphere. Despite the abundance of fossil remains, many questions related to their genetic diversity, population structure and dispersal route are still elusive. Here, we present both near-complete and partial mitochondrial genomes, as well as a partial nuclear genome from fossil bison samples excavated from Late Pleistocene strata in northeastern China. Maximum-likelihood and Bayesian trees both suggest the bison clade are divided into three maternal haplogroups (A, B and C), and Chinese individuals fall in two of them. Bayesian analysis shows that the split between haplogroup C and the ancestor of haplogroups A and B dates at 326 ky BP (95% HPD: 397-264 ky BP). In addition, our nuclear phylogenomic tree also supports a basal position for the individual carrying haplogroup C. Admixture analyses suggest that CADG467 (haplogroup C) has a similar genetic structure to steppe bison from Siberia (haplogroup B). Our new findings indicate that the genetic diversity of Pleistocene bison was probably even higher than previously thought and that northeastern Chinese populations of several mammalian species, including Pleistocene bison, were genetically distinct.

Introduction
During the Pleistocene, the genus Bison was widely distributed in Holarctic regions and evolved into different species, adapted to various environmental conditions [1][2][3][4]. Nowadays, there are only two extant species, the European bison (Bison bonasus) and the American bison (Bison bison), confined to relatively limited regions [5][6][7]. Among the recently extinct species, the steppe bison (Bison priscus) was the most wide-ranging one, inhabiting the Mid-Late Pleistocene steppes of the northern hemisphere. Its habitat ranged from Europe across Asia to North America, which is also called the great Pleistocene bison belt [8][9][10][11][12]. However, despite its abundant remains, this extinct species is still relatively poorly known regarding its origin, dispersal and population divergence [1,12,13].
In recent years, many genetic studies have been performed on steppe bison remains from Europe, Siberia and North America [3,10,[14][15][16][17]. Earlier maternal phylogenetic analyses indicated that the analyzed specimens represented a relatively homogeneous population, that is, all samples analyzed clustered into the main bison mitochondrial clade [10,18]. This situation changed with the discovery of a frozen early Holocene bison carcass (9497 ± 92 14C yr BP) from the mouth of the Rauchua River, Chukotka, Russia (69 • 30 N, 166 • 49 E) [13,19]. Analysis of its mitochondrial DNA revealed that this individual fell outside the known genetic diversity of steppe bison, forming a deeply divergent sister lineage to all other bison mitochondrial sequences [19], hereafter referred to as bison haplogroup C. Surprisingly, further investigation of mitogenomes from 26 additional bison samples collected from northeastern Siberia, including sites at Rauchua River, Middle Indigirka River and the Bilibinsky region, showed that all the newly analyzed samples fell into the main Eurasian steppe bison clade rather than haplogroup C [17]. Hence, the Pleistocene to Holocene distribution of bison haplogroup C has remained unknown and its geographical origin is still unclear. It has been suggested that bison carrying haplogroup C likely dispersed northward from their original habitat and escaped from the Late Pleistocene extinction event in the northern Siberian refuge [17], maybe similar to the Holocene survival of the giant deer (Megaloceros giganteus B.) [20,21]. However, so far it remains an open question where this lineage originated from and how and when its distribution changed. Situated to the south of Rauchua, northeastern China is an obvious candidate region for the geographical origin of this enigmatic bison lineage.
Northeastern China is one of the classic regions for the Late Pleistocene "Mammuthus-Coelodonta Fauna", in which members of the family Bovidae were common [22][23][24][25]. According to the fossil record, at least three genera of Bovidae (i.e., Bubalus, Bos and Bison) coexisted in this region during that time, including several named species, e.g., Bubalus wansjocki, Bubalus teilbardi, Bubalus brevicornnis, Bos primigenius, B. priscus, Bison exiguus and several others [22,23,[26][27][28]. Especially the members of the genus Bison dominated the local fauna during the Late Pleistocene. Fossil bison in northeastern China were initially identified as steppe bison (B. priscus) based on morphological analyses [27,29], while other scholars argued that it should be assigned to a separate species, B. exiguus, which was further divided into three subspecies, i.e., Bison exiguus exiguus, Bison exiguus curvicornis and Bison exiguus harbinensis [22,23,28]. So far, the taxonomic assignment of northeastern Chinese bison has remained controversial. Ancient DNA provides a powerful tool for revealing the evolution of extinct populations [30][31][32]. Unfortunately, to our knowledge, only two partial D-loop sequences were retrieved from fossil bison excavated from China so far [10]. Thus, almost no genetic information on Chinese fossil bison is available.
In this study, we generated near-complete and partial mitochondrial genomes from Late Pleistocene bison collected from northeastern China. A low coverage nuclear genome was also retrieved from our best-preserved bison fossil. Using the obtained sequences as well as previously published data, we performed phylogenetic analyses, principal component analyses (PCA) and admixture analyses to explore their relationships with their counterparts from Eurasia and America. We also estimated the divergence time and genetic distance of different bison lineages. Our results suggest a complex history of extinct bison lineages in northeastern China and contribute to filling the gap in understanding the evolution of Bison during the Late Pleistocene and early Holocene epochs in Asia.

Sampling Information
In total, 21 fossil bison specimens from northeastern China were included in this study; 19 samples were excavated from Heilongjiang province (Haerbin, Zhaodong, Zhaoyuan, Qinggang and Daqing sites), and two samples from Dalian, Liaoning province ( Figure 1, Table S1). Accelerator mass spectrometry (AMS)-radiocarbon dating was performed on two of these specimens (CADG456 and CADG467) at the β Analytic inc. in the USA (Lab Nos: β-609039 and β-583723). Detailed information on the samples is listed in

Sampling Information
In total, 21 fossil bison specimens from northeastern China were included in this study; 19 samples were excavated from Heilongjiang province (Haerbin, Zhaodong, Zhaoyuan, Qinggang and Daqing sites), and two samples from Dalian, Liaoning province ( Figure 1, Table S1). Accelerator mass spectrometry (AMS)-radiocarbon dating was performed on two of these specimens (CADG456 and CADG467) at the β Analytic inc. in the USA (Lab Nos: β-609039 and β-583723). Detailed information on the samples is listed in Table S1. The base map of the global world (chart no. GS (2021) 5443) was downloaded from the National Administration of Surveying, Mapping and Geoinformation of China (http://bzdt.ch.mnr.gov.cn, accessed on 22 March 2022).

DNA Extraction, Double-Stranded Library Construction and Sequencing
We performed DNA extraction and library preparation on 19 bison samples (Table  S1) in a dedicated ancient DNA laboratory at the China University of Geosciences (Wuhan). DNA extraction was performed following the protocol described in Hu et al. [33]. Each specimen was ground into about 200 mg fine powder with a mortar and pestle. The bone powder was digested in 4.5 mL of EDTA (0.5 M, pH = 8) and 0.06 mL of Proteinase K (20 mg/mL) by incubation in a hybridization oven at 37 °C under constant agitation for 16 h. After centrifugation at 7000 rpm for 10 min, we recovered the supernatant with an ultrafiltration tube (Millipore, Germany) and concentrated it to 100 µL at 7000 rpm for 35 min. Finally, we purified and eluted DNA into 80 µL EB buffer using the QIAquick PCR Purification Kit (Qiagen, Germany) following the manufacturer's instructions. To monitor potential contamination, we set up blank controls in both experimental steps.
Multiple double-stranded libraries for Illumina high-throughput sequencing were constructed from 20 µL DNA extract according to the protocol described by Meyer and Kircher [34]. In the blunt-end repair step, DNA templates were mixed with BSA (New England Biolabs, Ipswich, MA, USA) and NEB buffer 2 (New England Biolabs, Ipswich,

DNA Extraction, Double-Stranded Library Construction and Sequencing
We performed DNA extraction and library preparation on 19 bison samples (Table S1) in a dedicated ancient DNA laboratory at the China University of Geosciences (Wuhan). DNA extraction was performed following the protocol described in Hu et al. [33]. Each specimen was ground into about 200 mg fine powder with a mortar and pestle. The bone powder was digested in 4.5 mL of EDTA (0.5 M, pH = 8) and 0.06 mL of Proteinase K (20 mg/mL) by incubation in a hybridization oven at 37 • C under constant agitation for 16 h. After centrifugation at 7000 rpm for 10 min, we recovered the supernatant with an ultrafiltration tube (Millipore, Germany) and concentrated it to 100 µL at 7000 rpm for 35 min. Finally, we purified and eluted DNA into 80 µL EB buffer using the QIAquick PCR Purification Kit (Qiagen, Germany) following the manufacturer's instructions. To monitor potential contamination, we set up blank controls in both experimental steps.
Multiple double-stranded libraries for Illumina high-throughput sequencing were constructed from 20 µL DNA extract according to the protocol described by Meyer and Kircher [34]. In the blunt-end repair step, DNA templates were mixed with BSA (New England Biolabs, Ipswich, MA, USA) and NEB buffer 2 (New England Biolabs, Ipswich, MA, USA). After that, we ligated a 1:20 adapter dilution (Sangon Biotech, Shanghai, China) and the ends of DNA fragments using Quick Ligase Buffer (New England Biolabs, Ipswich, MA, USA) and Quick Ligase (New England Biolabs, Ipswich, MA, USA). The ligation reaction was then added to an adapter fill-in reaction. Indexing PCR amplifications were carried out with Q5 Hot Start High-Fidelity 2 × Master Mix (New England Biolabs, Ipswich, MA, USA) under the following conditions: 98 • C for 30 s and 17 cycles of 98 • C for 10 s, 60 • C for 75 s, and 60 • C for 6 min. Quantitative analyses of the libraries were conducted with Qubit 4.0 (Invitrogen, Carlsbad, CA, USA) and TapeStation 4150 (Agilent, Santa Clara, CA, USA). After pooling 10 or 40 libraries in equimolar ratios, calculated from concentration and DNA length, the libraries were finally sequenced on an Illumina HiSeq×10 platform at Annoroad Inc., Beijing, China.

Single-Stranded Library Construction and Hybridization Capture
In addition, we retrieved mitogenomes from two bison samples (CADG456 and CADG471) (Table S1) in a dedicated DNA laboratory at the University of Potsdam. DNA extraction was performed as previously published by Sheng et al. [35]. Then, we prepared single-stranded DNA libraries using 20 µL DNA extract of each sample following the protocol described by Gansauge and Meyer [36] with slight modifications as in Yuan et al. [37]. Hybridization capture of the complete mitochondrial genome was carried out using the procedures described in Gonzalez-Fortes and Paijmans [38]. We prepared baits using a modern cattle sample with the following steps: Firstly, four overlapping long-range PCR (LR-PCR) primer pairs (Table S2) were used to amplify the mitochondrial genome. Secondly, the amplified modern cattle mitochondrial DNA fragments were sheared, bluntend repaired and ligated to biotinylated adapters. Thirdly, two rounds of hybridization capture were carried out to improve the enrichment rate by the protocol of Yuan et al. [37]. Finally, the enriched libraries were pooled into a single pool in equimolar ratios and sequenced on the Illumina NextSeq 500 sequencing platform in 75 bp single-end runs, following Paijmans et al. [39].

Sequence Assembly
Cutadapt v1.4.2 [40] was used to trim adapters. After discarding reads shorter than 30 bp, the overlapping read pairs were merged with Flash v1.2.11 [41]. The trimmed reads were mapped to four bison complete mitochondrial sequences (GenBank Nos. KR350472, KM593920, KX269145 and GU946980) and the nuclear genome assembly of a modern American bison (GenBank No. GCA_018282365.1) processed with the "aln" and "samse" algorithms in Burrows-Wheeler Aligner v0.6.2 [42] with default parameters. Sequences with a map quality score less than 30 were removed by using "view" and the alignment was sorted on the reference genome by a 5 mapping position using "sort" in SAMtools v0.1.19 [43]. After removing all potential PCR duplicates by using "rmdup", the resulting bam files for each library were merged into a single bam file using "merge" from SAMtools. Then the merged bam file was imported to Qualimap v2.2.1 [44] and mapDamage v2.0 [45] to calculate read coverage and test the authenticity of ancient DNA, respectively (Figures S1-S5, Tables S1, S3 and S4). Finally, the mitochondrial consensus sequence was called with a minimum coverage of 2× and a base agreement greater than 75% and exported to Geneious v.10 (http://www.geneious.com/, accessed on 7 January 2022).
Identical methods were used to assemble mitochondrial genomes or nuclear genome data from short read data of two ancient steppe bison (875 and 3133), three aurochs, six American bison, ten European bison and two yaks downloaded from ENA (the European Nucleotide Archive) (Table S4).

Mitochondrial Phylogenetic Analyses
In order to investigate the phylogenetic status of the Chinese bison mitochondrial haplotypes, we performed a maximum-likelihood (ML) phylogenetic analysis with RAxML-HPC v8 [46] on the CIPRES server [47]. Three newly obtained mitochondrial genomes (CADG456, CADG465 and CADG467) were used to reconstruct the phylogenetic tree, which was aligned with 59 steppe bison, ten American bison and four yak sequences downloaded from GenBank and ENA (Table S5) using MAFFT v7.471 [48]. After removing the 891 bp long D-loop region that cannot be aligned properly, a 15,431 bp alignment was generated based on those positions where bases could be called with confidence. Considering the heterogeneity of evolutionary rates among partitions, one GTR, one GTR + G and two GTR + I nucleotide substitution models suitable for four partitions were generated by PartitionFinder v2.1.1 [49] with the greedy search algorithm, linked branch lengths, and the Bayesian Information Criterion. Using 500 bootstrap replicates, node support was also calculated (Figure 2a).
RAxML-HPC v8 [46] on the CIPRES server [47]. Three newly obtained mitochondrial genomes (CADG456, CADG465 and CADG467) were used to reconstruct the phylogenetic tree, which was aligned with 59 steppe bison, ten American bison and four yak sequences downloaded from GenBank and ENA (Table S5) using MAFFT v7.471 [48]. After removing the 891 bp long D-loop region that cannot be aligned properly, a 15,431 bp alignment was generated based on those positions where bases could be called with confidence. Considering the heterogeneity of evolutionary rates among partitions, one GTR, one GTR + G and two GTR + I nucleotide substitution models suitable for four partitions were generated by PartitionFinder v2.1.1 [49] with the greedy search algorithm, linked branch lengths, and the Bayesian Information Criterion. Using 500 bootstrap replicates, node support was also calculated (Figure 2a). To infer the divergence times of the different bison haplogroups, we constructed a time-tree with BEAST v1.8.3 [50] using complete mitochondrial genomes. We discarded sequences lacking accurate dating information and performed a Bayesian analysis including only one of our newly obtained mitochondrial genomes (CADG467) plus 32 steppe To infer the divergence times of the different bison haplogroups, we constructed a time-tree with BEAST v1.8.3 [50] using complete mitochondrial genomes. We discarded sequences lacking accurate dating information and performed a Bayesian analysis including only one of our newly obtained mitochondrial genomes (CADG467) plus 32 steppe bison and ten American bison (Table S5). We used the time of the most recent common ancestor (TMRCA) of bison from America (165 ky BP, 95% HPD: 195-135 ky BP) [14] to calibrate the phylogeny. Tip dates of the samples were also used as calibration points using either the median calibrated radiocarbon ages or stratigraphic ages. The best fitting nucleotide substitution model (HKY + I) was inferred by jModelTest v2.1.6 [51]. We ran the Markov Chain-Monte Carlo (MCMC) chain for 60 million iterations, sampling priors and trees every 1000 iterations. Then, the log file was visually inspected using Tracer v1.7 [52] to ensure all parameters had reached effective sample sizes (ESS) above 200. A maximum clade credibility tree was calculated on the tree file using TreeAnnotator v1.8.0 [50] after discarding the first 10% of sampled states as burn-in. The final tree was visualized and annotated in FigTree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree, accessed on 8 January 2022).
Following the phylogenetic analyses by Shapiro et al. [10], we also used a dataset that only contains 639 bp of the D-loop region and constructed another phylogenetic tree by Bayesian inference in MRBAYES v3.2.7 [53]. The alignment comprises our two newly obtained sequences (CADG456 and CADG467), 95 partial D-loop sequences described by Shapiro et al. [10] (Table S6) and 67 complete mitochondrial genomes downloaded from GenBank (Table S5). The obtained D-loop sequence of the sample CADG 465 is very short, hence it was excluded in this analysis. All other parameters were the same as in the previous Bayesian analysis.

Nuclear Analyses
A total number of 24 individuals were used in principal component analysis (PCA) and admixture proportion estimation (NGSadmix), which contains our best-preserved sample (CADG467) and 23 bovid genomes downloaded from ENA (Table S4). We only considered the autosomal chromosomes to avoid different rates of evolution and inheritance. Genotype likelihood estimations (GL) were obtained using angsd v0.921 [54] with the following flags: minimum base quality (-minQ 30), minimum mapping quality (-minMapQ 30), infer major and minor alleles from GL (-doMajorMinor 1), remove all reads showing multiple hits (-uniqueOnly 1), only consider SNPs with p-value below (-SNP_pval 1e-6), remove transitions (-rmtrans 1), calculate site allele frequencies (-doMaf 1), minimum minor allele frequency (-minmaf 0.0476), skip triallelic sites (-skiptriallelic 1), only consider sites where all steppe bison have coverage (-minind 21), only include scaffolds above 1 MB (-rf). The downstream PCA and NGSadmix analyses are based on the GL estimations. PCAngsd v0.98 [55] was used to construct a covariance matrix on the obtained GL beagle file. Individual admixture proportions were calculated using NGSadmix v32 [56], specifying K = 2-6. To evaluate convergence in our results for each K, we ran each calculation 100 times with an independent seed value. If the best likelihood could be repeated twice, then the result was considered to have converged and thus be reliable. In addition, we later performed another GL estimation for PCA with only nine individuals, i.e., our one sample CADG467, two steppe bison and six American bison.
In order to investigate the nuclear phylogenetic position of our sample, we constructed a rooted nuclear neighbor-joining (NJ) tree on the same data set that was used for the first GL estimation which is based on a total of 24 individuals (Table S4). The distance matrix was calculated by using the same flags as above in the GL estimation except for adding the filter -doIBS 1 -makeMatrix. Then the NJ tree was calculated using the R package ape [57].

Result
In total, 21 bison individuals excavated from northeastern China were included in this study, DNA was obtained from six individuals (CADG456, CADG465, CADG467, CADG471, CADG605 and CADG627) ( Table S1). Five of these six samples, excluding CADG465, show a higher coverage when mapping against the mitogenome of the Rauchua bison (GenBank No. KR350472) than those of other steppe bison (GenBank No. KM593920 and KX269145) and American bison (GenBank No. GU946980) (Table S2). Considering coverage issues, only three samples (CADG456, CADG465 and CADG467), which show a mean coverage of 9.1-fold, 2.2-fold and 8.3-fold, respectively, when mapping against the consensus mitogenome of the Rauchua bison (GenBank No. KR350472), were included in downstream analyses. The mitochondrial DNA fragments are short (average length: 51 bp for CADG456, 47 bp for CADG465 and 63 bp for CADG467), and CADG456 and CADG465 show elevated C to T substitution rates at the 5 ends and elevated G to A substitution rates at the 3 ends, in agreement with the damage characteristics expected for double stranded libraries constructed from ancient DNA (Figures S1 and S2).
For the analysis of the nuclear sequences, we used one sample (CADG467) with sufficient coverage and 23 previously published bovid genomes, which were all mapped to the American bison genome (GenBank No. GCA_018282365.1). This analysis resulted in 3,274,372 mapped reads with 0.0759× average genome-wide coverage for CADG467, and 2,365,117 to 749,141,984 mapped reads with 0.2071-35.3602× average genome-wide coverage across the downloaded individuals from ENA (Table S4). After filtering, we retained 6,000,104 nuclear SNPs for further analyses.
Three newly obtained mitochondrial genomes (i.e., CADG456, CADG465 and CADG467), together with 69 published sequences of the family Bovidae were used to infer a mitochondrial ML phylogenetic tree (Figure 2a). We obtained a topology consistent with previous studies [19,58], in which extant American bison fall within the variability of the extinct steppe bison clade. Similar to the result by Vershinina et al. [17], the bison clade is further divided into three haplogroups, i.e., A, B and C (Figure 2a). Haplogroup A mainly contains extant American bison and extinct steppe bison individuals from North America, whereas haplogroups B and C are extinct lineages without surviving descendants. Haplogroup B was the main lineage of steppe bison, once widely distributed in Eurasia during the Late Pleistocene. Notably, in the ML phylogenetic tree, our three samples do not form a monophyletic clade and fall within two recognized maternal haplogroups. Two of them (CADG456 and CADG467) cluster with high bootstrap support value (100%) with the Holocene Rauchua bison and form the distinct haplogroup C, which is the first lineage that splits from the bison clade. The other sample (CADG465) falls within the mitochondrial diversity of steppe bison haplogroup B (Figure 2a).
To estimate the coalescence time of nodes in the bison mitochondrial phylogeny, we carried out a Bayesian analysis in BEAST based on near-complete mitogenomes, using both root and tip-dating for temporal calibrations. We obtained a topology similar to the ML tree, in which haplogroups A and B diverged at 195 ky BP (95% HPD: 225-171 ky BP), while the split between haplogroup C and the lineage leading to haplogroups A and B was dated at 326 ky BP (95% HDP: 397-264 ky BP). The TMRCA of haplogroups A, B and C are 167 ky BP (95% HDP: 186-153 ky BP), 151 ky BP (95% HDP: 175-131 ky BP) and 130 ky BP (95% HDP: 172-94 ky BP), respectively (Figure 2b).
To expand the number of available bison individuals, we conducted a phylogenetic analysis using a 639 bp D-loop sequence dataset based on Shapiro et al. [10]. The Bayesian phylogenetic tree divides the bison samples into clades 1, 2, 3, 4 and 5 ( Figure 3). The number of clades is different from the three main branch structures suggested by complete mitochondrial genomes, which indicate only three haplogroups in Figure 2. By comparing the two trees (Figures 2a and 3), we found that clades 1, 2 and 4 correspond to haplogroup A. Clade 3 corresponds to haplogroup B, and clade 5 corresponds to haplogroup C. Although too few complete mitochondrial genomes may lead to incomplete details of branches, our two samples (CADG456 and CADG467) together with the Rauchua bison (KR350472) keep clustering together and fall in the newly named clade 5, for which there were no representative samples in Shapiro et al. [10], forming a basal branch in all phylogenetic analyses (Figures 2 and 3). Genes 2022, 13, 1684 8 of 16 To investigate the nuclear phylogenetic position of bison carrying mitochondrial haplogroup C, we used the 6,000,104 Single Nucleotide Polymorphisms (SNPs) in the dataset, which contains our single sample (CADG467), as well as one Siberian and one Yukon steppe bison, two yaks, three aurochs, six American bison and ten European bison. The nuclear phylogenetic analysis also supports a basal position of the haplogroup C individual within the American/steppe bison clade (Figure 4). To investigate the nuclear phylogenetic position of bison carrying mitochondrial haplogroup C, we used the 6,000,104 Single Nucleotide Polymorphisms (SNPs) in the dataset, which contains our single sample (CADG467), as well as one Siberian and one Yukon steppe bison, two yaks, three aurochs, six American bison and ten European bison. The nuclear phylogenetic analysis also supports a basal position of the haplogroup C individual within the American/steppe bison clade (Figure 4). Genes 2022, 13, 1684 9 of 16 To further explore the genetic affinities of haplogroup C within the extinct steppe bison and extant American bison, we performed a PCA analysis using our available lowcoverage nuclear sequence (CADG467) together with nuclear genomes from two steppe bison and six American bison individuals. The result is shown in Figure 5 with 22.5% and 13.2% of the variance explained by PC1 and PC2, respectively. It indicates that the three analyzed steppe bison samples clustered together. However, compared to one steppe bison individual from Yukon (haplogroup A), our sample CADG467 (haplogroup C) shows a closer relationship to a Siberian steppe bison carrying a haplogroup B mitochondrial sequence, indicating a certain phylogeographic structure (Figures 5 and S6). To further explore the genetic affinities of haplogroup C within the extinct steppe bison and extant American bison, we performed a PCA analysis using our available low-coverage nuclear sequence (CADG467) together with nuclear genomes from two steppe bison and six American bison individuals. The result is shown in Figure 5 with 22.5% and 13.2% of the variance explained by PC1 and PC2, respectively. It indicates that the three analyzed steppe bison samples clustered together. However, compared to one steppe bison individual from Yukon (haplogroup A), our sample CADG467 (haplogroup C) shows a closer relationship to a Siberian steppe bison carrying a haplogroup B mitochondrial sequence, indicating a certain phylogeographic structure ( Figure 5 and Figure S6). Additionally, we also performed admixture analyses to investigate the genetic structure of steppe bison, American bison, aurochs, European bison and yak ( Figure 5). When K = 3 and 4, five groups were identified, i.e., steppe bison from Siberia, steppe bison and American bison from North America, aurochs, European bison and yak ( Figure 6). Additionally, When K = 2-5, a similar genetic structure is found among the steppe bison and American bison individuals, although the proportion of ancestry varies between the samples. CADG467 (bison haplogroup C) shares a higher proportion of ancestry with the bison individual from Siberia (bison haplogroup B) than with the steppe bison from Yukon (bison haplogroup A), while the latter has the same population background as American bison. Additionally, we also performed admixture analyses to investigate the genetic structure of steppe bison, American bison, aurochs, European bison and yak ( Figure 5). When K = 3 and 4, five groups were identified, i.e., steppe bison from Siberia, steppe bison and American bison from North America, aurochs, European bison and yak ( Figure 6). Additionally, When K = 2-5, a similar genetic structure is found among the steppe bison and American bison individuals, although the proportion of ancestry varies between the samples. CADG467 (bison haplogroup C) shares a higher proportion of ancestry with the bison individual from Siberia (bison haplogroup B) than with the steppe bison from Yukon (bison haplogroup A), while the latter has the same population background as American bison.

Discussion
ML and Bayesian phylogenetic tree reconstructions using complete mitochondrial genomes indicate that Late Pleistocene fossil bison from Europe, Asia and America group into three major clades (i.e., haplogroups A, B and C), and samples from northeastern China fall into at least two bison lineages. One of our samples (CADG465) clusters into the main Eurasian bison mitochondrial clade (haplogroup B) (Figure 2a), which is consistent with the findings of Shapiro et al. [10]. Haplogroup C, a deeply divergent maternal bison lineage, previously represented by only a single sample from Rauchua, Far Eastern Siberia [14,17,19], is confirmed by two samples (CADG456, CADG467) in this study (Figures 2a and 3). Considering the overlapping radiocarbon dates of the available samples from northeastern China (Figures 2-4, Table S1), the results indicate that the identified two bison evolutionary haplogroups (i.e., haplogroups B and C) coexisted in northeastern China during the Late Pleistocene. Our findings thus indicate a more complex evolutionary scenario of fossil bison in northeastern China than previously thought.
The confirmation of the basal mitochondrial lineage C in this study (Figures 2-4) is strikingly similar to the aurochs from northeastern China, which were assigned as an unidentified haplogroup, with a basal position in the auroch's tree [59]. In addition, it is also true for extinct cave hyenas excavated from the same region, which form a deeply diverging mitochondrial haplogroup of Crocuta [33]. It should be noted though, that for nuclear DNA, Chinese cave hyenas form a clade with other Eurasian cave hyenas to the exclusion of African spotted hyenas, whereas in the bison, the clade C individual, for which we

Discussion
ML and Bayesian phylogenetic tree reconstructions using complete mitochondrial genomes indicate that Late Pleistocene fossil bison from Europe, Asia and America group into three major clades (i.e., haplogroups A, B and C), and samples from northeastern China fall into at least two bison lineages. One of our samples (CADG465) clusters into the main Eurasian bison mitochondrial clade (haplogroup B) (Figure 2a), which is consistent with the findings of Shapiro et al. [10]. Haplogroup C, a deeply divergent maternal bison lineage, previously represented by only a single sample from Rauchua, Far Eastern Siberia [14,17,19], is confirmed by two samples (CADG456, CADG467) in this study (Figures 2a and 3). Considering the overlapping radiocarbon dates of the available samples from northeastern China (Figures 2-4, Table S1), the results indicate that the identified two bison evolutionary haplogroups (i.e., haplogroups B and C) coexisted in northeastern China during the Late Pleistocene. Our findings thus indicate a more complex evolutionary scenario of fossil bison in northeastern China than previously thought.
The confirmation of the basal mitochondrial lineage C in this study (Figures 2-4) is strikingly similar to the aurochs from northeastern China, which were assigned as an unidentified haplogroup, with a basal position in the auroch's tree [59]. In addition, it is also true for extinct cave hyenas excavated from the same region, which form a deeply diverging mitochondrial haplogroup of Crocuta [33]. It should be noted though, that for nuclear DNA, Chinese cave hyenas form a clade with other Eurasian cave hyenas to the exclusion of African spotted hyenas, whereas in the bison, the clade C individual, for which we could obtain sufficient nuclear data also forms a basal lineage in the nuclear phylogenetic reconstruction (Figure 4). These results suggest that at least some members of Pleistocene fauna in northeastern China are genetically different from their counterparts from other regions and may indicate a long-term ecological refugium for these species in this region. In future studies, mammalian material from this region should receive particular attention to better understand the evolutionary histories of Pleistocene fauna.
Until now, fossil bison carrying haplogroup C were only identified from two isolated sites. Since it was first identified by an Early Holocene sample from Rauchua [19], Vershinina et al. [17] hypothesized that populations carrying haplogroup C possibly persisted south of Rauchua during the Pleistocene. In this study, we detected a Late Pleistocene presence of this haplogroup in northeastern China (Figures 2 and 3), which is indeed south of the location where it was first discovered. In addition, investigation of numerous bison samples from Eurasia and the American continent did not detect additional haplogroup C individuals [10,14,15,18,60,61]. Therefore, haplogroup C might represent a relatively isolated lineage compared to haplogroup B, but its past geographical distribution range still remains unclear due to the very limited number of fossils that were so far found to carry haplogroup C.
The origin of bison haplogroup C has attracted attention since it was identified by ancient DNA analyses [19]. It has been argued that the last surviving Holocene steppe bison population, represented by the Rauchua bison (haplogroup C) "was not a direct descendant of the major population that dominated Siberia during MIS3 and MIS2" [1]. Our phylogenetic time-tree analysis revealed that haplogroup C separated from the main bison clade first, at around 326 ky BP (95% HPD: 397-264 ky BP) (Figure 2b), and the TMRCA of haplogroup C was estimated at 130 ky BP (95% HDP: 172-94 ky BP), indicating that this haplogroup represents an ancient bison lineage, which appeared no later than the late Middle Pleistocene. Given this time scale, it remains an open question not only where it originated, but also, where the population carrying it survived for several hundreds of thousands of years.
However, our admixture analyses, albeit based on limited data, suggest a relatively high similarity between Siberian steppe bison carrying haplogroup B mitogenomes and our Chinese haplogroup C specimen. The PCA also suggests that the analyzed samples of the genus Bison show at least some phylogeographical structure (Figure 6), since the two steppe bison samples from Siberia (haplogroup B) and northeastern China (haplogroup C) fall closely together, whereas one steppe bison from Yukon (haplogroup A), for which enough nuclear DNA data are available, is genetically more distant from the analyzed steppe bison samples from Asia.
This finding not only contradicts the results of the maternal phylogenetic trees (Figures 2 and 3), it also questions the interpretation that haplogroup C represents a distinct population as previously suggested [1]. Indeed, it has been pointed out before that caution is warranted in the interpretation of single genetic loci, as the stochastic nature of the genealogical process can affect inference from such loci [62]. Clearly, high-quality whole paleogenome sequences from a sufficient number of geographically and temporally well-spread steppe bison samples will be needed to get a better understanding of the evolutionary history of this widespread Pleistocene species.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/genes13101684/s1, Figure S1: Estimated mitochondrial endogenous fragment length distributions for the 3 bison samples analyzed in this study, Figure S2: Mitochondrial cytosine deamination frequency inferred from the 3 bison samples analyzed in this study. X axis represents the position from 5 (left) and 3 (right) read ends. Red line corresponds to C to T substitutions, blue line to G to A substitutions. Double stranded libraries were constructed from CADG456 and CADG465, while a single stranded library was obtained from CADG465, Figure S3: Coverage plots of the mitochondrial genomes obtained in this study. The three black lines represent the mean coverage of the three mitochondrial genomes, Figure S4: Estimated nuclear endogenous fragment length distributions for bison sample CADG467, Figure S5: Nuclear cytosine deamination frequency inferred from bison sample CADG467. The X axis represents the position from 5 (left) and 3 (right) read ends. Red line corresponds to C to T substitutions, blue line to G to A substitutions, Figure S6: Principal component analysis (PCA) computed using genotype likelihoods with sample CADG467, two steppe bison, two yak, three aurochs, six American bison and ten European bison based on 6,000,104 SNPs, Table S1: The information and mitochondrial mapping results of bison samples analyzed in this study, Table S2: Information of the four overlapping long range PCR (LR-PCR) primer pairs used in hybridization capture library construction, Table S3: Mapping results against different bison mitochondrial genomes as reference, Table S4: Information of the nuclear genomes used in this study, Table S5: Complete mitochondrial genomes and ages of steppe bison, American bison and yak downloaded from NCBI, Table S6: D-loop sequences of steppe bison used in this study, Table S7: nucleotide substitutions in the mitochondrial genome that were unique to the individuals CADG456, CADG467 and the mitochondrial reference KR350472, Table S8: Mitochondrial genomes of three steppe bison haplogroups. References [63][64][65][66][67]   Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Two ancient mitochondrial genomes can be found under the Gen-Bank Accession codes OL741537 (CADG456) and OL741538 (CADG467), The address is as follows: GenBank www.ncbi.nlm.nih.gov/genbank, accessed on 10 January 2022.


A distinct role of STING in regulating glucose homeostasis through insulin sensitivity and insulin secretion

Significance The role of STING in maintaining glucose homeostasis remains unknown. Herein, using global and β-cell–specific STING knockout mouse models, we revealed a distinct role of STING in the regulation of glucose homeostasis through β-cells and peripheral tissues. Specially, while global STING knockout beneficially alleviated insulin resistance and glucose intolerance induced by high-fat diet, it surprisingly impaired islet glucose-stimulated insulin secretion (GSIS). Further analyses revealed that STING deficiency down-regulated expression of β-cell key transcription factor Pax6, impairing Pax6 nuclear localization and binding activity to the promoters of its target genes, including Glut2 and Abcc8, causing impaired GSIS. These data highlight pathophysiological significance of fine-tuned STING signaling in β-cells and insulin target tissues for maintaining glucose homeostasis.

Supplemental Materials and Methods
Intraperitoneal Glucose and Insulin Tolerance Tests. The intraperitoneal glucose tolerance testing (IPGTT) was performed by intraperitoneally injection of 2 g/kg body weight of glucose after 6-hour fasting. Tail vein blood glucose was monitored at 0, 15, 30, 60, 90, and 120 min. Blood glucose levels were determined using an automatic glucometer (Roche, ACCU-CHEK, Germany) and glucose test results higher than 33.3 mmol/L were recorded as 33.3 mmol/L. The area under the curve (AUC) was calculated using the trapezoidal rule. At the same time, blood samples were collected at individual time points and serum insulin levels were determined with insulin ELISA test kits (EZassay, China) according to the manufacturer's protocol. For intraperitoneal insulin tolerance test (IPITT), the mice were intraperitoneally injected with 0.75 U/kg body weight of insulin after a 4-hour fasting. All experiments described above were performed on age-and sex-matched cohorts.
Islet perifusion assays. For islet perifusion assays, islets were isolated and incubated overnight in RPMI medium containing 10% FBS, 50 μg/mL streptomycin, and 50 U/mL penicillin at 37°C in a 5% CO2 humidified environment. For each assay, 150 islet equivalents [IEQ] (approximately 100 size-matched islets) were used. Islets loaded in the perifusion chamber were first perifused with KREB's buffer containing 2.8 mmol/L glucose for 30 min, followed by perifusion with KREB's buffers containing 16.7 mmol/L glucose, and 2.8 mmol/L glucose. Flow rate of fluid input (2.8 mmol/L glucose and 16.8 mmol/L glucose in KREB's buffer) was 1 ml/min and fractions were collected every minute. Insulin secretion at each time point was quantified using HTRF insulin assay kit (62INSPEC, Cisbio) according to manufacturer's instructions, and normalized to the insulin content.
Whole islet calcium imaging. After being washed twice with KRBH buffer plus 2 mmol/L glucose, islets were loaded with calcium indicator Fluo-4 AM (Thermo Fisher Scientific, F14201) in KRBH buffer containing 5.6 µmol/L Fluo-4 AM, 0.4‰ Pluronic F127 (Sigma-Aldrich, P2443), and 2 mmol/L glucose for 45 min. Islets were then washed three times with KRBH buffer plus 2 mmol/L glucose, and seeded on 35 mm Glass-Bottom Dish (MatTek) for 15 min. The dish was placed into a 37°C humidified cell culture chamber containing 5% CO2 and mounted on a Zeiss LSM 800 inverted confocal microscope. A solid-state laser provided wavelengths of 488 nm to excite Fluo-4. Images were acquired every 10 seconds before and after adding glucose to the final concentration of 16.8 mmol/L to monitor the dynamic changes of intracellular calcium levels, and analyzed with Image J 1.52q software (NIH).
Biochemical Assays. Blood samples were obtained from mice, and serum was separated by centrifugation. Serum total cholesterol (TC), triglyceride (TG), low-density lipoprotein-cholesterol (LDL-C) and high-density lipoprotein-cholesterol (HDL-C) were determined by assay kits (Nanjing Jiancheng Bioengineering Institute, China) using microplate reader (Thermo Scientific). Serum glucagon levels were determined using glucagon ELISA test kits (Ruixin Biotechnology Co., Ltd., China) according to the manufacturer's protocol.
Immunofluorescence and Immunohistochemical Staining. Paraffin-embedded pancreases were cut into 5 μm sections. Immunohistochemical staining were performed using appropriate antibodies as previously described (25) Transmission electron microscopy (TEM). Fresh pancreases from 8-week-old mice were fixed with 2.5% glutaraldehyde. The pancreases were cut into pieces and then fixed with 2.5% glutaraldehyde in 0.1 M phosphate buffer (pH 7.4) for 16 h and subsequently treated with 1% osmium tetroxide at 4°C for 2 h. The specimens were dehydrated in a graded ethanol series and then embedded in EMBed 812. Ultrathin sections were obtained using the ultra microtome (Leica EM UC7, Germany) and stained with 2% uranium acetate saturated alcohol solution and 2.6% Lead citrate and observed by TEM (HT7700, Japan).

RNA Extraction, cDNA Synthesis, and Real-Time Quantitative PCR
(RT-QPCR). Total RNA of islets was extracted using TRIzol reagent (Invitrogen, USA) and cDNA was synthesized using a PrimeScriptTMRT reagent Kit with gDNA Eraser (Takara, Japan, Cat: RR047A) following the manufacturer's instructions.

RT-QPCR was conducted with the TB Green Premix Ex TaqTM (Tli RNaseH Plus)
(Takara, Japan, Cat: RR420A). GAPDH was used as standards. We normalized the data to Control = 1 for fold change which were calculated by the 2 -△△ Ct method.
Primers were designed and chemically synthesized by TSINGKE Biological Technology (Beijing, China) Co., Ltd. The forward and reverse primer sequences used in this study were listed in supplemental Table S1. pmol siRNA mixed with transfection reagent. After 6 h incubation, the supernatant was discarded and was replaced with normal medium. The efficiency of RNA knockdown was evaluated by RT-QPCR and western blotting. All RNAi experiments were repeated at least 3 times.

Protein Extraction and Western Blot
Statistics. All of the experiments above were performed independently at least three times. All data are presented as the means ± SEM. Graph-Pad Prism 8 software was used for all statistical analyses. We analyzed data using unpaired Student's t test and used the traditional threshold P < 0.05 to declare statistical significance.

Supplemental figures
Supplemental Figure S1. Global deficiency of STING alleviated HFD-induced hyperlipidemia. Blood samples were used to measure the levels of lipids, including triglyceride (TG) (A), total cholesterol (TC) (B), Low-Density Lipoprotein Cholesterol (LDL-C) (C), and High-Density Lipoprotein Cholesterol (HDL-C) (D). n ≥ 5 mice/group. Values are shown as mean ± SEM. WT + ND vs WT + HFD, * P < 0.05, ** P < 0.01, WT + HFD vs STING -/-+ HFD, # P < 0.05, ## P < 0.01, ns, not significant. Figure S2. The expression of insulin gene and amount of insulin content in islets of global STING knockout mice. (A) mRNA levels of Ins1 and Ins2 in islets from 16-week-old male mice fed with ND or HFD were measured by RT-QPCR. (B) Freshly isolated islets from 16-week-old male mice were directly lysed. Western blots were performed with anti-insulin antibody. (C) The quantification of insulin content from at least three independent experiments shown in Fig. S2B. The relative level of insulin content in islets of control mice fed with ND were set as 1. n ≥ 3 mice/group. Values were shown as the mean ± SEM. WT + ND vs WT + HFD, *** P < 0.001, WT + HFD vs STING -/-+ HFD, # P < 0.05, ns, not significant.