High throughput SARS-CoV-2 variant analysis using molecular barcodes coupled with next generation sequencing

The identification of SARS-CoV-2 variants across the globe and their implications on the outspread of the pandemic, infection potential and resistance to vaccination, requires modification of the current diagnostic methods to map out viral mutations rapidly and reliably. Here, we demonstrate that integrating DNA barcoding technology, sample pooling and Next Generation Sequencing (NGS) provide an applicable solution for large-population viral screening combined with specific variant analysis. Our solution allows high throughput testing by barcoding each sample, followed by pooling of test samples using a multi-step procedure. First, patient-specific barcodes are added to the primers used in a one-step RT-PCR reaction, amplifying three different viral genes and one human housekeeping gene (as internal control). Then, samples are pooled, purified and finally, the generated sequences are read using an Illumina NGS system to identify the positive samples with a sensitivity of 82.5% and a specificity of 97.3%. Using this solution, we were able to identify six known and one unknown SARS-CoV-2 variants in a screen of 960 samples out of which 258 (27%) were positive for the virus. Thus, our diagnostic solution integrates the benefits of large population and epidemiological screening together with sensitive and specific identification of positive samples including variant analysis at a single nucleotide resolution.

Introduction
The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is an RNA virus that causes the coronavirus disease 2019 . The virus was first identified in Wuhan (China) in December 2019, spread rapidly and was declared a pandemic in March 2020 by the World Health Organization. Today, more than two years later, COVID 19 has already resulted in over 400 million confirmed infection cases and over 5.7 million deaths across the world. COVID-19 symptoms are varied, ranging from none to severe according to the Centers for Disease Control and Prevention (CDC). Many patients develop moderate pneumonia and in the most critical cases respiratory failure and multiorgan disfunction occur. The world average mortality is 2%.
Prevention of viral transmission is essential as COVID-19 spreads through the respiratory route. Transmission prevention is achieved through social distancing, masks and hand washing in addition to early identification of potential infected individuals by different diagnostic methods.
SARS-CoV-2, a positive-sense single-stranded RNA virus is a highly pathogenic member of the coronavirus family [1]. The sequence of its whole genome was published (GenBank no. MN908947), encoding 9860 amino acids. It is composed of genes that express both structural and nonstructural proteins.
Several methods have been developed in order to assess the infection status of individuals. Most of them rely on the detection of viral RNA.
The gold-standard method for diagnosis of COVID-19 is real-time reverse transcription polymerase chain reaction (one step rRT-PCR) [3,4], which detects specifically the presence of a number of viral RNA fragments. This test is performed on nasopharyngeal swabs or saliva samples [5] and is widely applied and considered very robust. Results are generally delivered within 24 hours.
While, real time PCR combined with robotic systems, can be used as a high throughput method, the number of amplicons that can be detected are limited, and fluorescent probes and special reagents are needed in order to increase the sensitivity of the method. Moreover, as the global spread of the virus continues despite the use of preventative measures and positive effects of vaccination [6], SARS-CoV-2 mutations generate different viral variants, some with the potential to become more virulent and contagious or resistant to vaccination. The formation and spread of new variants, requires the development and implementation of cost-effective diagnostic methods that can both detect infection and new single nucleotide variants (SNVs) in real time to warn authorities prior to their spread.
Next-generation diagnostic solutions, which are based on the current PCR test coupled with next generation sequencing (NGS) offer major advantages, such as significantly larger scale testing (several hundred-fold scale up versus the currently available test), reduced cost per test, reduced amount of reagents per test and potential readout of viral sequence variation.
In NGS-based methods, individual samples can be uniquely labeled with molecular barcodes and multiple fragments can be amplified in parallel and pooled together for high throughput detection and processing of thousands of specimens together. A few different protocols have been developed [7,8] and one of them has even received an Emergency Use Authorization from the FDA [9]. As next-generation sequencers are now widely available, including inside hospitals, this platform offers a key solution for COVID-19 mass detection worldwide.
Here we propose a solution that allows high throughput testing by barcoding and pooling test samples, up to hundreds at a time, in one well (as illustrated in Fig 1). The protocol combines multiplexed PCR, patient specific molecular barcoding, pooling, NGS, bioinformatics and machine learning analysis, for identification of positive individuals by detection of three different viral genes and one human gene (as internal control) from a pool of thousands of specimens, enabling sensitive and accurate molecular-based diagnosis of COVID-19. Using this method, we have identified six known and one new viral-variants in a mass screening analysis of 960 samples. The test can prove to be vital in the identification of various viral variants as part of population and epidemiological screening.

Barcodes
A bank of 2133 unique 10-bases and about 21,000 unique 12-bases barcodes was generated.
For the purpose of this research, the most distinguishable 96 barcodes were selected for the experiments (see Methods).

Pilot
Positive and negative SARS-CoV-2 standards (Bio-Rad) were run as a pilot with 17 barcodes and a multiplex of two viral genes (N1 and E) and one human internal control (RNaseP). Following multiplexing, amplicons were pooled. A single library was prepared and ran on a MiSeq (Illumina) instrument. Following demultiplexing, number of reads for each barcode was obtained. RNase P reads were used for internal control and were required to be positive in order to consider the result as valid.
The viral genes were used for the diagnosis and detection of SARS-CoV-2 specifically in each specimen. For each barcode, the associated viral reads were counted and the results are presented in Fig 2a. A clear distinction was observed between positive and negative samples.

Barcode stability
Following this experiment, the stability of 96 barcodes was tested using positive SARS-CoV-2 standard as an identical template for each barcode. Furthermore, an additional viral region (N2) was amplified and added to the multiplex. The stability of 96 barcodes is presented in Fig  2b and 2c. No specific influence of the barcode sequence was noticed.
Moreover, the N2 viral fragment added additional information regarding the viral identification and can provide an improvement in diagnosis of the specimens.

New genes
In order to improve specificity and sensitivity, additional viral regions for barcoded amplification were examined. Each new set of primers was first tested on positive and negative SARS--CoV-2 standards. The count of reads is presented in Fig 2d. Human RNase P served as the internal control gene. A threshold was set to determine the genes that should be part of the multiplexes. This threshold represents the minimal number of reads of positive samples for a specific target in order to include this specific amplified amplicon in the multiplex assay. As the SC2 assay failed to yield any reads for both positive and negative templates, it was not comprised in the further multiplexes. Additionally, we observed a weaker signal for ORF8 in comparison to ORF1. Thus, N1, N2, E, and ORF1 were selected as the multiplex target fragments for the next step.

Multiplex optimization
In order to choose the most effective combination of genes targeted, the following multiplexes were tested on 96 clinical RNA specimens for optimization (Fig 3). Multiplex 1 included viral genes N1, N2 and E; multiplex 2 included viral genes N1, N2 and ORF1; and multiplex 3 included viral genes N1, E and ORF1. Human RnaseP was included in all three mixes as internal control.
Our results demonstrate that when genes N1 and N2 are multiplexed together, a non-specific DNA fragment is generated due to the physical proximity of the two amplicons (unshown data). This 944bp DNA fragment was formed as the elongation product of the N1 forward primer and N2 reverse primer. As this fragment is added to the library preparation and run together with all the amplicons, it can lead to a competition in the NGS run and to a lower number of reads that are obtained from the fragments of-interest.

Detection of SARS-CoV-2 in a 960 samples pool
In order to scale up our technology, 960 clinical RNA specimens were tested using multiplex 3. Amplicons were pooled (each plate separately), libraries were prepared for each plate with a different index and all ten libraries were sequenced together using the Next-Seq instrument (illumina). Following demultiplexing, obtained sequences were analyzed and positive samples were identified using the unique and specific barcodes (Fig 4a). The number of viral reads in our pool was in negative correlation with the Ct number obtained by regular real-time PCR. When the Ct was high, the number of viral reads in the pool were low (Fig 4b). Our results indicate that at a Ct higher than 30 the number of viral reads is similar between positive and negative specimens. Three different multiplexes including three viral genes among N1, N2, E and ORF1a genes, and including one human internal control gene RNAseP were tested separately with the same clinical specimens in order to compare between them and determine the optimal combination to be use for the test. Using a Gradient Boosting machine learning classification model, results were obtained with an accuracy of 93.3%, precision of 91.7%, sensitivity of 82.5% and specificity of 97.3% (Fig 4c) in a validation set. We further trained a model defining Ct<30 (N-gene) specimens as positive. The sensitivity was 100% with 98.5% specificity and a positive predictive value (PPV) of 94.7% in a validation set.

Variant detection
In addition to the detection of the virus with high specificity and sensitivity (especially for Ct <30), SARS-CoV-2 barcoded pooling provides a platform for detection of single nucleotide mutations (Fig 4d) in pre-determined gene regions. We have sequenced 3 viral fragments for detection of the virus in all 960 specimens. Following NGS analysis of the multiplexes, 7 single nucleotide variants were identified. By comparing the observed mutations to SARS_CoV-2 databases using Blastn in betacoronavirus genbank-https://blast.ncbi.nlm.nih.gov/) we have identified six known and one newly found missense mutation. Among the 258 SARS-CoV-2 PCR-positive samples, 30 carried the same N-gene missense mutation and 6 carried an additional synonymous substitution in ORF1a (Fig 4d). Eight PCR-positive and two PCR-negative samples were found with the five remaining single nucleotide variants, one variant per sample.

