Uncovering strain- and age-dependent differences in innate immune response 1 to SARS-CoV-2 infection in nasal epithelia using 10X single-cell sequencing

Assessing the impact of SARS-CoV-2 variants on the host is crucial with continuous emergence of new 29 variants. We employed single-cell sequencing to investigate host transcriptomic response to ancestral and 30 Alpha-strain SARS-CoV-2 infections within air-liquid-interface human nasal epithelial cells from adults 31 and adolescents. Strong innate immune responses were observed across lowly-infected and bystander cell- 32 types, and heightened in Alpha-infection. Contrastingly, the innate immune response of highly-infected 33 cells was like mock-control cells. Alpha highly-infected cells showed increased expression of protein 34 refolding genes compared with ancestral-strain-infected adolescent cells. Oxidative phosphorylation- and 35 translation-related genes were down-regulated in bystander cells versus infected and mock-control cells, 36 suggesting that the down-regulation is protective and up-regulation supports viral activity. Infected adult 37 cells revealed up-regulation of these pathways compared with infected adolescents, implying enhanced pro- 38 viral states in infected adults. Overall, this highlights the complexity of cell-type-, age- and viral-strain- 39 dependent host epithelial responses to SARS-CoV-2 and the value of air-liquid-interface cultures. 40 41 42


119
Comparative loss of cilia observed in Alpha-infected cultures compared with WT in adults, but intact cilia in Alpha-infected adolescent cultures.

120
Cellular differences also observed in child 3 (donor 6). Stained for α-tubulin (AcTub, green), nucleoprotein (NP, red) and nuclei (DAPI, blue). Both The traditional landscape of the human nasal epithelium is mainly composed of ciliated, basal and secretory 129 cells. 21 Additionally, rarer cell-types such as ionocytes, 14,22 deuterosomal cells 23,24 and transitional cell-types 130 with cell signatures from more than one cell-type may be present. 25 In our data, we observed ciliated, goblet, 131 basal, suprabasal, secretory, cycling-basal, brush/tuft, deuterosomal and ionocyte cells -consistent with 132 other scRNA-seq data of the human airway epithelium 23,24 (Figures 2a & Table S2). We additionally found 133 the presence of cell-types which were unable to be clearly classified into any one cell-type. These cells 134 were assigned as transitional cell-types, which included Secretory-Ciliated and Goblet-Ionocyte cells 135 (Figure 3a). Secretory-ciliated cells have been identified previously in human airway epithelial cells 136 (HAECs). 25 Additionally, we observed a small cluster of Goblet-Ionocyte cells with high expression of 137 markers of Goblet (MUC5AC) and Ionocytes (CFTR), which have not been previously described to our 138 knowledge. 23

148
We describe cells with fewer than 10 viral unique molecular identifiers (UMI) counts as being uninfected, 149 and cells with at least 10, 100, 1,000 and 10,000 viral UMI counts as having a low, medium, high or very-150 high level of infection respectively (Figure 3b). Of all the detected cell-types, the largest group of cells 151 were suprabasal cells with >20,000 cells overall, which were mostly lowly infected or uninfected. Ciliated 152 and Secretory-Ciliated cell-types had the highest proportion of medium, high or very-high levels of 153 infection. Within these groups we observed a clear separation of cells with mostly high-very-high viral 154 loads and have labelled these sub-clusters as "Ciliated+SC2" and "Secretory-Ciliated+SC2". These 155 subclusters showed 50.7% and 34.8% of cells with very-high level of infection, respectively (Figure 3c). 156 The results are consistent with the SARS-CoV-2 cellular tropism shown in the literature, as ciliated cells 157 are most susceptible to the SARS-CoV-2 virus in the HAECs. 14,16,25 Although the susceptibility of secretory-158 ciliated cells to SARS-CoV-2 has been noted in HAECs, 19,25 the exceptionally high rate of infection of 159 infected secretory-ciliated cells has not been identified previously, to our knowledge. We also found two 160 subsets of goblet and secretory cells -"Goblet+IFN-stim", "Secretory+IFN-stim/Secretory-3" exhibiting 161 high interferon (IFN) stimulation, with elevated levels of IFN-stimulated genes (ISGs), which increased 162 proportionally to infection level ( Figure S3 & Data S1-2). While most infected cells were classified as 163 lowly infected, we identified cells with high and very-high levels of infection within Secretory-Ciliated 164 (5.59%), Secretory+IFN-stim (4.59%), Goblet+IFN-stim (3.40%), Deuterosomal (2.34%), Ciliated-1 165 (1.40%), Goblet-Ionocyte (1.30%) and Ciliated-2 cells (1.06%) (Figure 3b). This highlights the SARS-166 CoV-2 susceptibility of secretory/goblet cells in addition to ciliated and transitional cell-types, which again 167 is consistent with the literature. 26 Alpha-infected datasets showed the highest proportion of infected cells, 168 followed by WT 72 hpi and then WT 48 hpi-infected cells (Figures S4a-e). We note very few cells (0.33% 169 of all WT 48 hpi data from adults and adolescents) were infected in WT-strain-infected cells harvested at 170 48 hpi, in comparison to 72 hpi with both the WT and Alpha strains (Figure 3c). 171 172 Additionally, the changes in cell-type distributions upon infection against the mock-control datasets were 173 compared within the same age-groups. As expected, the increase in Ciliated+SC2, Secretory-Ciliated+SC2, 174 Firstly, we assessed the general levels of entry-associated genes -ACE2, TMPRSS2 and FURIN in all 190 clusters which involved an assortment of one or more of the following categories -uninfected, bystander 191 and infected cells. ACE2 mRNA expression levels were found to be low across all cell-types in our data 192 regardless of infection or treatment status, consistent with previous studies 16 (Figures 4a-c). Furthermore,193 TMPRSS2 and FURIN expression was higher than ACE2 levels across many different cell-types, with 194 TMPRSS2 expression being the lowest in basal and FURIN lacking in ciliated cell populations (Figure 4a-195 c). We noticed comparatively higher ACE2 expression in a subset of secretory cells - 196 and Secretory+IFN-stim (Figure 4a). These cells showed high-levels of IFN responses with elevated gene 197 expression of ISGs and were only robustly present in conditions usually associated with higher levels of 198

