Spatiotemporal quantification of cell dynamics in the lung following influenza virus infection

Abstract. Lung injury caused by influenza virus infection is widespread. Understanding lung damage and repair progression post infection requires quantitative spatiotemporal information on various cell types mapping into the tissue structure. Based on high content images acquired from an automatic slide scanner, we have developed algorithms to quantify cell infiltration in the lung, loss and recovery of Clara cells in the damaged bronchioles and alveolar type II cells (AT2s) in the damaged alveolar areas, and induction of pro-surfactant protein C (pro-SPC)-expressing bronchiolar epithelial cells (SBECs). These quantitative analyses reveal: prolonged immune cell infiltration into the lung that persisted long after the influenza virus was cleared and paralleled with Clara cell recovery; more rapid loss and recovery of Clara cells as compared to AT2s; and two stages of SBECs from Scgb1a1+ to Scgb1a1−. These results provide evidence supporting a new mechanism of alveolar repair where Clara cells give rise to AT2s through the SBEC intermediates and shed light on the understanding of the lung damage and repair process. The approach and algorithms in quantifying cell-level changes in the tissue context (cell-based tissue informatics) to gain mechanistic insights into the damage and repair process can be expanded and adapted in studying other disease models.


Introduction
The pulmonary airway successively branches from trachea to bronchi; the bronchi continue to divide within the lung to bronchioles and to alveoli. In a mouse, bronchioles are covered with Clara cells whereas alveoli are covered with alveolar type I cells (AT1s) interspersed with alveolar type II cells (AT2s). The lung is constantly exposed to environmental assaults, such as chemicals and infectious agents, leading to significant tissue damages. Studies with bleomycin-and naphthalene-induced lung injury have shown that Clara cells, which express secretoglobin family 1A member 1 [Scgb1a1 or Clara cell secretory protein (CCSP)], can repair damaged bronchiolar epithelium by selfrenewal and give rise to bronchiolar ciliated cells. [1][2][3] AT2s, which express pro-surfactant protein C (pro-SPC), can selfrenew as well as differentiate into AT1s in response to alveolar damage. [4][5][6] Recently, we have established a lung damage and repair model in mice with the influenza virus infection. We showed that Clara cells can give rise to AT2s during the repair of damaged alveoli through a novel intermediate cell type that exhibits the Clara cell morphology, but expresses the AT2 maker pro-SPC in the bronchioles (Zheng et al. submitted). We referred these pro-SPC þ bronchiolar epithelial cells as SBECs and further showed that SBECs progress from Scgb1a1 þ to Scgb1a1 − stage during the alveolar damage repair. The Scgb1a1 þ to Scgb1a1 − transition represents shutting off bronchiolar marker and turning on alveolar marker in the bronchiolar progenitor cells, which is the key phenotype of bronchiolar to alveolar trans-differentiation.
Studying lung damage and repair requires quantification of tissue damage and repair. Various methods have been developed to characterize specific lung tissue damages. Based on the sickness scores of inflammation, edema, or fibrosis, a semi-quantitative grading system has been used to classify different levels of lung damage. [7][8][9][10][11][12][13][14] The average size of the alveolar sac has been used to quantify hyperoxia-induced lung injury, where alveolar networks are destroyed. 15,16 The tissue thickness has been measured from hematoxylin and eosin (H&E) stained images to quantify bleomycin-induced lung injury leading to fibrosis. 17 Various automated image analysis methods have been reported. Aerated lung volume has been calculated from CT images in the study of bleomycin-induced lung fibrosis and emphysema. 18 Cell infiltration in the lung has been quantified by texture analysis in study of allergic asthmas, which distinguishes infiltrated parenchyma from the rest of the structures. 19 These automated methods were validated with quantification methods, such as a histology score with some degrees of consistency and operational ease; however, these quantification methods are only at the tissue level. Tissue level quantification alone without cellular specificity is not sufficient to reveal mechanistic insights from the cellular changes in the tissue damage and repair process. Other studies applying image analysis with specific molecular markers of the cells quantify cell numbers or molecular information, but not in the spatiotemporal context of the damage and repair process in the tissue environment. [20][21][22][23] Neither of the quantification approaches is able to provide an understanding of the cellular mechanism of lung damage and repair.
Here, we aim to quantify large numbers of various types of lung cells spatially mapping into the entire lung tissue structure over the course of lung damage and repair post infection by acquiring high content images using an automatic digital slide scanner. [24][25][26][27] We have developed automated algorithms, using high content images obtained from the automatic slide scanner, to quantify cellular dynamics in lung damage and repair following the influenza virus infection to gain mechanistic insights into the lung damage and repair process. The cell-specific tissue informatics approach and the algorithms developed can also be expanded and adapted to study the disease, injury, and repair processes in other animal models of various organs.