PLOS ONE
SARS-CoV-2 variant analysis using molecular barcodes Discussion SARS-CoV-2 (severe acute respiratory syndrome coronavirus 2), a beta coronavirus, is the novel coronavirus that caused the COVID-19 outbreak originating in Wuhan, China in 2019. This virus has become a worldwide pandemic requiring immense and rapid viral detection methods throughout the world. The classic method for detection of infection involves sampling nasopharyngeal swabs from patients in order to detect viral infection. The swab is placed in transfer buffer and then lysis buffer is added to lyse the potentially infected cells. RNA is extracted from the solution then subjected to a one-step RT-PCR reaction. The RNA is first converted to cDNA by reverse transcription which is immediately followed by Multiplex quantitative Real time PCR (qPCR) with primers and a probe specific to the virus to detect viral presence.
Generally, an internal control is added in order to eliminate the possibility of false negative results and includes detection of a human gene as a proxy for host RNA content. However, some assays only amplify an exogenic internal control gene (added during the extraction step) which provides evidence for the validity of the qPCR reaction yet, is not reliable enough to exclude false negatives.
Here we describe a method we developed for combined screening of an infectious agent and simultaneously determine the presence of a genomic variant in a population, by using a DNA-barcoded sample from each subject in a population. The method includes performing one-step reverse transcription and amplification of RNA extracted from the samples with barcoded primers, pooling and sequencing the obtained amplicons, and identifying the infected individual out of all tested samples in the pool via the subject-specific barcode. The different primers target three viral genes as well as one human endogenous gene. The detection of an endogenic human internal control assay that is specific to each specimen, even though they are sequenced as a pool, is of extreme importance and represents one of the major advantages of this technique.
As the virus has spread globally, it has naturally acquired thousands of mutations. Most of these are deleterious or neutral, however, some confer selective advantage and have given rise to many viral variants that became the dominant viral strain in many countries. Some of these variants as they have been shown to increase infection rate, cause higher mortality and most importantly may evade the immune response obtained by vaccination [10]. It is therefore essential that health care authorities collect evidence on the extent and spread of such viral variants within the population. However, currently the RT-PCR method used to identify SARS--CoV-2 infection, does not enable the identification of sequence variations in specific genomic locations. Additionally, while NGS enables variant analysis and lineage tracking, it is not yet set as a standard mass-detection method which can be efficiently used for population screening.
We provide a method that can be used to identify and screen specific loci on the viral genome which are of concern at the population level. For SARS-CoV-2, one of the loci can be the entire 'Spike' sequence, while other areas can potentially be identified as more variants and lineages evolve. We suggest that once new SNVs are identified in this screen, specific samples can be identified based on their unique molecular barcode followed by full genome sequencing to provide a comprehensive sequence-analysis for lineage tracking.
Our method utilizes amplified cDNA produced from the infectious agent and at least one SNV in a sequencing read, to determine via the subject-specific barcode whether the subject is infected and simultaneously enables the identification of a genomic variant. The method presented here can potentially be modified using primers designed upstream of the variable S gene to enable the detection of the known variants (UK-B.1.1.7 lineage, South Africa-B.1.351 lineage, Brazil-P.1 lineage, NY-B.1.526 lineage) whose sequences are mutated in this specific locus (Spike gene). Our aim was to integrate sequence-analysis which can provide information at a population level by searching for SNVs on multiple short sequences from hundreds of test samples and by positively identifying SARS-CoV-2 variants at an individual level. Finally, sensitivity and specificity can be improved by applying machine learning algorithms. This integrative approach can be used to identify SARS-CoV-2 as a stand-alone method or, to be used as an additional information layer which screens the RT-PCR positive samples to identify the viral variants.
In order to generate a SARS-CoV-2 whole genome read with high sensitivity and specificity (98% and 97%), at least 1,000 copies of the viral genome per milliliter are required [9]. While full-genome analysis provides information on potential changes in the entire viral RNA sequence, this approach limits the potential amount of individual screening tests that can be generated in a population. Thus, the number of samples providing information from the whole viral genome limits its usage as a mass-screening tool in this current pandemic.
To overcome this limitation, additional NGS-based mass screening approaches have provided proof that such methodologies can be used to identify new variants, detect genetic epidemiology and follow viral lineages [11]. However, for low viral counts, the sensitivity and specificity of this system are not clear. The method we described here, further expands this high throughput approach through a unique individualized barcoding system, and a machine learning module to further increase sensitivity and specificity.
Based on our findings, we propose coupling individualized sample identification and pooling with machine learning approaches to optimize current NGS resources to quickly identify viral variants in the population. This approach combines the benefits of high throughput screenings by focusing on specific areas-of-interest in order to increase the number of individual screening tests and employs machine learning algorithms to increase diagnostic accuracy.
Population mass-screening diagnostic methods are required to be cost-effective. Based on a calculation of reagents and labor, we have generated a cost estimate to analyze each sample using the suggested method for SARS-CoV-2 variant detection. Given the large-capacity and pooling nature, which stand at the basis of this diagnostic method, it is estimated that an individual test would roughly cost 1-2 USD (starting from the extracted viral RNA). This estimate provides organizations which currently run large PCR-based SARS-CoV-2 screening assays, a viable cost-effective option to add an important variant-analysis layer to their data.
In conclusion, in comparison to the gold standard real time PCR method, that identifies amplified cDNA produced from the infectious agent and is only able to provide a binary answer (SARS-CoV-2 positive or negative), the method presented here, enables specific diagnosis of SARS-CoV-2 variants within a population of thousands. As this solution can either provide a stand-alone identification method for SARS-CoV-2 or be integrated with the current RT-PCR viral test, we propose adopting this technique to large population and epidemiological screening practice.

Clinical specimens
Oro-and nasopharyngeal swabs were collected from patients by trained health workers with personal protective equipment and directly placed into e-swab tubes (COPAN) for conservation and transport. The study number 0196-20-EMC was granted Ethics approval by the head

RNA extraction
Following sampling and storage of the oro-and nasopharyngeal swabs, nucleic acids were isolated and purified from swab specimens using STARMag 96 X 4 Universal Cartridge Kit (Seegene) on Star or Starlet instruments (Hamilton) following the manufacturer's instructions.

Primer sequences
The targets analyzed were amplified using previously reported sets of known primers, while unique barcodes were added to the forward primer to enable pooling of the products before sequencing. These primers sequences are summarized in Table 1.
The complete table of primers including barcodes appears in the supplementary data of this publication.

Singleplex RT-PCR
The purified nucleic acids were reverse transcribed into cDNA and amplified using Qscript XLT one-step RT-PCR (Quantabio) and specific primers for each viral and human amplicons. Each forward primer included a 10-nucleotide barcode at the 5'-end, which was unique for each sample tested and enabled pooling of the PCR products prior to sequencing. The total volume of reaction mix (25 μl) contained the following: 12.5 μl of 2X concentrated ToughMix; 0.4 μM reverse primer, 0.4 μM barcoded forward primer; 1 μl of 25X qScript XLT and 5 μl of SARS-CoV-2 positive and negative standard (BioRad). The one step RT-PCR was conducted in a C1000 touch thermocycler (BioRad) as per the following: cDNA Synthesis (48˚C, 20 min), Initial denaturation (94˚C, 3 min), PCR cycling -25 cycles (Denaturation: 94˚C, 15s; Annealing: 60˚C, 40s; Extension: 72˚C, 30s).

Multiplex RT-PCR
The purified nucleic acids were reverse transcribed into cDNA and amplified using Qscript XLT one-step RT-PCR (Quantabio) and specific primers for 3 viral amplicons (N1, N2, E in multiplex 1; N1, N2 and ORF1 in multiplex 2; and N1, E and ORF1 in multiplex 3) and one human (RNaseP as internal control). Each forward primer included a 10-nucleotide barcode at the 5'-end, which is unique for each individual tested (the barcode generation is described in the next paragraph). The total volume of reaction mix (25 μl) contained the following: 12.5 μl of 2X concentrated ToughMix; 0.4 μM each reverse and barcoded forward primer; 1 μl of 25X qScript XLT and 8 μl of the RNA sample. The one step RT-PCR was conducted in a

Barcode generation
Sequences of unique ten-base barcodes were created using DNABarcodes R package [13]. Barcodes were generated in such a manner that they were at least three Sequence-Levenshtein distance apart. In addition, to meet the required chemical properties, sequences that contained triplets that were self-complementary or that showed an unbalanced GC ratio were filtered out. Based on distance metrics between each pair of barcodes, most distinguishable 96 barcodes were selected for the experiments based on the distance between every possible pair among the 96 barcodes. A pool of 960 clinical samples was tested using the technology. Ten 96 plates with 96 different specimens (about 25% positive) were amplified using the barcoding-multiplex PCR. Wells from each plate were pooled together and a library was prepared using unique index. All ten libraries were sequenced together using the Next-Seq instrument (illumina).

