A tool for mapping microglial morphology, morphOMICs, reveals brain-region and sex-dependent phenotypes

Animals

C57BL/6J (no. 000664) and B6.129P-Cx3cr1tm1Litt/J (no. 005582, named here Cx3cr1GFP/−, only heterozygous were used) mice were purchased from The Jackson Laboratory. All animals were housed in the ISTA Preclinical Facility, under a 12-h light–dark cycle, with food and water provided ad libitum. Animals from both sexes were used. The number of animals used for each condition is detailed in Supplementary Table 5. All animal procedures were approved by the Bundesministerium für Wissenschaft, Forschung und Wirtschaft (bmwfw) Tierversuchsgesetz 2012, BGBI (I Nr. 114/2012, idF BGBI. I Nr. 31/2018 under the nos. 66.018/0005-WF/V/3b/2016, 66.018/0010-WF/V/3b/2017, 66.018/0025-WF/V/3b/2017, 66.018/0001_V/3b/2019 and 2020-0.272.234).

5xFAD and CK-p25 mice were obtained from the Tsai laboratory at Massachusetts Institute of Technology (MIT). All animal work was approved by the Committee for Animal Care of the Division of Comparative Medicine at MIT. 5xFAD mice (B6SJL-Tg(APPSwFlLon,PSEN1*M146L*L286V)6799Vas/Mmjax, stock no. 34840-JAX) were obtained from The Jackson Laboratory. CK-p25 mice69 were generated by breeding CaMKIIα promoter-tTA mice (CK controls) (B6;CBA-Tg(Camk2a-tTA)1Mmay/J, Jackson Laboratory, stock no. 003010) with tetO-CDK5R1/GFP mice (C57BL/6-Tg(tetO-CDK5R1/GFP)337Lht/J, Jackson Laboratory, stock no. 005706). CK-p25 mice were conceived and raised in the presence of doxycycline-containing food to repress p25 transgene expression. To induce p25 transgene expression, mice were fed a normal rodent diet. p25 transgene expression was induced in adult mice at the age of 3 months. For MorphOMICs, we compared CK-p25 brains upon drug withdrawal with our reference C57BL/6J adult microglia population. Notably, doxycycline withdrawal might affect the gut microbiome70, which can influence the microglial population in the brain and could cause some variability.

Mice were housed in groups of three to five on a standard 12-h light/12-h dark cycle, and all experiments were performed during the light cycle. Food and water were provided ad libitum.

Brain samples and analyzed brain regions

We analyzed brains of both sexes from C57BL/6J adult mice (8–12 weeks) exposed to 1×, 2× or 3× KXA (100 mg per kg body weight ketamine, MSD Animal Health, A137A01; 10 mg per kg body weight xylazine, AniMedica, 7630517; 3 mg per kg body weight acepromazine, VANA, 18F211) and recovered 3 d, 1 week and 2 weeks after 3× KXA71; Cx3cr1GFP/− mice at P7, P15 and P21; 5xFAD mice after 3 and 6 months; and CK-p25 mice 1, 2 and 6 weeks after doxycycline withdrawal. We focused on the following brain regions: the glomerular layer of the olfactory bulb (OB), cortical layer III–V of the frontal cortex (FC) and the primary somatosensory cortex (S1), the dentate gyrus of the hippocampus (DG), the substantia nigra (SN), the cochlear nucleus (CN) and the third lobe of the cerebellum (CB). The sagittal view of the brain sections analyzed (Figs. 15 and Extended Data Fig. 5a) was taken from the Allen Developing Mouse Brain Atlas-Sagittal and modified to show brain regions of interest72.

Ovariectomy

Adolescent C57BL/6J females at P20 were anesthetized with 5% isoflurane in 0.5 l min−1 O2 during the anesthesia induction and 2% isoflurane in 0.5 l min−1 O2 during the maintenance phase. Using an electric razor, the fur was shaved to expose the skin over the lumbar spine and the region was sterilized with 70% (vol/vol) ethanol. A midline incision of approximately 1 cm was made on the skin in the lower back, below the chest. The subcutaneous tissue was gently dissected to expose the muscular fascia, and the ovarian fat pad was identified under the muscular layer. The peritoneal cavity was cut with a 0.5-cm incision. The Fallopian tube was exposed, and the ovary identified and cut at the level of the oviduct. The blood vessels were cauterized to prevent bleeding. The remaining part of the Fallopian tube was placed back in the peritoneal cavity, and the muscular fascia was sutured. The same protocol was repeated for the contralateral ovary. At the end, the skin was sutured. The animals received metamizole (Sanofi Aventis, no. Ay005, 200 mg per kg body weight during surgery) and meloxicam (Boehringer-Ingelheim, no. KPOEH3R, 5 mg per kg body weight subcutaneously after surgery every 24 h for 3 consecutive days), 2 mg per kg body weight subcutaneously after surgery. Animals were euthanized at P60.

