Longitudinal host transcriptional responses to SARS-CoV-2 infection in adults with extremely high viral load

Current understanding of viral dynamics of SARS-CoV-2 and host responses driving the pathogenic mechanisms in COVID-19 is rapidly evolving. Here, we conducted a longitudinal study to investigate gene expression patterns during acute SARS-CoV-2 illness. Cases included SARS-CoV-2 infected individuals with extremely high viral loads early in their illness, individuals having low SARS-CoV-2 viral loads early in their infection, and individuals testing negative for SARS-CoV-2. We could identify widespread transcriptional host responses to SARS-CoV-2 infection that were initially most strongly manifested in patients with extremely high initial viral loads, then attenuating within the patient over time as viral loads decreased. Genes correlated with SARS-CoV-2 viral load over time were similarly differentially expressed across independent datasets of SARS-CoV-2 infected lung and upper airway cells, from both in vitro systems and patient samples. We also generated expression data on the human nose organoid model during SARS-CoV-2 infection. The human nose organoid-generated host transcriptional response captured many aspects of responses observed in the above patient samples, while suggesting the existence of distinct host responses to SARS-CoV-2 depending on the cellular context, involving both epithelial and cellular immune responses. Our findings provide a catalog of SARS-CoV-2 host response genes changing over time.


Introduction
Severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) is the etiologic agent of the coronavirus disease 2019 (COVID- 19) pandemic. The clinical spectrum of COVID-19, caused by SARS-CoV-2, is wide, ranging from asymptomatic infection to fatal disease. Risk factors for severe illness and death include age, sex, smoking, and comorbidities, such as obesity, hypertension, diabetes, and cardiovascular disease. Studies suggested that SARS-CoV-2 viral load can predict the likelihood of disease spread and severity [1][2][3] . A higher detectable SARS-CoV-2 plasma viral load was associated with worse respiratory disease severity 4 . Conversely, robust immune responses putatively mediate non-severe illness, in part, by controlling the replication of SARS-CoV-2 5,6 . Emerging evidence indicates that age and sex differences in the innate and adaptive immune response can explain the higher risks observed in older adults and male cases 7,8 .
Initial site of SARS-CoV-2 replication is the upper respiratory tract, and replication usually peaks within the rst week of infection 6 . The amount of virus produced at the respiratory epithelium is considered to be a critical element in determining SARS-CoV-2 transmissibility, duration of illness or severity, although it is not the only factor 9,10 . Higher viral loads have been observed in hospitalized patients with severe disease, have been attributed to high transmission and superspreading events, and have resulted in prolonged viral RNA shedding 1,[11][12][13][14][15] . Speci c anatomic site or host cell type where viral replication occurs, can also determine the course of infection. For example, angiotensin-converting enzyme 2 (ACE-2) and transmembrane serine protease 2 (TMPRSS2) receptors expression is highest in the upper respiratory tract and decreases in the distal or lower respiratory tract, incidentally SARS-CoV-2 infection mirrored this pattern, with high replication in proximal (nasal) versus distal pulmonary (alveolar) epithelial cells 16 . Control of viral replication and resolution of the in ammatory response is believed to be dependent, in part, on viral load and route of infection as well as the host immune response 17 . The early host immune response is regulated closely by the epithelial cell cytokine signaling in response to active viral replication 18 . Rapid and robust activation of the antiviral innate immune response at the site of viral replication is required to control and clear the virus. A delayed cytokine response can result in prolonged viral replication and worst clinical outcome as seen for other respiratory viruses 19 Our understanding of the viral dynamics of SARS-CoV-2 and host responses driving the pathogenic mechanisms in COVID-19 is evolving rapidly. Multiple studies have reported various characteristics of immune/in ammatory responses to SARS-CoV-2. Cytokine or chemokines-related host in ammatory responses such as CCL2/MCP-1, CXCL10/IP-10, CCL3/MIP-1A, and CCL4/MIP-1B were detected in bronchoalveolar lavage samples of SARS-CoV-2 infected adults while activation of apoptosis and the P53 signaling pathway were observed in lymphocytes 20 . In ammatory cytokine such as IL-1, IL-18, and IL-33 were enriched in the airways of COVID-19 patients 21 . In addition, a shotgun host transcriptomic analysis on nasopharyngeal samples revealed a wide range of antiviral responses. These included gamma and alpha interferon responses, elevated levels of ACE-2, interferon stimulated genes (ISGs), and interferon inducible (IFI) genes 22 . Very few studies have demonstrated the temporal correlation between viral load and host gene expression. Variation in viral load was associated with the SARS-CoV-2 disease and the host response dynamics via innate and adaptive immunity (To et al., 2020). Another study revealed that expression of interferon-responsive genes, including ACE-2, increased as a function of viral load, while transcripts for B cell-speci c proteins and neutrophil chemokines were elevated in patients with lower viral load 23 . Rouchka et.al. reported that cellular antiviral responses strongly correlated with viral loads. However, COVID-19 patients who experienced mild symptoms had a higher viral load than those with severe complications 6 .
We previously reported on a small group of adults with extremely high SARS-CoV-2 viral load, who had the potential to be super spreaders and a large group of adults with low SARS-CoV-2 viral load, both groups had mild illness 14 . Here, we wanted to determine the host response in relation to the viral load early during infection. We conducted a longitudinal study to investigate gene expression patterns detected in the secretion of the nasal epithelium during the acute phase of SARS-CoV-2 infection. The cases included SARS-CoV-2 infected individuals with an extremely high viral load early in their illness matched to individuals who either had a low SARS-CoV-2 viral load early in their infection or were otherwise stable patients who tested negative for SARS-CoV-2 prior to their outpatient surgical or aerosol generating procedure. We also determined the transcriptional response of a human nose organoid (HNO) line infected with SARS-CoV-2 and compared it to transcriptomic pro les generated from the upper respiratory tract secretion collected by nasal swabs from SARS-CoV-2 infected individuals.