Next generation sequence analysis
Reads processing pipeline. First, a raw FASTQ file was demultiplexed by DNA-barcodes using FASTX barcode splitter, i.e., splitting into numerous files based on DNA-barcode matching, one FASTQ per barcode (FASTX-toolkit version 0.0.14, http://hannonlab.cshl.edu/ fastx_toolkit/). One mismatch [13] between the known and the observed sequence of a barcode was allowed. Then, each barcode specific FASTQ file was trimmed using FASTX trimmer, leaving 41 bases of each read, from the 11th base to the 51st base, representing a part of the designated primers. The trimmed reads were aligned to N (nucleocapsid phosphoprotein), E (envelope protein) and ORF1ab (ORF1a; ORF1b polyproteins) SARS-CoV-2 genes and human RNASE P (from NCBI Genbank) using BBMap version 38.18, which is a short-reads aligner for DNA or RNA-seq data (http://sourceforge.net/projects/bbmap/). Finally, the set of reads were quantified using Salmon version 0.14.2 [14].
Variant calling pipeline. For the purpose of variant calling, the barcodes sequences were trimmed at the 11th base, leaving the entire reads. Quality control and preprocessing of FASTQ files were performed using fastp version 0.20.1 [15]. Each FASTQ per sample was mapped to SARS-CoV-2 reference genome (GenBank no. MN908947) using bwa version 0.7.17-r1188 [16]. BAM files were sorted by SAMtools version 1.9 [17], and variants were called using bcftools version 1.9 [18], filtered by quality-measure of at least 20 and read depth of at least 50.
Machine learning models. Classification models were implemented using scikit-learn module in Python, version 0.23.2 [19]. For the initial model, the samples were randomly divided into 60% train set, 20% validation set and 20% test set. SARS-COV-2 genes counts, normalized for total read number per sample, were selected as features, while the label of a positive/negative diagnosis was based on traditional RT-PCR. Hyperparameters were selected after an exhaustive grid search, with 5-fold cross validation. Selected hyperparameters were those which provided the maximal Area Under the Receiver Operating Characteristic (ROC) curve (AUC) scores. Then, the model was evaluated by predicting the diagnosis in the test set. Additional models, i.e., based on different gene combinations, were methodologically similar.

Code availability
The R code used to generate results presented in Figs 2-4 and Python workflow used to build the machine learning models are available on https://github.com/Adarya/SARS-CoV-2_NGS_ Barcodes.


Functional analysis of the nonstructural protein NSs of tomato zonate spot virus

Tomato zonate spot virus (TZSV), a member of the genus orthotospovirus, causes severe damage to vegetables and ornamental crops in southwest China. The NSs protein is an RNA silencing suppressor in various orthotospovirus like TZSV, but its mechanism and role in virus infection are poorly understood. Here, we observed that an NSs-GFP fusion protein was transiently expressed on the plasma membrane and Golgi bodies in Nicotiana benthamiana plants. The TZSV NSs gene was silenced and infiltrated into N. benthamiana and N. tabacum cv. K326. RT-qPCR and Indirect enzyme-linked immunosorbent assay (ID-ELISA) showed that the transcription and the protein expression of the NSs gene were inhibited by more than 90.00%, and the symptoms on silenced plants were alleviated. We also found that the expression of the Zingipain-2-like gene significantly decreased when the NSs gene was silenced, resulting in co-localization of the NSs-GFP and the Zingipain-2-like-mCherry fusion protein. The findings of this study provide new insights into the mechanism of silencing suppression by NSs, as well as its effect on systemic virus infection, and also support the theory of disease resistance breeding and control and prevention of TZSV in the field.

If the data are held or will be held in a public repository, include URLs, accession numbers or DOIs. If this information will only be available after acceptance, indicate this by ticking the box below. For example: All XXX files are available from the XXX database (accession number(s) XXX, XXX.).
• If the data are all contained within the manuscript and/or Supporting Information files, enter the following: All relevant data are within the manuscript and its Supporting Information files.
• If neither of these applies but you are able to provide details of access elsewhere, with or without limitations, please do so. For example: Data cannot be shared publicly because of [XXX]. Data are available from the XXX Institutional Data Access / Ethics Committee (contact via XXX) for researchers who meet the criteria for access to confidential data.
The data underlying the results presented in the study are available from (include the name of the third party • Introduction 28 Manucript without track changes Click here to access/download;Manuscript;Manucript -2021.11.14.docx Tomato zonate spot virus (TZSV), the dominant species, which belongs to the 29 orthotospovirus genus of the Bunyaviridae family in Yunnan Province, is transmitted 30 by thrips [1]. In recent years, TZSV has been prevalent in Yunnan, Guizhou, Guangxi, 31 and other regions of China and Southeast Asia. TZSV infection is often associated with 32 severe disease symptoms, including concentric rings and spots on fruits and necrosis of 33 leaves. TZSV has a very broad host range, infecting more than 20 types of economically 34 important crops and weed species, which belong to 7 families [2,3]. This prevalence 35 has not only led to production losses and quality problems for important vegetable and 36 ornamental plants in Southwest China but also seriously threatened economically 37 important crops for local farmers [4]. 38 TZSV consists of spherical, enveloped particles that become distributed in the 39 cytoplasm and the endoplasmic reticulum in mesophyll cells [1,5]. Like other members  PCR (RT-PCR) and the special primers (Table S1). The PCR fragments of NSs and Zingipain-2-like-mCherry, respectively. The Golgi marker ST52-mCherry was 134 amplified from the total RNA isolated from both tobacco and Arabidopsis [21].

136
The pTRV-pTV00-NSs constrcture was transformed into Agrobacterium tumefaciens  Confocal laser scanning microscopy and co-localization 155 The leaf epidermis was dissected from the areas of the agroinfiltrated N. 156 benthamiana leaves and placed in water between two coverslips. The confocal images 157 were captured with the inverted TCS SP8 and 10× water immersion objective lenses.

158
GFP was excited at a wavelength of 488 nm, and emissions were captured at 497-520  buffer was used as a blank control.

185
NSs was localized to the PM and Golgi bodies 186 Previous studies have shown that confocal laser scanning microscopy (CLSM) 187 was used to analyze the protein localized in living cells [22]. N. benthamiana is also a  be used as a calibration product to determine the copy numbers of the gene.

215
The PDS gene was used as a positive control to ensure the success of silencing. hypersensitive response [25,26]. To confirm whether the host factor for the Zingipain- The importance of NSs for tospoviral infection in plants was first discovered in 263 the early 1990s. A higher virulence of TSWV isolates and more severe symptoms were 264 observed with the elevated levels of NSs expression [27]. The TSWV NSs protein acts 265 as a suppressor of RNA silencing through binding small and long dsRNA and 266 suppresses short and long-distance viral accumulation and movement [28]. It also 267 represents the Avr factor of the Tsw resistance gene in pepper [29,30]. NSs of TSWV plants [37]. In this experiment, the temperature was strictly controlled, and thus, the 307 gene was silenced at relatively high levels, and the duration was relatively long.

Materials
The pTV00, pBINTRA, and pTV00-PDS vectors were provided by the Professor Jianqiang Wu's's laboratory, at the Kunming Institute of Botany, Chinese Academy of Sciences (KIB CAS) [20]. The vector of pCAMBIA-GFP, pBI121-mCherry were obtained from our lab in the institute of biotechnology and germplasm resources. The TZSV YN-Chili isolate was collected from the infected tomato field tomato plants in Yuanmou, Yunnan Province, China, and was maintained onin N. benthamiana [1]. N.  (Table S1).

TZSV inoculation
The frictionalfriction-based inoculation method was used to artificially infect N. tabacum cv. K326 and N. benthamiana plants at the 6 six-leaf stage. TZSV-infected N.
benthamiana leaves were homogenized in the PBS buffer (100 mg/mL) (containing 137 mM NaCl, 1 mM KH2PO4, 8 mMNa2HPO4·12H2O, and 3 mM KCl), and there afterafter, which the homogenate was applied uniformly to 3 leaves per plant (1 mL). Ten minutes after the TZSV inoculation, the inoculated plant leaves were rinsed with ddH2O.
PlantsPlant leaves inoculated with PBS buffer were used as controls, and five replicates were used for each sample. AfterFive days after inoculation for 5 days, when symptoms appeared, five replicates were used for each sample. The amplified products were subsequently analyzedimaged byusing a UVP gel-imaging system.

Plasmid constructs
(ⅰ) The fragments of the NSs gene to-be-silenced NSs gene were cloned in theinto pTV00 vectors. PCR was used to amplify the desired fragments with specific primers (Table S1)  (ⅱ) The sequences of NSs and Zingipain-2-like genes sequence was were amplified from the total RNA isolated from tobacco plants infected withby TZSV using reverse transcription-PCR (RT-PCR) and the special primers NSs-G-F/NSs-G-R (Table S1). The NSs-N and NSs-C PCR fFragments of NSs and Zingipain-2-like genes were digested with endonuclease Nco1 and Spe1 and inserted into vector of pCAMBIA-GFP and pBI121-mCherry using the same restriction sites to obtain pCAMBIA-NSs-N-GFP and pBI121-Zingipain-2-like-mCherrypCAMBIA-NSs-C-GFP, respectively. The Golgi markermarkers ST52-mCherry was amplified from the total RNA isolated from both tobacco and Arabidopsis [21].

Infiltration of the VIGs vector in Tobacco Leaves
The pTRV-pTV00-NSs vector constrcture was transformed into

Confocal laser scanning microscopy and co-localization assays
The leafLeaf epidermis waswere dissected from the areas of the agroinfiltrated leaf area of N.N. benthamiana leaves and placedmounted in water between two coverslips. The confocalConfocal images were captured with thean inverted TCS SP8 and 10× water immersion objective lenses. GFP was excited at a wavelength ofwith 488 nm, wavelength and emissions were captured at 497-520 nm captured. Moreover, mCherry waswere excited atwith a wavelength of 561 nm wavelength and emissions were captured at 585-615 nm captured. Images were processed using athe TCS SP8 and Adobe (San Jose, CA, USA) Photoshop.

ID-ELISA
The contentscontent of the genes proteins werewas tested by ID-ELISA according to the antibody instructions for antibodies to determine the antiviral activity of the VIGS. Leaves (0.2 g) were homogenized usingwith a mortar and pestle and diluted 1: 3 in a 1/3 dilution with PBS buffer. Crude extracts (100 µL) were added intoto ELISA plate wells and incubated at 37 °C for 2 h. The plate was then washed with PBST buffer.
TZSV NSs rabbit antibodies were diluted in a conjugationconjugate buffer. Afterward, and afterward, 100 µL of AP-conjugated goat anti-rabbit IgG-AP conjugate (Sigma, USA) was added to each well. The color-developing solution was dissolved in pnitrophenyl phosphate disodium hexahy drate (Sigma-Aldrich) in substrate buffer to obtain a final concentration of 1 mg/mL. The absorbance was determined at 405 nm was measure dusing athe ELx 808 microplate ELISA reader (Bio-Tek, USA). Healthy leaves were used as negative controls, PBS buffer was used as a blank control, and TZSV-infected leaves were used as positive controls. The PBS buffer was used as a blank control.

NSs wasis localized to the PM and Golgi bodies
Previous studies have shown that confocal laser scanning microscopy (CLSM) wasalways used to analyzeanalysis the protein localized in living cells [22]. N.N. benthamiana by agroinfiltration. We detected NSs-GFP was then detected to be associated within a structure that was the plasma membrane structures by CLSM (Fig   1 A-C). To determine whether the NSs-GFP bodies are colocalized with the Golgi stacks, we also checked the localization of NSs-GFP for Golgi bodies using the marker ST52-mCherry [23]. At 48 h afterAfter agroinfiltration onof N.N. benthamiana 48 h, we found that the NSs-GFP bodies are colocalized with the Golgi stacks (Fig 1 D-F), suggesting that the NSs-GFP protein was targetedtargets to the Golgi apparatus.

Silencing of the NSs gene
TRV vector wasis widely used to study the interactions between viruses and hosts and the functions of plant genes [24]. In the present study, Specificspecific primers (TZVNSsF1/TZVNSsR1) containing BamHI and HindIII restriction enzyme recognition sites were used to amplify the NSs gene fragments (Fig S1), and the DNA was inserted into athe pEASY-Blunt-Zero vector (TransGen, Beijing) for sequencing to ensure that the base sequences were not mutated.
The copies of the NSs gene were analyzedwas determined by RT-qPCR. The with the pTRV-pTV00-NSs vector before inoculation with TZSV (Fig 2A). Ssevere leaf shrinkage also occurred in positive control plants (inoculated only with only TZSV) exhibited severe shrinkage symptoms (Fig 2B)., and The shrinkage also appeared on the leaves of N. benthamiana plants injected with the pTRV-pTV00-NSs vector before inoculation with TZSV began to display light shrinkage symptoms (Fig 2A) ( Fig 3A, B). Taken together, theThe results indicateindicating that the pTRV-pTV00-NSs VIGS vector was successfullysuccessful constructed, and the NSs gene might be associated with TZSV infection.

Inhibition of the NSs protein expression
The NSs protein levelscontents waswere measured by ID-ELISA at 3 d, 5 d, 7d, and 9 ddays post-inoculation with TZSV,. and found to be significantly decreased The

The dependence of NSs gene silencing depending on the zingipain-2like gene
The Zingipain-2-like gene was a homocysteine protease, which possesses cysteinetype endopeptidase activity and participates in the regulation of plant-type hypersensitive response [25,26]. To confirm whether the host factor offor the  (Fig 4 A-C). Our results clarifiedclarify that the RNA silencing suppressor by (NSs) made closelyhad a close relationship with Zingipain-2-like when TZSV infected the plants.

Discussion
The importance of NSs for tospoviral infection in plants was first discovered in plants in the early 1990s., aA higher virulence of TSWV isolates and more severe symptoms werewas observed whenwith the elevated levels of NSs expression [27]. The TSWV NSs protein which acts as a suppressor of RNA silencing, through binding small and long dsRNA and suppressing suppresses local short and long-distance systemic silencing to enhanced viral us accumulation and movement [28]. It also represents the Avr factor of the Tsw resistance gene in pepper [29,30]. NSs of TSWV directly interacts with MYC2, a regulator of the JA signaling, to disable JA-mediated activation of terpene synthase genes and attenuates attenuate host defenses, increasingincreases the attractivenessattraction of the plants to thrips, toand thus transmittingtransmit the disease [31]. So far, the ability of NSs of TSWV NSs to counteract defense mechanisms mediated bythe RNA silencing defense mechanism in plants has been demonstrated was illuminated,; however, the mechanism of RNA silencing suppression by NSs in TZSV including microbe-associated (MAMP) or pathogen-associated molecular pattern (PAMP)-triggered immunity (PTI) and effector-triggered immunity (ETI). PTI is mediated by the corresponding membrane-anchored pattern recognition receptors (PRRs) inof plants, which serve as the first line of defense against the pathogen. Many plant viruses and their encoded proteins that could inhibit PTI were reported, such as like TSWV the NSs protein of TSWV that suppressed the production of reactive oxygen species (ROS), have been reported and PTI was inhibited [30,32,33]. In this study, we In this study, TRV vectors were used to construct VIGS vectors of the NSs gene tofor the analysis of their the functionsfunction. The results of RT-qPCR, andas well as the plant disease symptoms, showed that the gene replication was inhibited to up to 90%. ID-ELISA showed that the protein contents also significantly decreased. The high efficiency of gene silencing can be verifiedmeasured by sampling and testing immediately atwhen the onset of the disease appears on the plants. At the same time, temperature also hadhas effects on the silencing phenotypes inof plants [37]. In this experiment, the temperature was strictly controlled;, and thus, the gene was silenced at relatively high levels, and the duration was relatively long. TRV has a wide range of hosts, withand there is a significant difference in sensitivity to TRV between species and cultivars [38,39]. For instance, TRV sensitivity testing was carried out on 21 Gerbera gerbera cultivars, and the results revealed that only 5 cultivars showed photo bleached PDS-silencing symptoms on newly developed leaves [40]. In this study, the N. benthamiana leaves were bleached, but there were no phenotypic changes in the leaves of N. tabacum cv. K326, despite the fact that both species belongbelonging to the same genus, indicating that the TRV-VIGS vector exhibited differentdiffering sensitivities to different host species. TRV has a wide range of hosts, with a significant difference in sensitivity to TRV between species and cultivars [38,39]. For instance, TRV sensitivity testing was carried out on 21 gerbera cultivars, and the results revealed that only 5 cultivars showed photobleached PDS-silencing symptoms on newly developed leaves [40]. The VIGS method can be used for reverse genetics studiesresearch and the analysis ofto determine the functions of unknown genes.
In summary, our results presented here revealed that NSs, a suppressor of RNA silencing NSs ofin TZSV, was localized to the PM and Golgi bodies and might be alsoalso be associated with Zingipain-2-like to activate PTI-like responses taking advantage ofusing VIGS and subcellular localization predictionmethods. pTRV-pTV00-NSs +TZSV: means infiltrated pTRV-pTV00-NSs VIGS vectorconstruct and inoculated with TZSV; TZSV: positive control and means only inoculated with TZSV; pTV00-TZSV: means infiltrated pTRV-pTV00 VIGS vector and inoculated TZSV; CK: negative control and means the leaves with no treatment. All values are means ± SE. Means in a column followed by different letters are significantly different at P ≤ 0.05.

Fig 4 The NSs protein in TZSV NSs is co-localized with the Zingipain-2-like protein
NSs-G-F/ NSs-G-R (Table 1). The NSs-N and NSs-C PCR Fragments were digested with Nco1and Spe1 and inserted into pCAMBIA-GFP using the same restriction sites to obtain pCAMBIA-NSs-N-GFP and pCAMBIA-NSs-C-GFP, respectively. The Golgi markers ST52-mCherry was amplified from total RNA isolated from tobacco Arabidopsis [21]" has been changed into "The sequences of NSs and Zingipain-2-like genes were amplified from the total RNA isolated from tobacco plants infected with TZSV using reverse transcription-PCR (RT-PCR) and the special primers (Table S1).
The PCR fragments of NSs and Zingipain-2-like genes were digested with endonuclease and inserted into vector of pCAMBIA-GFP and pBI121-mCherry to obtain pCAMBIA-NSs -GFP and pBI121-Zingipain-2-like-mCherry, respectively. The Golgi marker ST52-mCherry was amplified from the total RNA isolated from both tobacco and Arabidopsis [21]." Answer: the sentence "The fragments of the to-be-silenced NSs gene were cloned in the pTV00 vectors." has been transferred to the location before the sentence "The vectors were then transformed into E. coli DH5α competent cells.".
Answer: The vector of pCAMBIA-GFP was obtained from our lab in the institute of biotechnology and germplasm resources, and insert this sentence before the sentence "The TZSV YN-Chili isolate was collected from infected field tomato plants in Yuanmou".
2) Vector and construct are different. pTRV-pTV00-NSs is a construct, not a vector.
Answer: Inoculation was performed using a 1-mL syringe, and the inoculation solution The reviewer asked where was pCAMBIA-GFP originally obtained from? Please also indicate where in the revised ms such information is provided.
Where is this information in the revised ms?
Where is the information in the revised ms? 2) Line 212-217: the sentence needs to be reorganized, Fig.2A goes before Fig. 2B Answer: Page 8, line 206-209: the sentence "TZSV was inoculated after 5 days, the positive control (inoculated with only TZSV) exhibited severe shrinkage symptoms ( Fig 2B), and the leaves of N. benthamiana plants injected with the pTRV-pTV00-NSs vector before inoculation with TZSV began to display light shrinkage symptoms (Fig   2A)," has been changed into "After 5 days of inoculation with TZSV, the shrinkage appeared on the leaves of N. benthamiana plants injected with the pTRV-pTV00-NSs vector before inoculation with TZSV (Fig 2A). Severe leaf shrinkage also occurred in positive control plants (inoculated only with TZSV) (Fig 2B);".