SARS-CoV-2 infection (Alpha-infected datasets in adults & adolescents and WT 72 hpi datasets in 199
adolescents) (Data S1). The higher levels of ACE2 in secretory/goblet cells has also been shown previously 200 in the human nasal epithelia. 27,28 Interestingly, despite being associated with higher viral-load datasets, these 201 clusters only involved low-levels of viral RNA (Figure 4b). Also, ACE2 was minimally expressed in 202 highly-infected cell-types such as Secretory-Ciliated+SC2 and Ciliated+SC2 (Figure 4a). We then 203 performed differential expression (DE) analysis between infected and mock-control datasets to investigate 204 whether SARS-CoV-2 infection caused up-regulation of ACE2. Significant up-regulation of ACE2 in only 205 the same high IFN-stim clusters (Goblet+IFN-stim, and Secretory+IFN-stim) and Alpha-infected Ciliated-206 1 cells were observed, which were also largely lowly infected (Figures 4d-k & 4b & Table S3, padj < 207 0.05). Furthermore, other infected cell-clusters with robust IFN-responses were not found to significantly 208 up-regulate ACE2 (Figure 6a). These results revealed that while SARS-CoV-2 infection generally causes 209 an up-regulation of ISGs when compared with control cells, not all cells are able to elicit a strong interferon 210 stimulated gene (ISG) induction, especially in the high-viral load cell-types. Even within the group of cell-211 types with robust ISG induction upon infection, only three cell-types with elevated ISG induction showed 212 up-regulation of ACE2 (Figures 4d-k). This suggests that SARS-CoV-2 infection alone does not lead to 213 increased ACE2 expression but requires higher IFN-responses compared with other ISGs for its up-214 regulation. To confirm these results, the level of ACE2 was compared in bystander cells compared with 215 mock-control cells. Bystander cells are cells which are exposed to the pathogen but have not been identified 216 to be infected, 14 but should be exposed to IFN through paracrine activity from neighboring infected cells. 217 Following the assumption that ACE2 requires high IFN-stimulation compared with other ISGs to be up-218 regulated, bystander cells did not show a significant increase in expression level of ACE2 compared to 219 mock-control, again, despite up-regulation of other ISGs such as IFITM3 (Figures S5a-c, 5b). Furthermore, 220 ACE2 expression correlated positively with levels of STAT1 (Figures 4d-k). These results are in line with 221 evidence showing that ACE2 is stimulated by IFNs and has expression correlation with STAT1. 27 In the 222 study by Ziegler et al.,27 the promoter of ACE2 was found to contain two STAT1 binding sites, revealing 223 the importance of this relationship.

Cell-types with low levels of infection show increased innate immune responses compared with cell-237 types with high levels of infection 238
We investigated whether the SARS-CoV-2-infected cells showed any differences in host immune response 239 compared with mock-control cells on a cell-type basis. Within both adults and adolescents, Alpha-infected 240 cells showed strong enrichment of infection-specific GO terms such as defense response to virus, response 241 to virus, type I interferon signaling pathway, innate immune response and the negative regulation of viral 242 genome replication (Figure 5a). Similarly, in WT-infected adolescents, IFN-response related terms were 243 enriched, in similar groups of cells. The reactome pathway analysis largely agreed with the GO biological 244 term analysis ( Figure S6a). Comparison of expression between Secretory-Ciliated+SC2 cluster in infected 245 samples with the Secretory-Ciliated cluster from mock-control revealed a marked reduction of significant 246 enrichment of these GO terms, regardless of age or strain. Furthermore, the absence of GO terms enriched 247 in Ciliated+SC2 cells versus Ciliated-1 in mock-control was observed. These results suggest that these cells 248 were either unable to produce a robust IFN response or showed viral suppression of the IFN response, 249 leading to higher viral loads compared with other cell-types (Figures 5a). Additionally, for Alpha-250 infections in both age groups and WT-infected children we note these highly infected Secretory-251 Ciliated+SC2 cell clusters showed down-regulation of genes related to cilium movement and cilium 252 assembly in comparison to mock-control Secretory-Ciliated cells, consistent with loss of cilia observed 253 from microscopy with Alpha-infected adult cells (Figure 2). In adolescent samples, we found enrichment 254 of these processes in other cell-types such as Suprabasal cells in infected vs mock-control comparisons, 255 showing the amplified enrichment in adolescents compared with adults upon infection. However, 256 microscopy results showed comparatively intact cilia in Alpha-infected adolescents compared with Alpha-257 infected adults, suggesting the disconnection between protein and mRNA levels. Furthermore, in contrast 258 to our data, Ravindra et al. 14 showed up-regulation of cilia-related genes in infected human bronchial 259 epithelial cells (HBECs) compared with control (Table S1). Alpha-associated Basal-1, Ciliated-3, Secretory-2, Suprabasal in adults and Basal-1, Ciliated-3, Secretory-278 2 in adolescents compared with control. Also, similar up-regulation was observed in Basal-1, Ciliated-1, 279 Ciliated-3, Cycling-basal, Secretory+IFN-stim and Suprabasal cells in WT 72 hpi-associated adolescent 280

