SARS-CoV-2 infection produces chronic pulmonary epithelial and immune cell dysfunction with fibrosis in mice

A subset of individuals who recover from coronavirus disease 2019 (COVID-19) develop post-acute sequelae of SARS-CoV-2 (PASC), but the mechanistic basis of PASC-associated lung abnormalities suffers from a lack of longitudinal tissue samples. The mouse-adapted severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) strain MA10 produces an acute respiratory distress syndrome (ARDS) in mice similar to humans. To investigate PASC pathogenesis, studies of MA10-infected mice were extended from acute to clinical recovery phases. At 15 to 120 days post-virus clearance, pulmonary histologic findings included subpleural lesions composed of collagen, proliferative fibroblasts, and chronic inflammation, including tertiary lymphoid structures. Longitudinal spatial transcriptional profiling identified global reparative and fibrotic pathways dysregulated in diseased regions, similar to human COVID-19. Populations of alveolar intermediate cells, coupled with focal up-regulation of pro-fibrotic markers, were identified in persistently diseased regions. Early intervention with antiviral EIDD-2801 reduced chronic disease, and early anti-fibrotic agent (nintedanib) intervention modified early disease severity. This murine model provides opportunities to identify pathways associated with persistent SARS-CoV-2 pulmonary disease and test countermeasures to ameliorate PASC.

COVID-19 is generally characterized as biphasic, with an acute phase dominated by active SARS-CoV-2 infection and a postviral clearance phase dominated by host reparative and immunologic processes (10). Human autopsy samples highlight the lung disease manifestations in patients who succumbed to COVID-19 (11,12), with broad features of chronic active pneumonia (CAP), alveolar architectural destruction, dense cellularity, and pulmonary fibrosis (PF) with myofibroblast proliferation and collagen deposition (13)(14)(15)(16)(17)(18)(19). Survivors of previous emerging coronavirus infections reported severe postinfectious fibrotic lung sequelae long after virus clearance, and autopsy data suggest that similar late sequelae will  Picrosirius Red staining (bright pink-red) highlights collagen fibers. Scale bars, 1000 m (low magnification) and 100 m (×400 images). (H) Disease incidence scoring is shown for indicated time points: 0 = 0% of total area of examined section, 1 = less than 5%, 2 = 6 to 10%, 3 = 11 to 50%, 4 = 51 to 95%, and 5 = greater than 95%. Graphs represent individuals necropsied at each time point (C and H), with the average value for each treatment and error bars representing SEM. Mock- fig. S5, A and B). IHC was used to quantitate the kinetics of CD4 + and CD8 + T cells ( fig. S5, C and D). Increased CD4 + T cells appeared as early as 2 dpi, peaked at 7 to 15 dpi, and persisted through 120 dpi ( fig. S5A). CD8 + T cell accumulation peaked at 15 dpi and remained at lower frequencies through 120 dpi ( fig. S5, A and D). B220 + B cell accumulation was observed at 7 dpi and sustained thereafter. CD68 + macrophages were increased at 7 dpi and remained elevated at 120 dpi in dense cellular regions, whereas inducible nitric oxide synthasepositive M1 and Arginase + M2 macrophages peaked at 2 and 7 dpi, respectively, and remained elevated at lower frequencies thereafter, suggesting the involvement of multiple subsets of macrophages in inflammatory and reparative process with different kinetics.
Flow cytometry at 30 dpi revealed that total cells, CD45 + immune, and CD31 + endothelial cells were increased (fig. S5, E and F), consistent with IHC results (fig. S5, A and B). CD4 + T and CD19 + B cells were increased in infected mice ( fig. S5G), consistent with prolonged inflammatory immune responses in pulmonary fibrotic diseases (36). Within the monocyte/macrophage lineage, interstitial macrophages were elevated in infected mice at 30 dpi ( fig. S5H), consistent with a documented role that macrophages play in lung remodeling in PF (37).