s2fieldsofstudy: Biology Medicine

text:
Theaflavin Attenuates TBHP-Induced Endothelial Cells Oxidative Stress by Activating PI3K/AKT/Nrf2 and Accelerates Wound Healing in Rats

The treatment of wounds remains a clinical challenge because of poor angiogenesis under the wound bed, and increasingly, the patients’ need for functional and aesthetically pleasing scars. Previous reports have shown that Theaflavin can induce angiogenesis and terminate the progression of ischemic cardiovascular disease, but limited therapy is available for the management of cutaneous wounds. In this study, our in vitro work discovered that human umbilical vein endothelial cells (HUVECs) exposed to Theaflavin can alleviate apoptosis and cell dysfunction induced by tert-butyl hydroperoxide (TBHP). The cellular activity of HUVECs were assessed by cell tube formation, migration and adhesion. Mechanistically, Theaflavin protected HUVECs from TBHP-stimulated cell apoptosis through the activation of the phosphatidylinositol-3-kinase (PI3K)/protein kinase B (AKT)/nuclear factor (erythroid-derived 2)-like 2 (Nrf2) axis, so Nrf2 silencing can partly eliminate the cytoprotective effect of Theaflavin treatment. In in vivo experiments, administering Theaflavin orally can enhance vascularization in regenerated tissues and accelerate wound healing. In summary, our data served as a novel evidence for the wound healing treatment with Theaflavin, and certified the potential mechanism of Theaflavin, which can be used as a potential agent for cutaneous wound therapy.