Transcardiac perfusion

For histological analysis, animals were quickly anesthetized with isoflurane (Zoetis, no. 6089373) and secured to the perfusion plate. The chest was open to expose the heart. The left ventricle was cannulated and the inferior vena cava cut. Animals were initially perfused with 30 ml of PBS with heparin (100 mg l−1; Sigma, no. H0878), followed by 30 ml of 4% (wt/vol) paraformaldehyde (Sigma, no. P6148) in PBS using a peristaltic pump (Behr, no. PLP 380, speed of 25 r.p.m.). Animals were decapitated, the brain explanted, fixed in 4% (wt/vol) paraformaldehyde for 30 min and post-fixed in 4% (wt/vol) PBS overnight (16 h). Then the tissues were washed in PBS and stored at 4 °C with 0.025% (wt/vol) sodium azide (VWR, no. 786-299). For cryoprotection, the tissue was transferred to 30% (wt/vol) sucrose (Sigma, no. 84097) in PBS and incubated overnight at 4 °C. To increase antibody permeability, the brain slices were frozen over dry ice and thawed at room temperature for three cycles.

Vibratome sections

Cryoprotected samples were embedded into 3% (wt/vol) agarose/PBS to obtain coronal brain sections. The brain was sliced in 100-µm coronal sections on a vibratome (Leica VT 1200S).

Immunofluorescence staining

The brain slices were incubated in blocking solution containing 1% (wt/vol) BSA (Sigma, A9418), 5% (vol/vol) Triton X-100 (Sigma, T8787), 0.5% (wt/vol) sodium azide (VWR, 786-299) and 10% (vol/vol) serum (either goat, Millipore, no. S26, or donkey, Millipore, no. S30) for 1 h at room temperature on a shaker. Afterwards, the samples were immunostained with primary antibodies diluted in antibody solution containing 1% (wt/vol) BSA, 5% (vol/vol) Triton X-100, 0.5% (vol/vol) sodium azide, 3% (vol/vol) goat or donkey serum, and incubated for 48 h on a shaker at room temperature. The following primary antibodies were used: rat α-CD68 (AbD Serotec, no. MCA1957, clone FA-11; 1807, 1:250 dilution), goat α-Iba1 (Abcam, ab5076, FR3288145-1; 1:250 dilution) and rabbit anti-Iba1 (GeneTex, no. GTX100042, no. 41556, 1 vol/vol 750). The slices were then washed three times with PBS and incubated for 2 h at room temperature on a shaker protected from light, with the secondary antibodies diluted in antibody solution. The secondary antibodies raised in goat or donkey were purchased from Thermo Fisher Scientific (Alexa Fluor 488 goat anti-rabbit IgG no. A11034, Alexa Fluor 647 goat anti-rat IgG 21247; 1:2,000 dilution). The slices were washed three times with PBS. The nuclei were labeled with Hoechst 33342 (Thermo Fisher Scientific, no. H3570; 1:5,000 dilution) diluted in PBS for 15 min. The slices were mounted on microscope glass slides (Assistant, no. 42406020) with coverslips (Menzel-Glaser no. 0) using an antifade solution (10% (vol/vol) Mowiol (Sigma no. 81381), 26% (vol/vol) glycerol (Sigma, no. G7757), 0.2 M Tris buffer (pH 8) and 2.5% (wt/vol) Dabco (Sigma, no. D27802)).

Confocal microscopy

Images were acquired with a Zeiss LSM880 upright Airy scan or with a Zeiss LSM700 upright using a Plan-Apochromat ×40 oil-immersion objective 1.4 NA. Then, 2 × 2 z-stack tail images were acquired with a resolution of 1,024 × 1,024 pixels.

Image processing