cells, with an absence of enrichment in negative regulation of viral genome replication in Suprabasal cells. 281
Similarly, reactome pathways such as interferon alpha/beta signaling and interferon signaling were 282 enriched in these datasets in the same cell-clusters ( Figure S6b). This reciprocates the results observed in 283 Alpha-strain-infected cells as described above. These results suggest that consistent with existing evidence Next, accounting for the baseline expression in the respective control cells and genetic variability between 297 all donors, we compared the differences between infected adults and adolescents. Enriched GO biological 298 terms included oxidative phosphorylation, mitochondrial electron transport, NADH to ubiquinone, 299 cytoplasmic translation (Figure 6a). Secretory and ciliated cell-types were especially involved in 300 significant gene set enrichment in both Alpha and WT infections (Figure 6a) and cilium movement and 301 cilium assembly were largely enriched in Alpha-infected datasets. Similarly, enriched reactome pathways 302 included respiratory electron transport, translation and influenza infection ( Figure S6c). We note that our 303 bulk RNA-seq study interrogating differentially expressed genes in WT SARS-CoV-2 infections within 304 continuous epithelial cell lines compared with control showed enrichment of other viral-infected 305 pathways/terms. 30 Therefore, the enrichment of terms such as influenza infection is unlikely to be due to an 306 asymptomatic influenza infection of these ALI cultures, but rather a cross-over from a respiratory viral 307 infection pathway. Furthermore, this has been partially verified via metagenomics analysis showing no 308 clear evidence of RNA from other pathogens in the mock-control datasets from donor 6 (Data S3). Overall, 309 adolescents showed overall down-regulation of gene expression compared with adults in these enriched 310 terms/pathways. This suggests either the lower requirement of these genes/pathways or the increased viral 311 suppression of these pathways upon infection in adolescents and absent or diminished suppression in adults, 312 supporting the idea of an age-dependent response to SARS-CoV-2.   Table S1.

Alpha-variant induces increased protein folding and innate immune responses compared to WT 324 strain 325
We next explored differences in host responses to the Alpha variant compared with the WT-strain of SARS-326 CoV-2. Anti-viral terms such as response to virus, negative regulation of viral genome replication, innate 327 immune response and defense response to virus were enriched in both adults (Ciliated-1, Secretory-1) and 328 adolescents (Basal-1, Deuterosomal, Secretory-2) (Figure 6b). In terms of reactome pathways, interferon 329 signaling, interferon alpha/beta signaling, antiviral mechanism by IFN-stimulated genes were enriched in 330 the same cells ( Figure S6d). The genes involved in these processes were up-regulated in the Alpha-variant 331 infections compared with WT. These results highlight the heightened anti-viral host responses in Alpha vs 332 WT infections. Furthermore, these results showed that although similar processes are elicited in the two 333 age-groups, there is a divergence within the groups of cells which are involved. Additionally, we observed 334 similar up-regulation of genes involved in protein refolding and response to topologically incorrect protein 335 GO biological terms in Ciliated+SC2 and Secretory-Ciliated+SC2 datasets in adolescents (Figure 6b). 336 Therefore, this provides evidence that the Alpha variant elicits a greater post-translational activity related 337 to refolding aberrantly folded/unfolded proteins in the most infected cluster of cells, at least within 338 adolescents. We note that we did not have enough WT-infected cells in those clusters in adults to compare 339 with adolescents. 340 341

Translation and oxidative phosphorylation up-regulation in infected cells vs bystander cells 342
Infected cells were then compared with bystander cells to understand differences in host response to 343 infection by SARS-CoV-2 compared with IFN stimulation. Firstly, we noted that Ciliated-3 cells had most 344 enrichment of GO biological terms out of all datasets (Figure 7a). In Ciliated-3 cells, we observed an up-345 regulation of genes involved in translation and oxidative phosphorylation compared with bystander cells 346 in both Alpha-infected adults and adolescents, and oxidative phosphorylation in WT 72 hpi adolescents 347 ( Figure 7a). Additionally, we also observed the enrichment of oxidative phosphorylation in Basal-1 cells 348 and Secretory-2 cells in adults but only additionally in Secretory-2 cells in adolescents, revealing age-349 dependent responses. Like the infected vs control datasets, the enriched GO terms overwhelmingly involved 350

up-regulated genes in infected cells compared with the bystander cells. Enriched reactome pathways 351
included the citric acid (TCA) cycle and respiratory electron transport, respiratory electron transport, 352 adenosine triphosphate (ATP) synthesis by chemiosmotic coupling and heat production by uncoupling 353

