Identification of proliferative and non-proliferative subpopulations of leukemic cells in CLL

Pathogenesis in chronic lymphocytic leukemia (CLL) is strongly linked to the potential for leukemic cells to migrate to and proliferate within lymph-nodes. Previous in vivo studies suggest that all leukemic cells participate in cycles of migration and proliferation. In vitro studies, however, have shown heterogeneous migration patterns. To investigate tumor subpopulation kinetics, we performed in vivo isotope-labeling studies in ten patients with IgVH-mutated CLL (M-CLL). Using deuterium-labeled glucose, we investigated proliferation in sub-populations defined by CXCR4/CD5 and surface (sIgM) expression. Mathematical modeling was performed to test the likelihood that leukemic cells exist as distinct sub-populations or as a single population with the same proliferative capacity. Further labeling studies in two patients with M-CLL commencing idelalisib investigated the effect of B-cell receptor (BCR) antagonists on sub-population kinetics. Modeling revealed that data were more consistent with a model comprising distinct sub-populations (p = 0.008) with contrasting, characteristic kinetics. Following idelalisib therapy, similar labeling suppression across all sub-populations suggested that the most proliferative subset is the most sensitive to treatment. As the quiescent sub-population precedes treatment, selection likely explains the persistence of such residual non-proliferating populations during BCR-antagonist therapy. These findings have clinical implications for discontinuation of long-term BCR-antagonist treatment in selected patients.

INTRODUCTION
The clinical course and outcome of chronic lymphocytic leukemia (CLL) varies markedly between patients. Investigation of the biological basis of this heterogeneity has contributed much to our understanding of disease pathogenesis. As for many tumors, genomic abnormalities are key determinants of outcome, however, functional properties of the tumor, especially its capacity to migrate into lymphoid tissues and proliferate in response to B-cell receptor (BCR) derived signals, are also important.
In addition to differences between patients with CLL, there is also evidence for heterogeneity within the leukemic clone within individual patients. Large scale sequencing studies have revealed significant subclonal genomic diversity [1,2] that evolves over time and in response to therapy [3]. Whether this is accompanied by functional intraclonal diversity, however, remains an open question since phenotype, capacity to migrate, signal and proliferate may all change in response to anatomic location and factors in the tumor microenvironment. For example, leukemic cells within the lymph node (LN) express higher levels of the prognostic marker CD38 than those in the peripheral blood (PB) [4], whilst the chemokine receptor CXCR4 is rapidly downregulated following B-cell receptor signaling, and upregulated in its absence [5].
Intraclonal diversity has been directly studied in vivo using the non-radioactive isotope deuterium, administered as heavy water, to label the DNA of dividing cells in patients with CLL [6,7]. This technique enables the proliferation rate, lifespan and location of populations and subpopulations of tumor cells to be determined from measurements of isotope enrichment in sorted subsets over time. Labeled cells took several days to appear in the PB, compatible with proliferation within and subsequent release from lymphoid tissues, a finding later confirmed by tissue direct sampling [8]. Analysis of sorted subfractions showed that the CXCR4 lo /CD5 hi subset appeared in the blood first and had the highest deuterium enrichment whilst CD5 int /CXCR4 int and CD5 lo /CXCR4 hi subsets show progressively lower and later labeling. In the light of these observations, it was proposed that CXCR4 lo /CD5 hi cells represent the most recently proliferated migrants from tissue to blood, and that, in the PB, these cells quiesce, re-express CXCR4, downregulate CD5 and then re-enter lymphoid tissue where further proliferation may occur. Since deuterium was eventually found in all sorted subsets, these data suggest that the CLL clone behaves in a homogeneous manner and that all cells are similarly able to participate in cycles of migration and proliferation.
Although this model of CLL biology, in which the leukemic clone is functionally homogeneous in respect of its capacity to migrate and proliferate, fits much of the available data, several observations indicate that function may be heterogeneous within the clone. First, although inhibitors of BCR signalling (BCRi) undoubtedly have a major impact on the size of the leukemic clone, in many patients, small numbers of leukemic cells persist [9] suggesting that not all tumor cells are equally dependent on this pathway for survival. Second, CXCR4 lo /CD5 hi peripheral blood CLL cells expressing a similar phenotype to LN CLL cells show preferential migration across vascular endothelium [10], suggesting the existence of a subpopulation primed for LN entry and proliferation. Finally, long term follow-up of patients who have undergone whole exome sequencing has shown a clear relationship between subclonal complexity and tumor growth characteristics in vivo [11].
In view of the above, we re-examined whether the leukemic clone in CLL is always homogeneous in respect of proliferative potential. The deuterium labeling studies described above used heavy water ( 2 H 2 0) which, being administered over several weeks, labels multiple generations of cells. Such studies therefore have a limited capacity to resolve subpopulation evolution over time. To address this issue, we opted to use deuterated-glucose ( 2 H 2 -glucose) as a labeling agent. Because it has a very short half-life it labels a cohort of cells dividing over a period of less than a day and allows the progress of that population of cells to be tracked through subsequent phenotypic transitions and migration. Although CLL cells divide slowly, the high precursor enrichments achieved with deuterated-glucose (about one log higher than 2 H 2 O) allow measurable labeling after a short-pulse label. We have previously used this approach to estimate the lifespan of CLL cells, generating average results very similar to those estimated using 2 H 2 O [12].
In the present study, we focused our investigation on patients with M-CLL, the subgroup of patients most likely to exhibit persistent lymphocytosis during BCR inhibitor therapy [9], and therefore most likely to have more readily detectable functional subclonal heterogeneity. Mathematical modeling techniques were then used to determine whether subsets defined by phenotypic markers associated with migration and BCR signalling capacity (CXCR4/CD5, surface IgM [13] and BCR internalization [14]) are more likely to be homogenous, deriving from each other (Model A, Fig. 1), or instead behave as separate subpopulations that remain distinct over time (Model B, Fig. 1). We also report the impact of the BCR antagonist idelalisib on the release and survival of CLL subpopulations labeled prior to and again during therapy.

Patients and sampling
Two cohorts of subjects with CLL with mutated immunoglobulin heavy chain genes (IGHV) were studied (Table 1). In the first, 10 treatment-naïve patients with non-progressive Binet Stage A or B M-CLL were labeled with oral deuterium-glucose; PB was drawn and CLL cells isolated at set intervals after labeling (Fig. 2). Duplicate studies were performed 8 weeks later to estimate reproducibility.
The second cohort comprised two individuals participating in the CALiBRe study (EudraCT 2012-003631-36) which investigated idelalisib, a phosphoinositide-3 kinase (PI-3K) inhibitor, in patients with previouslyuntreated symptomatic CLL [15]; deuterium labeling was included as an optional sub-study. CALiBRe closed early due to safety concerns with idelalisib, hence recruitment to the labeling sub-study was limited to two patients. In these subjects, deuterium-glucose administration was timed to begin 24 h before commencing idelalisib treatment (150 mg bd); labeling and sampling matched those in the untreated cohort (Fig. 2). All subjects gave written, informed consent; studies were approved by the UK National Research Ethics Service (13/LO/0434 and 15/YH/0020).