Confocal tile images were stitched using the software Imaris Stitcher 9.3.1.v. Then, the confocal images were loaded in Fiji 1.52e (http://imagej.net/Fiji). To remove the background, the rolling ball radius was set to 35 pixels, and images were filtered using a median 3D filter with x, y and z radii set at 3. Image stacks were exported as .tif files, converted to .ims files using the Imaris converter and imported into Imaris 8.4.2.v. (Bitplane Imaris).

Quantification of CD68 volume within cells

Surface renderings were generated on microglia and CD68 z-stacks using the surface-rendering module of Imaris 9.2.v Surfaces were generated with the surface detail set to 0.2 µm. To determine the CD68 surface within microglia, the surface–surface coloc plugin was used. This analysis was performed on the entire image. The total ratio of CD68 volume within microglial volume (CD68-to-microglial volume) was calculated per image. To compute the CD68 fold change, the total CD68-to-microglial volume from each condition (sex/time point) was scaled to the CD68-to-microglial volume ratio from the respective adult control brain region. CD68 fold change > 1 means an increase in CD68 volume, while CD68 < 1 means a decrease in CD68 volume. CD68 fold change = 1 denotes no change in CD68 volume.

Quantification of microglia density and statistical analysis

The spot-function plugin of Imaris 9.2.v was used to count the number of cells, that is, the soma of iba1-positive microglia within every confocal image. Microglial cell density was estimated as the total number of cells obtained in this way, divided by the size of the imaged sample in mm2.

Reconstruction of 3D-traced microglia

After filtering and background subtraction, images were imported in Imaris 9.2.v (Bitplane Imaris). Microglial processes were traced in 3D with the filament-tracing plugin. Because the filament-tracing plugin provides a semiautomated reconstruction, this eliminates the need for a user-blind approach for selecting representative microglia. New starting points were detected when the largest diameter was set to 12 µm and with seeding points of 1 µm. Disconnected segments were removed with a filtering smoothness of 0.6 µm. After the tracing, we manually removed cells that were sitting at the border of the image and were only partially traced so that these cells would not be analyzed. The generated skeleton images were converted from .ims format (Imaris) to .swc format73 by first obtaining the 3D positions (x, y and z) and the diameter of each traced microglial process using the ImarisReader toolbox for MATLAB (https://github.com/PeterBeemiller/ImarisReader/) and then exporting for format standardization using the NL Morphology Converter (http://neuroland.org/). Artifacts from the 3D reconstructions automatically failed to be converted into .swc format.

Analysis of morphometric features

Classic morphometric features were calculated from the .swc files using the functions Length (for total process length), N_branch (for number of branches), N_bifs (for number of branching points) and N_tips (for number of terminal points) from L-measure74 (http://cng.gmu.edu:8080/Lm/).

Sholl analysis

Sholl curves were calculated from the .swc files using the sholl_crossings function of the NeuroM Python toolkit (https://github.com/BlueBrain/NeuroM). In brief, concentric Sholl spheres centered on the soma of a given traced microglia are constructed with a given step size radius. The number of microglial processes that intersect each Sholl sphere are determined. This step is performed for each traced microglia in the data. From this, Sholl curves of a microglial population are then calculated as the average number of intersections across the population.

Topological morphology descriptor

A topological data analysis algorithm, the TMD, was used to extract topological phenotypes, called persistence barcodes, from 3D morphological structures (https://github.com/BlueBrain/TMD/;14). In brief, the 3D-reconstructed microglia is represented as a tree T rooted in its soma. The TMD summarizes this tree by calculating ‘persistence barcodes’, where each bar represents a persistent microglial process with respect to a filtering function, that is, the radial distance from the soma. Note that the persistence barcode that the TMD associates with T under this filtering function is invariant under rotations about the root and rigid translations of T in R3.

Each bar is described by two numbers: the radial distance, di, at which a process originates; and the distance, bi, when it merges with a larger, more persistent process or with the soma. A bar can be equivalently represented as a point (di, bi) in a ‘persistence diagram’. We could therefore convolve each point in the persistence diagram with a Gaussian kernel and discretize it to generate a matrix of pixel values, encoding the persistence diagram in a vector, called the ‘persistence image’.

Average and bootstrapped persistence images

To construct the ‘average persistence image’ of a given condition, all the persistence barcodes of microglia from the same condition are combined before Gaussian convolution and discretization are performed. We also constructed average persistence images by performing first the Gaussian convolution and discretization of individual microglia persistence barcodes before taking the pixel-wise average. This produced qualitatively similar results.

The bootstrapping method subsamples the microglial population within a given condition, thereby introducing variations around the average persistence image. Starting from the population of all microglia from the same condition, called the ‘starting population’ of size n (Supplementary Table 4), the persistence barcodes of a predefined number of unique microglia, called the ‘bootstrap size’, are combined to calculate the ‘bootstrapped persistence image’. We iterated this process a predefined number of times, nsamples, with replacement to obtain the ‘bootstrap sample’.

Subtraction images and topological morphology descriptor distance

The subtraction image is the pixel-wise difference between two given persistence images. From the subtraction image, the TMD distance can be computed as the sum of the absolute pixel-wise difference. For stability of the TMD distance, we refer the reader to Kanari et al.14.

Hierarchical clustering

Hierarchical clustering allowed us to find similarities between microglia across several conditions. Hierarchical clustering was done on the basis of the average persistence images. Clusters were then identified hierarchically using the average linkage criterion with the TMD distance metric and was implemented using cluster.hierarchy.linkage from SciPy v1.6.2 (https://www.scipy.org/). Dendrograms were generated using cluster.hierarchy.dendrogram to visualize the arrangement of the resulting cluster.

Dimensionality reductionUniform manifold approximation and projection

A fast, nonlinear dimensionality reduction algorithm, UMAP75, was applied to visualize the high-dimensional pixel space of bootstrapped persistence images using a 2D representation while preserving local and global structures in the bootstrap samples (https://github.com/lmcinnes/umap/)75. Given a bootstrap sample containing multiple conditions, a TMD distance matrix containing pairwise distances between bootstrapped persistence images in the bootstrap sample is calculated. Principal components are then obtained using a singular value decomposition of the TMD distance matrix. The first seven principal components, where the elbow in the singular values is located, were used as input to UMAP with n_neighbors = 50, min_dist = 1.0 and spread = 3.0. Note that we tested for a wide range of parameter values that did not qualitatively change any of the aforementioned observations (Extended Data Fig. 2f).

t-SNE

An alternative dimensionality reduction algorithm is t-SNE (https://github.com/DmitryUlyanov/Multicore-TSNE/), which finds a dimensionality-reduced representation where similar points are pulled closer together while dissimilar points are pushed farther apart with high probability. The first seven principal components were taken as an input to run t-SNE with perplexity = 50.

Pseudo-temporal ordering

The concept of morphological phenotypes as encoded in the persistence images can be likened to transcriptional phenotypes in single-cell RNA-sequencing studies. Bootstrapped persistence images, which encapsulate morphological phenotypes of microglial populations from similar conditions, are comparable. Furthermore, it is reasonable to assume that morphological changes in bootstrapped microglial populations from control-to-disease conditions occur with incremental differences in the persistence images. This conceptual similarity allowed us to use the pseudo-temporal trajectory-inference algorithms that are well used in the single-cell RNA-sequencing community to study the morphological progression during microglial development and degeneration.

Palantir

Palantir76 uses principles from graph theory and Markov processes to calculate the pseudo-time and the probability of a cell reaching each of the terminal conditions in the sample (https://github.com/dpeerlab/Palantir/). First, the principal components of the bootstrapped persistence images were obtained using palantir.utils.run_pca with n_components = 100 and use_hvg = false. The diffusion maps were then calculated from the PCA projections using palantir.utils.run_diffusion_maps with n_components = 10 and knn = 20 which outputs the Palantir pseudo-times. Harmony77 was then used to construct an augmented affinity matrix from the Palantir pseudo-times to connect the Palantir pseudo-times and construct a trajectory using a force-directed graph (https://github.com/dpeerlab/Harmony/).

Monocle

To corroborate the Palantir trajectories, an alternative pseudo-temporal trajectory-inference algorithm called Monocle was used. Monocle78 uses reversed graph embedding, which learns a principal graph that approximates a lower-dimensional manifold to construct a pseudo-time trajectory (https://github.com/cole-trapnell-lab/monocle3/)78. Similar to Palantir implementation, the principal components of the bootstrapped persistence images were first obtained using preprocess_cds with num_dim = 100. A 2D UMAP representation was then obtained using reduce_dimension with umap.metric = ‘manhattan’, umap.min_dist = 1.0, and clusters were identified using cluster_cells with cluster_method = ‘leiden’. Finally, the pseudo-temporal trajectory was then obtained using learn_graph with use_partition = FALSE and close_loop = FALSE.

Stable ranks analysis

An alternative representation of the persistence barcodes is through stable ranks79. Stable ranks are functional summaries of persistence that depend on pseudometrics to compare persistence barcodes. Given a pseudometric d, the stable rank \(\widehat _d}\left( X \right)\left( t \right)\) of a persistence barcode X is a function that assigns to t the number:

$$\widehat }_d\left( X \right)\left( t \right) = }\left\\left( Y \right)|d\left( \right) \le t} \right\}.$$

whereby ‘rank(Y)’ denotes the number of bars of the persistence barcode Y. The stable rank \(\widehat }_d\left( X \right)\left( t \right)\) associates to a persistence barcode a non-increasing and piece-wise constant function with values in [0, ∞). An important property is that this mapping is continuous with respect to the chosen pseudometric d and the Lp metric on the space ℳ of measurable functions.

A class of pseudometrics on persistence barcodes can be constructed from density functions79, which intuitively are used to vary the weight along the filtration scale parametrizing a barcode. With such pseudometrics, the stable rank is a bar count based on length of bars as scaled by the density. The standard stable rank is defined by a density function with constant value one. In this case, \(\widehat }_d\left( X \right)\left( t \right)\) is the number of bars in X with length greater than or equal to t, that is, all filtration scales are weighted equally.

Stable ranks can be used in place of persistence images in the MorphOMICs pipeline. Similarly to MorphOMICs, the persistence barcode X of a given microglia is calculated using the TMD algorithm. To obtain ‘bootstrapped standard stable ranks’, we combined the persistence barcodes of a predefined number of microglia and computed their standard stable ranks. Dimensionality reduction was then implemented similar to the methods above (‘Dimensionality reduction’).

Classification accuracy using stable ranks

To support and quantify the impact of bootstrapping on the regional segregation visualized in the reduced UMAP space (Fig. 1f), we performed a classification task for microglia morphologies represented by their standard stable rank and labeled by brain region. We used an SVM with a specific kernel based on stable ranks80,81 for the classification. For persistence barcodes X and Y, the stable rank kernel with respect to a pseudometric d is given by

$$K_d\left( \right) = \mathop \limits_0^\infty _d} \left( X \right)\left( t \right)\widehat }_d\left( Y \right)\left( t \right)}t.$$