proteins, infectious disease, influenza infection, interleukin-1 signaling, viral mRNA translation and 354
translation in Ciliated-3 in Alpha-associated adults and adolescents ( Figure S6e). Additionally, the citric 355 acid (TCA) cycle and respiratory electron transport, respiratory electron transport, ATP synthesis by 356 chemiosmotic coupling and heat production by uncoupling proteins were enriched in Basal-1 and 357 Secretory-2 cells in Alpha-associated adults, Secretory-2 in Alpha-associated adolescents and Ciliated-3 358 cells in WT 72-associated adolescents. Furthermore, we then combined all cells from Ciliated-cell clusters 359 (Ciliated 1-3 & + SC2) and performed DE analysis between infected and bystander cells. From these results, 360 we noted the consistent up-regulation of genes NFKBIA, JUN and SOX4 in Alpha-infected cells when 361 compared with bystander cells in both age-groups as well as NFKBIA and JUN in WT-infected cells vs 362 mock-control cells in adolescents (Figures 7b-e). This is generally consistent with results from the study 363 by Ravindra et al. 14 However, in contrast to the study by Ravindra et al.,despite similar strain (USA 364 ancestral) and MOIs used, there was no significant up-regulation overall as well as these three specific 365 genes (NFKBIA, JUN, SOX4) in the WT adult data (Figure 7d) and nor did we find many significantly 366 down-regulated genes in all comparisons (Figures 7b-e). These differences may be attributed to the tissue 367 type (nasal vs bronchial cells), method for DE (pseudo-bulk vs non-pseudo-bulk) or strain-differences 368     Table S1.

Donor variability in ALI-HNECs affects cellular diversity and infection 401
In our results, one child donor (donor 6) showed notable differences to other donors. In uninfected 402 cultures the clustering data showed the expansion of goblet cells, and a separated group of basal (Basal-403 2), secretory (Secretory-2) and ciliated (Ciliated-2) cells ( Figure S3a). This was also largely observed in 404 the other treatment datasets (WT 48 WT 72 and Alpha 72 hpi) (Figures S3b-d), which suggests that these 405 cells clusters were present in donor 6 datasets prior to infection in the natural state and persists even with 406 infection. The infection level of this donor was also reduced in comparison to the other adolescent donors 407 ( Figures S4b&d). Physically, the cells were larger and more oval in shape compared with the other two 408 child donors, and cell recovery and viability were slightly reduced ( Table S4). This could have 409 contributed to a lower cell sample size, leading to a certain bias in sampling, causing a different mix of 410 cells to be sampled from the other donors. To understand the reasoning behind the differences observed in 411 this donor at a deeper level, we applied various tests. Firstly, epithelial-mesenchymal transition (EMT) 412 can occur in ALI-culturing, 31 which can decrease the level of infectable cells. When the mesenchymal 413 marker Vimentin (VIM) in mock-control cells of donor 6 was compared with mock-control samples in all 414 other donors, VIM was up-regulated (padj < 0.05, Figure S7a). However, no significant DE was observed 415 when compared with only the other child donors (padj ≥ 0.05, Figure S7b). Then, to rule out any 416 asymptomatic co-infections within this donor, we applied a metagenomics pipeline Kraken2Uniq to 417 search for evidence of reads mapping to any common pathogens within mock-control cells (see 418 Methods). However, we were unable to find any clear indications of co-infections with most reads 419 mapping to human and unclassified categories (Data S3). Next, we searched for any immune-state 420 differences between mock-control cells. When donor 6 was compared with all other donors, we also did 421 not find a strong enrichment of immune-related GO terms ( Figure S7c). However, we note that when 422 donor 6 was also compared with the other two adolescent donors, we observed a decrease in IFI27L1 423 (log2FC=-1.14). This gene is part of an IFN-inducible protein gene family of IFI27, which has been 424 shown to be increased in the blood of severe COVID patients. 32 Therefore, perhaps the down-regulation 425 of IFI27L1 could contribute to protective anti-viral functions in donor 6. However, overall, these results  (Figure 3b). The infection of such transitional cells have been noted before. 25 We speculate 438 that this phenomenon can be attributed to 1) secretory cells being precursors to ciliated cells, 34 and therefore The main objective of this study was to understand the age and strain-dependent responses to SARS-CoV-461 2. At the 72 hpi time point, the TCID50 (Figures 1a&b) and short-read sequencing data (Figure 4c) showed 462 that the Alpha variant infected adults and adolescents similarly and the WT strain resulted in much higher 463 viral load in adolescents compared with adults. Also, generally the Alpha-variant infections yielded higher 464 viral titers (Figures 1c&d) and viral reads (Figure 3c) than WT strain-infections. The elevated viral titers 465 with Alpha-infections compared with WT was expected, due to the increased transmissibility in the variant 466 via enhanced receptor binding affinity due to the N501Y mutation. 6,13 However, the distinctly increased 467 viral titers and reads in WT-infected adolescents compared with adults in all replicates were unexpected. 468 This is because it has been thought that children, especially with WT-infections are less susceptible to 469 SARS-CoV-2. 36 One potential reason for this is that the lower abundance of ACE2 receptors in the upper 470 airways of children compared with adults. However, there are mixed reports regarding age-dependent ACE2 471 mRNA levels, showing either that ACE2 mRNA expression can increase with age in the nasal epithelia, 37 472 or there are no age-dependent effects. 38 In this study, we did not observe significant differences in ACE2 473 mRNA expression level between infected adolescents vs adults as well as mock-control adolescents and 474 adults in all tested cell-clusters (padj ≥ 0.05), perhaps due to this study involving a single-cell method, 475 preserving the heterogeneity in the sample. Therefore, our data appeared to be in line with previous studies 476 suggesting that ACE2 transcriptional levels do not correlate with susceptibility to SARS-CoV-2. 16 Despite 477 the lack of differences between ACE2 at the mRNA level, we note that the protein level may be contrasting, 478 as staining of ACE2 protein has been shown previously in ALI-human nasal epithelial cells (HNECs), 39 479 which could cause these differences in viral load. Zhu et al. 40 used immunofluorescence staining to show 480 lower surface levels of ACE2 in children compared with adults in ALI-HNECs. However, the authors noted 481 no quantitative protein level differences in ACE2 or TMPRSS2 between adults and children via western 482 blotting. Interestingly, in the same study, the authors revealed that the WT-strain infects less in children vs 483 adults, which is directly in contradiction to our results 40 and align more with the stain results than the 484 western blot results. We note that these contrasts between our study and the study by Zhu et al. may be due 485 to the differences in ages in the child donors, where our study has focused on adolescents of ages 12-14, 486 whereas this contrasting study involved younger children with ages under 12. However, consistent with our 487 results, another study showed lower viral titers in older adult HNECs compared with children and younger 488 adults, 41 though in our study the ages of adults match more with the younger adult group in the study by 489 Capraro et al. 41 Although unlikely, there is also a possibility that the MOI used for infecting the adolescents 490 may have been unintentionally increased compared with the adults with the WT infections, e.g. due to 491 overestimation of the number of host cells. Upon investigating the number of cells counted for the 492 uninfected cultures for each donors, the results showed donor-donor variability and surprisingly showed 493 lower MOI in child donors compared with adult donors and higher MOI in donor 6 ( Table S4). This 494 contrasts with the idea of MOI overestimation and further work will be required to deconvolute the 495 relationship between the age-dependent host responses to SARS-CoV-2. 36 496 497 ACE2 mRNA expression was low across all cell-types (Figures 4a&c). However, we noticed that the gene 498