Host transcriptional profiles are spatially and temporally altered in response to SARS-CoV-2 infection
GeoMx DSP was used to interrogate viral and mouse RNA expression in pulmonary lesions from a subset of mock versus infected 1-year-old mice at 2, 15, and 30 dpi ( Fig. 2A). GeoMx DSP allows for quantitative analyses of RNA transcripts within targeted regions of interest (ROIs) using barcoded antisense oligos hybridized to more than 19,000 host and viral transcripts on formalin-fixed paraffin-embedded tissue sections. Because SARS-CoV-2 MA10 primarily infects alveolar type II (AT2) cells and terminal bronchiolar secretory club cells (28), we focused on these two regions. At 2 dpi, alveolar ROIs were selected on the basis of the presence of SARS-CoV-2 MA10 RNA-positive cells. Bronchiolar ROIs at 2 dpi were selected to represent a range of SARS-CoV-2 MA10 infection. At later time points (15 and 30 dpi), the heterogeneity of alveolar lung infection/responses was sampled by obtaining ROIs from morphologically "diseased" regions with hypercellularity versus morphologically "intact" regions. All distal airways appeared normal at 15 and 30 dpi with ROIs defined as intact. After data quality control and normalization, 60 alveolar and 36 bronchiolar epithelial ROIs from SARS-CoV-2 MA10-infected or mock mice were sampled at acute (2 dpi) and late (15 and 30 dpi) time points (Fig. 2, B and C, and data file S3). Quantification of viral RNAs demonstrated clearance of viral RNAs from intact and diseased alveolar ROIs by 15 dpi (Fig. 2D), concordant with clearance of infectious virus (Fig. 1C).
Principal components analysis (PCA) of expressed genes identified time-, region-, and virus-dependent effects (Fig. 2, E and F). High-virus transcript-positive regions at 2 dpi clustered away from mock in both distal airway and alveolar regions. Furthermore, the alveolar ROIs selected from diseased regions of infected mice at 15 and 30 dpi separated from mock, suggesting persistent alterations of host transcriptomes (Fig. 2F). In contrast, the ROIs selected from intact airway and alveolar regions at 15 and 30 dpi clustered near mock healthy ROIs, suggesting recovery (Fig. 2, E and F).
Consistent with PCA, viral infection induced major changes in transcriptome profiles in infected mouse lungs (Fig. 3, A and B, and data files S2 to S4). In both alveoli and bronchioles, virally infected disease ROIs at 2 dpi were characterized by a broad and robust up-regulation of viral infection-induced acute inflammatory genes, represented by enrichment of IFN, IL-1, and nuclear factor B (NF-B) signaling pathways (Fig. 3, A to C, and data files S2 and S5). Up-regulated IFN-stimulated genes (ISGs) were consistent with ISGs reported in human cells after emerging coronavirus infection ( fig. S6, A to C, and data file S2) (38,39), suggesting that common antiviral pathways are activated in human and mouse pulmonary cells. As noted in other human lung cell types after coronavirus infection (40), ISG expression patterns in airway and alveolar ROIs were not identical, with some ISGs more robustly up-regulated in airway epithelial (Ifitm1, Lap3, and DSP pathway analyses revealed down-regulation of biological oxidation (bronchiolar and alveoli) and surfactant metabolism (alveoli) in infected mice at 2 dpi (Fig. 3, A and B), associated with loss of secretory club (Cyp2f2, Scgb1a1, and Scgb3a2) and AT2 (Sftpc, Lamp3, and Abca3) cell markers ( fig. S7A). RNA-ISH confirmed that SARS-CoV-2 MA10 RNA was predominantly localized in Scgb1a1 + secretory club cells and Sftpc + AT2 cells at 1 dpi in bronchioles and alveoli, respectively ( fig. S7, B and C). Loss of club (Scgb1a1) and AT2 (Sftpc) cell marker expression accompanied SARS-CoV-2 MA10 infection at 1 and 2 dpi, followed by restoration to baseline values by 15 dpi (fig. S7, A to E). The early loss of Scgb1a1 and surfactant protein genes is consistent with reported human COVID-19 autopsy data (41). Ciliated (Foxj1, Dnah5, and Rsph1) and AT1 (Ager, Hopx, and Cav1) cell markers were minimally affected by MA10 infection at any time point ( fig. S7, A to C and F).
The transcriptomic analyses also revealed notable temporal differences in gene expression in alveolar versus bronchiolar regions (Fig. 3, A to C). Consistent with the failure of diseased alveolar regions to return to histologically intact-like states, pathway analyses at 30 dpi revealed persistently up-regulated cellular senescence, hypoxia signaling, complement activation, p53 damage responses, signaling by the transforming growth factor- (TGF-) receptor complex, collagen formation, and extracellular matrix organization pathways, unique to diseased alveolar regions. The difference in postinfection recovery between the bronchiolar (rapid and complete) and alveolar regions (slow and incomplete) was also notable. Because apoptosis is reported to be less inflammatory than necrotic cell death (42), we investigated whether apoptotic cellular responses to infection were different between the two regions ( fig. S7G). At 2 dpi, SARS-CoV-2 MA10-infected bronchiolar epithelial cells expressed evidence of activated apoptotic pathways (cleaved caspase-3). In contrast, alveolar regions were characterized by widespread infection but little cleaved caspase-3. These differences in apoptotic activity are consistent with reports that murine airway epithelial cells are more primed for apoptosis than alveolar epithelial cells in basal states (43).
Using ADI/DATP/PATS signature genes reported from mouse acute lung injury (ALI) models (44)(45)(46), the SARS-CoV-2 MA10 DSP data demonstrated enrichment of ADI/DATP/PATS signatures in diseased alveolar ROIs at 2, 15, and 30 dpi (Fig. 4A). The ADI/ DATP/PATS signature genes were categorized into three expression clusters ( Fig. 4B and data file S2). The first cluster (Cdkn1a/F3/Timp1) was enriched in diseased ROIs at 2 dpi and decreased after 15 dpi, suggesting that these genes may play a role in AT2 cell transdifferentiation into ADI/DATP/PATS cells. The second cluster (Krt8/Cxcl16/Cstb) exhibited increased expression at 2 to 30 dpi. The third gene cluster (Clu/Eef1a1), including a variety of ribosomal protein genes, exhibited increased expression at 15 dpi and later. To further characterize the relationships between ADI/DATP/ PATS cells and disease, combined RNA-ISH and DSP analyses of reported transitional ADI/DATP/PATS cell markers (Cdkn1a and Krt8) (47,48) were serially performed after infection (Fig. 4, C and D). DSP data demonstrated that (i) Cdkn1a was up-regulated at 2 dpi and waned at late time points and (ii) Krt8 was also up-regulated at 2 dpi and in diseased ROIs at 15 dpi (Fig. 4C). Although Krt8 + /Cdkn1a + RNA-ISH signals were not detectable in alveolar regions in mock mice, increased numbers of dual Krt8 + and Cdkn1a + cells were observed by RNA-ISH in SARS-CoV-2-infected alveolar regions at 1 and 2 dpi (Fig. 4D), consistent with the DSP data (Fig. 4, B and C). The murine DSP gene signatures exhibited features similar to ADI/ DATP/PATS signature genes identified in human COVID-19 autopsy lungs (47) (Fig. 5A and data file S2), including p53, apoptosis, and hypoxia pathways (Fig. 3, B and C).

Persistent inflammation and fibrosis are a chronic manifestation in SARS-CoV-2 MA10-infected mice
In diseased alveolar ROIs at 15 and 30 dpi, multiple genes involved in adaptive immune signaling and extracellular matrix deposition were highly up-regulated, consistent with a wound repair/profibrotic environment (Fig. 6, A and B). Recent human COVID-19 autopsy and transplant lung studies identified abundant interstitial profibrotic monocyte-derived macrophages characterized by increased expression of SPP1, MMP9, and CTSZ (16,47,49). These macrophage features, coupled with up-regulated extracellular matrix remodeling (SPARC and CTSK) and M-CSF signaling genes (CSF1 and CSF1R), defined a profibrotic macrophage archetype in human IPF samples (50). Our DSP analyses identified features associated with this profibrotic macrophage archetype in diseased alveolar ROIs at 15 and 30 dpi, including increased Spp1, Sparc, and Csf1r expression (Fig. 6C). RNA-ISH confirmed a persistent increase in Spp1 expression in SARS-CoV-2 MA10infected mice after 7 dpi (Fig. 6, D and E). These chronic fibrotic manifestations were consistent with IHC and flow cytometry data demonstrating increased interstitial macrophage populations during chronic SARS-CoV-2 MA10 infection ( fig. S5H). In addition, adaptive immune cell signatures, such as immunoglobulin (Igha, Igkc, and J chain) and major histocompatibility complex II (H2-Ea, H2-Eb1, and H2-Ab1) genes, were up-regulated in diseased alveolar ROIs at 30 dpi (Fig. 6B), consistent with the accumulation of interstitial macrophages and CD19 + B cells observed by IHC and flow cytometry ( fig. S5, A, G, and H).

Mouse SARS-CoV-2 MA10 infection recapitulates features of human lungs from fatal COVID-19 cases
We next compared mouse and published human data to a human COVID-19 autopsy cohort. Analyses of human COVID-19 autopsy by DSP, histology scoring, and IHC revealed biological  and TP63/MUC5AC networks were enriched in some COVID-19 lungs, which are consistent with histopathologic IPF features that exhibit infiltration of fibrotic alveoli with airway basal cells and "honeycombing cysts" lined by mucus-producing ciliated epithelia (51,53). The absence of this finding in the mouse may reflect a dearth of basal cells in the bronchiolar region of mice or unknown preexisting lung disease in patients with COVID-19 (31,53,54).
EIDD-2801 reduces chronic pulmonary lesions, and nintedanib decreases peak fibrotic disease EIDD-2801 (molnupiravir) is a Food and Drug Administration (FDA)-approved direct-acting antiviral that rapidly clears SARS-CoV-2 infection in mice and humans (55,56). We treated infected 1-year-old female BALB/c mice with EIDD-2801 or vehicle twice daily from 12 hours after infection to 5 dpi and followed survivors through 30 dpi. As reported (55), EIDD-2801 administration reduced weight loss, mortality, virus titers, gross lung congestion, diffuse alveolar damage (DAD), and ALI during the acute phase of infection (Fig. 7, A to F). At 30 days, profibrotic disease prevalence was reduced compared to vehicle controls (Fig. 7, G and H).
Nintedanib is an FDA-approved antifibrotic therapeutic agent that prevents IPF progression in humans (57,58). Nintedanib inhibits platelet-derived growth factor, fibroblast growth factor, and vascular endothelial growth factor receptors and interferes with fibroblast proliferation, migration, differentiation, and secretion of extracellular matrices (59). Older BALB/c mice that received nintedanib continuously from 7 dpi showed no differences in weight loss or recovery compared to vehicle-treated mice through 30 dpi (Fig. 7I). Nintedanib treatment beginning at 7 dpi did not affect the clearance of infectious virus by 15 dpi (Fig. 7J). Nintedanib treatment, however, decreased gross tissue congestion scores, fibrotic prevalence scores, and collagen deposition, at 15 dpi compared to controls (Fig. 7, J to M). Serum nintedanib concentrations were confirmed by ultrahigh-performance liquid chromatography time-of-flight mass spectrometry (UHPLC-TOF MS) to be within a range previously reported in mice (Fig. 7N) (60).

DISCUSSION
SARS-CoV-2 infection causes ALI, ARDS, and PASC. Although PASC encompasses nonrespiratory sequelae, including cardiovascular and neurologic disease (61), pulmonary manifestations are especially common, including CAP and PF (62,63). CT scans reveal chronic COVID-19 pulmonary findings as evidenced by ground glass opacities (44%) and fibrosis (21%) after acute COVID-19 infection (64) and fibrotic-like changes (35%) 6 months after severe human COVID-19 pneumonia (65). Pathology studies of COVID-19 lungs obtained at autopsy reveal similar late findings, such as CAP and PF (53,66,67). Accordingly, we focused our studies of PASC in the SARS-CoV-2 MA10 mouse model of COVID-19 and specifically on the pulmonary features of PASC. Currently, our understanding of PASC and COVID-19-induced CAP and PF is poor, and countermeasures are limited due to the wide spectrum of potential disease pathophysiologies. Although better studied, these limitations are also observed in infections with other viral respiratory pathogens such as influenza (68)(69)(70)(71). Human and animal model data comparing the chronic  sequelae of a spectrum of respiratory viruses should, therefore, identify unique versus shared disease manifestations and mechanisms of disease that will contribute not only to the knowledge of COVID-19 PASC but also to virus infection-mediated sequelae of the lung in general (72)(73)(74).
Recently, a chronic (30-dpi) SARS-CoV-2 infection model was reported in immunosuppressed, humanized mice characterized by persistent virus replication and chronic inflammation with fibrotic markers, typical of rare infections seen in immunosuppressed humans who cannot clear virus (75). In contrast, we report a mouse model of long-term pulmonary sequelae of SARS-CoV-2 infection that persisted after virus clearance and was more characteristic of disease outcomes seen in the general patient population. In the SARS-CoV-2 MA10 model, surviving older mice cleared infection by 15 dpi but exhibited damaged pulmonary epithelia accompanied by secretion of a spectrum of proinflammatory and profibrotic cytokines often up-regulated in fibrotic disease in humans, including IL-1, tumor necrosis factor-, granulocyte-macrophage CSF, TGF-, 15 30 IL-33, and IL-17A (76). Similar to humans, surviving SARS-CoV-2infected mice developed heterogeneous, persistent pulmonary lesions of varying severity by 30 to 120 dpi (77-79), presenting with abnormally repairing AT2 cells, interstitial macrophage and lymphoid cell accumulation, myofibroblast proliferation, and interstitial collagen deposition, particularly in subpleural regions. Micro-CT detected heterogeneous subpleural opacities and fibrosis in surviving mice, similar to human studies (64). Although most acute cytokine concentrations returned to normal values by 30 dpi, DSP and RNA-ISH data revealed focally prolonged up-regulation of cytokine signaling, including TGF-, in subpleural fibrotic regions. Similar heterogeneous cellular and fibrotic features in subpleural regions are also evident in patients with late-stage COVID-19 (80).
SARS-CoV-2 MA10 infection principally caused acute loss of distal airway club cell (Scgb1a1) and alveoli AT2 cell (Sftpc) marker expression, phenotypes consistent with SARS-CoV-2 cellular tropism in humans (81). The expression of club and AT2 cell genes was variably restored by 15 dpi, as demonstrated by DSP and RNA-ISH data. We speculate that a key variable determining the ability of the alveolar region to repair, or not, reflects the capacity of surviving or residual AT2 cells to regenerate an intact alveolar epithelium. The failure of AT2 cells to replenish themselves or AT1 cells and repair alveolar surfaces in subpleural regions may reflect the intensity of SARS-CoV-2 infection. On the basis of data from COVID-19 autopsy lungs, an accumulation of replication-defective and proinflammatory (ADI/DATP/PATS) transitional cells emerges early after SARS-CoV-2 infection and may persist, associated with continued inflammation and failure of repair (47,48). Our longitudinal mouse model data support this notion as evidenced by the observation that ADI/ DATP/PATS cells were detected at 2 dpi and persisted through 30 dpi in diseased, but not morphologically intact, alveolar regions. These ADI/DATP/PATS cells were notable for up-regulation of pathways associated with senescence, Hif1, and proinflammatory cytokines such as IL-1, consistent with low cycling rates, a failure to replenish AT2 and AT1 cells, and a proinflammatory phenotype (45). However, as evidenced by the return of Sftpc expression by 15 dpi in intact alveolar regions, a fraction of the ADI/DATP/PATS cells likely regenerated mature Sftpc-expressing AT2 cells. Our longitudinal studies revealed that the gene expression profiles of ADI/ DATP/PATS cells are dynamic over the evolution of lung disease.
As reported in humans, CD4 + and CD8 + T cell populations increased in SARS-CoV-2-diseased areas of mouse lungs, and peripheral lymphoid aggregations were a feature of chronic disease. These features were consistent across all analyses, including IHC, DSP, and flow cytometry data. A notable macrophage feature, identified by DSP and flow cytometry data, was the expansion of the interstitial macrophage population, consistent with human data (16). The subpleural regions exhibited the most notable histologic evidence of immunologic cell recruitment and activation of adaptive immune, hypoxia, fibrotic, and extracellular matrix pathways in association with ADI/DATP/PATS cells.
Final clues to the etiology of the late-stage alveolar CAP/PF response emerged from comparisons to infection in bronchioles. The alveolar regions exhibited persistent CAP/PF disease, particularly in subpleural regions. This finding is consistent with proposed relationships between the maximal pulmonary mechanical stretch imposed on the subpleural region during tidal breathing and activation of stretch-induced fibrotic pathways, such as TGF--mediated signaling, during periods of injury (82,83). In contrast, despite similar infection, bronchioles repaired without evidence of organizing or fibrotic sequelae. Bronchioles may be protected from this adverse fate by tissue-specific ISG responses to control the duration or severity of infection. In this context, several ISGs, including Ifitm1 and Ifitm2, exhibited clear differences in tissue-specific expression or persistence through 30 dpi. Other possible relevant variables that may favor bronchiolar repair include (i) more "controlled" cell death, such as apoptosis, (ii) a less damaged basement membrane architecture, and (iii) inability of club cells to enter an intermediate, ADI/DATP/PATS cell equivalent.
Mouse models of acute and chronic viral disease are critical also for countermeasure development. Molnupiravir is one of three FDA-approved direct-acting antivirals that clear virus and reduce morbidity, mortality, and time to recovery (55,84). Early molnupiravir treatment attenuated chronic PASC in the SARS-CoV-2 MA10 mouse model. Although speculative, early direct-acting antiviral treatment may forestall chronic lung and other organ PASC manifestations. On the basis of preclinical studies of antifibrotic agents in reducing the severity of PF responses to chemical agents, we tested the concept that early intervention with an antifibrotic agent may reduce the severity of PF after SARS-CoV-2 infection (59). Nintedanib administered from 7 dpi blunted maximal fibrotic responses to virus at 15 dpi, supporting the concept that early intervention with antifibrotic agents may attenuate post-SARS-CoV-2 severe disease trajectories. This suggests that early administration of direct-acting antivirals or antifibrotic drugs may help reduce human PF, and combination therapies may further increase efficacy and prognoses. Additional studies of other antifibrotic candidates and host immune modulators will be important in continuing to develop PASC treatments. COVID-19 in mice and humans represents key findings that may prove translatable to other future emerging coronavirus disease pathologies. Moreover, comparative models of viral-induced chronic lung disease are needed to identify common and unique pathways associated with virus-induced CAP and are key for the development of new therapeutic options for treating PASC.
With respect for study limitations, the SARS-CoV-2 MA10 model was developed using the ancestral clinical SARS-CoV-2/USA-WA1 strain that was unable to use murine angiotensin-converting enzyme 2 (ACE2) as a receptor and was serially passed to acquire increased virulence. Recently, variant-of-concern (VOC) human strains have been isolated subsequent to the USA-WA1 strain, which contains spike receptor binding domain substitutions that permit direct infection of standard laboratory mice. However, similar to the parental strain of MA10, strain SARS-CoV-2 MA (85), these VOC isolates do not cause substantial acute disease in mice, and persistent postviral-phase disease has not been reported. Our mouse-adapted SARS-CoV model has been used to evaluate human neutralizing antibodies, antiviral drugs, and vaccines that have progressed through phase 3 human trials, resulting in FDA-approved products for human use (86)(87)(88)(89)(90)(91)(92). Here, we show that acute and chronic disease phases in SARS-CoV-2 MA10-infected mice strongly recapitulate the pulmonary pathology observed in patients with COVID-19 and provide an excellent model for studies of pathogenesis and selected countermeasures. Transgenic and vectored expressions of human ACE2 mouse models are also commonly used to understand SARS-CoV-2 pathogenesis (93)(94)(95), but studies investigating the long-term pulmonary effects of infection in these models have not been reported. In addition, this study was limited to the long-term consequences of SARS-CoV-2 MA10 infection only in female mice due to housing constraints for long-term studies within our biosafety level 3 (BSL3) facility. It remains unclear whether sex-related effects account for long-term disease progression and recovery after SARS-CoV-2 infection in mice.
The current study provides new data for modeling chronic SARS-CoV-2 and other respiratory viral pathogens. By extending the studies of SARS-CoV-2 MA10 sequelae in mice to 120 days in 1-year-old BALB/c mice, we observed that many chronic phenotypes first observed at 15 dpi were maintained for the entire 120-day observational period, extending from acute to ongoing to chronic COVID-19-defined disease classifications used in human populations (29). The observation that fibrotic pulmonary disease in young BALB/c and aged C57BL/6J mice peaked at 15 dpi, but waned by 30 dpi in many animals compared to aged BALB/c mice, suggests that multiple time points for extended time intervals will be required for informative studies across viral pathogens. In summary, the SARS-CoV-2 MA10 mouse model provides opportunities to longitudinally study the molecular mechanisms and pathways mediating long-term COVID-19 pulmonary sequelae as it relates to human PASC and to evaluate treatments. Although future studies will be needed to determine whether other chronic, extrapulmonary organ sequelae develop after acute infection, the current model supports high-priority research directions that include SARS-CoV-2 infection of transgenic lineage tracing reporter mice to define longitudinally the fates of infected club and AT2 cells, ADI/DATP/PATS cell transitions, mechanisms of cell death, and epithelial cell regeneration and repopulation after infection. This study also provides the foundation for understanding the role of sex, host genetics, and immunological contributions, through knockout and Collaborative Cross studies (96)(97)(98), in defining PASC outcomes. Additional studies should be feasible in this model to investigate the nonpulmonary sequelae of PASC, including cardiovascular, neurological, and behavior manifestations in mice. With respect to countermeasures, 1-year-long clinical trials are required to assess the therapeutic benefit for lung fibrosis, emphasizing the utility of the SARS-CoV-2 MA10 model to rapidly test agents that may counter the pulmonary CAP/PF effects of COVID-19 (57,99). Thus, the murine SARS-CoV-2 MA10 model permits longitudinal selection and validation of therapeutic targets, accelerated timelines, and controlled experimental settings for testing of additional therapeutic agents.

Study design
The goal of these studies was to determine the prolonged pulmonary manifestations and underlying mechanisms of SARS-CoV-2induced lung injury in mice, with intent to highlight host pathways as possible therapeutic targets to treat PASC in patients surviving COVID-19. Animal experiments were restricted to only female mice because of long-term housing constrains within our BSL3 facility. Animal experimental cohorts were designed to be sufficiently large for at least three to five mice to be under each experimental condition and time point, accounting for survival rate after SARS-CoV-2 MA10 infection. Animals were randomized into experimental groups with predetermined endpoints at the start of the experiment. Mice for experimental replicates were age-matched. Experimental endpoints were chosen on the basis of previous work understanding acute SARS-CoV-2 MA10 infection and time points chosen for other respiratory infections in mice. For DSP of whole-transcriptome analyses (WTAs), three representative animals per group (mock, 2 dpi, 15 dpi, and 30 dpi) were selected on the basis of limitations of glass slide occupancy of the GeoMx instrument. ROIs were selected by an unblinded veterinary pathologist to ensure that all tissue types of interest were included in the analysis. All animals and samples were included in downstream data analysis. RNA-ISH and IHC were used to validate the data obtained from DSP analyses. Representative pathology images were selected from animals that represent the mean score after blinded quantification and were taken by a veterinary pathologist. Findings from mouse DSP data were validated on samples from matched and independent samples from repeated experimental cohorts. Human autopsy samples were limited by the number of samples we were able to obtain and used to corroborate findings in mice.

Ethics and biosafety
The generation of SARS-CoV-2 MA10 was approved for use under BSL3 conditions by the University of North Carolina at Chapel Hill Institutional Biosafety Committee (UNC-CH IBC) and by a Potential Pandemic Pathogen Care and Oversight committee at the National Institute of Allergy and Infectious Diseases. All animal work was approved by the Institutional Animal Care and Use Committee (IACUC) at the University of North Carolina at Chapel Hill, according to guidelines outlined by the Association for the Assessment and Accreditation of Laboratory Animal Care and the U.S. Department of Agriculture under UNC IACUC protocol 21-122. All work was performed with approved standard operating procedures and safety conditions for SARS-CoV-2, including all virologic work that was performed in a high-containment BSL3 facility, and personnel wore powered air-purifying respirators and Tyvek suits and were double-gloved. Our institutional BSL3 facilities have been designed to conform to the safety requirements recommended by Biosafety in Microbiological and Biomedical Laboratories, the U.S. Department of Health and Human Services, the Public Health Service, the Centers for Disease Control and Prevention (CDC), and the National Institutes of Health (NIH). Laboratory safety plans have been submitted, and the facility has been approved for use by the UNC Department of Environmental Health and Safety and the CDC.

Viruses and cells
Serial in vivo passaging of parental SARS-CoV-2 MA virus (85) in mice led to the plaque purification of a passage 10 clonal isolate (SARS-CoV-2 MA10) (28). A large working stock of SARS-CoV-2 MA10 was generated by passaging the plaque-purified clonal isolate sequentially on Vero E6 cells at 37°C (passage 3, SARS-CoV-2 P3). SARS-CoV-2 MA10 P3 was used for all in vivo experiments. Vero E6 cells were cultured in Dulbecco's modified Eagle's medium (Gibco) with the addition of 5% FetalClone II serum (HyClone) and 1× antibiotic/antimycotic (Gibco). Working stock titers were determined using plaque assay by adding serially diluted virus to Vero E6 cell monolayers. After incubation, monolayers were overlayed with media containing 0.8% agarose. After 72 hours, Neutral Red dye was used to visualize plaques.

In vivo infection
BALB/c mice used in this study were purchased from Envigo (BALB/cAnNHsd; strain 047) and C57BL/6J mice were purchased from Jackson Laboratories (strain 000664). Because of BSL3 animal housing constraints magnified by the duration of these studies, only female mice were studied. For intranasal infection, mice were anesthetized using a mixture of ketamine and xylazine. SARS-CoV-2 MA10 (10 4 PFU) diluted in phosphate-buffered saline (PBS) was used for inoculation of 10-week-old BALB/c and 1-year-old C57BL/6J mice. One-year-old BALB/c mice were infected with 10 4 PFU of SARS-CoV-2 MA10. Weight loss and morbidity were monitored daily as clinical signs of disease, whereas lung function was assessed at indicated time points using whole-body plethysmography (WBP; DSI Buxco respiratory solutions, DSI Inc.). Lung function data were acquired as previously described (100) by allowing mice to acclimate in WBP chambers for 30 min with a data acquisition time of 5 min. Data were analyzed using FinePointe software.
At indicated harvest time points, randomly assigned animals were euthanized by an overdose of isoflurane, and samples for analyses of titer (caudal right lung lobe) and histopathology (left lung lobe) were collected. Animals recorded as "dead" on nonharvest days were either found dead in the cage or were approaching 70% of their starting body weight, which resembles the criteria for humane euthanasia defined by respective animal protocols. Viral titers in lungs were determined by plaque assay for which caudal right lung lobes were homogenized in 1 ml of PBS and glass beads, monolayers of Vero E6 cells inoculated, and 72 hours after incubation stained with Neutral Red dye for visualization of plaques as described above.

Disease incidence scoring
Profibrotic disease incidence was scored by a blinded veterinary pathologist using serial hematoxylin and eosin (H&E)-and Picrosirius Red-stained slides. Ordinal scoring was defined by the percent of total parenchyma affected on the sampled section: 0 = 0% of total parenchyma, 1 = less than 5%, 2 = 6 to 10%, 3 = 11 to 50%, 4 = 51 to 95%, and 5 = greater than 95%. Instances of rare and isolated alveolar septa with gentle fibrotic changes were excluded from scoring.

Chemokine and cytokine analysis
Chemokine and cytokine profiles of serum and lung samples were assessed using Immune Monitoring 48-plex mouse ProcartaPlex Panel kits (Invitrogen). Briefly, 50 l of either a 1:4 dilution of serum or 50 l of straight clarified lung homogenate was incubated with magnetic capture beads containing analyte-specific antibodies. After washing, 96-well plates containing samples and magnetic beads were incubated with detection antibodies and streptavidin conjugated to phycoerythrin. Results were collected using a MAGPIX machine (Luminex), and quantification was achieved by comparing to a standard curve; both were done in xPONENT software. Values below the limit of detection (LOD) were set to LOD, and hierarchical clustering heatmaps were generated with the Bioconductor R package ComplexHeatmap after scaling the values across samples.

Preparation of lung cell suspensions for flow cytometric analysis
Enzymatic digestion of lung tissue was performed by intratracheal instillation using a 20-gauge catheter of 1 ml of collagenase I (5 mg/ml; Worthington Biochemical Corp.) and deoxyribonuclease I (0.25 mg/ml; Sigma-Aldrich) prepared in RPMI 1640 media (Life Technologies) before instilling 0.5 ml of 1% (w/v) low-melting agarose (Amresco), similar to previous protocols (101). Lungs were then incubated at 37°C for 30 min. Lungs were then minced and triturated through a 5-ml syringe. Cell suspensions were then filtered through a 50-ml conical 100-m filter (Thermo Fisher Scientific) before red blood cell lysis and stained as previously described (101).

Multicolor flow cytometry
The prepared lung cells were suspended in about 1 ml of PBS buffer supplemented with 1.5% (w/v) bovine serum albumin (Sigma-Aldrich) and 2 mM EDTA (Sigma-Aldrich). The total cell count was determined by a hemocytometer with trypan blue (VWR). For each sample, 1.5 × 10 6 cells first underwent Fc receptor blockade with rat anti-mouse FcRIII/II receptor (CD16/32; BD Biosciences). After Fc receptor blocking for 5 min on ice, cells were surface-stained using antibodies listed in table S1 and as previously described (101). For intracellular staining, the cells underwent fixation and permeabilization with the Foxp3/Transcription Factor Staining Buffer Set (eBioscience). Fixed and permeabilized single-cell suspensions were subsequently stained with intracellular antibodies to characterize differences in specific populations. Neutrophils and macrophage subpopulations were identified through gating, as demonstrated in prior reports (101,102) and adapted from previously published methods (103).
Flow cytometry was performed using a Cytoflex flow cytometer (Beckman Coulter) and analyzed using CytExpert (Beckman Coulter) software. To determine the total number of a specific population in the lung, we first calculated the population's percentage with respect to the total live single-cell population. Next, we multiplied this percentage to the total cell count as determined by hemocytometer measurements to calculate the specific population's total number per mouse lung.

Specimen CT imaging
Phosphotungstic acid staining was performed to increase soft tissue conspicuity for specimen CT imaging. Lungs were inflated and fixed with 10% formalin at 20 cm of H 2 O pressure for 7 days. Samples were initially washed three times in 70% ethanol (EtOH) in 50-ml nonreactive tubes before staining. Each lung was then immersed in 0.3% (w/v) phosphotungstic acid hydrate (Sigma-Aldrich, P4006) in 70% EtOH for 7 days on an oscillating table. They were subsequently air-dried before imaging.
Specimen CT scanning of the dried lungs was performed on a Sanco CT 40 (ScanCo Medical AG). Imaging was performed at 70 kilovolt peak at 114 A of current and 200 ms of integration time. Images were reconstructed using a conebeam algorithm at 16 m of voxel size in a DICOM file format. Images were viewed with ImageJ and analyzed by blinded radiologists (M.K.S. and Y.Z.L.).

RNA-ISH, IHC, and quantification
For histopathological analyses on mouse lung tissue sections, left lung lobes were stored in 10% phosphate-buffered formalin for at least 7 days before transferring out of the BSL3 facility for further processing. Histopathological scoring was performed after tissue samples were embedded with paraffin, sectioned, and stained. IHC was performed on paraffin-embedded lung tissues that were sectioned at 5 m. IHC was carried out using the Leica Bond III Autostainer system. Slides were dewaxed in Bond Dewax solution (AR9222) and hydrated in Bond Wash solution (AR9590). Heat-induced antigen retrieval was performed for 20 min at 100°C in Bond-Epitope Retrieval solution 2 (pH 9.0) (AR9640). After pretreatment, slides were incubated with primary antibodies (table S1) for 1 hour followed with Novolink Polymer Detection System (Leica, RE7260-K) secondary antibody. Antibody detection with 3,3′-diaminobenzidine (DAB) was performed using the Bond Intense R detection system (DS9263). Stained slides were dehydrated and coverslipped with Cytoseal 60 (Thermo Fisher Scientific, 8310-4).
RNA-ISH was performed on paraffin-embedded 5-m tissue sections using the RNAscope Multiplex Fluorescent Assay v2 or RNAscope 2.5 HD Reagent Kit according to the manufacturer's instructions (Advanced Cell Diagnostics). Briefly, tissue sections were deparaffinized with xylene and 100% EtOH twice for 5 and 1 min, respectively, incubated with hydrogen peroxide for 10 min and in boiling Target Retrieval Reagent (Advanced Cell Diagnostics) for 15 min, and then incubated with Protease Plus (Advanced Cell Diagnostics) for 15 min at 40°C. Slides were hybridized with custom probes at 40°C for 2 hours, and signals were amplified according to the manufacturer's instructions.
Stained mouse tissue sections were scanned and digitized using an Olympus VS200 slide scanner. Images were imported into Visiopharm software (version 2020.09.0.8195) for quantification. Lung tissue and probe signals for targeted genes detected by RNA-ISH were quantified using a customized analysis protocol package to (i) detect lung tissue using a decision forest classifier and (ii) detect the probe signal based on the intensity of the signal in the channel corresponding to the relevant probe. The same methodology was applied to quantify CD4 + and CD8 + cells identified by IHC. Positive signals for CD4 + cells were determined using a contrast of red-blue channels at a determined threshold to exclude background; similarly, CD8 + cells were determined using a contrast of green-blue channels. All slides were analyzed under the same conditions. Results were expressed as the area of the probe relative to total lung tissue area.
Paraffin-embedded mouse and human tissue sections (5 m) were used for fluorescent IHC staining. According to the previously described protocol (104), sections were baked at 60°C for 2 to 4 hours, followed by a deparaffinization step including xylene and graded EtOH. Antigen retrieval was achieved after rehydration by boiling slides in 0.1 M sodium citrate (pH 6.0) in a microwave. Slides were allowed to cool down and rinsed with distilled water before quenching of endogenous peroxidase was performed with 0.5% hydrogen peroxide in methanol for 15 min. After a PBS wash, slides were blocked with 4% normal donkey serum for 60 min at room temperature, followed by incubation with primary antibodies (table S1), where antibodies were diluted in 4% normal donkey serum in PBS with 0.1% Tween 20 (PBST), at 4°C overnight. Isotype control (species-matched gamma globulin) was diluted in the same manner as the primary antibody. Slides were incubated for 60 min at room temperature with secondary antibodies after being washed in PBST. Reduction of background staining was achieved by utilization of the Vector TrueVIEW Autofluorescence Quenching Kit (Vector Laboratories). Tissue sections were covered in glass coverslips by adding ProLong Gold Antifade Reagent with 4′,6-diamidino-2phenylindole (DAPI; Invitrogen). Stained human tissue sections were scanned and digitized using an Olympus VS200 slide scanner.

GeoMx DSP
Five-micrometer-thick formalin-fixed, paraffin-embedded sections were prepared using the RNAscope and DSP combined slide prep protocol from NanoString Technologies. Before imaging, mouse tissue morphology was visualized by IHC for CD45 and by RNAscope for SARS-CoV-2 RNA, and DNA was visualized with 500 nM Syto83. Human tissue morphology was visualized by IHC for the immune cell marker CD45 and epithelial cell marker panCK/ Syto83 and for KRT5 (IHC) and SARS-CoV-2 (RNA) on serial sections. Mouse or Human Whole Transcriptome Atlas probes targeting more than 19,000 targets were hybridized, and slides were washed twice in fresh 2× SSC then loaded on the GeoMx Digital Spatial Profiler (DSP). Briefly, entire slides were imaged at ×20 magnification, and 6 to 10 ROIs were selected per sample. ROIs were chosen on the basis of serial H&E-stained sections and morphology markers (mouse: DNA/CD45 IHC/SARS-CoV-2 RNA; human: CD45/PanCK/Syto83 IHC and SARS-CoV-2 RNA/KRT5/ DAPI IHC) on serial sections by an unblinded veterinary pathologist (S.A.M.). The GeoMx then exposed ROIs to 385-nm light (ultraviolet), releasing the indexing oligos and collecting them with a microcapillary. Indexing oligos were then deposited in a 96-well plate for subsequent processing. The indexing oligos were dried down overnight and resuspended in 10 l of diethyl pyrocarbonatetreated water.
Sequencing libraries were generated by polymerase chain reaction (PCR) from the photo-released indexing oligos and ROI-specific Illumina adapter sequences, and unique i5 and i7 sample indices were added. Each PCR used 4 l of indexing oligos, 4 l of indexing primer mix, and 2 l of NanoString 5× PCR Master Mix. Thermocycling conditions were 37°C for 30 min, 50°C for 10 min, and 95°C for 3 min; 18 cycles of 95°C for 15 s, 65°C for 1 min, and 68°C for 30 s; and 68°C for 5 min. PCRs were pooled and purified twice using AMPure XP beads (Beckman Coulter, A63881) according to the manufacturer's protocol. Pooled libraries were sequenced at 2 × 27 base pairs and with the dual-indexing workflow on an Illumina NovaSeq.

Analysis of mouse GeoMx DSP data
For mouse samples, raw count, third quartile (Q3)-normalized count data of target genes from ROIs were provided by the vendor, which were used as input to downstream analyses (data file S3). Mouse Q3-normalized data were used for PCA using the R package ade4 and visualized using factoextra package. Raw count data were used for differential expression analysis using the Bioconductor R package variancePartition (105), with transformation of raw counts by voom method (106). The dream function from variancePartition allows fitting of mixed-effects models to account for ROIs obtained from the same animal and assay slides as random-effect factors. Differentially expressed genes (DEGs) were defined as genes that passed the filters of Benjamini-Hochberg-adjusted P value of <0.05 and absolute log 2 fold change of >1. Preranked gene set enrichment analysis was performed using the Bioconductor R package fgsea (107), with gene set collections obtained from Gene Ontology (GO) Biological Process (108), and Reactome pathways (109). Various gene lists of interests were curated manually from published literature, and human gene symbols from references were converted into homologous mouse genes using bioDBnet (https://biodbnet-abcc. ncifcrf.gov/). Plots and hierarchical clustering heatmaps were generated using the R package ggplot2 (110) and Complex-Heatmap (111).
For the human samples, WTA and COVID-19 spike-in gene targets were assayed. FASTQ data were first converted to digital counts conversion format. Probe outlier tests were performed on each set of negative probes (one set of negative probes for the WTA panel and one set for the COVID-19 spike-in panel). Specifically, for a given negative probe pool, the geometric mean of all counts (across all probes and all samples) was computed. A probe was identified as a low-count outlier if its probe-specific geometric mean divided by the grand geometric mean was less than the threshold of 0.1. From the remaining probes, the Rosner test was used to detect local outliers on a sample-specific case using the R package EnvStats (112) with parameters k equal to 20% of the number of negative probes and alpha equal to 0.01. A negative probe was considered a global outlier if it was found to be a local outlier in more than 20% of samples and was discarded from downstream analysis. For each panel pool, the negative probe geometric mean and geometric SD were computed. The sample-specific limit of quantification (LOQ) was estimated from these moments by multiplying the geometric mean by the geometric SD and then squaring that quantity. Gene targets in the COVID-19 spike-in, which contains multiple probes per target, were collapsed to a single floating point value using the geometric mean. After outlier filtering, the sequencing saturation for each sample was computed as one minus the number of deduplicated reads divided by the number of aligned reads. One sample yielded a sequencing saturation below the 0.67 cutoff (range of other samples: 85.9 to 96.8) and was removed. In addition, one sample had an LOQ of more than 2.7 SD from the mean in the WTA panel and 4.2 SD from the mean in the COVID-19 spike-in pool and was removed from the analysis. Filtering gene targets was also performed. If a gene target was below LOQ in more than 10% of samples, then it was filtered out. After the above probe, sample, and target filtering steps, the data matrix was normalized using the Q3 method (see above).
Preliminary analysis of the log 2 -transformed and scaled Q3normalized data identified a putative batch effect between two runs, as identified using the PCA in the R package FactoMineR. The following batch correction algorithm was used before downstream data analysis. We first ensured that the batching factor was not itself confounded with group (healthy or COVID-19) or region (alveolar, bronchiolar, and disorganized). This was done by creating a design matrix and checking for any linearly dependent terms using the core R package stats (113). No factors were correlated with batch using a correlation threshold of 0.3. Batch correction was performed for each gene target by modeling its log 2 Q3 expression (dependent variable) in a mixed-effects model that included a random intercept for the fixed portion and batch as a random effect with random intercept. Modeling was done in the R package lme4. For each model, the residuals of the model were extracted and converted back to the linear scale. These residuals were then multiplied by the model's estimated intercept (also linear scale) to shift the values to an intensity similar to the original Q3 data. To evaluate how well the above approach removed the batch effect, we regressed the first five PC scores against batch for both the Q3 and the batch-corrected data using a series of analyses of variance (ANOVAs). Of the five PC axes, only the first was associated with the batching factor (P < 4 × 10 −36 ; all others, P > 0.23) in the Q3 data. After correction, no axes were associated with batch (all P > 0.80).

Analysis of human GeoMx DSP data
For human samples, raw count and Q3 + batch-corrected count data of target genes from ROIs were provided by the vendor (table S2  and data file S6). Before downstream analysis, Q3 + batch-corrected data were log 2 -normalized. PCAs were performed on the top 1000 highly variable genes on the log-normalized data. Coexpression network analysis was performed on 11,556 expressed genes using the Weighted Gene Coexpression Network Analysis R package (114). Differential gene and network expression between groups was evaluated under a linear mixed-model approach accounting for multiple ROIs per donor using R package Ime4. Statistical significance of the estimates was evaluated with R package lmerTest (115), using the Satterthwaite's degrees of freedom method. Sets of DEGs were tested for overrepresentation of the genes in the databases (GO: Biological Process, GO: Molecular Function, GO: Cellular Components, Kyoto Encyclopedia of Genes and Genomes, and Reactome) using R package enrichR (116). For each network, genes were selected on the basis of the degree of correlation with the network eigengene. To cluster ROIs obtained from healthy and COVID-19 donors, hierarchical clustering was performed on the basis of the 50 most correlated network genes from each of the seven identified networks using ward.D2 agglomeration method. As a result, healthy ROIs were separated from COVID-19 ROIs, and COVID-19 ROIs were segregated into three subtypes, including COVID 1, COVID 2, and COVID 3. Various plots and heatmaps were generated using the R packages ggplot2 (110) and heatmap3 (117).

Histological scoring of human COVID-19 lung tissue
The H&E-stained ROIs were scored by a blinded pulmonary pathologist (S.D.G.) grading each section on a semiquantitative scale between zero and three, with zero representing a normal human lung section and three representing the most severe histologic change encountered in clinical practice. The features scored in each ROI are interstitial inflammation, airspace fibrin exudates (acute phase of lung injury), the fibroblastic/organizing phase of lung injury, and mature fibrosis. Human donor information can be found in table S3.

Human lung tissue and quantification of Picrosirius Red and SMA signals
Control lungs were obtained from lung transplant donors without any history of pulmonary disease, whose lungs were unsuitable for transplant because of size mismatch provided by the UNC Tissue Procurement and Cell Culture Core (institutional review boardapproved protocol #03-1396). COVID-19 autopsy lung tissue sections were obtained from R.  table S3. Early-and late-phase specimens were defined as autopsy tissues obtained ≤20 and >20 days after an onset of symptoms, respectively.
Stained areas of Picrosirius Red and SMA detected by IHC in the alveolar regions were quantitated using Fiji software. Alveolar regions were randomly selected and cropped from the field. Optimized threshold value was determined by adjusting the threshold accurately representing the original images. The optimized threshold values were applied to identify Picrosirius Red or SMA signals. The Picrosirius Red-or SMA-stained areas were measured and normalized to alveolar areas.
In vivo drug treatment EIDD-2801 (Emory Institute of Drug Design) was dissolved in a solution of 2.5% cremaphor (Sigma-Aldrich), 10% polyethylene glycol 400 (PEG 400; Fisher Chemical), and 87.5% molecular biology-grade water (HyClone) by bath sonication at 37°C for 10 min, as described previously (55). Drug solution was made, at a concentration of 62.5 mg/ml, fresh daily for a final dose of 250 mg/kg per mouse (500 mg/kg twice per day). Mice were dosed by oral gavage with 100 l of vehicle [2.5% cremaphor (Sigma-Aldrich), 10% PEG 400 (Fisher Chemical), and 87.5% molecular biology-grade water (HyClone)] or EIDD-2801 solution twice daily, beginning at 12 hours after infection, and were dosed every 12 hours until 120 hours after infection.
Nintedanib (MedChemExpress) suspension was made in molecular biology-grade water (HyClone) with 1% Tween 80 (Sigma-Aldrich), fresh daily at a concentration of 15 mg/ml for a final dose of 60 mg/kg per mouse (118,119). Mice were dosed once daily by oral gavage with either 100 l of nintedanib suspension or vehicle [molecular biology-grade water (HyClone) with 1% Tween 80 (Sigma-Aldrich)] starting at 7 dpi until final harvest at either 15 or 30 dpi. Mouse serum was harvested at indicated time points after nintedanib administration, inactivated for BSL3 removal with 0.05% Triton X-100 and heating at 56°C, and was analyzed using UHPLC-TOF MS. Samples were prepared by precipitating protein with acetonitrile (Sigma-Aldrich) containing diazepam (Cerilliant) as an internal standard. The supernatant was separated using a Flexar FX-20 UHPLC system (PerkinElmer) with a Kinetex C18 biphenyl column (2.6 m, 50 mm by 3 mm Phenomenex) at 45°C with 98% MS-grade water (Sigma-Aldrich), 10 mM ammonium acetate (Hagn Scientific), and 98% methanol (Sigma-Aldrich) and 0.1% formic acid (Hagn Scientific) gradient elution at a flow rate of 0.6 ml/min. The PerkinElmer Axion2 TOF mass spectrometer operated in positive-ion electrospray ionization mode was used to detect accurate mass spectra of nintedanib at 540.2605 [M + H]+. The method was linear from 1 to 500 ng/ml with a lower LOD of 1 ng/ml. The results for nintedanib concentration in mouse serum for this study were in agreement with the serum concentrations reported previously (60).

Statistical analysis
Raw data are available in data file S7. Wilcoxon rank sum test was used to test the difference in CD4 + or CD8 + T cells ( fig. S5, C and D), as well as Picrosirius Red-or SMA-stained areas (fig. S8, C and D), identified by IHC between two groups. Flow cytometry data were analyzed by Wilcoxon rank sum test ( fig. S5E) or ANOVA, followed by Sidak's multiple comparisons test ( fig. S5, F to H). The difference in DSP Q3-normalized counts for targeted genes in ROIs between each condition and time point was statistically tested using a linear mixed-effects model using the R package Ime4 (120), with condition and time point as fixed-effects factor and replicate mice as random-effects factor ( Fig. 4C and fig. S6, D and E). Statistical significance was evaluated with the R lmerTest package (115), using the Satterthwarte's degrees of freedom method. Multiple post hoc comparisons of subgroups were performed using the R multcomp package (121). P < 0.05 was considered statistically significant.