where we used the pseudometric induced by the constant function with value one.

We performed pairwise classifications. For each pair of brain regions, we constructed a dataset consisting of 400 bootstrap samples, that is, 200 from each region and bootstrap sizes of 10, 20 or 50 (the results are reported separately for these three values). We randomly partitioned the dataset for cross-validation wherein 240 samples were used for SVM training (training set) and 160 samples for validation (test set). We reported the average accuracy over ten repeated cross-validations on the test set. The SVM was trained using the implementation in the Python library sklearn (https://scikit-learn.org/stable/) with default settings except for the usage of the stable rank kernel.

Bootstrapped morphometric features and bootstrapped Sholl curves

To understand whether classical morphology analysis pipelines are able to recapitulate the microglial dynamics recovered by MorphOMICs, a similar bootstrapping analysis was also done where we pooled a predefined number of microglia. Each morphometric quantity in the extended list enumerated in Supplementary Table 3 was then averaged to obtain a 27-dimensional vector, with each dimension corresponding to a morphometric feature, called the ‘bootstrapped morphometric features’. On the other hand, Sholl curves averaged across the pooled microglia to obtain the ‘bootstrapped Sholl curves’. Dimensionality reduction was then implemented similarly to the methods above (‘Dimensionality reduction’).

Mapping morphologies onto the reference atlas

We have generated a larger reference atlas with nsamples = 2,000 bootstrapped persistence images for each condition to construct the Palantir trajectories. Palantir coordinates (x, y) were rescaled to (0, 1). We took and filtered the bootstrapped persistence images keeping only the 500 most highly variable pixels across the images in all conditions. We used linear regression, one for each axis, to learn the mapping from the filtered bootstrapped persistence image to the rescaled Palantir coordinates. Given a novel condition, we generated the bootstrapped persistence images and filtered them with the 500 most highly variable pixels identified earlier. We used the trained regression model to infer the locations of each image in the reference atlas. Then, we calculated the mean position denoting the center of the inferred locations and indicated the spread using the standard deviation.

Statistics and reproducibility

Each experiment was repeated independently with similar results. Supplementary Table 4 provides the numbers of 3D traces obtained per condition, sex and brain region. Supplementary Table 5 describes the number of animals per condition, their sex and brain region.

In Extended Data Figs. 1a and 10a, statistical analysis was performed using scipy.stats (v1.6.2) and scikit-posthocs (v0.6.7). These morphometric features were first tested for normality using the Kolmogorov–Smirnov test (scipy.stats.kstest). After determining the non-normal distribution of the features, we performed non-parametric pairwise tests for independence between measurements from two brain regions using the Kruskal–Wallis test (scipy.stats.kruskal). We used Bonferroni-corrected P values, calculated using Dunn’s test via scikit_posthocs.posthoc_dunn (Supplementary Table 1).

In Extended Data Fig. 4a,b, statistical analysis was performed using R (v3.4.4). Normality of the microglial average densities was tested using Shapiro–Wilk’s test (shapiro.test()). Differences in densities were compared with two-sided t-test. When comparing between brain regions, Tukey’s post hoc test (TukeyHSD()) was used for multiple comparisons (Supplementary Table 2).

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

留言 (0)

沒有登入
gif