Mice and Tissue Sample Preparation
Female C57BL/6 mice were purchased from the Centre for Animal Resources, Singapore and housed in specific pathogen-free BSL2 facilities at the National University of Singapore (NUS). Mice were anesthetized with ketamine (100 mg∕kg body weight) and infected with a sub-lethal dose (100 PFU/mouse) of the influenza virus A/Puerto Rico/ 8/34 (H1N1) by intratracheal instillation. Five mice were sacrificed at each time point: before infection (control) and 1,3,5,7,9,11,13,15,17,19,21, and 28 days post infection (dpi). Lung tissues were collected. The left lobes were fixed in 10% neutral buffered formalin solution (Sigma-Aldrich, HT501128-4L) overnight and processed with Tissue Processor (Leica, TP1020, Shanghai, China) and embedded in paraffin blocks. For each lobe, 15 transverse sections were cut with microtome (Leica, RM2165, Nussloch, Germany) from the middle part of the lobe and with 50 μm in between sections. Sections were mounted on polylysine-coated slides (Thermal Fisher Scientific, J2800AMNZ). Ten sections were used for H&E and five sections for immunofluorescent staining. The right lobes were used for lavage, measuring virus titer and other analysis. All animal protocols were approved by the Institutional Animal Care and Use Committee at NUS and Massachusetts Institute of Technology.

Histological Staining
Paraffin sections were de-waxed in Histo-Clear solutions twice (National Diagnostics, HS-200) and rehydrated first in absolute ethanol three times, and then once each in 90% ethanol, 70% ethanol, and 50% ethanol. H&E staining was processed according to standard protocol. For immunofluorescent staining, antigen retrieval was performed by incubating the lung section with proteinase K solution (Sigma-Aldrich, 20 mg∕ml, in 50 mM Tris-Cl, 1 mM EDTA, pH 8.0) at 37°C for 30 min. Sections were then blocked in blocking buffer (3% BSA, 0.2% Triton X-100 in PBS) for 1 hr. Polyclonal rabbit anti-CCSP (also known as Scgb1a1) antibody (USbiology, C5828) was used in 1∶200 dilution, and goat anti-pro-SPC (Santa Cruz sc-7706) antibody was used in 1∶50 dilution. Incubation was performed at 4°C overnight in blocking buffer. Secondary Alexa Fluor 488-labeled donkey anti-rabbit antibody (Invitrogen, A21206) or Alexa Fluor 546-labeled donkey anti-goat antibody (Invitrogen, A11056) were used in 1∶200 dilution. Incubation was performed at 4°C for 1 h. Cover slips were mounted on stained section with antifade reagent containing DAPI (Invitrogen, S36939).

Image Acquisition
A high-resolution MIRAX MIDI system (Carl Zeiss, Jena, Germany) equipped with bright field and fluorescence illumination was used to scan the stained lung sections. Images were captured with Axiocam MR(m) (Carl Zeiss, Thronwood, New York), and converted to TIFF format with Miraxviewer software (Carl Zeiss). The original resolution of all images is 1.86 × 1.86 μm 2 per pixel. The images used for quantification were down-sampled eight times from the original images to reduce the computation time and the use of memory space.

Histology Analysis
Digital images of H&E stained lung samples from different time points post infection were evaluated blindly by three independent researchers. A semi-quantitative grading system was applied as previously reported, 9 in which subjective scores from 0 to 3 were given to reflect absent, mild, moderate, and extensive cell infiltration, respectively. Scores from the three researchers were averaged to obtain the average score of cell infiltration at each time point post infection.

Feature Extraction and Quantification
Automated algorithms were developed to compute four indices for quantitative evaluation of lung damage and repair based on image analysis. All image processing and computation algorithms were implemented with Matlab with the image processing toolbox (The Math Works Inc., Natick, Massachusetts). The Matlab codes are available on request.

Infiltration index
The infiltration index was defined as the percentage of infiltrated areas to total alveolar area in H&E stained lung section images. The infiltration index was calculated by infiltration algorithm as follows. The original H&E stained RGB image was first converted to greyscale. The image pyramid of greyscale image was then created by subsampling the image by a factor of 1, 2, 4, and 8, respectively, along each spatial dimension. 30 A threshold was manually adjusted and applied to identify low intensity pixels in each level of the image pyramid. A pixel was identified within the infiltrated area if it was selected as low intensity pixels in all levels of the image pyramid. The mask of the infiltrated area was generated by dilating all pixels identified within the infiltrated area to connect small regions.
The dilation was performed with a "disk" operator of 8 pixels in size. The total alveolar area was segmented from the original image by excluding large empty blood vessels and bronchioles with a threshold of 100 pixels in size. The percentage of infiltrated area to total alveolar area was then computed as the infiltration index and used as an index of general lung damage and repair at tissue level ( Fig. 1).