was particularly up-regulated in IFN-stimulated populations of cells such as Secretory+IFN-stim and 499
Goblet+IFN-stim cells and Ciliated-1 cells (Figures 4d-k). This was in line with evidence that ACE2 is an 500 ISG. 27 Interestingly, a similar study involving ALI-HAECs infected with SARS-CoV-2 did not show 501 increased levels of ACE2 mRNA after infection 14 and similarly between HNECs derived from COVID-19 502 patients vs healthy controls. 16 We speculate that perhaps if an IFN-stimulated cluster was also separated 503 from the main body of cells, these datasets may also show similar results. Furthermore, ACE2 was not found 504 to be up-regulated in bystander cells when compared with mock-control cells, although induction of ISGs 505 occurred in these cells (Figures S5a-c) showed enrichment of terms and pathways related to innate immune responses, and an absence of these 513 responses were observed in highly-infected cell-clusters (Figure 5a). We speculate that these cells did not 514 mount a robust IFN-response during early infection, leading to a higher viral load. Otherwise, multiple 515 virion particles may have infected these cells. This would lead to increased antagonization of IFN-responses 516 and therefore high viral loads as certain SARS-CoV-2 open reading frames have antagonistic properties to 517 IFN-responses. 43 Finally, due to the low MOI applied, multiple rounds of viral infection may occur, which 518 would lead to some cells infecting earlier than others. Hence, these high viral-load cells may be part of this 519 early-infection group. The challenges of identifying the source of cell-to-cell heterogeneity in virus-infected 520 scRNA-seq datasets have been reviewed by Suomalainen and Greber (2021). 44 521 522 We also assessed the differences in responses between infected adolescents and infected adults. Loske et 523 al. 17 showed the pre-primed immunity to SARS-CoV-2 as well as increased ISG induction in children 524 compared with adults. Contrary to this finding, we did not observe strong enrichment of IFN-response-525 related GO terms when infected adolescent ALIs were compared with infected adult ALIs (Figure 6a). 526 However, we did find the enrichment of translational, oxidative phosphorylation and cilia-related terms, 527 where most genes were found to be down-regulated in the infected adolescents compared with the infected 528 adults. Similarly, with reactome pathways, the level of expression of genes involved in translation, 529 respiratory electron chain and influenza infection pathways were decreased compared with infected adults 530 ( Figure S6c). This suggested that in the infected state, adults require higher levels of these genes compared 531 with adolescents. Considering that cilia-related GO terms was found to be enriched with the significant 532 genes when expression levels of COVID-19 airways was compared with healthy controls, 19 potentially the 533 higher requirement of these genes in adults may be indicative of a greater damage caused by SARS-CoV-534 2 in adults compared with adolescents. However, the direction of magnitude of DE was opposite in our 535 ALI-HNEC data when compared with the ALI-HBEC datasets ( Table S1). 14 In addition, when we 536 compared mock-control cells between adolescents and adults, we did not observe significant up-regulation 537 of viral-sensing related genes -IFIH1, DDX58, DHX58 -which were found to be up-regulated in healthy 538 child nasal airway cells compared to adults. 17 Overall, these results contrast with the evidence of a stronger 539 IFN response gene expression in children compared with adults. 14,17,19 540 541

The Alpha-variant infected cells showed increased expression of genes involved in protein refolding and 542 response to topologically incorrect protein, compared with WT-infected cells in adolescents, in the clusters 543
with the highest level of infection (i.e. Ciliated+SC2 and Secretory-Ciliated+SC2) (Figure 6b). Under 544 normal conditions, the protein refolding response is not activated, and is switched on after an accumulation 545 of unfolded/misfolded proteins occurs under endoplasmic reticulum (ER) stress. 45 The aggregation of 546 unfolded proteins may occur as a host defense mechanism, or viral manipulation to increase replication or 547 viral immune evasion. 46,47 The induction of the unfolded protein response (UPR) due to ER stress has been 548 documented with the spike protein of SARS-CoV-1 48 and also SARS-CoV-2. 49 Overall, this increased 549 activity in adolescents may be in part due to the increased transmissibility of the Alpha strain increasing 550 the build-up of misfolded/unfolded proteins. This may be also due to the accumulation of excess viral 551 proteins during infection, overwhelming the cellular architecture and therefore negatively affecting both 552 host and viral post-translational modifications and proper protein folding. We note that the matching data 553 with Ciliated+SC2 and Secretory-Ciliated+SC2 clusters in adults was unavailable. While we cannot 554 comment on an age-dependent/independent effect, we speculate these processes could have been 555 reciprocated in the adult datasets had we been able to analyze these datasets due to strain-dependent viral-556 load being observed between WT and Alpha-infections in adults. 557