Study cohort
Ten SARS-CoV-2 cases were randomly selected from our population of adults with extremely high viral load (Ct <16 to N1 target) at the time of their rst RT-PCR positive test. For each high viral load case, two additional human subjects were matched based on gender, week of rst SARS-CoV-2 RT-PCR test, age, and home zip code. These additional subjects consisted of either 1) SARS-CoV-2 infected adults with low viral load (Ct 31-<40) (SARS-CoV-2 low viral load case) or 2) stable adults who were SARS-CoV-2 RT-PCR negative (SARS-CoV-2 negative control) for their out-patient surgical or aerosol generating procedure. Each high viral load case had two to three subsequent SARS-CoV-2 RT-PCR positive mid-turbinate (MT) swab samples collected over a 4-week period. Each SARS-CoV-2 low viral load case had similarly spaced SARS-CoV-2 positive MT swab samples matched to its respective extremely high viral load case. On the other hand, the SARS-CoV-2 negative control only had one MT-swab sample collected with no longitudinal follow-up and was used to establish the transcriptomic baseline in the respiratory epithelium during the time the extremely high and low viral load matched cases were identi ed. The demographic and visit characteristics for the cohort is presented in Table 1. In general, age, gender, race, ethnicity, and zip code were comparable between the extremely high viral load, low viral load, and SARS-CoV-2 negative adults. The adults in the SARS-CoV-2 negative group were mostly asymptomatic at the time of testing, although their demographic information was not signi cantly different compared to that of both the SARS-CoV-2 extremely high and low viral load groups. The median Ct value difference between the extremely high and low viral load groups were 794,672-fold (19.6 Ct difference) and 724-fold (9.5 Ct difference) different at Visit 1 and Visit 2, respectively. At Visit 3, approximately 14 to 17 days after their Visit 1, the Ct values were comparable between the two groups.
RNA sequencing of serially collected specimens.
Of the 73 MT swab samples from the extremely high and low viral load SARS-CoV-2 groups with longitudinal follow-up and SARS-CoV-2 negative controls, only 44 (60.3%) MT swab samples from 20 (66.7%) individuals were of good quality to generate RNA-sequence data to study the host response to SARS-CoV-2 infection over time (Table 2). Demographic factors such as age, gender, race, ethnicity, zip code, disease severity and co-morbid conditions were comparable between the extremely high viral load, low viral load groups, and SARS-CoV-2 negative control group. Host response data were available on eight cases (extremely high viral load) with 23 samples. Six of the 8 extremely high viral load cases had gene expression data for Visits 1, 2, and 3, and two others for Visit 1 and 3. On the other hand, eight low viral load cases had 17 samples with gene expression data. Only two of the low viral load cases had gene expression data for Visit 1, 2, and 3. Another two low viral load cases had gene expression data at Visit 1, Visit 2 or 3, and Visit 4. The remaining four low viral load cases had gene expression data at Visit 1 only (n=1), Visit 2 only (n=2) or Visit 1 and 2 (n=1). Only 4 of the 10 SARS-CoV-2 negative control adults had gene expression data. All together 44 MT swab samples were sequenced for RNA to observe gene expression changes in the host response of the cases with extremely high viral load over time, as compared to the SARS-CoV-2 low viral load matched cases and the negative controls.