Clara cell coverage index
The Clara cell coverage index reflected the coverage density of Clara cells on the interior wall of bronchioles. Bronchioles were automatically recognized from immunofluorescent images as empty spaces with no nuclei inside and surrounded by Scgb1a1-expressing Clara cells [ Fig. 2(a)]. Nuclei were segmented from the original image by manual thresholding DAPI fluorescence intensities [ Fig. 2(b)]. Centers of each nucleus were located and denoted as nodes. A delaunay triangulation diagram was generated by introducing all nodes to connect neighboring nodes while none was inside the circumcircle of any triangle formed 31 [ Fig. 2(c)]. Nodes at the boundary of bronchioles (bronchiolar epithelium) tended to have longer average edge length than others, and nonboundary nodes with shorter average edge length were neglected accordingly. Boundary nodes belonging to the same bronchiole were then clustered together based on their connectivity in the original triangulation. Centers of each boundary node cluster were located and region growing 32 was performed from the centers to determine the exact boundaries of potential bronchioles [ Fig. 2(d)]. To further recognize bronchioles from all candidates, Clara cells were identified from the original image by manual thresholding Scgb1a1 fluorescence intensities. The Clara cell coverage ratios around the boundary area were calculated for each candidate, and a minimum threshold (1%) was applied to exclude false ones with completely no Clara cells surrounded [ Fig. 2(e)]. The remaining true bronchioles were then clustered into the high-Clara-cell-coverage ones (HC) and low-Clara-cellcoverage ones (LC) using the Otsu method 33 based on the Clara cells coverage ratio calculated. The Clara cell coverage index was defined as the ratio of HC frequency to LC frequency and used for the evaluation of the lung bronchiolar damage and repair.

Relative density of AT2
The two-dimension density of AT2s in lung tissue sections changes over the course of damage and repair post infection. Since AT2s were depleted in damaged areas, they were identified in healthy areas only. It was observed from the image that healthy lung tissue expressed dim auto-fluorescence in pro-SPC channel, while damaged areas emitted much brighter auto-fluorescence  Fig. 2(g)]. Accordingly, a k-means clustering algorithm was performed to group the pixels into four categories based on pro-SPC fluorescent intensities. 34 The brightest cluster with the highest pixel intensities represented SPC expressing cells, including AT2s in healthy areas and SBECs in damaged areas. The second and third brightest clusters referred to areas with bright autofluorescence and dim auto-fluorescence signals, respectively. The last cluster with minimum fluorescence intensity represented empty spaces in the lung tissue such as bronchioles, blood vessels, and alveolar sacs. To create the mask of damaged areas, we dilated the pixels in the group representing areas emitting bright auto-fluorescence with a "disk" operator of 3 pixels in size, so that most of them were connected to form blocks of areas. Then the small areas with less than 100 pixels were ignored and the remaining were defined as masks of damaged areas [ Fig. 2(h)]. Damaged areas were then subtracted from the total tissue section to form the mask of healthy areas. Pixels from the first cluster were then selected and overlapped with the mask of healthy areas to create mask of AT2s. The watershed method 35 was performed on the AT2 cell mask to separate cojoint cells and the number of AT2s was counted [ Fig. 2(f)]. The relative AT2 density was estimated by computing the number of AT2s per unit tissue area, and used as an indicator for lung alveolar damage and repair.

Frequency of SBECs
The frequency of the SBECs in bronchioles was defined as the percentage of SBECs-containing bronchioles to total bronchioles. Bronchioles were identified as described in Sec. 2.6.2. Within these bronchioles, the presence of Scgb1a1 þ and Scgb1a1 − SBEC were determined by checking the colocalization of Scgb1a1 and pro-SPC signals (at least 1%) at bronchioles in the damaged areas, which were segmented as described in Sec. 2.6.3 [ Fig. 2(i)]. The frequencies of Scgb1a1 þ and Scgb1a1 − SBEC-containing bronchioles were quantified by calculating the percentage of Scgb1a1 þ and Scgb1a1 − SBEC-containing bronchioles to total bronchioles, respectively [ Fig. 2(i)].