558
We demonstrate the importance of oxidative phosphorylation in infection. Oxidative phosphorylation is an 559 integral part of the energy production process, 50 which produces reactive oxygen species (ROS) and an 560 excessive level can lead to oxidative stress and inflammation. We observed down-regulation of oxidative 561 phosphorylation in bystander cells versus both infected and mock-control cells (Figures 5b & 7a), 562 indicating that down-regulation may occur as a result of IFN-stimulation, which is then reverted to original 563 levels by viral factors. Down-regulation oxidative phosphorylation may lead to the reduction in ATP in the 564 host cell, which is required for viral processes for positive strand viruses, 51  We have also shown the importance of oxidative phosphorylation between age groups. We found that 575 infected adolescents had lower levels of expression oxidative phosphorylation genes than infected adults 576 (Figure 6a), which may provide a more anti-viral state in adolescents. These results match some results 577 from a recently published study, which showed elevated levels of oxidative phosphorylation in 578 nasopharyngeal samples from adults when compared with pediatric patients, although these changes were 579 non-significant after FDR-correction. 59 580 581 Similarly, the GO term translation was down-regulated in bystander cells versus both infected and mock-582 control cells (Figures 5b & 7a) but not between infected and mock-control (Figure 5a), and also down-583 regulated in adolescents versus adults (Figure 6a). This coupling of translation and oxidative 584 phosphorylation was not surprising due to ROS being found to increase phosphorylation of eukaryotic 585 translation initiation factor 2a (eiF2a) with enterovirus infection. 60 The down-regulation of translation thus 586 appears to be driven by host-factors, while its up-regulation driven by viral factors, and down-regulation in 587 adolescents indicates a stronger anti-viral state in this group. 588 589