Gene expression changes by viral load
From our RNA-seq dataset, we could identify widespread gene expression changes from the nasal epithelium attributable to transcriptional host responses to SARS-CoV-2 infection. By comparing the expression levels of each gene with the sample viral load (representing the inverse correlation with Ct value) across the 44 MT swab samples, 425 genes were statistically correlated at p<0.01 signi cance level and 112 genes at p<0.001 (Figure 1a, Pearson's correlation). A stricter statistical cutoff would involve fewer expected false positive genes from multiple testing. However, the above 425 genes with p<0.01 would still be highly enriched for true positives, as revealed by integrating these genes with information from external databases, as described below. We also compared the expression levels of genes at individual time points during infection of both the extremely high viral load and low viral load groups with the SARS-CoV-2 negative control group (Figure 1b). Comparing Visit 1 MT swab samples from the extremely high viral load cases (n=8 samples from eight subjects) with the MT swab samples in the SARS-CoV-2 negative control group (n=4) yielded the highest number of genes with statistically signi cant correlated expression, as opposed to comparisons involving later times for the extremely high viral load group or involving the low viral load group. The gene expression from the extremely high viral load cases at Visit 1 highly overlapped with the differentially expressed genes of the low viral load group at Visit 1 ( Figure 1c) and remained highly correlated throughout their last visit. Interestingly, genes from the extremely high viral load group that did not overlap with the low viral load group did not show signi cant overlap with information from external databases.
To further delineate the differences in host gene expression between extremely high and low SARS-CoV-2 viral load groups, we performed an upset plot analysis to identify unique and common intersecting genes between the samples ( Figure S1). Among all the differentially expressed genes (DEGs) in the samples, 614 DEGs were unique to the subjects in the extremely high viral load group at visit 1 ( rst visit) and 226 genes were unique for the extremely high viral load at the last visit. The low viral load subjects on the rst and last visit showed 157 and 93 unique DEGs respectively. There were 31 DEGS that were common between all the groups. We performed the Gene ontology (GO) analysis of the unique and overlapping DEGs sets, and we found signi cant enrichment (FDR <0.05, count =3) of the biological processes including defense response to virus, negative regulation of viral genome replication, innate immune response, response to virus ( Figure S2) that were uniquely expressed in the extremely high viral load group at visit 1. SARS-CoV-2 infection in the low viral load group at either the early or later phase of the infection and the extremely high viral load group at the last visit did not show statistically signi cant enrichment of GO biological process. These ndings indicate that subjects with extremely high viral load at their initial visit were responding to the infection with increased immune responses, and thus preventing prolonged viral infection with a poor prognosis.
Differentially expressed gene in respiratory samples from extremely high viral load adults Focusing on the 112 top gene expression correlates of viral load across the 44 MT swab samples (p<0.001, Pearson's), 108 of these genes were higher in the SARS-CoV-2 infected adults with extremely high viral load. When visualizing the differential expression patterns of these 108 genes by heat map (Figure 2a), the genes were highest at Visit 1 of the extremely high viral load group, then decreased in expression with subsequent time points, tracking with the decrease in viral load (i.e., increase in Ct value). The 108 genes showed intermediate relative expression levels in the low viral load group and low expression in the SARS-CoV-2 negative control group. The 367 genes increased with extremely high viral load at p<0.01 were highly enriched for functional gene categories, as de ned by GO annotation terms.
Enriched GO terms (Figure 2b, p<=3E-5, one-sided Fisher's exact test) included 'immune system process', 'response to virus', 'type I interferon signaling pathway', 'cytokine-mediated signaling pathway', 'response to stress', 'regulation of viral life cycle', 'immune response', 'response to cytokine', 'innate immune response', 'response to interferon-gamma', 'regulation of I-kappaB kinase/NF-kappaB signaling', 'JAK-STAT cascade', 'protein ubiquitination', 'regulation of cell death', 'T cell activation', 'vesicle-mediated transport', and 'complement activation'. Of the 17 functional gene categories, there were ve gene categories -'response to virus', 'type 1 interferon signaling', 'regulation of viral life cycle', 'response to interferon-gamma', and 'JAK-STAT cascade' -where approximately 20% or higher of the genes were over expressed for that pathway. Overall, the above gene categories were highly indicative as representing a host immune response to an acute viral infection. Some of the genes that were upregulated were EIF2AK2 (eukaryotic translation initiation factor 2 alpha kinase 2), and ZC3HAV1 (zinc nger CCCH-type containing, antiviral 1), which have anti-viral activity. Other genes like IFIT2 and IFIT3 (interferon induced protein with tetratricopeptide repeats) aid in apoptosis. Chemokine genes like CXCL9 and CXCL10 that are involved in T-cell tra cking were also highly expressed. In contrast, the 62 genes decreased with high viral load at p<0.01 were not highly enriched for GO terms. Some of the genes that were downregulated included OR4A16 and OR10X1, involved in olfactory responses; SALL3 and MAGB6, which aid in downregulation of transcription; and TUBA3E, MLN, and ISTN1, affecting tubulin functions.  (Figure 3c). This in part may re ect the potential differences in the respective SARS-CoV-2 variants causing the infection or differences in the illness severity of the host.
Comparison of our respiratory sample gene sets to the transcriptional response of the human nose organoid infected with SARS-CoV-2 As another means to identify host transcriptional responses to SARS-CoV-2 infection, we generated RNAseq data on the human nose organoid model HNO 26 . We sampled HNO cells infected with SARS-CoV-2 and mock control cells at 6hrs, 72hrs, and 6 days post-infection, and we pro led these samples for gene expression. In the HNO204 RNA-seq dataset, 1760 genes were statistically signi cant at p<0.05 signi cance level and 341 genes, at p<0.01, exceeding chance expected. The top 867 genes over- In our study samples, the numbers of genes that were upregulated were much higher compared to down regulated genes (367 vs 62). Some of genes that were downregulated included those which operate olfactory functions (OR4A16 and OR10X1), downregulation of transcription (SALL3 and MAGB6), and tubulin functions (TUBA3E and MLN, and ISTN1). Previous studies have reported larger numbers of down regulated host response genes especially involving olfactory receptor pathway, neutrophil degranulation, and vesicle formation-indicating the role of these genes in loss of olfactory function in SARS-CoV-2 infections as well as the viral control of host-cell machinery 20,22,23 . One other study also showed very low number of downregulated genes with SARS-CoV-2 infection 29 , one reason for the low number of downregulated genes observed in our longitudinal study could conceivably relate to the mild illness experienced by both the extremely high and low viral load groups, in addition to the timing of sample collection as compared to other studies as well as the SARS-CoV-2 variants respectively involved.
Remarkably, the highest number of signi cant expressed genes were driven by the extremely high viral load group at Visit 1 ( rst visit). Also, all the genes that were upregulated with the low viral load group at Visit 1 completely overlapped with the extremely high viral load group at Visit 1 except for one gene, -CNN2, which plays a role in cell adhesion and muscle contraction. The predominant sets of genes involved in defense response to virus, type I interferon signaling pathway, cytokine-mediated signaling pathway-such as CXCL10, TGFB, IFIT2, IFIT3, OAS1, and IRF1-were not found signi cantly upregulated in the low viral load group. Consonant with this, Rouchka et. al. also observed that subjects with high viral loads had robust interferon and cellular anti-viral response and even exhibited strong inverse correlation with disease severity 6 . We previously noted that some SARS-CoV-2 infected adults with low viral load experienced prolonged viral shedding and low uctuation in viral load over time 14 . Absence or low expression of the anti-viral response in the low viral load group strengthens our observation of prolonged shedding in adults with a low viral load early in infection.
In our longitudinal study, the up-regulated host response genes that correlated with SARS-CoV-2 viral load over time in the respiratory secretion collected by the MT swabs were similarly differentially expressed across independent data sets of SARS-CoV-2 infected lung and upper airway cells 24 . About 170 of the differentially expressed genes observed in our study overlapped with SARS-COV-2 infected Calu-3 lung adenocarcinoma cell line but not with A549 cells. The observed difference across the cell lines could possibly be attributed to A549 cells not supporting robust replication of SARS-CoV-2 due to the low expression of ACE-2 30 . Similarly, 207 up-regulated genes from our longitudinal study overlapped with nasopharyngeal swabs from SARS-CoV-2 infected patients (3). Genes involved in cytokines and in ammatory response pathways were the ones that overlapped the most, demonstrating that anti-viral innate immune responses are common with SARS-CoV-2 infections. In addition, the up-regulation of differentially expressed genes related to an in ammatory response in COVID-19 patients can result in the induction of interleukin-6 (IL-6), CXCL10 (IP-10), and TNF-α with hyperactivation of Th1/Th17 responses that results the recruitment and activation of pro-in ammatory neutrophils and macrophages into the airways 31 . This has been proposed as the prime reason for failure to resolve in ammation in severely symptomatic patients 31,32 .
To better understand the contribution of epithelial cellular responses to SARS-CoV-2, we compared differentially expressed genes in the respiratory secretion of adults infected with SARS-CoV-2 to those that were expressed in HNO infected with SARS-CoV-2. . This could re ect the difference in cellular complexity between the cell lines and greater diversity of the HNO epithelium resulting in fewer overlapping up-regulated genes. HNO204 is a complex pseudostrati ed epithelium composed of at least 9 different cell types including ciliated, goblet, secretory and basal cells 33,34 . In contrast, the Calu-3 cell line, was generated from a bronchial adenocarcinoma, a submucosal gland cell line of a single cell type 35 Previous studies have demonstrated high expression of ACE2 in SARS-CoV-2 infected nasopharyngeal samples and these were greatly elevated in high viral load subjects, suggesting that higher replication occurs with increased receptor expression 22 . In our cohort we did not observe a statistically signi cant increase in ACE2 expression in both extremely high and low viral load groups. However, the expression of ACE2 was elevated in our HNO infected with SARS-CoV-2 but not TMPRSS2, which has increased expression in nasal airway epithelial brushings 36 .
In summary, our longitudinal study investigated gene expression patterns in SARS-CoV-2 infected individuals with an extremely high viral load displayed strong immune responses that decreased over time, and eventually became comparable to those with low viral loads. We detected hundreds of upregulated genes that were highly correlated to the SARS-CoV-2 viral load. Enriched cellular pathways involved in the innate immune response, antiviral interferon responses were observed in other cohorts of SARS-CoV-2 infected adults. A limited but highly signi cant up-regulated gene response overlapped with our human nose organoid line, a complex pseudostrati ed ciliated epithelium, suggesting that the gene expression pro le detected in SARS-CoV-2 infected adults is generated from both the epithelial and cellular immune responses. In conclusion, high SARS-CoV-2 viral loads primarily elicit a heightened host immune response for the control of viral replication and clearance.

Materials and Methods
Ethical Approval RT-PCR testing was performed as a service to BCM, the collection of metadata was performed under an Institutional Review Board of Baylor College of Medicine approved protocol (H-47423) with waiver of consent. All methods were performed in accordance with the relevant guidelines and regulations.

Study cohort
Ten extremely high, viral load SARS-CoV-2 positive cases were matched to 10 low viral load SARS-CoV-2 positive adults, and 10 stable adults (SARS-CoV-2 negative controls) who were cleared for having an outpatient surgical or aerosol generating procedure. The cases and controls were selected from our population of 17,644 adults (24,822 samples) evaluated in the outpatient clinics at Baylor College of Medicine (BCM) and their a liate institutions from March 18, 2020, through January 16, 2021, as previously described 14 . Three distinct adult populations were tested: 1) symptomatic employees utilizing occupational health services, 2) patients evaluated at medical and surgical clinics, and 3) patients who required clearance for an out-patient surgical or aerosol generating procedure. Serial samples were obtained from individuals who came back to be tested for evidence that the virus was cleared or were enrolled as sub-study to determine the viral shedding kinetics. Testing for SARS-CoV-2 was performed in our Clinical Laboratory Improvement Amendments (CLIA) Certi ed Respiratory Virus Diagnostic Laboratory (ID#: 45D0919666).
The extremely high viral load cases consisted of adults with an extremely high viral load (Ct <16) for the N1 target on their rst mid-turbinate (MT) sample and had at least two subsequent positive MT samples 14 . Of the 104 individuals with an extremely high viral load in their rst test, 30 individuals met the criteria for multiple positive samples over the ensuing 4 weeks. Adults from two other groups were matched to each extremely high viral load case: a low viral load (Ct 31-<40) SARS-CoV-2 positive adult (SARS-CoV-2 low viral load) and an otherwise stable control who tested negative for SARS-CoV-2 (SARS-CoV-2 negative control) and was cleared for an out-patient surgical or aerosol generating procedure. Of the 453 individuals with a low viral load in their rst test, 126 individuals met the criteria for multiple positive samples over the ensuing 4 weeks. The extremely high viral load cases were matched to the other two groups by gender, week of rst test (+ 1 week), age (+ 1 year) and zip code (5 digits). If a match could not be found the range of the factors were expanded to + 3 weeks of rst test, + 10 years and 3 digits for the zip code. The ten extremely high viral load cases were randomly selected from our pool of 30 individuals with an extremely high viral load with multiple positive MT samples. The best matched SARS-CoV-2 low viral load case and negative control were then selected for each extremely high viral load case.

RNA extraction, library preparation and sequencing
Samples were extracted using the Qiagen RNeasy mini kit (#74104 rev. 10/19) following the manufacturer's protocol for samples <5e6 cells. Samples were eluted in 50ul RNase-free water. RNA quality and quantity were estimated using Agilent Bioanalyzer OR Caliper GX. To monitor sample and process consistency, 1 µl of the 1:50 diluted synthetic RNA designed by External RNA Controls Consortium (ERCC) (4456740, ThermoFisher) was added. Whole transcriptome sequencing (total RNAseq) data was generated using the Illumina TruSeq Stranded Total RNA with Ribo-Zero Globin kit (20020612, Illumina Inc.) cDNA was prepared following rRNA and Globin mRNA depletion, and paired-end libraries were prepared on Beckman BioMek FXp liquid handlers. For this, cDNA was A-tailed followed by ligation of the TruSeq UD Indexes (Cat # 20022370) and ampli ed for 15 PCR cycles following manufacturer's recommendation. AMPure XP beads (A63882, Beckman Coulter) were used for library puri cation. Libraries were quanti ed using a Fragment Analyzer (Agilent Technologies, Inc) electrophoresis system and pooled in equimolar ratios. This pool was quanti ed using qPCR to determine loading concentration for sequencing. Sequencing was performed on the NovaSeq 6000 instrument using the S4 reagent kit (300 cycles) to generate 2x150bp paired end reads.

Primary Analysis for Total RNASeq
The RNA-Seq analysis pipeline cleans and processes raw RNA sequencing data (FASTQs), providing robust QC metrics and has the exibility to map the reads to GRCh38 reference genome (after excluding the alternate contigs). The latest versions of software for sequence alignment (STAR v2.7.3a), for marking of duplicate reads (Picard v2.

Statistical analysis
All p-values were two-sided unless otherwise speci ed. We performed all tests using log2-transformed gene expression values. False Discovery Rates (FDRs) due to multiple testing of genes were estimated using the method of Storey and Tibshirini 40 . Even in instances of nominally signi cant genes only moderately exceeding chance expectations by FDR, the nominally signi cant genes were found in downstream enrichment analyses (involving functional gene sets and results of external SARS-CoV-2related RNA-seq datasets) to contain molecular information representing real biological differences. We evaluated enrichment of GO annotation terms 41 within sets of differentially expressed genes using SigTerms software 42 and one-sided Fisher's exact tests. Visualization using heat maps was performed using JavaTreeview (version 1.1.6r4) 43,44 . Gene ontology (GO) analysis of DEGs used in the upset plot ( Figure S1) was performed using the web-based Database for Annotation, Visualization, and Integrated Discovery (DAVID; version -v2023q1) 45,46 Data Availability The RNA-seq dataset of serially collected samples and of nose organoids will be deposited at Gene Expression Omnibus (GEO) (GEO accession number pending). In terms of previously published data, we