Statistical Analysis
Pearson correlation coefficients and p-values have been calculated to measure the consistency between: (1) cell infiltration kinetics obtained by automated algorithm and manual grading; (2) cell infiltration kinetics and Clara cell recovery kinetics from 5 to 14 dpi, and (3) Clara cell recovery kinetics and SBEC induction kinetics from 5 to 14 dpi. All statistics were computed with Matlab (The Math Works Inc., Natick, Massachusetts). The Matlab codes are available on request.

Quantification of Cell Infiltration in the Infected Lung
The influenza virus infection induces immune responses in the infected mice, leading to immune cell infiltration into the lung. We first quantify the cell infiltration as a frame of reference for the damage and repair process. The infiltrating cells appeared as densely packed nuclei in H&E staining [ Fig. 3(a)]. To quantify the cell infiltration in the infected lung, we developed an infiltration algorithm by segmenting dark infiltrated areas from the healthy areas in images of H&E stained lung sections, and computing the infiltration index by averaging the percentage of infiltrated area to total alveolar area in each lung sections (Fig. 1). By calculating the infiltration index over the 28 days following infection, we obtained a dynamic picture of cell infiltration [ Fig. 3(b)]. The infiltration was clearly detected 5 to 7 days post infection (dpi), corresponding to the length of time required to induce a T cell response. 36 The level of infiltration increased continuously and peaked around 14 dpi. By 21 and 28 dpi, most of the infiltrating cells were gone, corresponding to recovery from infection. The levels of cell infiltration as calculated by the automated algorithm matched well with those obtained from a semi-quantitative grading system, 9 in terms of overall kinetics (Pearson correlation coefficient equals to 0.964, p ¼ 1 × 10 −7 ) and peak days [ Fig. 3(c)]. Several imaging-based quantitative analyses of H&E stained lung sections have previously been reported. In the hyperoxiainduced lung injury where the network of alveolar sacs was destroyed, the quantification was based on the size of the alveolar sac. 15,16 In bleomycin-induced lung injury the quantification was based on the tissue thickness because of lung fibrosis. 17 Following the influenza virus infection, cell infiltration into the lung is prominent; therefore we used the levels of cell infiltration as a measure of the lung damage and repair. Compared to the semi-quantitative grading, 9 our automated algorithm has similar dynamic range and accuracy.

Quantification of Clara Cells, AT2s, and SBECs
In healthy lung, Scgb1a1-expressing Clara cells cover the bronchiolar epithelium whereas pro-SPC-expressing AT2s are dispersed in the alveolar epithelium [ Fig. 4(a) and 4(b)]. Influenza viruses infect Clara cells and AT2s directly 37 (Zheng et al. submitted), leading to gradual loss of both cell types. Loss of Clara cells was observed as early as 1 dpi (data not shown) and reached the peak level 5 dpi when Clara cells sloughed off from bronchiolar epithelium We have developed an automated algorithm to quantify the loss and recovery of Clara cells using immunofluorescent images of lung sections. Bronchioles were first identified in the images as closed empty structures with some Scgb1a1 signals (at least 1%) at boundaries (bronchiolar epithelium), and the coverage of Clara cell was then computed by calculating the ratio of bronchiolar epithelium length covered by Scgb1a1 þ cells to total bronchiolar epithelium length for each bronchiole [ Fig. 2(a) to 2(e)]. The bronchioles were divided into the high-Clara-cell-coverage (HC) and low-Clara-cell-coverage (LC) bronchioles; and the ratio of HC to LC bronchioles was defined as the Clara cell coverage index. Calculation of the Clara cell coverage index over the 21 dpi revealed that Clara cells were lost from bronchioles immediately after infection, reached the peak loss 5 dpi, and then gradually recovered to normal level by 14 dpi [Fig. 5(a)]. Similarly, we quantified the loss and recovery of AT2s using immunofluorescent images of lung sections. Healthy areas in the lung without apparent cell infiltration were first distinguished from damaged areas by k-means clustering-based segmentation. 34 AT2s were then identified as pro-SPC-expressing cells in the healthy regions of the lung. The relative density of AT2s was calculated as the number of AT2s per unit tissue section area. Quantification of AT2 density over 21 days post infection showed that AT2s were gradually lost after infection and the most severe loss was around 7 dpi [ Fig. 5(b)]. By 21 dpi, most of AT2s were regenerated in the damaged parenchyma.
SBECs were first detected around 7 dpi and became abundant by 14 dpi [ Fig. 4(g) and 4(h)]. Some SBECs were Scgb1a1 þ whereas others were Scgb1a1 − (Fig. 6). We have quantified the frequencies of Scgb1a1 þ and Scgb1a1 − SBECs based on immunofluorescent images. Damaged areas in the lung with severe cell infiltration were first identified using k-means clustering-based segmentation. Bronchioles in the damaged areas were then identified using the same way in the calculation of the Clara cell coverage index. SBEC-containing bronchioles were identified as bronchioles containing pro-SPC signals in the damaged areas. Frequencies of Scgb1a1 þ and Scgb1a1 − SBEC-containing bronchioles were quantified by examining the colocalization of Scgb1a1 and pro-SPC signals in the identified bronchioles in the damaged areas. As shown in Fig. 5(c), the Scgb1a1 þ SBECs first appeared around 7 dpi, peaked around 14 dpi, and decreased by 21 dpi. Scgb1a1 − SBECs were not observed until 9 dpi and then followed the similar dynamic profile as Scgb1a1 þ SBECs, suggesting that Scgb1a1 − SBECs come from Scgb1a1 þ SBECs.