Limitations 590
Whilst we have utilized a pseudo-stratified airway model approach for our study which is superior 591 compared with continuous clonal cell lines, we acknowledge that these results may not directly translate to 592 in vivo situations. Particularly, the absence of immune cells in the model may distort the relevance of these 593 results within this study. However, via this model, we were able to determine the epithelial immune 594 responses to SARS-CoV-2 without confounding effects from immune-epithelial cell interactions. 595 Additionally, we have also used a low MOI of < 0.02, which may be more clinically relevant as low numbers 596 of virions will initiate an in vivo infection, but this means that the infection stage of each cell will be 597 temporally asynchronized. Finally, one donor (donor 6) showed some cellular and infection differences 598 compared with other donors. However, most cell-types aligned with other donors, as shown by the 599 clustering analysis and this effect has been minimized by applying filters such as minimum number of donor 600 replicates and minimum number of cells for pseudo-bulking. Furthermore, due to ambient RNA being able 601 to be encapsulated into 10X droplets, there is a potential overestimation of true viral RNA load in each cell. 602 We have applied a threshold of 10 UMI viral counts per cell to be deemed as infected to account for the 603 potential contamination, according to the empirical threshold found by Ravindra et al. 14  Melbourne (HREC/2057111). 33 Written consent was obtained from all participants (or participant's 635 guardian) prior to collection of biospecimens. All samples were de-identified before tissue processing. 636 Human nasal epithelial cell culture methods have been described previously. 33,61-64 Briefly, three healthy 637 adult (PDI-5 (Male/32Y), PDI-1 (Female/32Y), and PDI-4 (Male/26Y)) and child/adolescent (PDI-8 638 (Female/12Y), PDI-9 (Female/13Y), and PDI-10 (Male/14Y)) biobanked cells were utilized. Nasal 639 turbinate brush samples were taken before the COVID-19 pandemic, ensuring no subject encountered 640 SARS-CoV-2 exposure prior to in vitro infection. To initiate differentiation of ALI cultures, cryo-preserved 641 cells were thawed and seeded on to 6.5 mm Transwell inserts (Corning) pre-coated with PureCol-S collagen 642 type I (Advanced BioMatrix). The cells were incubated at 37°C and 5% v/v CO2 until confluency in 643 PneumaCult TM -ExPlus media (STEMCELL Technologies) for 4-7 days before being switched to ALI 644 culture conditions by removing the apical media and feeding the basal side with PneumaCult TM ALI 645 medium (STEMCELL Technologies). The cultures were incubated 3-4 weeks to achieve mucocilliary 646 differentiation evidence by the presence of mucus and beating cilia. 647 648

Immunofluorescence and confocal microscopy 667
Immunofluorescence and confocal microscopy imaging was performed as previously described. 33 In brief, 668 at the time of harvest the cells were washed three times with PBS++. Cells were then fixed with 4% 669 paraformaldehyde (#15710, Electron Microscopy Sciences, USA) for half an hour at room temperature. 670 The fixative was removed and replaced with 100 mM glycine in PBS++ for 10 minutes to neutralize the 671 remaining fixative. Cells were permeabilized with 0.5% Triton-X in PBS++ for half an hour on ice, before 672 being washed 3 times with PBS++ at room temperature. At this stage, the membranes were carefully excised 673 from the Transwell inserts, cut into half, one for test antibodies and the other for control antibodies, and 674 blocked for 90 minutes at room temperature in immunofluorescence (IF) buffer (PBS++ with 0.1% bovine 675 serum albumin, 0.2% Triton, 0.05% Tween 20) supplemented with 10% goat serum. After this, the block 676 buffer was replaced by block buffer containing the primary antibodies, anti-acetylated α-tubulin (Sigma-677 Aldrich #T7451, diluted at 1:250) and anti-SARS Nucleocapsid Protein (Novus Biologicals #NB100-678 56683, diluted at 1:200). After incubation for 48 hours at 4 o C, the primary antibodies were washed off with 679 IF buffer 3 times; then, fluorophore conjugated secondary antibodies, goat-anti-mouse Alexa Fluor 488 680 (Invitrogen #A11001) and goat-anti-rabbit Alexa Fluor 647 (Invitrogen #21244) and Hoechst, were added 681 and incubated for 3 hours at room temperature in the dark. Secondary antibodies were then washed off 5 682 times with IF buffer. The membranes were incubated with DAPI for half an hour, washed once with PBS++ 683 and transferred to slides where they were mounted in FluoroSave Reagent (#345789 EMD Millipore). The 684 confocal microscopy imaging was acquired on the Zeiss LSM 780 system. The acquired images were 685 processed using ImageJ software. 686

10X Genomics single-cell preparation 688
The ALI-cultured HNECs were prepared for the 10X Chromium step according to the Single Cell Protocols 689 Cell Preparation Guide General Sample Preparation RevC (10X Genomics). Briefly, cells were dissociated 690 using trypsin and filtered through a 40 µm strainer and pipette-mixed to ensure a single-cell suspension. 691 The cells were washed with PBS with 0.04% BSA. Once cells were counted, they were harvested for mock-692 The cycling parameters were as follows: Read 1 -150 bp, Index 1-10 bp, Index 2 -10 bp,  The sequencing was carried out by Ramaciotti Centre for Genomics at the University of New South Wales 705 (UNSW). A total of ~7.95 billion reads were acquired. 706

Genotyping 707
Genotyping for each of the donors was required to accurately demultiplex the mixed population of cells 708 used as inputs into the 10X Chromium preparation. DNA extraction for genotyping was carried out with 709 the DNeasy Blood and Tissue Kit (Qiagen) according to the manufacturer's guidelines (Purification of 710 Total DNA from Animal Blood or Cells (Spin-Column Protocol)) with minor modifications. Briefly, ALI 711 cell culture membranes were excised from the inserts and placed into tubes containing PBS and proteinase 712 K. Once Buffer AL was added, the sample was briefly vortexed and incubated at 56°C for 10 minutes with 713 a Thermomixer C (Eppendorf) at 1000 rpm. 100% ethanol was added to the reaction and tubes briefly 714 vortexed. The reaction was loaded on to a spin-column and all subsequent spins were carried out at 12,000 715 rpm except for the step 6 of the protocol, where after adding AW2 buffer, the columns were spun at 14,000 716 rpm. 200 µL of buffer AE was used to elute the DNA and passed through the column in total of three times 717 to concentrate the sample. Quality control of DNA was carried out using Qubit 4.0 Fluorometer via the 718 Qubit 1X dsDNA HS Assay Kit, BioAnalyzer 2100 using the High Sensitivity DNA Assay and NanoDrop 719 2100 Spectrophotometer (ThermoFisher Scientific). Genotyping was carried out using DNA derived from 720 each individual donor using the Infinium Global Screening Array (GSA) v2.0 BeadChip (Illumina) and 721 performed by Macrogen (Korea). The reference annotation used was GRCh37. 722 723

Data analysis 724
Genotyping analysis 725 PLINK v1.9 66 was utilized to convert the output of the genotyping data (from section Genotyping) to VCF 726 files. Firstly, the sex of the samples was checked, and this information was incorporated into the data. 727 Heterozygous haploid hardcalls, all female chrY calls were erased from the data. The resulting file was 728 converted to VCF file with '-recode'. Variants with one or more multi-character allele codes and single- Genomics) with the 'mkfastq' function. Count files were produced with Cellranger 'count' using the 738 reference package 'refdata-gex-GRCh38-2020-A'. BAM files generated from GRCh38/hg38 were lifted 739 with liftover tool CrossMap v0.5.4 67 to GRCh37/hg19 to enable incorporation of the donor genotype 740 information which was analyzed using GRCh37. Similarly, a Cellranger reference was created for SARS-741 CoV-2 reference genome using Cellranger 'mkref' using the Ensembl reference ASM985889v3, INSDC 742 Assembly GCA_009858895.3, Jan 2020 and a custom GTF file by setting the whole genome as an exon. 743 Viral counts were determined separately but similarly to host counts using Cellranger 'count'. For the viral 744 counts, raw matrices instead of filtered matrices were utilized for downstream analysis as the effect of 745 filtered matrices (i.e. filtering of artifactual cells) was not applicable to the viral counts. The VCF files and 746 the sorted GRCh37 BAM files were used as inputs to Demuxlet v11022021. 68 This enabled the donor 747 assignment to cell barcodes and estimated the number of doublets in the data. 748 749

Data filtering and unsupervised clustering 750
For downstream analysis of Illumina datasets, Seurat v4.0.5 was implemented. Seurat objects were created 751 separately for both viral and host counts data. The demultiplexing information from Demuxlet was 752 incorporated into the Seurat object using importDemux from dittoSeq v1.0.2. Firstly, the viral data was 753 separated based on infection tier as follows; uninfected < 10, low = 10-99, medium =100-999, high = 1,000-754 9,999, very-high ≥ 10,000. This information was incorporated into host meta data per cell-barcode. 755 Uninfected cells (<10 UMI counts) with exposure to virus were assigned as 'bystander' cells. Then, the 756 data was filtered for singlets, according to the Demuxlet results. Cells which had <20% mitochondrial RNA, 757 >5% ribosomal RNA were kept for analysis. Also, cells with greater than 200 and less than 9000 detected 758 genes were kept, and genes expressed in at least 3 cells were kept for analysis. Each of the 16 libraries (i.e. 759 8 x main and 8 subsample) were analyzed separately. Data was normalized and scaled using scTransform 760 v0.3.3 (using the original method) and differences in cell cycle were regressed out by the alternate workflow 761 (regressing out the G2M -S phase scores). Then, these datasets were merged using Seurat 'merge'. 762 Utilizing cell markers from scRNA-seq datasets and within the literature 14

Differential gene expression (DGE) analysis via pseudo-bulking 779
After identifying the different cell-types, DGE analysis was performed on the host via a pseudo-bulking 780 method. Briefly, the data was separated, and the counts were aggregated by a unique combination of cell-781 type, treatment, age-group, infection status and donor information. variability between donors were regressed out blocking the 'batch' (i.e. each donor) factor variable as a 791 random effect. The effect of sex was added as a fixed effect in the design matrix. We note that any 792 comparisons between Ciliated+SC2, Secretory-Ciliated+SC2, Secretory+IFN, Goblet+IFN and mock-793 control dataset was compared with controls from Ciliated-1, Secretory-Ciliated, Secretory-1&2 and Goblet 794 clusters, respectively, due to lack of control cells in the clusters. The populations of matching control cells 795 were determined by closest cell states. Gene Ontology (GO) biological and reactome pathways were 796 visualized using an in-house visualization tool -multiGO ( Table S1). The parameters used were 797 pv_thresh=0.05, enrichment pv_thresh=0.005 and logFC_thresh=1. Background lists of genes were curated 798 from all genes tested for DGE in all groups which were displayed in each multiGO analysis via setting the 799 'Background set for DE' parameter as 'gene_list' (Table S1). 800 801

Immune profiles & test for EMT 802
As donor 6 showed lower viral load (Figures S1b&d) and showed some differences in cellular 803 composition (Figures S2f & 3a-d) to other donors, we investigated the differences in immune profiles 804 and the expression of epithelial-mesenchymal transition (EMT) marker VIM for donor 6 without 805 infection. DE was performed using a similar pseudo-bulking approach as the main DE analysis. Firstly, 806 all mock-control cells from donor 6 were compared against the other two child donors, and then between 807 all other donors as a bulk analysis. This was carried out via the limma-voom method through limma 808 v3.44.3 73 and edgeR v3.30.3. 72 As with the limma-voom analysis in the main text, the differences in 809 genetic variability between donors were regressed out blocking the 'batch' (i.e. each donor) factor 810 variable as a random effect. The effect of sex was added as a fixed effect in the design matrix. Gene 811 Ontology (GO) biological and reactome pathways were visualized using an in-house visualization tool -812 multiGO. The parameters used were pv_thresh=0.05, enrichment pv_thresh=0.005 and logFC_thresh=1. 813 Background lists of genes were curated from all genes tested for DE in all groups which were displayed in 814 each multiGO analysis via setting the 'Background set for DE' parameter as 'gene_list'. 815

Asymptomatic co-infection testing 816
Following the reasons described in the previous section (Immune profiles & test for EMT), we next 817 wondered whether this was due to an asymptomatic co-infection. This was carried out using a 818 metagenomic testing approach. The output BAM files from mock-infected child donor datasets from the 819 larger group (85% cells, 'Short_read_uninfected_child') which were mapped to the human genome from 820 Cellranger were demultiplexed with data from Demuxlet. This was carried out via isolating singlet cell 821 barcodes matching to each of the child donors and extracting the data using Samtools v1.9. The 822 demultiplexed files were then filtered for unmapped reads using Samtools 'view -b -f 4' to deplete already 823 human-mapped reads. The unmapped reads were converted back into paired-end FASTQ format using 824 Cellranger's 'bamtofastq' function. Reads [from the depletion step] were classified using Kraken2Uniq 825 protocol for pathogen detection with a minimum hit groups setting of 3 74,75 and the PlusPF database based 826 on RefSeq (2022-09-08, https://benlangmead.github.io/aws-indexes/k2). Taxonomic classification reports 827 were inspected manually for the presence of viral, bacterial and eukaryotic taxa that may cause co-828 infections, evaluating abundance, number of reads classified and the number of distinct minimizers in 829 relation to the number of reads 74 . Kraken2 v2.1.2 was utilized for this analysis with the parameters '--db 830 31angmead_pluspf_64GB/ --minimum-hit-groups 3 --report-minimizer-data --threads 32--output 831 ${name}.kraken2uniq --report ${name}.kraken2uniq.report'. 832 833

Data availability 834
All code is available on Github: https://github.com/cjy-23/ALI_scRNA_seq_SC2. Raw sequencing data 835 will be released upon publication and will be available at BioProject PRJNA956316. DNA genotyping data 836 will not be released due to ethical restrictions.      Table S1. multiGO links to GO biology terms/reactome pathway enrichment analysis for DGE results. See also Figures 5-7   884 Table S4. Cell viability for each mock-control ALI-HNEC used for 10X Chromium preparation. See also Figures 1-3, S1-