INTRODUCTION
The skin can block the body from direct external stimulation, so it is regarded as the most susceptible damaged tissue (Sen et al., 2016). Skin wounds can be induced by mechanical, chemical, electrical and thermal wounds or persistent illnesses. Skin wounds pose a widespread global threat to the quality of life and socio-economic conditions (Boer et al., 2016;Gallo, 2017). These wounds not only cause pain, but also become infectious, requiring amputations among a large group of patients worldwide. Therefore, it is of great clinical significance to take novel strategies to accelerate the healing process of wound.
Wound healing is a complicated procedure that encompasses the collaboration between keratinocytes, fibroblasts, endothelial cells, macrophages, and platelets (Broughton et al., 2006;Velnar et al., 2009). Angiogenesis is thought to play a substantial role during the healing of wounds, which is primarily done by endothelial cells (Tonnesen et al., 2000;Veith et al., 2019). During wound healing, the microvasculature axis containing angiogenic capillaries can provide regenerated tissues with nutrients and oxygen (Davidson and Mustoe, 2001;Tejada et al., 2019). But when oxygen is used to generate energy through oxidative phosphorylation, reactive oxygen species (ROS) are synthesized, thereby increasing oxidative stress at the wound site and causing additional endothelial damage (Liu et al., 2019). Hence, it is probably a potential therapeutic target for protecting endothelial cells from damage and improving wound healing.
Heme oxygenase (HO)1 is activated under the condition of oxidative stress. It is a key molecule for the abrogation of oxidative stress (Zhang H. et al., 2018;Xiao et al., 2018). In addition, in the absence of HO1, endothelial cells become highly susceptible to stressrelated injury. Nevertheless, the stress-related injury can be alleviated by activating HO1 (Kinderlerer et al., 2009;Böckmann and Hinz, 2020). The upregulation of HO1 is related to the transcriptional agent nuclear factor-erythroid 2 associated factor 2 (Nrf2) (Sudan et al., 2019;Dong et al., 2020). The activated Nrf2 becomes nuclear-bound where it binds to the antioxidant response element, so as to promote balance in the ratio of oxidants/antioxidants after oxidative stress damage (Ma, 2013;Shaw and Chattopadhyay, 2020). Therefore, we hypothesized that stimulation of the Nrf2/HO-1 signaling axis was probably a protective treatment for endothelial cells against oxidative stress-driven destruction.
In recent years, a variety of natural products possessing antioxidant effects have been attracting significant attention as therapeutic agents for the prevention of oxidative stress-related diseases Xu et al., 2021). Among them, theaflavin (TF), a principal constituent of black tea, owns multiple health benefits and is well known for its antioxidant properties owing to its effect on superoxide anion scavenging (Ilacqua et al., 2017;Li et al., 2021). A recent study has clarified that TF ameliorates ionizing radiation-induced hematopoietic stem cell damage mainly by decreasing oxidative stress through activating the Nrf2/HO-1 pathway (Han et al., 2017). It has also been found that TF could attenuate cerebral ischemia/reperfusion injury by abolishing Nrf2 inhibition and reducing oxidative stress (Li et al., 2019). However, their role in wound healing therapy remains unknown.
Previous research has reported that growth factors such as bFGF, PDGF, and VEGF could promote angiogenesis, and several signal proteins including phosphatidylinositol-3-kinase (PI3K), protein kinase B (AkT) and mTOR participate in the wound healing process (Gupta and Qin, 2003;Muñoz-Chápuli et al., 2004;Zhang R. et al., 2018). Interestingly, a research study confirmed that activation of PI3K and AKT can promote the cells proliferation, migration, differentiation, angiogenesis . However, the potential angio-modulatory and wound healing roles of TF remain poorly understood.
Consequently, in this study, we assumed that TF could promote angiogenesis and accelerate wound healing through activating the PI3K/AKT/Nrf2 signaling pathway. The proangiogenic effect and the underlying mechanism of TF were investigated in vitro.

Cell Culture and Treatment Protocols
HUVECs (ATCC, Manassas, VA, United States) were grown in DMEM/F12 ((Gibco, Invitrogen, Grand Island, NY)) containing 10% thermal-inactivated FBS and 1% penicillin and streptomycin in an incubator at 37°C with 5% CO 2 . Briefly, fresh cells were seeded at a density of 2,500-3,000 cells/cm 2 in T-75 flasks (Falcon). The medium was changed every 48 h. By light microscopic examination we observed that cultures reached confluence after 6-7 days.
According to previous group experiments (Carracedo et al., 2012), cells passaged 2-9 times were defined as young, cells passaged 15-25 times were defined as intermediate and those passaged >30 times were considered as senescent.
Each huvec type (young, intermediate or senescent) was stimulated for 24 h with different concentrations of Theaflavin (0, 25 and 50 µM) with or without TBHP (500 µM) to determine the apoptosis, oxidative stress and cell proliferation. All assays were quantified and analyzed by flow cytometry on a FACSCalibur cytometer (Becton Dickinson Biosciences (BD); San Jose, CA, United States) equipped with standard CellQuest software. Data were collected for 10,000 cells per sample.

Cell Viability Assay
Cell Counting Kit-8 (CCK-8) assay (MedChemExpress LLC; Monmouth Junction, NJ, United States) was used to determine cell viability. TBHP can stably cause oxidative stress in endothelial cells, thus it was used as an in-vitro stimulus to simulate the oxidative stress process in the survival of wound healing (Jiang et al., 2021). Before being co-treated with Theaflavin (0, 25 and 50 µM) with or without TBHP (500 µM), second-generation HUVECs (1×10 4 cells/well) were grown on a 96-well plate and incubated in DMEM/F12 medium at 37°C. 24 h later, the cells were rinsed with PBS and then DMEM/F12 serum-free culture containing 10 μL CCK8 was added to each well to incubate for another 2 h. Eventually, the absorbance of cells at 450 nm in each well was recorded by a microwell Plate reader (Thermo Fisher).

Western Blotting Assay
Western blotting assay was performed as per the conventional protocol (Gu et al., 2021). HUVECs were lysed in 1 mM PMSF (phenylmethanesulfonyl fluoride) radioimmunoprecipitation analysis buffer to extract protein and then quantified with BCA protein assay kit (Beyotime). By following the manufacturer's guidelines, we used a commercial kit to isolate the nuclear and cytoplasmic protein fractions. Equal amount of protein was isolated by sodium dodecyl sulfate-polyacrylamide gel electrophoresis, transferred onto a polyvinylidene fluoride membrane (BioRad, United States), and then blocked with 5% skim milk at room temperature for 2 h, followed by exposure to the following primary antibodies at 4°C overnight: Cleaved-caspase3 (C-caspase3) (1:1000), Bax (1:1000), Bcl-2 (1:1000), Cytochrome C (1:1000), Nrf2 (1:1000), HO1 (1:1000), or Lamin B (1:800). After that, they were exposed to corresponding HRP-conjugated secondary antibodies for 2 h at room temperature. The washed bands were visualized by using a chemical XRS + imaging system (Biolard, United States) and Image LabV 3.0 (Bio-Rad, United States) was used to quantify protein expression.

Tunel Assay
In Situ Cell Death Detection Kit (Roche, SouthSanFrancisco, CA, United States) was used to measure apoptotic cells, according to the manufacturer's guidelines. To put simply, the treated HUVECs (1×10 4 cells/mL) were rinsed with PBS for three times and then fixed in paraformaldehyde (4%) for 20 min at room temperature. Later, the cells were exposed to 3% hydrogen peroxide and freshly-made 0.1% Triton X-100 for 10 min. For each step, the HUVECs were rinsed with PBS for three times and then exposed to TUNEL reagent. We quantified the apoptotic cells in 3 separate arbitrary fields under a fluorescence microscope (Olympus Inc., Tokyo, Japan).

EdU Staining Assay
According to the manufacturer's steps, huvec proliferation was assessed by using the intake of 5-acetyl-2′-deoxyuridine (EDU) in DNA with the Clickit EdU Microplate Test Kit (Invitrogen). The treated HUVECs (1 × 10 4 cells/mL) in each group were incubated with Oregon Green azide conjugated EdU. And then the cells were permeated and reacted with HRP-conjugated anti-Oregon green antibody and Amplex infrared. The stained HUVECs were observed under a fluorescence microscope (Olympus Inc., Tokyo, Japan).

Mitochondrial Function Determination
HUVECs (1×10 4 cells/mL) were treated with different concentration TF (0, 25 and 50 μM) for 24 h with or without 500 μM TBHP for 2 h. To detect the changes in mitochondrial membrane potential (MMP), the treated HUVECs were incubated in JC-1 solution (T3168, Invitrogen) at a concentration of 10 mg/L for 20 min at 37°C and then washed 3 times with PBS to remove excess JC-1 solution. Then, the stained samples images were immediately captured under a confocal microscope (Olympus, Tokyo, Japan). Lastly, the relative MMP was analysed with the ImageJ software. The relative MMP was calculated as a ratio of the mean fluorescence intensity of red fluorescence (excitation wavelength, 525 nm; emission wavelength, 590 nm) to green fluorescence (excitation wavelength, 490 nm; emission wavelength, 530 nm). To detect the changes in mitochondrial ROS generation, the treated HUVECs were stained in 5 nM MitoSOX Red Mitochondrial Superoxide Indicator solution (M36008, Invitrogen) for 30 min at 37°C. The fluorescence images were captured via a microscope (Olympus, Tokyo, Japan) under 510 nm excitation wavelength and 580 nm emission wavelength. Quantitation of mean fluorescence intensity by ImageJ was used to compare the ROS changes of the mitochondria.

Tube Formation Assay
Huvec tube formation was performed on a chamber glass slide coating by means of matrix coagulation. An ECMatrix gel solution was applied on the u-slide plate, and then placed it in a 37°C incubator for 1 h to solidify the matrix. The pretreated HUVECs (1 × 10 4 cells/mL) were harvested by using trypsin/ ethylenediaminetetraacetic acid. The treated HUVECs were plated on the Matrigel processed earlier, and then incubated at 37°C for 6 h. The tube formation was observed under a phase contrast microscope (40 × ) and then estimated by quantifying the selected area of each well randomly.

Cell Migration Assay
Huvec migration was measured on an 8-μm-pore polycarbonate membrane Boyden chamber insert within a transwell system (Costar, Cambridge, MA, United States). The cells were exposed to Theaflavin (TF) and TBHP as previously described. After that, the cells were separated, centrifuged, and resuspended. 1 × 10 4 cells were placed in 200 μL of non-FBS DMEM/F12 medium inside a transwell device, and then 700 μL of medium with 1% FBS was introduced to the bottom chamber. After being placed in a 5% CO 2 incubator for 12h, the membrane was rinsed with PBS for three times, and then fixed with 4% paraformaldehyde. The transwell apparatus were stained with crystal violet, and cells on the upper surface were removed with swabs. Those that moved to the lower surface were quantified in 3 arbitrary fields (40×).

Cell-Matrix Adhesion Assay
An huvec adhesion assay was performed on a 6-well plate. HUVECs were exposed to TBHP and Theaflavin (TF) as previously described. The plate was pre-coated with fibronectin (5 μg/ml) at 37°C for 1 h. Equal numbers of harvested cells (1 × 10 4 cells/mL) were grown on each coated plate and placed in the incubator for 30 min. Later on, the nonadherent cells were rinsed with PBS and then fixed with 4% paraformaldehyde. The adherent cells were detected by DAPI staining. The quantity of cell adhesion was determined in three separate representative areas of each well.