Insights into Lung Tissue Damage and Repair
We compared the kinetics of virus titers, cell infiltration, Clara cell and AT2 loss and recovery, and induction of SBECs in the lung during the 21 days following the influenza infection. For easy comparison, the values were normalized and plotted over the same time course (Fig. 7). First, cell infiltration into the lung continued even when most of viruses were cleared. Second, the  continued cell infiltration preceded the Clara cell recovery. Third, the kinetics of SBEC induction starting at 7 dpi lagged behind but paralleled Clara cell recovery starting at 5 dpi (Pearson correlation coefficient equals to 0.9587, and p ¼ 0.0487 between the kinetics of Clara cell recovery and Scgb1a1 þ SBEC induction from 5 to 14 dpi. The Pearson correlation coefficient equals to 0.9763, and p ¼ 0.0237 between the kinetics of Clara cell recovery and Scgb1a1 − SBEC induction from 5 to 14 dpi). Finally, in contrast to Clara cells and SBECs that peaked at 14 dpi, recovery of AT2s did not reach the peak level even 21 dpi. These results are consistent with the notion that Clara cells give rise to Scgb1a1 þ SBECs and then Scgb1a1 − SBECs, which then differentiate into AT2s during the repair of damaged alveolar epithelium (Zheng et al. submitted). The continued cell infiltration into the lung when viruses were already cleared suggests that infiltrating cells likely contribute to tissue damage. The parallel kinetics between cell infiltration and recovery of Clara cells also raises the intriguing possibility that infiltrating immune cells may be involved in lung tissue repair. Thus, our systematic and quantitative analyses shed new light on lung tissue damage and repair following influenza virus infection.
In this study, we report new algorithms for quantifying cell infiltration, Clara cells, AT2s, and SBECs during the influenza virus-induced lung damage and repair. The methods are based on images acquired with an automatic slide scanner, which can scan a large area, such as the entire transverse lung section, with high speed and adequate resolution. As a result, we could obtain a comprehensive picture over the entire tissue section and specific spatially related features. For example, two stages of SBECs could be accurately identified and quantified in bronchioles in the damaged lung areas. To obtain the same information with microscopy modalities would be impractical because large numbers of small areas of the lung have to be analyzed and then combined. Thus, our quantification algorithms take advantage of the images with high content from a slide scanner, and extract specific features over large areas of tissue for quantification. Furthermore, our analyses extract features, such as nucleus density and spatial distribution of molecular markers; the same features could be used for studying cell dynamics in lung injury induced by other agents with no or minor modifications of the existing methods. Our spatiotemporal approach of image analysis could be potentially applied in study of cellular mechanisms of other diseases.

Conclusions
We have developed algorithms for quantifying cell infiltration, loss and recovery of Clara cells and AT2s, and induction of SBECs in the lung following the influenza virus infection. The methods were designed to use high content images of large tissue sections, so as to obtain spatiotemporal information in the entire tissue. Our analyses reveal a dynamic change of various cell types during the influenza virus-induced lung damage and repair, and provide insights into the relationships of these cell types in the lung repair process. Prolonged immune response was observed long after the virus has been cleared in the lung, which implies possible contribution of immune cells in lung repair. Clara cells were observed to experience more rapid loss and recovery than AT2, which sets the prerequisite for Clara-to-AT2 differentiation. Scgb1a1 þ SBECs were observed to be induced later than Scgb1a1 − SBECs; and the kinetics of both phases lagged behind and parallel with Clara cell recovery, which strongly support the notion that SBECs are originated from Clara cells. The cell-specific tissue informatics approach and algorithms developed in this study provide a useful way in investigating other disease progression and recovery models.