Fig. 1
Models for CLL cell proliferation and recirculation. a Model A represents the hypothesis that there is a single recirculating pool of cells. We assume all proliferation (p) occurs in a non-blood compartment, denoted "Tissue" and that the blood compartment can be considered as 3 phenotypically distinct subpopulations. Black arrows denote transition between compartments and grey arrows denote cell death or disappearance. b Model B represents the alternative hypothesis that there are two (semi-)independent cell populations. The PB2 compartment recirculates to the tissue compartment, leaving PB3 as a distinct non-recirculating pool of cells. Cell disappearance (db, disappearance from blood compartment; dl, disappearance from non-blood compartment is thought (but not assumed) to be primarily due to cell death. Version with all parameters marked is shown in Supplementary Fig. 1.

Deuterium-glucose labeling and analysis
In vivo labeling was performed by oral administration of 60 g deuteriumglucose (6,6-2 H 2-glucose, Cambridge Isotope Laboratories, USA) over 10 h, as previously described [16]. Since labeled glucose disappears rapidly (within hours), this approach labels cells dividing during an approximately 12-hour window. PBMC from follow-up samples (Fig. 2) were cryopreserved in 10% dimethyl sulphoxide prior to sorting (as below), then analysed for DNA deuterium enrichment by gas chromatography-mass spectrometry (GC-MS) as previously described [16,17]. Results were expressed as the percentage of newly-proliferated cells per day after normalizing for the level and duration of blood glucose labeling directly measured during the labeling phase, as previously described [16] (Figs. 3, 4, 5, 6).

Immunofluorescence cell sorting
Cryopreserved cells were thawed and stained with fluorochromeconjugated antibodies to CD5/CD19 either alone or in combination with either CXCR4, sIgM or IgM linked to pHrodo, a dye that fluoresces following internalization into acidic endosomes; for details see Supplementary  Table 1. Cells were then flow sorted into fractions representing the whole CLL clone (CD5+/CD19+) or subpopulations defined by expression of CD5 and CXCR4, or IgM at the cell surface (sIgM) or within acidified endosomes (eIgM). Each selected CLL subpopulation comprised 10-15% of the whole CLL cell population, which appeared as a continuum (Figs. 4a, 5a, c). Fluorochromes and gating strategies are detailed in Supplementary Methods.

Statistical analysis and mathematical modeling
We fitted mathematical models describing the kinetics of CLL cells to the glucose enrichment data. The models assumed that CLL cells undergo proliferation in a compartment "tissue" that is not in immediate dynamic equilibrium with blood, i.e. we assumed that cells neither divide in blood nor in a compartment that freely and immediately recirculates with blood; this assumption precludes a model in which PB1 immediately recirculates with "tissue". The precise anatomical correlate is immaterial to the modeling. In addition to tissue we included three PB subsets: PB 1 , PB 2 , and PB 3 representing the three sorted subpopulations: PB 1 : CXCR4 lo /CD5 hi , sIgM hi or eIgM hi ; PB 2 : CXCR4 int /CD5 int , sIgM int or eIgM int ; PB 3 : CXCR4 hi /CD5 lo , sIgM lo or eIgM lo . Since sIgM and eIgM levels in CLL are closely correlated [14], for the purposes of this study the two were considered together.
Two models of CLL kinetics were considered. In the first (Model A) there are no distinct subpopulations; instead, newly-proliferated CLL cells transit from the tissue compartment, through pools PB 1 , PB 2 and PB 3 then back to tissues. In the second (Model B), CLL cells transit from tissues to blood pools PB 1, then PB 2 , then back to tissues, whilst the third subpopulation, PB 3 (CXCR4 hi /CD5 lo , sIgM lo or eIgM lo cells), does not re-enter tissues but forms a stable non-proliferating subpopulation in blood. Both models were fitted to deuterium labeling data from bulk and sorted subsets. The small sample corrected Akaike Information Criterion (AICc) was used to determine which model best described the data (based on fit to data penalized by model complexity). Model parameters are fully defined and equations describing the movement between tissue and blood compartments provided in Supplementary Methods.

RESULTS
Proliferation, release, and disappearance rates of labeled CLL cells DNA deuterium labeling curves for bulk tumor cells from this cohort of patients with non-progressive M-CLL are shown in Fig. 3 (for detail of days 0-7 see Supplementary Fig. 2). The curves reflect three processes; (i) the first part of the curve, generally a delay before the appearance of labeled cells, represents the time taken for cells to transit between their site of division and blood, the site of sampling; (ii) the peak primarily reflects the proliferation ratemore proliferation resulting in greater uptake of deuterium during the short labeling phase; (iii) the subsequent disappearance of label from blood primarily Fig. 2 Study schemas. In the observational study, 10 patients underwent labeling with deuterium-labeled glucose (6,6-2 H 2-glucose), shown by the yellow bars, on two occasions (Cycles 1 and 2) separated by eight weeks, denoted day 0 and day 56 respectively. In the CALiBRe study, the labeling schedule was performed in an identical way but treatment with idelalisib was commenced on day 1, starting 24 h after commencement of the first labeling cycle, and continued for a minimum of 16 weeks (as indicated by the red arrow). Blood was drawn at regular intervals as shown by the blue arrows (Cycle 1: baseline, daily for the first 4 days, then weekly for 4 weeks; Cycle 2: baseline, day 4, weekly for 4 weeks, then on day 56). reflects the rate at which labeled cells die or migrate out of blood. Analysis of DNA deuterium enrichment curves with these processes in mind led us to four general observations. Firstly, we noted striking inter-individual heterogeneity in labeling patterns, particularly evident in the release rates of newly proliferated CLL cells into blood. While some subjects such as K8 and K10, demonstrated rapid release of highly-labeled cells, others, such as K5, K6, K9, showed slower release with a lower peak, and some (e.g. K1) showed very little labeling at all. The time to maximum PB deuterium enrichment ranged from 2 to 56 days (median 10.5 days) and was longest in patients with splenomegaly; four of the five patients with very delayed release (≥14 days) had splenomegaly (K4, K5, K6, K9).
Secondly, we observed generally low deuterium labeling rates in tumor cells, consistent with previously reported low proliferation rates in CLL cells [12,18,19]. In this study, taking peak labeling as a crude minimum estimate of proliferation (disregarding cell death), we estimated a median proliferation rate of 0.26 %/d (IQR = 0.20-0.32 using cycle 1 data, n = 10) equivalent to a median half-life of 271 days (0.29%/d if also including cycle 2 data, n = 19).
Thirdly, we found that labeled CLL cells disappeared remarkably slowly (Fig. 3); labeled cells remaining detectable at 56 days in 5/ 10 cases, suggesting that death rates of labeled cells are low (i.e. that CLL cells have a long lifespan), in keeping with previous observations [6]. Fourthly, although there was wide inter-patient variability, intra-patient reproducibility in repeat studies was generally good (Supplementary Fig. 3).
CXCR4/CD5 expression predicts proliferative subclonal behaviour and identifies a non-proliferative subpopulation Analysis of subfractions of PB CLL cells sorted into three noncontiguous populations according to CXCR4 and CD5 expression (Fig. 4) showed that the CXCR4 lo /CD5 hi subfraction appeared in the PB earliest, consistent with previous studies [7]. The intermediate population (CXCR4 int /CD5 int ) contained less deuterium and took longer to appear in PB, compatible with progressive acquisition of CXCR4 and loss of CD5 expression by the labeled cohort of cells (Fig. 4). The CXCR4 hi /CD5 lo subpopulation, however, exhibited minimal or undetectable labeling throughout the eight-week period of monitoring. Patient K5 proved to have two phenotypically distinct clones and was not included in subsequent analyses.

Expression of BCR identifies subclonal populations with distinct kinetic characteristics
In view of the previous observation that expression of sIgM is dynamic and that low expression is related to a reversible state of anergy, we went on to investigate proliferation in CLL subsets defined by sIgM expression [5,13,20] (Fig. 5a). Assessment of deuterium labeling revealed that CLL cells with the highest sIgM expression (sIgM hi ) represent the most recently proliferated population with progressively lower labeling in sIgM int and sIgM lo subpopulations (Fig. 5b). Interestingly, in some patients, no detectable label appeared in the sIgM lo cells over the entire two months of the study. We then went on to further investigate the functional relationship between subpopulations defined by the presence of IgM within acidified endosomes (eIgM) in three patients; this measure reflects BCR internalisation and correlates with surface expression [14]. We found that cells with high rates of BCR internalisation (eIgM hi ) had the highest rates of proliferation whilst the subpopulation with the lowest levels of eIgM remained unlabeled suggesting a fixed state of quiescence (Fig. 5c, d).
Mathematical modeling frequently predicts the presence of independent subpopulations We next assessed the relationship between subpopulations using mathematical modeling, testing two alternative hypotheses: (i) Fig. 3 Deuterium labeling in CLL cells in ten patients with non-progressive disease. DNA deuterium enrichments in purified CLL cells from 10 subjects (K1-K10) without intervention. Data are expressed as fraction of new cells produced (F, as a percentage) normalised by the glucose enrichment area under curve to the equivalent of one day's labeling, i.e. equivalent to %/day. Note different axes for K8 and K10. Cycle 2 results were similar ( Supplementary Fig. 2). Time represents days following oral labeling with deuterium-labeled glucose (labeling day = d0); error bars represent the standard deviation of at least four replicate measurements by GC-MS (error bars not shown where they are smaller than the symbols).
that low labeling in the CXCR4 hi /CD5 lo, sIgM lo, eIgM lo populations could be explained by dilution in sequentially labeled pools, or (ii) that CXCR4 hi /CD5 lo, sIgM lo, eIgM lo cells represent a distinct quiescent population. These hypotheses are represented by the two models A and B respectively (Fig. 1), to which the deuterium labeling data was fitted.
Overall, across all data sets (i.e. pooling datasets sorted by CXCR4/ CD5 with those sorted by sIgM and eIgM, n = 14) Model B provided     (Table 2). Of these 7 cases, the difference in AICc was large (>=3) for 5 datasets (Table 2). In the other 6 cases, the AICc did not discriminate between the models. Model fits are shown in Supplementary Material. We were therefore able to reject the null hypothesis that Model A always wins with certainty (p = 0.008 Wilcoxon signed rank test two sided). This demonstrates that it is unlikely that there is a continuous single population of leukemic cells and that, at least in some patients, recirculation to LNs most likely occurs at the intermediate stage (Model B).

Idelalisib has a heterogeneous effect on functional subclones
To address the effect of BCR inhibition on these functionally heterogeneous subpopulations, we performed in vivo labeling in two patients (CAL01 and CAL02) with symptomatic CLL participating in the CALiBRe clinical trial. Both had mutated IGHV genes and CD38 expression rates >7% (Table 1). Total lymphocyte count and expression of sIgM by CLL cells was recorded over time to assess the kinetics and phenotype of the leukemic clone following drug administration. As anticipated, there was a rapid rise in the lymphocyte count in the days following treatment initiation, peaking on approximately day 4 and falling thereafter (Fig. 6a). This correlated with an abrupt rise in sIgM expression, again maximal on day 2-4 before returning to baseline by day 7, falling to a mean of 43% of peak values by day 56 (Fig. 6b).
Labeling with deuterium-glucose started 24 h before the commencement of idelalisib and was repeated eight-weeks later. The first cycle of pre-treatment labeling was planned to label all cells actively dividing in situ in the absence of drug; idelalisib treatment would not therefore be expected to affect labeling rates in tissue as any deuterium-labeled glucose would have cleared well before idelalisib administration. According to the previous experiments in untreated patients above, such labeled cells would normally remain at their sites of replication for several days before sequential release into blood, resulting in a peak in PB labeling at 7-14 days (Fig. 4). By contrast, both study patients showed a prompt increase in labeling rates in PB CLL cells, apparent within four hours of idelalisib dosing ( Fig. 6 and in detail for days 0-7 in Supplementary Fig. 4). Labeling peaked on day 4 and was greatest in the sIgM hi compartment; by contrast, labeling in the sIgM lo compartment remained undetectable indicating that sIgMhi cells represent the most recently proliferated population.
Strikingly, when rates of labeling were re-evaluated eight weeks into idelalisib treatment, there was a dramatic, nearcomplete cessation of proliferation in all leukemic cell subsets (Fig. 6c), evident as the absence of an expected second peak in deuterium enrichment. In untreated individuals this peak occurred predictably 0-7 days after labeling (Fig. 5b), especially in the sIgM hi subpopulation and, to a lesser extent, in sIgM int cells; for the second (d56) labeling, this peak would be expected at 56-63 days post-study initiation. Comparing the red lines at d0-7 in Fig. 5b with those at d56-63 in Fig. 6c, the extent to which proliferation appears to be blocked can be appreciated. At this time-point, both treated patients still had a persisting lymphocytosis but median sIgM expression levels were reduced compared to pre-treatment levels and labeling rates were now virtually undetectable in all sorted subpopulations suggesting that CLL cells expressing high levels of sIgM are most sensitive to BCR-directed treatment.

DISCUSSION
In this study, we investigated whether the leukemic cells of CLL patients are homogeneous in their capacity to undergo proliferation or contain distinct subpopulations that differ in this respect. Our results convincingly show that CLL does not always behave as though the tumor is functionally homogeneous and that, at least in this population of patients with mutated IGHV, distinct subsets of cells with different proliferative capacities can be identified using phenotypic markers.
Uniquely, our study used in vivo pulse-chase labeling with deuterated-glucose to investigate this important question. To increase our capacity to detect non-proliferating subpopulations, we focussed on patients with mutated IGHV who are already known to have lower proliferation rates. As expected, we found very low estimated rates of in vivo tumor proliferation (median 0.26%/day), consistent with previous reports with deuterated glucose (0.57%/d) [12] or deuterium-labeled water (0.36%/d) [19], albeit in different patient populations. The delay to peak appearance of newlyproduced leukemic cells in blood is consistent with cell division in the non-blood compartment for all subclones. This lag was highly variable and longest in the four patients with marked splenomegaly, suggesting that spleen and blood together form recirculating pool; when this is large, as in splenomegaly, labeling is slow. In keeping with previous studies, we found that the rates of disappearance of CLL cells from the peripheral blood were much slower than those reported for healthy B lymphocytes [12]. This finding might result from impaired apoptosis in CLL cells [12], but could also reflect trapping of CLL cells in blood.
Analysis of the deuterium content of flow-sorted subpopulations clearly revealed sequential incorporation into CXCR4 lo /CD5 hi , then CXCR4 int /CD5 int , then CXCR4 hi /CD5 lo fractions, consistent with a model in which CXCR4 lo /CD5 hi cells represent the actively-dividing population. (The alternative hypothesis that all subpopulations divide but appear sequentially because they have different release dynamics is refuted by the observation that CXCR4 lo /CD5 hi cells comprise the population enriched for proliferating cells within the LN compartment [8]). However, our observations that deuterium enrichments were substantially lower in CXCR4 int/ CD5 int cells, and even lower (occasionally undetectable) in the CXCR4 hi /CD5l o fraction, suggests that exchange of deuterium between phenotypic groups over time may be limited, in contrast to the conclusions of previous work using deuterated water [7].
Similar patterns were observed when cells were sorted according to sIgM or eIgM levels, both markers of BCR signaling capacity. Deuterium labeling was greatest and most rapid in the sIgM hi/ eIgM hi fraction, consistent with high sIgM expression in LNs, the sites of antigen encounter and proliferation [14,21]. (Although it may seem counterintuitive that cells remain sIgM hi despite recent BCR ligation, levels may be maintained by microenvironmentrelated factors such as CD4 T-cell derived IL4 which enhances sIgM expression and BCR signaling whilst inducing reciprocal downregulation of CXCR4 [22,23].) The most striking finding, however, was that, as for CXCR4/CD5 sorting, a subpopulation defined by reduced surface or endosomal BCR expression remained almost entirely unlabeled throughout the period of the study.
When these patterns of labeling of CLL subpopulations were analysed using mathematical modeling we found very little support for sequential labeling model (Model A) and refuted the null hypothesis that 'Model A is true in all cases of M-CLL'. Conversely, our data suggest that, at least in some if not all patients, most CLL cells return to LN as CXCR4 int/ CD5 int /sIgM int , eIgM int cells, to re-emerge after proliferation as CXCR4 lo /CD5 hi / sIgM hi/ eIgM hi cells. Meanwhile, the CXCR4 hi /CD5l o /sIgM lo /eIgM lo subpopulation forms a quiescent long-lived pool of cells that either does not migrate and proliferate, or if it does, only at rate lower than the limit of detection of the method described, a Deuterium labeling data in CLL subpopulations (defined by CXCR4/CD5 expression, sIgM expression and eIgM expression) were fitted to two mathematical models: Model A, a single population of leukemic cells in peripheral blood with cells transitioning between phenotypic sub-groups, and Model B, with an independent subpopulation of non-proliferating cells.
Lower AIC C (I.e. more negative) indicate a better fit to the data. The 'winning' model is defined by a difference in the AICc of at least 1. Differences of <1 are denoted '-'. The null hypothesis that model A is correct can be rejected with certainty (p = 0.008; Wilcoxon signed rank test two sided).
possibility which is perhaps more likely as some input of new cells is required to balance cell death and preserve population size. In terms of phenotypic characterisation, previous studies have shown that the CXCR4 hi /CD5l o population is enriched for CCR7 and CTLA4 expression, in contrast to the CXCR4lo/CD5 fraction which is enriched for CD38 and CD62L [7]. A hypothetical model, building upon previously-published models [7] is shown in Fig. 7. This interpretation is consistent with work performed using an in vitro circulation system showing that the leukemic clone has heterogeneous kinetics and that migration is limited to a small sub-population characterised by an activated lymph-node phenotype [10,24]. It is also consistent with recent studies showing that subclonal diversity correlates with the in vivo growth characteristics of CLL [11]. Interestingly, in their study, logistic tumor growth pattern, with plateauing of cell numbers after an initial increase, was observed in patients with lower subclonal complexity and mutated IGHV, and it is plausible that this may be a consequence of the emergence of the nonproliferating subset described herein.
Finally, although necessarily limited to two patients by regulatory issues beyond our control, pre-treatment deuterium labeling showed that idelalisib therapy induces rapid release of recentlyproliferated sIgM hi cells into PB, consistent with the observation that Ki67 [25] and sIgM expression [26] in circulating CLL cells rapidly rises following initiation of treatment with ibrutinib, likely as a consequence of release from lymphoid tissues. After eight weeks of treatment with idelalisib, deuterium labeling data showed cessation of proliferation of the sIgM hi subpopulation, as expected from reports that CLL cells expressing high levels of sIgM (the main determinant of BCR signaling) are most sensitive to BCR inhibition [13,20]. The quiescent sIgM lo population meanwhile persisted, perhaps because it is less dependent on BCR signaling and so less sensitive to its inhibition ( Supplementary Fig. 6). Selection of the non-proliferating CXCR4 hi CD5 lo , sIgM lo subset may explain the persistent low-level disease observed during long term ibrutinib therapy [9]. This idea is supported by changes in the tumor immunophenotype over time during ibrutinib therapy as sIgM, eIgM and CD5 expression levels initially increase then fall in the longer term [14,26].
In terms of limitations, although we studied only 10 + 2 subjects, our conclusions are supported by replicate measurements for every sample/time-point and repeated measurements over time within subjects, giving an overall database of over 500 enrichment measurements (Supplementary Material). Inter-individual variation was marked, as would be expected in such a heterogenous disease, and the overall pattern of labelling varied widely between patients (Fig. 3), but despite this the hierarchy and sequential nature of labeling in defined subpopulations was remarkably similar between patients (Fig. 4 & 5). Our primary finding regarding model comparisons is highly unlikely to have arisen by chance (p = 0.008).
Our results raise important further questions. Firstly, how do these functionally distinct subsets arise? The link between anergy and lack of capacity to migrate into LN follicles has been previously described in normal murine B-cells [27]. Our hypothesis is that the nonmigratory/non-proliferating CLL subset we have identified resembles anergic normal B-cells in that both are excluded from entry into follicles but, whereas the latter undergo apoptosis, CLL cells are rescued by bcl-2 expression. Detailed genetic and biochemical studies were outside the scope of the present study, but would determine whether the previously identified clonal heterogeneity [1,11,28] is linked to the functional and kinetic differences we have described here. Secondly, if there is a non-proliferative component in some cases of CLL, is it necessary to eliminate every detectable leukemic cell to obtain lasting clinical benefit? Many patients with monoclonal B-lymphocytosis or stage A CLL remain well for prolonged periods; if it is possible to eliminate the more dangerous Fig. 7 Hypothetical model of leukemic cell kinetics in untreated IGHV mutated CLL. Proposed model of tumor kinetics in CLL that fits the deuterium labeling data in the present study. A proliferative subset of cells transits between PB and lymphoid tissues. These cells become labeled with deuterium, change their phenotype from CXCR4 lo CD5 hi to CXCR4 int CD5 int and express high levels of surface BCR. A largely distinct non-proliferating subset is also present. This subset does not label with deuterium, has a CXCR4 hi CD5 lo phenotype and low expression of surface BCR. These cells generally do not migrate into lymphoid tissues but, unlike normal anergic B cells, which undergo apoptosis when excluded from germinal centres, these cells express bcl-2 and so survive.
proliferative component of the disease with drugs such as ibrutinib, then perhaps therapy can be discontinued once only nonproliferative/anergic cells remain. The observation of prolonged disease stability in some patients who have discontinued ibrutinib because of side effects suggests that this may be the case [29].


Agronomic performance of soybeans with the presence of volunteer RR corn 1

The presence of volunteer corn plants among soybeans, both resistant to glyphosate, has increased in recent years in Brazil. The objective of this study was to quantify the impact on the agronomic performance in RR soybeans by the presence of volunteer RR corn in field conditions. The experiment was conducted in Tiradentes do Sul (RS), 2018/19 harvest. The experimental design used in the experiment was randomized blocks with four replications: control in the clean (only soybean) and eight densities of volunteer corn in soybean crop with 1, 2, 3, 4, 5, 6, 7 and 8 corn plants m -2 . In the physiological maturation of soybean, plant height, number of lateral branches, stem diameter, height of insertion of the first pod, number of pods plant -1 and grain yield were evaluated. All variables analyzed were negatively affected by the increase in the density of volunteer corn infestation, where the variables grain yield and number of pods plant -1 were the most affected, 63% and 81%, respectively. The presence of volunteer corn has a high capacity to reduce the yield of soybeans in the proportion of 500 kg ha -1 plant -1 of corn m -2 , and even at low densities, it needs control.

INTRODUCTION
A widely used cultivation alternative in Brazil, especially in intensive systems of commercial areas, is the succession of off-season soybeans after corn cultivation. However, with the emergence of technologies such as the use of Roundup Ready (RR) grains, the succession of crops results in a new problem for planting systems, such as weeds (Aguiar et al., 2017). For example, the presence of volunteer RR corn plants in coexistence with soybean has been increasing since the introduction of glyphosate-resistant corn varieties on the market (Braz et al., 2019).
Several weeds are commonly found in soybean crops in Brazil and their levels of economic damage are known. However, plants that were not considered harmful to crops now have this characteristic, as in the case of volunteer corn with RR technology that ends up competing with soybean in succession and has become a problem for post-emergence control in RR soy (Marquardt et al., 2012).
The so-called volunteer corn originates from the germination of grains lost by mechanized harvesting that becomes a weed in the soybean crop competing for water, light, nutrients and CO 2 and cannot be eliminated by the application of the glyphosate herbicide, since it ends up being resistant to this type of herbicide.
The potential for volunteer corn to cause damage to soybean crops is high when compared to other weeds. For example, Ipomoea grandifolia, Euphorbia heterophylla or Brachiaria plantaginea plants in the density of one plant m -2 have low potential for loss of soybean yield (Voll et al., 2002). On the other hand, the high density of Raphanus sativus (55 plants m -2 ), in addition to the reduction of plant height and branch length, caused losses in soybean yield by 3 to 15% (Bianchi et al., 2011). Volunteer corn plants have great capacity for competition for growth and development above the canopy of soybean plants (Page et al., 2010), and can also cause metabolic changes due to the occurrence of oxidative stress in soybean plants (Piasecki et al., 2018a), resulting in a reduction in plant growth and loss of grain yield. As a weed, volunteer corn reduces soybean yield by 10 to 22% plant -1 m -2 on average (Aguiar et al., 2017). However, this can be over 40% with high densities (Marquardt et al., 2012;Alms et al., 2016;Aguiar et al., 2017;Piasecki et al., 2018b). Thus, the continuous use of RR cultivars tends to reduce crop yields, especially for off-season soybean in succession to corn. According to Piasecki et al. (2018b), the incorrect handling of herbicides and mechanical losses due to the previous harvest are the two main factors that are related to the presence of volunteer corn among soybean crops.
According to Braz et al. (2019), the level of volunteer corn interference as a weed in soybean crops depends on several factors such as infestation density, location, length of stay and crop stage. Also, there are few site studies, specifically the response to interference imposed by corn plants due to the sowing of soybean cultivars with different growth habits, development cycle, volunteer corn infestation densities and type of specific local climates (Nordby et al., 2007;Braz et al. 2019). There are studies in the literature that assess the interference of volunteer corn in soybeans, but few verify the interference in conditions of local agroecosystems, short-cycle cultivars and low densities of volunteer corn. The hypothesis of this study is that, due to the competition capacity, low densities of volunteer corn plants are capable of promoting loss in the yield in soybean culture in the mesoclimate of the Northwest region of Rio Grande do Sul. Thus, the objective of this study was to quantify the impact on agronomic performance caused by the presence of volunteer RR corn in RR soybean under field conditions.

MATERIAL AND METHODS
The field experiment, in a crop area, was carried out in the 2018/19 harvest in the locality of Linha Formosa, in Tiradentes do Sul, State of Rio Grande do Sul (RS), latitude 27º23'51" south and longitude 54º05'02" west, 407 meters above sea level. According to the Köppen classification, the climate in the region where the experiment was conducted is humid subtropical (Cfa). The soil at the beginning of the experiment was evaluated for its physical-chemical composition ( Table 1). The soil fertility levels were corrected following the recommendations of the CQFS RS/SC (2016).
The experiment was carried out in random blocks with four replications of the following treatments: Controlarea in the clear (RR soybean only) and eight densities of RR corn in the soybean crop with 1, 2, 3, 4, 5, 6, 7 and 8 corn plants m -2 . Each experimental plot contained 16 lines of 6 meters in length with a spacing between lines of 45 cm, totaling 43.2 m -2 . The useful area considered for each plot was the eight central lines (3.6 m) with 3.0 m in length (10.8 m -2 ). Before the installation of the experiment, a dose of 2.5 L ha -1 of glyphosate (356 g L -1 ) and another subsequently after five days of 2.0 L ha -1 of paraquat dichloride (200 g L -1 ) were applied to total desiccation of the plants in the experimental area.
The sowing of the soybeans was carried out in succession to early corn crop on the same day as the application of paraquat dichloride. Before sowing, the seeds were treated with Bradyrhizobium japonicum, with peat-type inoculant, using 0.5 kg of inoculant for every 60 kg of soybean seeds. The soybean cultivar resistant to glyphosate (NA 5909 RG) of the undetermined-growth type and semi-early cycle, was sown with a mechanized seeder on the December 1 st 2018, with a density of 292,000 plants ha -1 . Next, F2 volunteer corn, from the cultivar Coodetec 3410 RR, was sown manually and randomly distributed alleatory in between rows of each soybean plot, according to the respective density of treatments volunteer corn. The presence of weeds in the soybean crop, except for corn, was controlled manually. The application of fungicides and insecticides was carried out according to the technical recommendation of the soybean crop. The experiment was carried out under natural conditions and without irrigation (Figure 1).
At maturity, stage R8, ten soybean plants from the useful area of each plot were randomly selected for evaluation of the plant height, number of lateral branches, stem diameter, height of insertion of the first pod and number of pods plant -1 . To determine the height of the plants, a measuring tape was used for measuring from the surface of the soil to the apex of the plant. Determination of the height of insertion of the first pod was also used to measure the distance from the soil surface to the first pod. The diameter of the stem at ground level was measured with a digital caliper. Regarding the number of lateral branches, all the branches that had pods were counted, starting from the main stem of each plant. The soybean grain yield was determined by manual harvesting of the four central lines of each plot. The grain weight was adjusted to 13% moisture and the yield expressed in kg ha -1 . The data were subjected to analysis of variance, and when significant, adjustments between models for corn populations were made using polynomial regression analysis. The computer program software Genes was used (Cruz, 2006).

RESULTS AND DISCUSSION
The presence of volunteer corn positively influenced the height of the soybean, the higher the density of corn, the greater the height of the soybean plants, in the proportion of 0.67 cm for each corn plant added in the area (Figure 2). Braz et al. (2019) also observed a positive correlation in the height of soybean plants as the density of corn plants. Aguiar et al. (2017), observed an increase of approximately 1 cm in the height of the soybean plants for each corn plant, 49% higher to the result obtained in our study. These changes in the morphology of the soybean plants, which are higher than the normal pattern, indicate that they are in a competitive environment with the weeds, although this competition is apparently not critical for the development of both of them. According to Page et al. (2010), volunteer corn plants have a large growth capacity above the canopy of soybean plants, resulting in a reduction in the growth of soybean plants. Radosevich et al. (1997), found that the weed, in the case also applied to volunteer corn, reflects light at a certain wavelength, which is recognized by the plants of the other crop, which causes soybean to increase its height to capture the maximum of solar radiation. Ballare & Casal, (2000) report that from that moment, the plant begins to adjust its growth according to the presence of weeds. Similar data from older studies recorded by Esechie (1992) had already registered a relationship between the height increase of the planted species in relation to the weed. Studies from the last 10 years (Silva et al., 2009;Aguiar et al., 2017;Braz et al., 2019), report that, in addition to the competition for light and increase in plant height, with the coexistence between the cycle of species there will also be competition for soil nutrients, as well as for water, which may cause a drop in the crop yield.
The insertion of the first pod in the soybean corresponds to the height measurement in relation to the soil up to the point where the insertion of the first pod with grain appears in the main stem of the plant. The higher  the density of volunteer corn in soybeans, the higher the pods insertion (Figure 3). The insertion of the first pod was 21.8 cm, whereas in a competition environment, it was 0.35 cm higher for each volunteer corn plant in the area. Pires et al. (2012), state that plants with lower first pod insertion tend to have higher yield in grain production. In this sense, for the cultivar used, the result obtained indicates a reduction in grain production in soybean culture. According to Merotto et al. (2002), the interference to the passage of sunlight caused by the weed directly affects the insertion of the first pod in the soybean plant, as it reduces its photosynthetic efficiency. In addition, it is possible to observe a variation in the morphological behavior such as horizontal growth to increase the leaf area of the plant.
The measure of the stem diameter, measured on the main stem before the first branching of the plant, showed that the reduction in diameter was proportional to the increase in the density of volunteer corn plants ( Figure   4). A density of three volunteer corn plants per square meter has already caused a decrease by 30% in the stem diameter of soybean plants. This decrease reached more than 50% with the insertion of eight corn plants. A study by Aguiar et al. (2017), showed an average reduction of 31% in the stem diameter of two soybean cultivars with corn density of 1 to 8 plants m -2 . This decrease in stem diameter is directly linked to the increase in the height of soybean plants, which increased by about 10 cm with the presence of eight corn plants m -2 (Figure 2). According to Braz et al. (2019), in general, etiolate plants tend to have a smaller stem diameter, which leaves them more exposed to lodging. Crotser et al. (2003), report that the elongation of the plant causes a decrease in the diameter of the stem, as weeds induce a change in the quality of light in the lower canopy, therefore, reducing the diameter of the stem.
The number of side branches also decreased with the presence of volunteer corn ( Figure 5). While the stem   (Figure 4), the number of branches was reduced by almost 60% with this same plant density. This result corroborates with Braz et al. (2019), who observed a reduction of 15, 41 and 60% when competing with corn at densities of 2, 4 and 8 plants m -2 , respectively. In this context, the development of the soybean crop is delayed with the presence of a greater number of corn plants. Carvalho & Veline (2001) found a delay in the lateral growth of soybean plants when they competed for solar radiation and nutrients with weeds. Bianchi et al. (2011), showed that soybean cultivars in competition with Raphanus sativus reduced plant height, average branch length and soybean yield.
The number of pods plant -1 was the component most affected by the increase in the density of volunteer corn plants ( Figure 6). Volunteer corn plants caused up to 81% reduction in the number of pods plant -1 when subjected to the densities of 8 plants m -2 . This result corroborates those obtained by Aguiar et al. (2017), showing that the volunteer corn plants caused a reduction of 24 to 74% in the number of pods plant -1 when submitted to densities of 1 to 8 plants m -2 , respectively. Silva et al. (2008), found that in the areas of low, medium and high infestation of several species of weeds, at the end of the cycle, there was a reduction of 58, 71 and 78%, respectively in the number of soybean pod plant -1 . Juan et al. (2003), observed that among the components of the soybean yield, the number of pods plant -1 was the most affected by the competition of Euphorbia dentata, with reductions of 40% in relation to the treatment without weeds. However, according to Lamego et al. (2004) this effect is related to the soybean cultivar, showing that cultivars with highly competitive ability, in addition to tolerating competition, preserving the yield potential, also suppress the grain production of competing plants. The space limitation, both aerial and underground, promoted by the volunteer corn plants, could also affect the formation, growth and development of the reproductive structures of the soy plants; however, it is not possible to state which of the factors is the most relevant. The number of pods in a plant is a prior indication of the yield of soybeans in a given area. In our study, the number of pods showed a high degree of correlation with grain yield (r 2 = 0.92; p < 0.05). Juan et al. (2003), also observed high levels of correlation between grain yield, number of pods plant -1 and weight of 1000 soybean seeds, 98 and 85%, respectively. Thus, according to Silva et al. (2008), the number of pod plant -1 is one of the yield components most affected by weeds in the soybean crop.
Our study found that the final consequence of the presence of volunteer corn plants is the yield of soybean grains (Figure 7). It can be seen that the soybean yield decreased by 63% in the presence of 8 plants m -2 of volunteer corn. Piasecki et al. (2018b) observed a maximum loss of 83% with 13 maize plants m -2 , which is 32% higher than that observed in our study. Aguiar et al. (2017), observed a reduction of up to 58% in soybean grain yield, also with 8 plants m -2 . Braz et al. (2019) found for each volunteer corn plant m -2 a reduction of 9.3 and 13.5% for soybean cultivars BMX Potência RR and M8210 IPRO, respectively, in Brazilian Cerrado conditions. Marquardt et al. (2012), concluded that the reductions in soybean yield varied from 10 to 41%, where the densities of earlyemerged volunteer corn varied from 0.5 to 16 plants m -2 , respectively. Alms et al. (2016), showed that 0.3 and 4.4 plants m -2 of volunteer corn reduced 9 and 51% the soybean grain yield, respectively.
The reduction in grain yield is one of the factors most impacted by the presence of volunteer corn in soybeans. As already observed in the number of pods, grain yield has also been greatly reduced with the increase in the density of volunteer corn. Also, according to study by Piasecki et al. (2018a), in addition to the competition for light, water and nutrients, the intense reduction in soybean yield under the influence of volunteer corn may be related to oxidative stress resulting from the greater production of reactive oxygen species (ROS), such as hydrogen peroxide (H 2 O 2 ). This data is worrying when considering production costs, since the yield becomes less than 40% of the expected when the density is 8 plants m -2 (Piasecki et al., 2018a). Thus, it is important to carry out the proper management of crops and techniques to control volunteer corn plants to improve the grain yield of the soybean crop. The use of soybean cultivars with a semi-early cycle and indeterminate growth may have influenced the reduction in the number of pods and grain yield, even in low densities of volunteer corn, as cultivars with a shorter cycle are generally more susceptible to competition and stress.
It was also possible to estimate with our study, in addition to the negative effect of the presence of volunteer corn plants on soybean grain yield, the economic loss. The presence of only one corn plant per square meter resulted in a reduction of 500 kg of soybeans ha -1 (Figure 7). The average price offered per bag of 60 kg of soybeans in March 2020 in Brazil was $ 15.86 (CONAB, 2020). In financial terms, there would be a loss of $ 132.16 plant -1 of volunteer corn m -2 . Aguiar et al. (2017), in a more detailed study, with the expectation of high grain harvest, 5000 kg ha -1 and $ 27.3 per bag, presented an economic loss with 0.96 plants volunteer corn m -2 . Although the price of soybean and the cost of control are factors that affect the control decision, in practice, with one plant per square meter, it could indicate the level of need for controlling volunteer corn. Finally, it should be noted that the results obtained in our study come from and reflect the interference of low density of volunteer corn in a short-cycle soybean cultivar in mesoclimate and agroecosystem in the northwest region in the State of RS.


PacBio long-read amplicon sequencing enables scalable high-resolution population allele typing of the complex CYP2D6 locus

The CYP2D6 enzyme is estimated to metabolize 25% of commonly used pharmaceuticals and is of intense pharmacogenetic interest due to the polymorphic nature of the CYP2D6 gene. Accurate allele typing of CYP2D6 has proved challenging due to frequent copy number variants (CNVs) and paralogous pseudogenes. SNP-arrays, qPCR and short-read sequencing have been employed to interrogate CYP2D6, however these technologies are unable to capture longer range information. Long-read sequencing using the PacBio Single Molecule Real Time (SMRT) sequencing platform has yielded promising results for CYP2D6 allele typing. However, previous studies have been limited in scale and have employed nascent data processing pipelines. We present a robust data processing pipeline “PLASTER” for accurate allele typing of SMRT sequenced amplicons. We demonstrate the pipeline by typing CYP2D6 alleles in a large cohort of 377 Solomon Islanders. This pharmacogenetic method will improve drug safety and efficacy through screening prior to drug administration.

The authors describe a pipeline for accurate allelotyping of the CYP2D6 gene using long read amplicons. CYP2D6 is an important locus due to its relevance in metabolism of drugs and the presence of a pseudogene complicates the process of sequencing this locus. There have been many published papers on genotyping CYP2D6 using long read sequencing technologies (as summarized by the authors in Table 1). The strengths of this paper -as summarized by the authors -are a robust data processing pipeline and the large number of samples that they genotyped. Overall, the paper is well written and easy to follow. Since this problem has been tackled by several groups, I think that the authors could have done additional analysis in order to distinguish their work from previously published work and demonstrate the broad utility and robustness of their pipeline.
1. There have been several pipelines developed for genotyping CYP2D6 from Illumina WGS data (e.g. Cyrius, https://www.nature.com/articles/s41397-020-00205-5). How does the targeted long read based sequencing compare to the Illumina based allelotyping? One advantage of the Illumina WGS based genotyping is that copy number does not need to measured separately using qPCR. Table 1, the authors list previous studies for CYP2D6 typing. The Buermans et al. study is listed as having the pipeline available. It would be useful to analyze the data generated in this study using this pipeline for direct comparison.

In
3. Along the same lines as (2) above, the authors could analyze previously published targeted long read data for CYP2D6 to demonstrate that their pipeline works for different types of long read amplicon data and sequencing technologies.
Reviewer #2 (Remarks to the Author): Dear authors, I wish to congratulate the authors on this interesting work. The manuscript is well written and the experiments well described and executed. The application of long-read sequencing to these challenging parts of the genome is highly suitable. I have only minor comments. Sincerely,

Wouter De Coster
Page 4, line 107 mentions the accuracy of long-read sequencing platforms. These error rates are, at least for ONT, highly outdated, and 2D consensus sequencing has also been deprecated since 2017. The company claims the modal accuracy of the latest base callers and chemistry is at 99%, but conservatively ~5% would be a reasonable estimate of currently produced data. Finding citable references for such a moving target is challenging, but https://www.nature.com/articles/s41588-021-00865-4 reports a median accuracy of 11.6%. Better results are obtained with the 1D^2 method: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-02255-1 On page 14, line 337 recommendations are made for the number of reads that are necessary for the identification of fusion alleles and variant phasing. It is however not clear to me how these numbers were determined, and if a downsampling analysis was performed to obtain these estimates.
The authors used GATK HaplotypeCaller for variant calling, which may be appropriate for highquality HiFi reads, but I wonder about the reason to not use a variant caller tailored to long-read sequencing, e.g. DeepVariant, Clair, or Longshot? While this does probably not affect the results in a major way, some of these tools additionally leverage the phasing of variants or are more appropriate for the error types in long-read sequencing, and I believe it would be relevant to comment on the choice of variant caller.
The authors describe a pipeline for accurate allelotyping of the CYP2D6 gene using long read amplicons. CYP2D6is an important locus due to its relevance in metabolism of drugs and the presence of a pseudo-gene complicates the process of sequencing this locus. There have been many published papers on genotyping CYP2D6 using long read sequencing technologies (as summarized by the authors in Table 1). The strengths of this paper -as summarized by the authors -are a robust data processing pipeline and the large number of samples that they genotyped. Overall, the paper is well written and easy to follow. Since this problem has been tackled by several groups, I think that the authors could have done additional analysis in order to distinguish their work from previously published work and demonstrate the broad utility and robustness of their pipeline.
1. There have been several pipelines developed for genotyping CYP2D6 from Illumina WGS data (e.g. Cyrius, https://www.nature.com/articles/s41397-020-00205-5). How does the targeted long read based sequencing compare to the Illumina based allelotyping? One advantage of the Illumina WGS based genotyping is that copy number does not need to measured separately using qPCR.  We thank the reviewer for this comment, and we agree that Cyrius in particular appears to be a good approach to genotyping CYP2D6 genotypes from short-read WGS. While Cyrius has the advantage of not requiring a separate qPCR assay, we note that Cyrius does not have any means of typing novel CYP2D6 alleles, and running short-read WGS in place of multiplexed long-read amplicon sequencing for CYP2D6 genotypes alone is not cost effective. We have added this information and the reference to Cyrius in the introduction (line 96-98). Table 1, the authors list previous studies for CYP2D6 typing. The Buermans et al. study is listed as having the pipeline available. It would be useful to analyze the data generated in this study using this pipeline for direct comparison.  We thank the reviewer for this comment. The "pipeline" used in the Buermans et al. study we noted as being available is PacBio's long amplicon analysis protocol, and only carries out a portion of the analysis steps that are part of our much more comprehensive pipeline (essentially clustering/phasing and chimera removal, but not allele-typing). We have updated Table 1 with the column "End-to-end Pipeline available" to make it clearer that no previous publications have released an end-to-end pipeline. We do agree that it would be useful to have a benchmark of a pipeline used by another study, so we have implemented as close as possible the complete pipeline used by Buermans et al (which is only a subset of our more complete end to end pipeline) and have tested our data using this pipeline. We have updated the results section at lines 261-264 reflecting this benchmark, as well as the methods section at lines 553-562.

In
3. Along the same lines as (2) above, the authors could analyze previously published targeted long read data for CYP2D6 to demonstrate that their pipeline works for different types of long read amplicon data and sequencing technologies.  We thank the reviewer for this comment. The pipeline as it exists now is specific to PacBio data. In particular the pre-processing stage uses PacBio specific tools (ccs, lima). However, the pipeline has been released as free and open-source software, and it would be feasible for a separate pre-processing module to be written for other sequencing technologies. Additionally, no other study to our knowledge has publicly released PacBio CYP2D6 sequencing data.


s2fieldsofstudy: Biology Psychology

text:
Parvalbumin-Positive Interneurons Regulate Cortical Sensory Plasticity in Adulthood and Development Through Shared Mechanisms

Parvalbumin-positive neurons are the largest class of GABAergic, inhibitory neurons in the central nervous system. In the cortex, these fast-spiking cells provide feedforward and feedback synaptic inhibition onto a diverse set of cell types, including pyramidal cells, other inhibitory interneurons, and themselves. Cortical inhibitory networks broadly, and cortical parvalbumin-expressing interneurons (cPVins) specifically, are crucial for regulating sensory plasticity during both development and adulthood. Here we review the functional properties of cPVins that enable plasticity in the cortex of adult mammals and the influence of cPVins on sensory activity at four spatiotemporal scales. First, cPVins regulate developmental critical periods and adult plasticity through molecular and structural interactions with the extracellular matrix. Second, they activate in precise sequence following feedforward excitation to enforce strict temporal limits in response to the presentation of sensory stimuli. Third, they implement gain control to normalize sensory inputs and compress the dynamic range of output. Fourth, they synchronize broad network activity patterns in response to behavioral events and state changes. Much of the evidence for the contribution of cPVins to plasticity comes from classic models that rely on sensory deprivation methods to probe experience-dependent changes in the brain. We support investigating naturally occurring, adaptive cortical plasticity to study cPVin circuits in an ethologically relevant framework, and discuss recent insights from our work on maternal experience-induced auditory cortical plasticity.

INTRODUCTION
The cerebral cortex has a remarkable capacity to remodel its synaptic structures and refine its neuronal activity in response to changing conditions throughout life. Such cortical plasticity interacts with a wide range of fundamental biological processes including aging (Freitas et al., 2011), stress (Takatsuru and Koibuchi, 2015;McGirr et al., 2020;Muehlhan et al., 2020), injury and recovery (Grasso et al., 2020;Moreno-López and Hollis, 2021), and sensory learning (LeMessurier and Feldman, 2018). A substantial body of evidence supports the crucial role of GABAergic inhibition in sensory cortical plasticity both during development and adulthood.
The neocortex is built around a "canonical circuit, " a set of vertically connected circuit motifs composed of diverse excitatory and inhibitory cell types that establish common hierarchical neuronal computations (Pluta et al., 2015). This arrangement amplifies thalamic input and integrates that input with projections from other regions to accomplish the complex tasks of sensory processing, learning, and memory. Within the cortex, there is a wide disparity between the numbers of excitatory and inhibitory neurons; glutamatergic excitatory neurons compose the overall majority, whereas GABAergic inhibitory interneurons account for only ∼10-20% of neurons (Markram et al., 2004;Petilla Interneuron Nomenclature Group et al., 2008;DeFelipe et al., 2013;Chu and Anderson, 2015). Despite this mismatch in relative densities, cortical inhibitory drive tends to match excitatory drive (Shadlen and Newsome, 1994;Rubenstein and Merzenich, 2003). This suggests that changes to the synaptic connectivity, morphology, intrinsic properties, and firing patterns of inhibitory neurons play a disproportionate role in controlling the timing and specificity of cortical responses (Hensch, 2004;Lu et al., 2007;Trachtenberg, 2015).
In the last two decades, the contributions of specific inhibitory interneuron subtypes to sensory learning, memory, and plasticity have been studied in the visual, auditory, somatosensory, and olfactory cortices. Here we review some of that literature, focusing on cortical parvalbumin-expressing interneurons (cPVins), which are the most abundant inhibitory cell type in the brain. We discuss how the physiological and molecular features of cPVins equip sensory networks for plasticity. We argue that cPVins across sensory cortices share core intrinsic features and mechanisms that enable periods of experiential learning. We further propose that approaches exploiting ethologically relevant behaviors are important for understanding how cPVin directed plasticity is engaged naturally in the adult brain. While there are many varieties and mechanisms of "plasticity, " our use of the term here denotes adjustments to intrinsic and extrinsic connectivity in sensory cortices that optimize representation of stimuli. We particularly focus on the mechanisms that open and close episodes of heightened sensitivity of cortex to stimuli in response to developmental programs or novel sensory experience.

INHIBITORY CONTROL OF CORTICAL PLASTICITY
Many of the fundamental features of cortical plasticity were first identified by studies examining the visual cortex during developmental critical periods (e.g., Wiesel and Hubel, 1963;Levay et al., 1978;Gordon and Stryker, 1996;Hensch et al., 1998;Huang et al., 1999;Sawtell et al., 2003). Much of this work observed the consequences of restricting sensory experience, as occurs in monocular deprivation when vision is restricted in one eye. When performed during a specific period early in life, this manipulation weakens responses corresponding to the deprived eye and strengthens responses to the remaining eye, a phenomenon known as ocular dominance plasticity Fagiolini et al., 1994;Gordon and Stryker, 1996).
The onset and offset of developmental critical periods are defined by changes in cortical GABAergic transmission (Fagiolini et al., 2004;Maffei et al., 2006Maffei et al., , 2010. Premature maturation of GABAergic inhibitory networks terminates the V1 developmental critical period precociously . Mice that are deficient in glutamate decarboxylase (GAD), the rate limiting enzyme involved in synthesis of GABA, have altered neural circuitry and are insensitive to monocular deprivation (Hensch et al., 1998;Fagiolini and Hensch, 2000;Chattopadhyaya et al., 2007). Further, pharmacological disruption of GABA transmission is sufficient to impair typical critical period opening (Jiang et al., 2005;Deidda et al., 2015). On the other hand, complementary approaches that facilitate the development of inhibitory networks accelerate the onset of the critical period (Hanover et al., 1999;Fagiolini and Hensch, 2000;Di Cristo et al., 2007;Sugiyama et al., 2008).
Plasticity continues to be mediated by GABAergic transmission in adulthood (Nahmani and Turrigiano, 2014). The mouse visual cortex retains some plasticity in response to monocular deprivation following critical period closure (Sawtell et al., 2003;Hofer et al., 2006;Fischer et al., 2007;Sato and Stryker, 2008). In the adult V1, cortical potentials corresponding to stimulation of the non-deprived eye are strengthened with experience (Sawtell et al., 2003;Sato and Stryker, 2008). This is associated with enhanced visual acuity in the spared eye (Iny et al., 2006). Manipulations of inhibitory signaling in the adult visual cortex can reactivate experience-dependent plasticity mechanisms. For example, restoring inhibitory function in GAD mutants rescues the capacity for plasticity in response to visual deprivation (Hensch et al., 1998;Fagiolini and Hensch, 2000;Iwai et al., 2003;Chattopadhyaya et al., 2007). Interestingly, transplant of inhibitory precursor cells enhances adult cortical plasticity (Southwell et al., 2010;Tang et al., 2014;Davis et al., 2015). On the other hand, Harauzov et al. (2010) reported that reducing intracortical inhibition reactivated visual plasticity, possibly due to a reduction in chondroitin sulfate proteoglycan structures in the extracellular matrix (ECM) called perineuronal nets (PNNs).
Several pieces of early circumstantial evidence implicated cPVins, among all cortical GABAergic populations, as having a central role in regulating visual cortical critical periods. Mutant mice that exhibit an early critical period also show early cPVin maturation (Hanover et al., 1999;Krishnan et al., 2015). Plasticity can be restored in GAD mutants by manipulating GABA A receptors, which are highly and specifically enriched on cPVins (Fagiolini et al., 2004). More directly, dissolution of cPVin-specific extracellular structures (perineuronal nets-PNNs) reinstates plasticity in adult mice (Pizzorusso et al., 2002). Evidence from Kuhlman et al. (2011Kuhlman et al. ( , 2013 established that cPVins play a permissive role in critical period plasticity by disinhibiting the excitatory network. However, attempts to reproduce this effect by chemogenetically and optogenetically suppressing cPVin inhibition following monocular deprivation in older mice has yielded mixed results (Kuhlman et al., 2013;Saiepour et al., 2015;Kaplan et al., 2016). This may be due to methodological differences in the parameters used for manipulating PV neurons such as the depth of optical penetration, magnitude, duration, and frequency (El-Boustani and Sur, 2014) or the timing of administration of chemogenetic tools. Interestingly, recent work has shown that monocular deprivation can drastically change the number of synaptic inputs onto V1 cPVins (Ribic, 2020). This suggests that the observed fluctuations in inhibitory drive may be a compensatory mechanism downstream of aberrant synaptic pruning of excitatory inputs onto cPVins (Severin et al., 2021).

MOLECULAR AND STRUCTURAL PLASTICITY OF CORTICAL PARVALBUMIN-EXPRESSING ENTERNEURONS
At the finest spatial scale, cPVin control of cortical plasticity is evident at the level of subcellular interactions with perineuronal structures. The extracellular matrix (ECM) is a complex structure composed of elastin, fibronectin, integrin, glycoprotein, and polysaccharide macromolecules that surrounds cells, and which accounts for 10-20% of the brain's volume (Cragg, 1979). The ECM supports a wide range of physiological processes including remyelination (You and Gupta, 2018), stem cell storage (Roll and Faissner, 2014), protection against reactive oxygen species (Cabungcal et al., 2013), synapse development , and stabilization of mature synapses (Frischknecht and Gundelfinger, 2012).
Extracellular matrix components are not uniformly distributed across brain regions or cell types (Dauth et al., 2016). For example, perineuronal nets (PNNs) are chondroitin sulfate glycosaminoglycan (CS-GAG)-based, net-like structures found within ECMs that preferentially surround cPVins and their proximal dendrites (Härtig et al., 1999;Wegner et al., 2003). A hyaluronan backbone acts as a scaffold for assembling CS-GAG, link proteins, and other anchoring proteins into these condensed mesh structures which attach to cell membranes (Brückner et al., 1993;Day and Prestwich, 2002). PNNs contain a number of specific CS-GAGs (e.g., aggrecan, neurocan, brevican) which represent diverse sidechains that attach to lectican protein cores (Giamanco et al., 2010;Frischknecht and Gundelfinger, 2012). How the specific combination of these molecules influences PNN properties is an active area of investigation (Galtrey et al., 2008;Miyata et al., 2018).
Expression of the parvalbumin (PV) protein itself also fluctuates according to cPVin activity (Kamphuis et al., 1989;Patz et al., 2004). Downregulation of PV has broad consequences for the physiological and synaptic properties of cPVins (Caillard et al., 2000;Schwaller, 2012), as well as for behavior (Wöhr et al., 2015). These effects are likely related to PV's role in rapidly buffering excess calcium (Aponte et al., 2008). Downregulation of PV expression may also reflect shifts in regional excitatoryinhibitory balance resulting from pathology, for example, as broadly observed in animals models and post-mortem tissue from patients with Autism Spectrum Disorders (Schwaller, 2012;Filice et al., 2016). Moreover, there is evidence of strong mutual influence between the composition and maturity of PNNs and the excitability and activity of cPVins (Balmer, 2016;Favuzzi et al., 2017;Chu et al., 2018;Gottschling et al., 2019;Devienne et al., 2021). The complex relationship between PV expression, cPVin physiology, and PNN assembly are matters of ongoing study, and will further inform our understanding how cPVins regulate cortical plasticity.
Well-developed PNNs appear to be a necessary component for termination of developmental critical periods. Manipulations that disrupt critical period closure (e.g., dark-rearing) delay disinhibition and slow the formation of PNNs (Lander et al., 1997;Pizzorusso et al., 2002). Moreover, genetic knockout of PNN components (Rowland et al., 2021), enzymes required for the synthesis of those components (Hou et al., 2017), or link proteins necessary to establish PNN architecture (Carulli et al., 2010;Ribic, 2020) also prevent critical period closure. Meanwhile, pharmacological inhibition of endogenous MMP activity disrupts ocular dominance plasticity in the adult visual cortex, i.e., a shift in V1 activity to preference the spared eye (Pielecka-Fortuna et al., 2015;Akol et al., 2022).
Conversely, plasticity can be reinstated by pharmacological dissolution of PNNs in the adult cortex, which allows for enhanced synaptic plasticity and physical restructuring of mature PV+ synapses (Magnowska et al., 2016;Ferrer-Ferrer and Dityatev, 2018). For example, in the adult visual cortex, removal of PNNs reinstates most aspects of ocular dominance plasticity, although it is "incomplete" compared to that instated during developmental critical periods (Pizzorusso et al., 2002). This includes disinhibition seen as a drop in V1 cPVin firing rates (Lensjø et al., 2017), and remobilization of dendritic spines, allowing them to lengthen and potentially change the architecture of established circuitry (de Vivo et al., 2013;Faini et al., 2021). Likewise, the effects of monocular deprivation during development on V1 can be partially reversed during adulthood via environmental enrichment, potentially due to a drop in PNNs and reduced inhibitory drive (Sale et al., 2007).
Perineuronal net-fluctuations as a marker of plasticity can also be observed outside of the context of monocular deprivation. For example, reintroduction of light in dark-exposed animals triggers V1 plasticity through the endogenous degradation of the ECM by MMP-9, and resulting disinhibition of fast-spiking interneurons (Murase et al., 2017). In the barrel cortex, whisker trimming leads to downregulation of PV-PNNs (McRae et al., 2007). Finally, evidence from our lab and others also suggests that aberrantly high or low numbers of PNNs during critical periods of sensory learning may impair A1 sound processing (Krishnan et al., 2017;Wen et al., 2018;Pirbhoy et al., 2020).

PHYSIOLOGICAL PROPERTIES OF CORTICAL PARVALBUMIN-EXPRESSING ENTERNEURONS
cPVins represent an estimated 5-20% of cells in the brain, and up 40-50% of the inhibitory population, based on developmental progenitor and tracing studies (Tamamaki et al., 2003;Butt et al., 2005;Rudy et al., 2011;Rodarie et al., 2021). Despite similar genetic profiles (Tasic, 2018), and gross numbers of afferent thalamic inputs (Pouchelon et al., 2021) across cortical IN subtypes, individual cPVins have unique cellular physiological characteristics that enable them to locally shape sensory responses.

Temporal Control of Sensory Responses
One such characteristic is their high spiking rates they exhibited spontaneously in vivo or in response to somatic current injection (Connors and Gutnick, 1990), leading to their original label of "fast spiking" cells (Kawaguchi and Hama, 1987;Kawaguchi and Kubota, 1993;Kawaguchi, 1995). "Fast spiking" cells are now widely acknowledged to correspond to PV expressing GABAergic basket cells and to a lesser extent, chandelier cells (Kawaguchi et al., 2019). The fast spiking attributes of cPVins result from their expression of rapidly repolarizing Kv3 potassium channels (Rudy and McBain, 2001) and Ca 2+ -permeable AMPA receptors which have rapid gating kinetics (Jonas et al., 1994;Geiger et al., 1997).
The brisk spiking output and high reliability of cPVin firing enables them to exert tight temporal control over cortical excitation. They convey precisely timed, feedforward inhibition that lags just behind feedforward excitatory input from the thalamus, thereby imposing strict temporal limits (Wehr and Zador, 2003;Zhang et al., 2003;Mittmann et al., 2005;Wilent and Contreras, 2005;Heiss et al., 2009). The potency of this mechanism is increased by cPVins' particular sensitivity to onsets and offsets of sensory stimuli (Wehr and Zador, 2003;Weible et al., 2014;Keller et al., 2018;Li et al., 2021). Over broader timescales, cPVins make important contributions to stimulus specific adaptation, which allows the brain to maintain sensitivity to novel or behaviorally relevant stimuli among repetitive distractors (Pérez-González and Malmierca, 2012;Natan et al., 2015).
Precise temporal regulation of cortical activity by cPVins is also relevant for more persistent changes in sensory responses through long-term potentiation (LTP) or depression (LTD). These mechanisms rely on coincidence detection of preand post-synaptic neuronal spiking. Spike timing-dependent plasticity (STDP), which modifies synapses based on the precise relative timing of individual pre-and post-synaptic spikes (Bi and Poo, 1998), has been widely observed in the cortex (Bender et al., 2006;Meliza and Dan, 2006;Feldman, 2012) and supported by computational models of cortical cell assembles (Wenisch et al., 2005;Klampfl and Maass, 2013). Only more recently has the importance of STDP at inhibitory output synapses, and at cPVin specific synapses, been appreciated (D'Amour and Froemke, 2015; Vickers et al., 2018). Computational modeling predicts that maturation of cPVins helps stabilize STDP such that the most temporally coherent inputs onto cortical pyramidal neurons are strengthened . Experimental data indicate that LTD of cPVin output onto principal neurons implements cortical disinhibition to regulate plasticity, and is necessary for remodeling the structure of sensory maps (Vickers et al., 2018).
Many studies have examined the specific operations that cPVins perform on sensory responses through optogenetic activation and/or inactivation of the population. Collectively, these studies have uncovered a panoply of conflicting transformations of sensory representations by cPVins. These include non-selective divisive gain control (Atallah et al., 2012;Wilson et al., 2012;Hamilton et al., 2013), receptive field sharpening Kaplan et al., 2016), temporal sharpening (Andermann et al., 2011;Vecchia et al., 2020), non-selective stimulus adaptation (Natan et al., 2015), and spatial constraint of neuronal ensemble activity (Agetsuma et al., 2018).
Broad inhibition of cortical output by cPVins derives from their ability to integrate excitatory activity across a wide area via their extensive dendritic fields (Poo and Isaacson, 2009;Sohal et al., 2009;Packer and Yuste, 2011;Karnani et al., 2014;Hage et al., 2022). This allows cPVins to exert gain control over cortical activity and to maximize the signal-to-noise ratio of information carried to downstream targets. Their intrinsic features also support their synchronized firing (Galarreta and Hestrin, 2001;Bartos et al., 2007;Ferguson and Gao, 2018). First, cPVins are highly sensitive to a large number of high amplitude, phase-locked inputs as a result of their Kv1 channel expression (Goldberg et al., 2008). Second, cPVins are coupled by gap-junctions (Fukuda and Kosaka, 2000;Tamás et al., 2000;Galarreta and Hestrin, 2001;Bartos et al., 2002). Third, cPVins maintain highly branched dendrites to access synaptic input over a sprawling area (Bartholome et al., 2020). In sum, these features facilitate the spread and synchronization of phasic, oscillatory activity across large brain regions.

cPVin CONTROL OF CORTICAL NETWORK ACTIVITY AND BEHAVIORAL STATE
At the broadest spatiotemporal scale, cPVin networks have a capacity to coordinate cortical output of large areas. On a population level, the rhythmic firing of cPVins establishes gamma oscillations, a high frequency, (30-80 Hz) periodic signal that results from synchrony of local field potentials and synaptic activity (Buzsáki and Draguhn, 2004). Genetic reduction of the number of cPVins attenuates gamma oscillations (Kalemaki et al., 2018), while optogenetically driving cPVin pools triggers gamma activity (Cardin et al., 2009;Sohal et al., 2009;Etter et al., 2019). Computational models suggest that gamma wave strength directly correlates with the expression of PV and GAD67 (Volman et al., 2011).
Gamma oscillations critically depend on NMDA receptor expression by cPVins (Korotkova et al., 2010;Gonzalez-Burgos and Lewis, 2012). In contrast to the fast spiking properties of cPVins, NMDA receptors activate synaptic currents with slow kinetics (Forsythe and Westbrook, 1988). The key contribution of NMDA receptors to the establishment of gamma rhythms appears to be their capacity to stabilize and coordinate recruitment of cPVins into synchronous ensembles (Carlén et al., 2012;Cornford et al., 2019).
The functional roles of gamma oscillations are less well understood. Oscillations reflect cPVin task-dependent activity and recruitment of greater numbers of cPVins during sensory learning (Gotts et al., 2012;Brunet et al., 2014;Ainsworth et al., 2016). Changes in gamma activity are critical for learning rules, such as during operant conditioning (Caroni, 2015;Lintas et al., 2021). Moreover, they are necessary for updating cortical responses when there is a mismatch between the history of stimuli-salience association and the new rules of the rules of sensory learning (Cho et al., 2020).
Further, the presence of gamma oscillations also correlates with the familiarity of sensory stimuli, which suggests they facilitate memory traces in sensory cortices (Womelsdorf et al., 2006;Headley and Weinberger, 2011;Weinberger et al., 2013;Cooke et al., 2015;Kissinger et al., 2018Kissinger et al., , 2020. In addition these oscillations may contribute to storing memories (Galuske et al., 2019). Gamma waves can also be triggered in response to meaningful behavioral outcomes (Fries et al., 2001;Kim et al., 2015;Ray and Maunsell, 2015;Cardin, 2016;Cho et al., 2020), which could imply their contribution to salience encoding.

cPVin CONTROL OF MATERNAL-EXPERIENCE INDUCED AUDITORY CORTICAL PLASTICITY
In mice, plasticity in the auditory cortex (ACtx) is triggered by maternal experience with pups. A broad base of evidence supports the conclusion that auditory plasticity in adult females helps to sharpen and amplify responses to behaviorally salient distress vocalizations from the pups. Here we review these data, highlighting work by our lab identifying cPVin-specific physiological mechanisms underlying "maternal, " experiencedependent auditory plasticity. These mechanisms share features of classical models of experience-dependent cortical plasticity, including those observed in developmental critical periods. Importantly, however, this plasticity is activated not by an exogenous experimental manipulation, but simply by social experience.
Our mechanistic understanding of maternal experienceinduced plasticity has been enhanced by our research with Mecp2 mutant mouse models. While wild-type (WT) maternal surrogates readily learn to perform pup retrieval, adult female mice that are missing one copy of the X chromosome -linked Mecp2 gene (Mecp2 HET ) fail to perform this task (Krishnan et al., 2017). We examined potential molecular mechanisms underlying this behavioral phenotype and found that exposure to pups initiated a transient increase in expression of GAD67 in the ACtx of both Mecp2 HET subjects and wild type littermates (Mecp2 WT ) (Krishnan et al., 2017). However, pup experience also leads to overexpression of parvalbumin protein (PV) and cPVinassociated PNNs in the ACtx of Mecp2 HET subjects, relative to naïve and Mecp2 WT controls (Krishnan et al., 2017). In the context of the literature discussed above, we hypothesized that the overexpression of PV and PNNs in the brains of Mecp2 HET surrogates inhibits the ACtx plasticity underlying retrieval learning, paralleling restriction of developmental critical periods. That is, while depression of these markers was not observed in WT subjects, their aberrant elevation in Mecp2 HET surrogates is sufficient for impairing behavioral responses to USVs. Therefore we surmise that a threshold level of PV and PNN expression, below which ACtx cells of Mecp2 WT controls maintain expression, allows for normal learning.
Further, we predicted that dissolution of PNNs might reinstate ACtx plasticity, as has been found in other sensory cortices. Indeed, PNN dissolution in the ACtx just prior to birth of the pups ameliorated overexpression of PNNs and facilitated retrieval performance in Mecp2 HET subjects (Krishnan et al., 2017). Again this suggests a cPVin specific mechanism or disruption that underlies the ability of Mecp2 HET subjects to respond behaviorally to USV cues. We further linked cPVins to maternal learning by showing that deletion of Mecp2 in cPVins was sufficient to disrupt early retrieval (Krishnan et al., 2017).
We next examined the effects of pup exposure on stimulusdriven activity in the ACtx of awake mice. Consistent with a central role of cPVin inhibitory drive in the pup-induced ACtx plasticity, we found that, in Mecp2 WT subjects, responses of cPVins to USVs were weaker after experience with pups (Lau et al., 2020). The diminished output of individual cPVins was mirrored by a complementary increase in the magnitude of responses from excitatory, putative pyramidal ACtx neurons (Lau et al., 2020). That is, maternal experience triggers a drop in cPVin responses to USV playback in the ACtx of Mecp2 WT . This is complemented by increased pyramidal neuron responses to USVs.
In contrast, these physiological changes were not observed in the ACtx of Mecp2 Het . More specifically, cPVin activity in response to USVs remained high despite experience with pups. These firing patterns matched the observed elevation in PV protein expression in Mecp2 Het , but not Mecp2 WT subjects. That is, given the activity-dependent expression of PV, the continued high stimulus-evoked firing rates of cPVins in Mecp2 Het links the cPVin-specific molecular expression to a functional lack of plasticity.
This pattern of results is also consistent with our understanding of sensory cortical plasticity in other systems, which is often initiated by disinhibition. In sum, inhibitory ACtx plasticity maximizes neural "contrast" and appears to align temporal dynamics of cortical response to USVs, enabling pup retrieval (Liu and Schreiner, 2007;Shepard et al., 2015). These findings are very much in keeping with the physiological properties of cPVin that allow them to rapidly regulate the timing of cortical output.
As others have suggested, studying sensory processing in natural and/or socially salient contexts may be particularly effective at triggering rapid or robust plasticity (Taborsky et al., 2015). Our work uses a maternal experience paradigm to evoke a natural and spontaneous form of cPVin-mediated sensory cortical plasticity in adult animals. However, there are many remaining questions to be answered. What does the population level structure of cPVin activity, including cortical gamma oscillatory activity, in response to USV calls look like, and how does experience with pups change those patterns? Do large scale oscillations reflect, or regulate, synchronous activity that favors synaptic modification in response to behavioral state? How does the activity of cPVins, and other inhibitory cortical populations, differ in freely behaving mice, as opposed to awake, head-fixed mice? What is the contribution, if any, of other cortical inhibitory subpopulations to maternal experience-induced auditory plasticity? How is cPVinmediated, sensory cortical plasticity modulated by multisensory integration? And how might it be modulated by behavioral states including attention, arousal, and emotional salience or valence?
Here we briefly summarize the shared features of plasticity during visual, developmental critical periods and during adult, auditory learning in natural contexts. First, both are marked by dynamic downregulation of cortical inhibition that results from decreased stimulus-evoked firing by cPVins, and leads to a concomitant increase in excitability of principal cortical neurons. This likely reflects changes to relative synaptic strengths from cPVins that lower the LTP threshold of pyramidal cells (Smith et al., 2009). Second, this cortical disinhibition is permissive for changes in connectivity, including the number of inputs by cPVins onto pyramidal cells. Third, periods of plasticity are delineated by fluctuations in the expression of parvalbumin protein and perineuronal nets, including prominent elevation of these markers as windows of heightened plasticity close. These molecular fluctuations may aid in mobility, or pruning, of synaptic inputs (Huang et al., 2015) that facilitate synaptic plasticity. Finally, in both cases, changes that reflect cortical plasticity do not preclude mechanisms of plasticity upstream of the cortex, e.g., changes to thalamic projections onto cPVins spurred by sensory experiences (Sommeijer et al., 2017). Nevertheless, the similarity in cortical plasticity features provides a general framework through which we can examine cPVin-specific contributions to sensory learning.

CONCLUSION AND FUTURE DIRECTIONS
Here we have focused on the known contributions of functional and structural cPVin plasticity for supporting sensory learning. However, the current understanding is far from complete. There are many areas of on-going investigation which promise to enhance this understanding including: contributions of cPVin to encoding natural or other complex stimuli (Zhu et al., 2015), real-time dynamics between cPVins and other cell types that make up canonical units (Rikhye et al., 2021), modulation of cPVins by top-down brain-state dynamics (Pakan et al., 2016), better dissection of transcriptomic (Fishell and Kepecs, 2020) and spatial (Large et al., 2018) profiles, and cPVin circuit specificity (Jiang et al., 2021).
We have also emphasized, by highlighting our own work, the need for further dissection of cell type specific plasticity in naturally occurring, socially salient contexts. Technical advancements, such as the development of less invasive approaches, are particularly important for implementing such investigations. For example, extracranial delivery of optogenetic  and ultrasonic stimulations (Estrada et al., 2021) as well as wireless recording capabilities (Shin et al., 2022) will make observation and manipulation of cell type specific populations during freely moving task engagement less cumbersome.
Advancing our understanding of cPVin plasticity has implications beyond scientific inquiry. Aberrant cPVin circuitry has been implicated in a diverse set of neurocognitive conditions (Huang, 2014;Ferguson and Gao, 2018;Filice et al., 2020). Therefore, a thorough understanding of cell type specific contributions to synaptic plasticity and network-wide oscillatory patterns has wide-spread implications for diagnostic and therapeutic advancement.

AUTHOR CONTRIBUTIONS
Both authors listed have made a substantial, direct, and intellectual contribution to the work, and approved it for publication.

FUNDING
This work was supported by a pre-doctoral fellowship from Autism Speaks (grant #11100) to DR and funding from NIMH (MH106656) and the Feil Foundation to SS. Funding agencies were not involved in the conceptualization, design, or writing of this article or in the decision to submit it for publication.