Small Interfering Ribonucleic Acid Incorporation
To silence human Nrf2 gene (RiboBio, Guangzhou, China), double-stranded siRNA was designed and synthesized, with the following sequence: sense strand 5′-GGTTGAGACTAC CATGGTT-3'. According to the manufacturer's guidelines, the cells (1×10 4 cells/mL) were mixed with 50 nM siRNA and Lipofectamine 2000 Reagent (Thermo Fisher, UT, United States) for 36 h, exposed to TF and TBHP as previously described, and then subjected to Western blot analysis.

Real-Time PCR
TRIzol (Invitrogen) was used to isolate the total celluar RNA from treated HUVECs according to the methods provided by the manufacturer. After that, 1 µg of total RNA in each group was reverse-transcribed with the cDNA synthesis kit (MBI Fermantas, Germany). The PrimeScript-RT reagent kit (TAKARA, Japan) and SYBR Premix Ex Taq (TAKARA) in a CFX96 Real-Time PCR System (Bio-Rad Laboratories, CA, United States) were used to analyze the quantitation of PCR. Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) expression was taken as an internal control by using the 2 −ΔΔCt method, to evaluate the expression of target gene. The primer sequences were designed in the laboratory and synthesized by TsingKe Biotech based on the mRNA sequences obtained from the NCBI database as follows: HO-1, 5′-AGAGTTTCTTCGCCAGAGG-3' (forward), and 5′-GAGTGTGAGGACCCATCG-3' (reverse).

Animals and Ethics Statement
Fifty-four healthy male Sprague-Dawley rats (250-300 g) were purchased from Wenzhou Medical University (license no. SYXK [ZJ] 2020-0014) and caged under routine environment (temperature: 23 ± 2°C, humidity: 50 ± 5%, 12 h light/dark cycle). All animal surgical interventions followed the Guide for Care and Use of Laboratory Animals of the China National Institutes of Health. The experiments involving animals received approval from the Wenzhou Medical University's Animal Research Committee (wydw 2021-0256). For this study, rats were arbitrarily separated into two groups: Control (n = 18), TF (n = 18).

Skin Wounding Model Establishment and Drug Administration
2% (w/v) pentobarbital (40 mg/kg) was administered intraperitoneally to anesthetize the rats, and the vital signs of the rats were observed closely during the re-anesthesia. After shaving and disinfecting, two circular wounds with a diameter of 2 cm were made on both sides of the dorsal torso using surgical scissors. To determine the appropriate oral dose of TF, thirty rats were randomly divided into the five groups after surgery (n = 6) and were received TF orally for seven consecutive days with doses of 0, 10, 20, 30, 40 mg/kg, respectively. After the determination of the appropriate TF dose, animals were randomly divided into two groups. The TF group was administered TF (20 mg/kg) intragastrically each day until the day of sacrifice. The control group was treated with same dose of normal saline. On day 7, day 14, and day 21 after wounding, the rats were euthanized with excessive pentobarbital sodium, and the wounds and adjacent tissues were sampled for histological evaluation.

Laser Doppler Blood Flow Imaging
For the purpose of visualizing circulation under the flap, the rats were anesthetized, made to lie prone and then measured with a laser doppler instrument (Moor Instruments, Axminster, United Kingdom) day 7, day 14, and day 21 after surgery. Moor LDI View (ver.6.2; Moor Instruments) was adopted to make a data analysis, and perfusion units (PU) were used to assess the intensity of blood flow. The scan was performed for three times and the mean of each animal was used for statistical analysis.

Wound Tensile Strength Measurement
The tensile strength (TS) testing device is based on a specially shaped horizontal arm pulling one side of a sample with the opposite side fixed to a measuring tip of a force meter unit (OMEGA Engineering, Inc., Stamford, CT, United States). The moving arm is driven by a high-precision stepper motor MDI-17 (Intelligent Motion Systems, Inc., Marlborough, CT, United States) through a linear slider.
The measurement technique was described previously (Gál et al., 2009). Briefly, two 1-cm-wide skin strips were removed from each incision and placed lengthwise between the clamps of the TS testing device. Pulling was performed perpendicularly to the original direction of the wound. The maximal breaking strength was recorded for each sample. The TS was then calculated by using the formula: TS = MRS/A (MBS = maximal rupture strength [g], A = wound area [mm 2 ]) and expressed in g/mm 2 .

Histological Analysis
In each group, six tissue samples were collected to perform Hematoxylin and eosin (H&E) staining. These specimens first received fixation in 4% (v/v) paraformaldehyde overnight and then were embedded in paraffin wax. After that, the specimens were cut into 5 µm thick slices for H&E staining, and then observed under an optical microscope to assess histological changes. From each section, the microvessels were measured on six fields that were selected randomly, and the number of microvessels per unit area (/mm 2 ) was counted to quantify the level of microcirculation.

Immunofluorescence
For the purpose of cell immunofluorescence, HUVECs were grown on a 6-well plate for 24 h, exposed to 0, 25 or 50 μM TF for 24 h, and then exposed to 500 μM TBHP for 2 h. After that, they were rinsed with PBS, fixed with 4% paraformaldehyde for 15 min, permeabilized with 0.1% TritonX-100 in PBS for 10 min, and blocked with 5% bovine serum albumin (BSA) for 1 h at RT. After that, they were exposed to primary antibodies against Nrf2 (1:200; Abcam) overnight at 4°C, followed by exposure to Alexa Fluor ® 488-goat or Alexa Fluor ® 594-goat anti-rabbit IgG (H + L) secondary antibodies (1:300; Jackson Immunoresearch, PA, United States) at 37°C for 1 h. Subsequently, they were stained with DAPI (Beyotime, China) for 5 min, and observed under a Nikon ECLIPSE Ti microscope (Nikon, Japan). For each slide, six arbitrary fields were analyzed. Later on, fluorescence quantification was done with Image-Pro Plus (Media Cybernetics, Rockville, MD, United States) and then analyzed by researchers who were ignorant of the experimental information.
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org March 2022 | Volume 10 | Article 830574 For tissue immunofluorescence, the methods were similar to those mentioned in the immunohistochemistry procedure. As mentioned above, the tissue samples were deparaffinized, rehydrated, and washed with PBS. Endogenous peroxidase was quenched by using 3% hydrogen peroxide and the tissue antigen was repaired with 10.2 mM sodium citrate buffer. The specimens were exposed to primary antibody against a-SMA (anti-alpha smooth muscle actin (αSMA), 1:200) overnight at 4°C. After that, tissue sections were treated with a Texas red-conjugated anti-IgG secondary antibody. The nuclei were stained with DAPI. Observed under a fluorescence microscope (Olympus Inc., Tokyo, Japan).

Immunohistochemistry
For immunohistochemistry analyses, the dewaxed and hydrated sections were placed in boiling buffer (citrate buffer; pH 6.0) for 10 min. Each section was then incubated in one drop of 3% H 2 O 2 for 10 min and in one drop of primary monoclonal antibodies of CD31 (1:100) or VEGF (1:100) at 4°C overnight. The sections were washed with PBS, incubated with HRP goat anti-rabbit (1: 400; Abcam) for 1 h, washed again with PBS and incubated with a drop of a DAB chromogen kit (ZSGB-BIO, Beijing, China), then counterstained with hematoxylin (Beyotime Institute of Biotechnology, China) for 5 min. Sections were then counterstained, differentiated and blued. Finally, the sections were dehydrated, dried, sealed, aired and photographed using a Nikon ECLPSE 80i (Nikon, Japan), and analyzed by Image-Pro Plus 6.0 software.

Molecular Modeling
We next assessed whether there was any affinity between TF and Nrf2 or upstream proteins in the Nrf2 pathway via a computational molecular docking analysis (Luo et al., 2020). For this analysis, we utilized the TF chemical structure shown FIGURE 1 | Effects of TF on cell viability and proliferation. (A) TF chemical structure. (B) Different concentrations of TF treated with HUVECs for 24 h using CCK-8 assay to assess viability. (C) HUVECs were exposed to different concentrations of TBHP for 24 h using CCK-8 assay to examine cellular survivability. (D) Representative images demonstrating cell morphology alterations in TF (0, 25, 50 μM) with or without TBHP (500 μM) co-treatment in HUVECs (scan bar, 100 μm). (E) Evaluation of HUVECs viability, using CCK-8 assay, after TF (0, 25, 50 μM) with or without TBHP (500 μM) co-treated in HUVECs. (F) Representative images show that EdU staining assay results in TF (0, 25, 50 μM) with or without TBHP (500 μM) co-treatment in HUVECs. Scale bar: 50 μm. (G) Histogram plots show that percentage of EdU positive cells in each group. Data presented as mean ± SD, *p <0.05, **p < 0.01, versus the indicated group, n = 3.
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org March 2022 | Volume 10 | Article 830574 in Figure 1A. After examination of all generated models, we found TF to clearly interact with and dock in the Nrf2 binding site ( Figure 4A), with macro-and local-level views of these interactions shown using a ribbon model. We additionally utilized a space-filling model to illustrate between TF and the Asp422, Val420, Val608, and Val561 residues of Nrf2. These results suggested that TF may function to inhibit the development of wound healing in part through its ability to interact with Nrf2 in a manner that promotes its nuclear translocation.

Statistical Analysis
Data are represented by means ± standard deviation. Statistical analyses were conducted via GraphPad Prism (United States) (one-way analysis of variance (ANOVA) and Tukey's post hoc test. Comparisons of mean values between two groups were performed using the independent-sample t test. p values < 0.05 was regarded as significant difference.

Effects of TF on Cell Viability and Proliferation
To estimate whether TF exerts cytotoxic effects on HUVECs, we treated HUVECs with different doses of TF. No damaged cells were found in HUVECs when the dose of TF was below 50 μM. According to our data, 50 μM TF was selected as the maximum non-toxic concentration for subsequent experiments (p = 0.03, Figure 1B). Then, after TBHP was administered, the cell viability decreased, and 63.0 ± 4.3% cell viability was observed at the concentration of 500 µM (p = 0.008, Figure 1C). In subsequent experiments, compared with the TBHP-treated group, TF protected the cells from TBHP-induced damage in a dosedependent manner (61.6 ± 2.8%, 70.0 ± 5.7%, 84.3 ± 5.2%, respectively, at the doses of 12.5, 25, and 50 µM, indicating that TF produced a cytoprotective effect on TBHP-mediated cytotoxicity ( Figures 1D,E). Futhermore, 5-ethynyl-2′deoxyuridine (EdU) fluorescence assays was used to detect TF induces proliferation in HUVECs under different conditions. The results indicated that HUVECs proliferation was dramatically reduced by TBHP but reversed by TF treatment in a dosedependent manner (62.0 ± 3.8%, 22.6 ± 4.8%, 37.2 ± 3.8%, 58.8 ± 5.2%,respectively) ( Figures 1F,G).

TF Ameliorated Mitochondrial Functional Damage and Apoptosis in TBHP-Induced HUVECs
Oxidative stress played a significant role in wound healing during the ischemia-reperfusion process. The high levels of oxidative stress and mitochondrial dysfunction can severely influence the functional state of HUVECs. For this reason, we evaluated the effects of TF in TBHP-treated huvec on mitochondrial dysfunction and antioxidant enzyme activity. The levels of mitochondrial ROS production were tested by MitoSox staining. According to Mitosox Red staining, compared with the untreated group, the mtROS of TBHP treatment group was 168.0 ± 7.1% (p = 0.0012). We also discovered that the TBHP-stimulated cells treated with TF had a lower MMP than that of the TBHP treatment group (p = 0.0021). Compared with the untreated cells, TBHP treatment facilitated the production of superoxide anion. These effects were restored to near physiological levels by TF treatment (Figures 2A,B). Futhermore, TUNEL assay was used to evaluate the TFinduced cytoprotection of TBHP-induced HUVECs. The expression of TUNEL + cells was upregulated after TBHP treatment and alleviated depending on the dose of TF (9.9 ± 3.2%, 75.3 ± 5.9%, 22.0 ± 4.8%, respectively, p = 0.003, p = 0.002, Figures 2D,E). According to the western blot data, the results of the expressions of apoptosis-related proteins showed that TBHPinduced mitochondria dysfunction in cells was restored by TF treatment in a dose-dependent manner ( Figures 2F-J).

TF Promoted Cell Function in TBHP-Treated HUVECs
To determine whether TF had therapeutic benefits on cell functions, we did several cell function experiments on the TBHP-treated HUVECs. Firstly, a transwell migration assay was performed to examine potential huvec migration mediated by TF. As shown in Figure 3, there were typical representations of hematoxylin-stained cells that migrated to the basolateral membrane. The statistical results showed that TBHP treatment significantly lowered the number of migratory cells, while TF reversed this effect (p = 0.001, p = 0.003, p = 0.006, respectively, Figures 3A,C). Next, to evaluate the role of TF on huvec adhesion, a fibronectin adhesion assay was conducted. The statistical analysis indicated that TBHP stimulation significantly reduced the number of adherent cells, but this effect was reversed by TF treatment (p = 0.002, p = 0.005, p = 0.007, respectively, Figures 3B,D). Lastly, a tube formation assay was used to examine huvec neovascularization after TF treatment, and then the number of capillary-like structures within each group was observed. Based on our data, TBHP markedly reduced huvec neovascularization, while TF pretreatment can protect from the TBHP treatment in a dose-proportional manner (p = 0.012, p = 0.004, p = 0.006, respectively, Figures 3C,F).

TF-Mediated Cytoprotection via Promoting Nucleus Translocation of Nrf2 and Activating the PI3K/AKT/Nrf2 Signaling Pathway
To determine whether TF and Nrf2 or other proteins in the Nrf2 axis exists any affinity, a artificial intelligence software was utilized to simulate molecular docking (Luo et al., 2020). We applied the TF chemical structure ( Figure 1A) for further analysis. Upon establishment of all possible models, TF does exist clearly intract and dock with the Nrf2 docking site ( Figure 4A). We adapted the ribbon model to display the macro-and local-level views of these interactions. Moreover, a space-filling model was performed to demonstrate this interaction. Finally, we found that there was a high-affinity (−9.8 kcal/mol) hydrogen binding events between TF and Nrf2 Frontiers in Bioengineering and Biotechnology | www.frontiersin.org March 2022 | Volume 10 | Article 830574 residues, incluing Asp422, Val420, Val608, and Val561. Results of software simulation indicated that TF may influence the wound healing process through its ability to interact with Nrf2, thereby promoting its nuclear translocation. Since Nrf2 is an important transcription factor in the redox reaction system (Kanlaya et al., 2021), we explored the potential role of Nrf2 played in TFinduced the cellular protective effect in HUVECs exposed to TBHP-induced oxidative stress. In order to determine whether TF promoted nuclear translocation of Nrf2 in the HUVECs, we carried out immunofluorescence staining on Nrf2, which suggested that TF induced the nuclear translocation of Nrf2 in the HUVECs ( Figure 4B). Furthermore, the western blotting further confirmed that compared with the control, there were higher levels of Nrf2 in the nuclei of TF-treated HUVECs (p = 0.037, p = 0.004, respectively, Figures 4E,F). Later, we analyzed the expression level of HO1, a target protein of Nrf2. According to the expression of HO1 mRNA and western blotting, compared with the control group, TF-treated HUVECs significantly facilitated the expression of HO1 (Figures 4C,D,G). In the end, we examined whether TF resulted in different expressions of Nrf2 upstream molecules, including Akt and PI3K. In keeping with previous results, TF significantly promoted phosphorylation of Akt and PI3K ( Figures 4H,I). These results implied that TF increased the nuclear transcription of Nrf2 in HUVECs by activating AkT and PI3K.

The Cytoprotective Effect of TF is Eliminated by the Nrf2 Knockdown
To explore whether TF induced cytoprotection by activating PI3K/AKT/Nrf2 signaling pathway, we silenced Nrf2 by using siRNAs and determined its effect in TF-treated HUVECs. As evidenced by the results of western blotting, the levels of Nrf2 and HO1 proteins in the control group were higher than those in si-Nrf2-transfected HUVECs (Figures 5A,B). Furthermore, as expected, Nrf2 siRNA reversed TF-induced alteration in apoptosis-related proteins. Compared with the control group, the western expression levels of C-caspase3, Bax, Bcl-2 and Cytochrome C significantly decreased in the Nrf2-silenced TFinduced HUVECs (Figures 5C,D). The TUNEL assay also confirmed that Nrf2 siRNA abrogated the TF-mediated protection of HUVECs ( Figures 5E,F). Futhermore, MitoSox Red staining were significantly increased in the Nrf2-silenced TFtreated HUVECs compared to controls ( Figures 5G,H). Taken together, these results showed that the protective effect of TF can be eliminated by Nrf2 knockdown.

TF Promotes Wound Healing Process in Rats
According to our in vitro experiments, TF played a role in promoting angiogenesis. We hypothesized that it can also contribute to angiogenesis and accelerate wound healing process. To this end, we established a rat model of fullthickness skin wound to assess the effect of TF on wound repair after surgery. The surgical process of wounding is shown in Figure 6A. First, thirty rats were administered different doses of TF to determine the appropriate dose for wound healing ( Figure 6C). As seen in Figure 6D, the wound closure rate increased with the increase of TF doses from 0 to 20 mg/kg, whereas it was no significance with administered dose of TF from 20 to 40 mg/kg. So 20 mg/kg was determined as the appropriate dose for wound healing in subsequent experiments. On day 7, the closure rate of wounds was 42.5 ± 2.3% in the TF treatment group and 21.3 ± 1.7% in the control group respectively. On day 14, the healing in the TF treatment group accelerated and the closure rate remained significantly higher than the control group. On day 21, the wounds in the TF treatment group were nearly closed, while part of the wounds in the control group remained unhealed (p = 0.001, p = 0.003, p = 0.042, respectively, Figures 6B,E), indicating that TF accelerated wound healing in vivo. The results of wound tensile strength measurements on day7 and day21 are shown in Figure 6F. On day 7 after wounding, the tensile strength (TS) of wound significantly increased in the TF group (TF = 17.2 ± 3.3 g/mm 2 ), but only 5.9 ± 2.5 g/mm 2 in the control group. On the other hand, on day21, compared with the control, the TS of TF-treated skin wounds was significantly higher (TF = 83.1 ± 12.2 g/mm 2 , control = 42.3 ± 8.4 g/mm 2 ).
Taken together, there was a significant improvement in TS after 3 weeks in the TF treatment group than in the untreated group. Wound closure rates in the control and TF (20 mg/kg) group at different times. Healing rates of full-thickness cutaneous wounds were significantly increased by the TF treatment. (C) Effects of different doses of TF on POD 7. (D) Wound closure rates in the different doses of TF. (E) Representative digital images of wound closure for control and TF (20 mg/kg) group on day 0, day 7, day 14, and day 21. Data represented as mean ± SD, *p <0.05 and **p < 0.01, n = 5. (F) Wound tensile strengths of wounds removed from untreated control, TF-treated (20 mg/kg) groups at day 7 and 21 post-surgery. The data are presented as the means ± SD of five independent experiments. *p < 0.05, **p < 0.01, versus the indicated group, n = 5.

TF Enhanced Angiogenesis in Wound Healing Process
To detect whether TF exerted a prominent effect on the microvascular network on the back of rats, we performed Laser Doppler to visualize the circulatory system. On day 7, the results suggested that compared with the control group, TF group had a markedly stronger signal intensity in the flaps (422.1 ± 8.64PU, 460.6 ± 16.3PU, respectively, p = 0.004, Figures 7A,C). On day 14, the control group was 533.3 ± 23.57PU and the TF group was 633.3 ± 12.47PU. Eventually, on day 21, the control group was 610.2 ± 16.32PU and the TF group was 683.3 ± 12.47PU. In addition, we carried out hematoxylin-eosin staining on day 7 and day 21, and the results indicated that compared with control rats, the vessel density in the TF rats was significantly improved (p = 0.0014, p = 0.003, respectively, Figure 7E). Compared with control rats, the immunofluorescenc staining of α-SMA in the TF rats also significantly enhanced, which coincided with the results of H&E staining (p = 0.005, Figure 7G). In immunofluorescence, a large number of Nrf2 positive cells were identified in the dermis of the TF group relative to the control group, as depicted in Figures 7H,I (p = 0.002). IHC was employed to detect the expression of CD31. Compared with the control group, the CD31 expression was considerably enhanced in the TF group, as shown in Figures 7J,K. Similarly, the quantity of VEGF-positive vessels in TF-treated animals were obviously higher than that of the controls (Figures 7L,M). Collectively, these data showed that TF faciliated angiogenesis and played a substantial role in improving wound healing.

Different Ligands Docked in the Nrf2 Binding Site
We also tested the affinity of other molecules known to bind Nrf2, as shown in the Table 1. This affinity is an important indicator to measure whether the ligand can effectively bind to the receptor molecule, and it is the core parameter of AutoDock Vina software. The lower the affinity value is, the closer the ligand binds to the receptor. From the Table 1, we found that the affinity of these known molecules was close to TF, indicating that TF had a high affinity with Nrf2.

DISCUSSION
In our study, we have firstly certified that Theaflavin (TF) exposure stimulated the antioxidant PI3K/AKT/Nrf2 signaling axis followed by attenuating cellular apoptosis and ameliorating the oxidative-induced celluar dysfunction on HUVECs. In our in vivo experiment, the gastric delivery of TF (20 mg kg −1 ·day −1 ) significantly accelerated the wound healing process, which was characterized by faster wound closure, increased new capillaries. As TF could enhance the formation of new capillaries at the wound site and promote wound healing in vivo, it can be considered as a promising drug for the treatment of ischemic wounds.
The healing process of wounds involves a series of intricate events. And the healing process of wound is mediated by plenty of individual growth factors, cytokines, and hormones (Barrientos et al., 2008;Behm et al., 2012;Ridiandries et al., 2018). In our study, the ability for HUVECs proliferation, adhesion, migration and tube formation was regulated by TF treatment, demonstrating that TF had proangiogenic effects. In general, the angiogenic process could be divided into a few stages, including endothelial cell proliferation, separation and migration, adherence to the ECM and differentiation (Lamalice et al., 2007). In our present study, TF could evidently boost HUVECs proliferation, which can be regarded as a signal of the beginning of angiogenesis. In addition, TF could regulate the adhesion of HUVECs, with decreased cell-cell adhesion and increased cell-matrix adhesion. Finally, endothelial cells differentiate into new capillary vessels, which are confirmed by the faster tube formation rate and tube numbers ( Figures 3C,F). Taken together, TF exerted effective in proangiogenesis in vitro.
In our study, we discovered that TF (0-50 μM) didn't produce a detrimental effect on cell viability. And TBHP was applied to induce an oxidative stress state in HUVECs. In conformity to other publications, we observed that at a concentration of 500 μM TBHP treatment significantly reduced the viability of HUVECs, but this effect was restored by TF pretreatment in a dosedependent manner. This implied that when TBHP was administered, TF produced a cytoprotective effect in HUVECs. What's more, we also found that TF pretreatment can downregulate the levels of apoptosis-related protein Bax, Cytochrome C and C-caspase3, and upregulate the level of the anti-apoptotic protein Bcl-2. In normal physiological environments, mitochondrial membrane potential equilibrium is strictly controlled. However, TBHP induces apoptosis by increasing the production of mtROS associated with the interruption of mitochondrial membrane potential, which facilitates the release of cytochrome C from mitochondria to cytoplasm by activating apoptotic protein Bax and inactivating the anti-apoptotic protein Bcl-2 (Ju et al., 2017). The results of this study suggested that TF protected HUVECs from TBHPinduced apoptosis by inhibiting the production of mtROS and inhibiting caspase, implying that TF was an effective antioxidant that can protect HUVECs from oxidative stress. Our results also demonstrated that TF-pretreatment was able to significantly improve cellular dysfunction induced by TBHP stimulation on HUVECs.
We also tested the affinity of other molecules known to bind Nrf2, as shown in the Table1. This affinity is an important indicator to measure whether the ligand can effectively bind to the receptor molecule, and it is the core parameter of AutoDock Vina software. The lower the affinity value is, the closer the ligand binds to the receptor. From the Table1, we found that the affinity of these known molecules was close to TF, indicating that TF had a high affinity with Nrf2.
HO1 offers cellular protection in the treatment of vascular disease (Durante, 2010). HO1 is ubiquitous in the heart and blood vessels, protects from vasculopathy and offers cytoprotection for endothelial cells (Khitan et al., 2014;Sandrim et al., 2019). Mice, whose HO-1 gene is knocked down, are very vulnerable to severe vascular diseases (Hull et al., 2013). Nrf2 is widely accepted as an essential transcriptional agent, and its stimulation can increase the levels of downstream antioxidant-related genes, including HO1 (Park et al., 2017;Huang et al., 2018;Liu et al., 2018). In a normal internal environment, Kelch-like ECH-associated protein-1 (Keap1) marks Nrf2 for ubiquitin-mediated proteasomal destruction (Suzuki and Yamamoto, 2015;Bellezza et al., 2018). In the presence of an external stimulus, Keap1 becomes inactivated and the ubiquitination-mediated removal of Nrf2 is weakened, leading to the upregulation of nuclear Nrf2 and combination with antioxidant response factors (Nezu et al., 2017;Guo and Mo, 2020). This process enhances the levels of antioxidants and second-stage detoxification enzymes like HO1. Oxidative stress can stimulate the Nrf2/HO1 axis, a key determinant of homeostasis maintenance in case of oxidative stress damage (Lv et al., 2019;Wu et al., 2019). However, after coadministration with TF, the HO1 levels and nuclear accumulation of Nrf2 were upregulated in TBHP-induced HUVECs, while the apoptosis and celluar dysfunction persisted. This was probably because TBHP treatment stimulated of the Nrf2/HO1 axis improperly, so as to circumvent cytotoxicity. Multiple studies have suggested that in specific cell types, TF alleviated the oxidative stress-driven injury through the Nrf2/HO1 axis (Cheng et al., 2013;Li et al., 2021). In fact, our data revealed that TF significantly upregulated the nuclear accumulation of Nrf2 and HO1 in TBHP-stimulated HUVECs. Nrf2 siRNA merger was performed to inhibit the expression of Nrf2, then decrease the levels of HO-1, and increase the levels of apoptosisassociated proteins. Based on these data, TF-mediated cytoprotection of TBHP-treated HUVECs took advantage of the stimulation of Nrf2/HO-1 axis.
PI3K/Akt is a classical signaling pathway that plays an important role in cell proliferation, migration, apoptosis, angiogenesis, and other biological processes. AKT is a serine/ threonine kinase and a downstream signal molecule of PI3K, activated by it (Li et al., 2017). Several studies have reported that the PI3K/Akt pathway is essential for regulating Nrf2/Ho-1 pathway activation and is thus involved in protection against oxidative stress and apoptosis in multiple cell types (Di Tu et al., 2020;Feng et al., 2021;Wu et al., 2021). Under oxidative stress, the PI3K/Akt signaling pathway leads to Nrf2-dependent transcription and overexpression of HO-1 protein (He et al., 2019). We found that phosphorylation of AKT protein increased in the presence of TF, which possibly upregulated TF-induced Nrf2/HO-1 protein. Under normal physiological conditions, Keap1 sequesters Nrf2 in the cytoplasm. After electrophilic agents or ROS oxidize cysteine residues within Keap1, Nrf2 translocates to the nucleus following dissociation from its cytoplasmic docking protein. In the nucleus, Nrf2 activates the transcription of several phase II detoxifying enzymes and antioxidant enzymes including HO-1 by binding to their promoter regions (Ryu et al., 2014). Thus, we speculate whether TF also promotes the angiogenesis by activating AKT phosphorylation. As we expected, the phosphorylated expression of AKT was enhanced and the expression of Nrf2 in the nucleus was activated under the effect of TF, indicating that AKT participated in the proangiogenic process of TF. In addition, a large number of studies have shown that Nrf2 is a redox transcription factor and a major participant in antioxidant response (Li et al., 2009). PI3K can also promote neuronal survival by activating Akt phosphorylation and Nrf2 nuclear translocation as reported (de Oliveira et al., 2016). In our study, the p-PI3K and p-AKT were increased after TF treatment, confirmed by western blot results, in which demonstrating that TF could promote angiogenesis via the PI3K/AKT/Nrf2 signaling pathway.
In an effort to evaluate the effect of TF in vivo, we established full-thickness cutaneous wound model, to assess the healing process of wounds on day0, day7, day14 and day21 after surgery. Our results indicated that compared with the control group, the TF treatment quickened wound healing. On day 21, the wounds in the TF treatment group were completely closed, while some wounds in the control group remained unhealed. Angiogenesis was critical for wound healing. To evaluate TF's activity on angiogenesis, IHC staining of CD31-positive vascular cells were carried out, and it was found that TF had increased the density of blood vessels in the wound tissue. Then whether TF modulated VEGF was investigated, which were found to have enhanced the angiogenesis. VEGF contributed to multiple processes of angiogenesis (particularly mitosis of vascular cells) promoting the formation and maturation of neovascularization. In brief, our work confirmed that TF probably served as an agent for wound healing.

CONCLUSION
In summary, according to the results of in vitro and in vivo experiments, we found that TF is possible target for wound healing therapy. In HUVECs, TF-mediated cytoprotection will attenuate oxidative stress-driven apoptosis and ameliorate cellular function. And the stimulation of PI3K/AKT/Nrf2 signaling axis is probably of great significance for the improvement of the healing process of wounds treated with TF.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

ETHICS STATEMENT
The animal study was reviewed and approved by Wenzhou Medical University's Animal Research Committee (wydw 2021-0256). Written informed consent was obtained from the individual(s), and minor(s)' legal guardian/next of kin, for the publication of any potentially identifiable images or data included in this article.

AUTHOR CONTRIBUTIONS
G-NH designed the study; DC, JJ performed experiments; DC wrote the manuscript; DC, ZW, L-NW prepared figures and statistic alanalysis; G-NH revised the manuscript. All the authors had read and approved the submitted the manuscript version.