The Subtype-Free Average Causal Effect for Heterogeneous Disease Etiology

Studies have shown that the effect an exposure may have on a disease can vary for different subtypes of the same disease. However, existing approaches to estimate and compare these effects largely overlook causality. In this paper, we study the effect smoking may have on having colorectal cancer subtypes defined by a trait known as microsatellite instability (MSI). We use principal stratification to propose an alternative causal estimand, the Subtype-Free Average Causal Effect (SF-ACE). The SF-ACE is the causal effect of the exposure among those who would be free from other disease subtypes under any exposure level. We study non-parametric identification of the SF-ACE, and discuss different monotonicity assumptions, which are more nuanced than in the standard setting. As is often the case with principal stratum effects, the assumptions underlying the identification of the SF-ACE from the data are untestable and can be too strong. Therefore, we also develop sensitivity analysis methods that relax these assumptions. We present three different estimators, including a doubly-robust estimator, for the SF-ACE. We implement our methodology for data from two large cohorts to study the heterogeneity in the causal effect of smoking on colorectal cancer with respect to MSI subtypes.

Introduction
In recent decades, it has become clear that many diseases that share general clinical characteristics evolve through a range of heterogeneous molecular pathologic processes, which may be affected differently by the same exposure. For this reason, classifying a disease into subtypes according to criteria based on molecular characteristics can improve our etiologic understanding of the disease (Ogino et al., 2016).
In this paper, we consider the causal effect smoking may have on having colorectal cancer (CRC). One well-established CRC subtype classification is based on microsatellite instability (MSI). The two mutually-exclusive subtypes are MSI-high (approximately 15% of the cases) and non-MSI-high (approximately 85%). Smoking has been found to be a strong risk factor for the MSI-high subtype, while weaker evidence was found that smoking is a risk factor for the non-MSI-high subtype (Carr et al., 2018;Amitay et al., 2020).
To study etiologic heterogeneity, researchers often use a multinomial regression model in which being healthy or having each of the disease subtypes form the possible values of the outcome. However, the multinomial regression parameters do not correspond to well-defined causal effects, because the multinomial regression parameters are equivalent to the parameters obtained from a series of logistics regressions, each comparing one disease subtype to the healthy controls, resulting in selection bias (Nevo et al., 2021).
This form of bias is not limited to subtype comparisons and is not unique to multinomial regression. More generally, competing events create a challenge in making causal statements. Section 4.2 reviews related approaches and considers their applicability in our setup (Young et al., 2020;Stensrud et al., 2020;Nevo et al., 2021).
In this paper, we use data from two large cohorts to study well-defined and relevant estimands representing the effect of smoking on the two CRC subtypes. We propose an alternative estimand, inspired by the Survivor Average Causal Effect (SACE) (Zhang and Rubin, 2003;Rubin, 2006;Zehavi and Nevo, 2021). The SACE is a causal effect typicallyused to overcome the problem of truncation by death, which occurs when an outcome of interest is undefined for individuals who died before their outcome could be measured. The SACE is defined using principal stratification (Robins, 1986;Frangakis and Rubin, 2002) as the average causal effect of a treatment among individuals that would have survived under both treatment/exposure values. The SACE is not identifiable from the data under standard assumptions. Approaches to identify and estimate SACE include relying on additional strong assumptions, performing sensitivity analyses, and obtaining bounds.
Our proposed approach is based on the same idea as the SACE. First, we define principal strata with respect to the potential outcomes of other disease subtypes. Then, we define the Subtype-Free Average Causal Effect (SF -ACE) to be the average causal effect of the exposure on one disease subtype among the individuals who would have been free of the other MSI disease subtype under either exposure level (smoking status).
One commonly-made assumption in truncation-by-death problems is Monotonicity, namely that treatment/exposure cannot hurt survival (Zhang and Rubin, 2003). In the setting of disease heterogeneity with two subtypes, the definition of Monotonicity is more nuanced for two reasons. First, the assumption needs to be made (or not) for each subtype separately. Second, the disease could occur under both exposure statuses, but the subtype may differ. We term this phenomenon exposure-induced subtype switching and show it plays a key role in the (lack of) identification of the SF -ACE. We further show that unlike the SACE setting, the SF -ACE is non-parametrically identifiable under certain Monotonicity assumptions. Under weaker Monotonicity assumptions, the SF -ACE is no longer identifiable. For this scenario, we develop a sensitivity analysis approach as a function of the unidentifiable exposure-induced subtype switching probabilities.
We consider three SF -ACE estimators. The first two, standardization-based estimator and an inverse probability of treatment weighting (IPTW) estimator, may rely on parametric modeling assumptions. We therefore present a doubly-robust (DR) estimator that will be consistent if at least one of two models, but not necessarily both, is correctly specified.
The rest of the paper is organized as follows. Section 2 describes the epidemiologic problem and the available data. Section 3 presents the notations and the causal framework. Section 4 reviews existing approaches for etiologic disease heterogeneity and discuss available causal estimands. Section 5 presents a detailed study of the proposed causal estimand, the SF -ACE, including its interpretation, assumptions needed for its identification, sensitivity analyses and estimation. Section 6 summarizes the results of our simulation studies. Section 7 presents various analyses of the data to study the causal effect of smoking on the two CRC subtypes. Final conclusions are given in Section 8. The R package TheSFACE is available from CRAN implements our methodology. Reproducibility materials are available from https://github.com/amitSasson/SFACE Reproduce.

Data and problem description
We use data from two large US cohorts. The Nurses Health Study (NHS)  CRC diagnoses were reported by study participants. Formalin-fixed paraffin-embedded tumor tissue specimens were retrieved from hospitals across the US, and a CRC diagnosis was confirmed by the study pathologist (S.O.), who also conducted analyses to determine the MSI status. Further details can be found in Ugai et al. (2022) and references therein.
We set the baseline to be age 60, and define the exposure to be ever versus never smoking by age 60. Such binary exposure definitions are often used to minimize recall bias and measurement error. The outcome was CRC by age 70, subtyped by MSI status. Age 70 was chosen to minimize the impact of death as a competing event.
We removed from the data people who were older than 60 when enrolled in the study, or people who died, had CRC, or were lost to follow up before age 60. As in Ugai et al. (2022), we also removed people with missing baseline smoking data. In total, we were left with 114,947 people at baseline. Then, we removed 5,132 people (4.4%) who have died, without being diagnosed, prior to age 70. We were left with a final sample of 109,815. In this sample, 17,854 people were CRC-free and younger than 70 when the last questionnaire was collected (out of which 94.8% were older than 65). We treated them as CRC-free at age 70, and examined how it may affected our results by considering alternative analyses.
In our sample, the number of exposed (ever smoking) was 58,432 (53.3%) and the number of CRC cases by age 70 was 961 (0.9%). Out of them, 358 were diagnosed with non-MSI-high CRC (37.3%) and 61 were diagnosed with MSI-high CRC (6.3%). The subtype was missing for the remaining 542 (56.4%) diagnosed participants, as is often the case in such studies (Nevo et al., 2018;Liu et al., 2018). To have a clear presentation of the issues at end, we develop our theory and methodology under the assumption that subtypes are known for all disease cases. We address the issue of missing subtypes in our analysis in Section 7 using inverse probability weighting (Liu et al., 2018).
The NHS and HPFS cohorts include a vast number of covariates collected through biennial questionnaires. Table C.7 in the Appendix provides information on key baseline confounders we adjusted for in our analyses: gender (male/female), regular aspirin use (yes/no), first-degree family history of CRC (yes/no), history of lower endoscopy (yes/no), physical activity level (mean metabolic equivalent task score hours per week), body mass index (categorical), alcohol intake (g/day) and total calorie intake (kcal/day).

Preliminaries
Using the potential outcomes framework, let Y i (a) denote the potential disease status had the exposure of individual i = 1, ..., n been set to A i = a for a = 0, 1. In our study, A i = 0 indicates no smoking by age 60 (as captured by the cohort data), while A i = 1 indicates ever smoking by age 60. If individual i would have been diagnosed with disease subtype k = 1, ..., K under exposure value A i = a, then Y i (a) = k, and if they would be disease-free Note that because the subtypes are i (a) ≤ 1. In our study, K = 2, subtype 1 is non-MSI-high CRC, and subtype 2 is MSI-high CRC. So, for example, Y i (0) = 1 if person i would have been diagnosed with non-MSI-high CRC by age 70 had they not been smoking by age 60.
Throughout the paper, we assume there is no interference between individuals and that there are no multiple versions of the exposure leading to different potential outcomes.
Under these assumptions, known collectively as the stable unit treatment value assumption (SUTVA), the observed disease status for person i, hereafter denoted by Y i , with observed exposure value A i , equals to the appropriate potential outcome for that person. That is, . The latter is also known as the consistency assumption.
As in many observational studies, we also assume weak ignorability, namely that rich enough data were collected as part of the NHS and HPFS cohorts such that for a vector of measured covariates X, we have A |= Y (a) | X for a = 0, 1. Another typically needed assumption for identification of causal effects is positivity, meaning that the exposure assignment probability e(x) = Pr(A = a|X = x) is positive for a = 0, 1 and for all x.
In summary, the observed data consists of n i.i.d copies sampled from the distribution of (X, A, Y (1) , ..., Y (K) ), that can be equivalently represented as (X, A, Y ).

Review of existing approaches
A principled approach to causal inference typically starts with a definition of causal effects of interest, followed by the assumptions needed for their identification from observed data, and only then statistical models are discussed. However, as the existing approaches employed in practice for studying etiologic heterogeneity are model-based and are not grounded in the causal inference paradigm, we start with reviewing these models (Section 4.1). Then, we discuss relevant causal effects arising in the presence of competing events (Section 4.2).

Statistical models for disease heterogeneity
Let π k (a, x) = Pr(Y = k|A = a, X = x), k = 0, ..., K, be the probability of subtype k disease (or being subtype-free for k = 0) among those with exposure A = a and confounders X = x. Let also π(a, x) = (π 0 (a, x), ..., π K (a, x)). A commonly-used approach for studying etiologic heterogeneity is built upon the multinomial regression model for k = 1, ..., K and π 0 (a, x) = 1 − K k=1 π k (a, x). Under this model, exp(β k ) is the odds ratio (OR) between the exposure A and the observed subtype-specific outcome Y (k) compared with not having the disease, that is, exp(β k ) = π k (1,x) π 0 (1,x) × π 0 (0,x) π k (0,x) . The parameter exp(β k ) also approximates the risk ratio (RR) when the disease is rare, as in our study.
Model (1) was also extended for studies within which subtypes were not a-priori known and may be defined using multiple, possibly mis-meausred, biomarkers (Chatterjee, 2004;Wang et al., 2015;Nevo et al., 2016). See also Wang et al. (2016); Zabor and Begg (2017) for a survey of and comparison between methods based on this and similar models.
Estimation of Model (1) parameters is typically carried out by standard maximum likelihood approach, with some exceptions when further complications arise (e.g., Chatterjee, 2004;Nevo et al., 2016). In practice, researchers are often interested in whether there is heterogeneity in the effects of A across disease subtypes. When the null hypothesis H 0 : β 1 = ... = β K is rejected, researchers conclude that the exposure A contributes to the risk of at least one subtype differently from its contribution to the risk of other subtypes.
A recently-proposed alternative approach postulates K logistic regressions, each comparing a subtype k disease to the alternative of being free of that subtype, that is, being either disease-free or diagnosed with the disease, but of a different subtype (Sun et al., 2017). Formally, these models can be expressed by for k = 1, ..., K. Under this model, exp( β k ) is the OR between the exposure A and the observed indicator of subtype-specific outcome Because fitting these models separately may result in K k=0 π k = 1, estimation of the parameters in (2) calls for constrained maximization of the models' likelihood. A more computationally attractive approach, developed by Sun et al. (2017), uses constrained Bayesian estimation, and limits the posterior distribution of the parameter vector to have positive mass only for parameter values such that K k=0 π k = 1.
Comparing models (1) and (2), it has been recently argued (Sun et al., 2017) that variation in β k across the subtypes (across k) is not a satisfactory measure of etiologic heterogeneity, because the parameter space of a multinomial regression model requires a risk factor for one subtype to be a risk factor for the other subtypes, and that β k may be preferred, although this issue was also debated (Begg et al., 2018;Sun et al., 2018).

Causal effects for competing risks
To fix ideas, we focus on the effect of A on Y (1) , the first subtype, and start with a naïve approach that does not correspond to a causal effect. Researchers may opt for a comparison of the subtype 1 proportion of cases at each exposure level only among individuals who are free of other disease subtypes. When researchers also adjust for the confounders X, their statistical estimand on the difference scale is and can be analogously defined on the RR or OR scales. Because model (1) is equivalent to K separate logistic regression models, each comparing one disease subtype to the disease-free individuals (Wang et al., 2016), then exp(β k ) is the OR analogue of (3), and approximates the RR as CRC is a rare disease. However, (3) does not correspond to a causal effect, because of selection bias due to common causes of the different subtypes (Nevo et al., 2021). These common causes cannot be classified as classical confounders, as they can be independent of the exposure. If these common causes are also mediators of an exposure-subtype mechanisms that are shared between more than one subtype, even had they been observed, adjusting for them would have eliminated part of the causal effect of the exposure on the subtype (Nevo et al., 2021). Importantly, because MSI-high and non-MSI-high are subtypes of the same disease, CRC, they are likely to have unknown and/or unmeasured common causes. Young et al. (2020) considered possible causal estimands in the presence of competing events. The total effect, T E (1) = E[Y 1 (1) − Y 1 (0)], contrasts the population-level subtype 1 risk had the entire population been exposed, with the population-level subtype 1 risk had the entire population were unexposed. In the presence of measured confounders, this effect is identified from the observed data under SUTVA, weak ignorability and positivity by Unlike (3) is sensitive to the phenomenon of exposure-induced subtype switching. For example, assume some individuals would have been diagnosed with subtype 1 under A = 0 and would have been diagnosed with subtype 2 under A = 1, while the rest of the population will be free of subtype 1 for all A values. Then, the total effect on subtype 1 will be negative, T E (1) < 0, and it would appear as if the exposure protects from subtype 1, even though the individuals "protected" by A would still have the disease, just of another subtype.
The direct effect (Young et al., 2020) is the causal effect comparing subtype 1 rates under two hypothetical joint interventions. The first sets the exposure to A = 1 and eliminates the possibility of disease subtypes other than subtype 1 disease, and the second sets the exposure to A = 0 and also eliminates risks of other disease subtypes. This approach is, however, irrelevant for our study, because an intervention eliminating the risk of a particular CRC subtype is not currently available and is not expected to be in the near future.
Another alternative considers the separable effects (Stensrud et al., 2020), which rely on a conceptual separation of the exposure into the components affecting each disease subtype.
Such separation requires profound understanding of the mechanisms leading to each CRC subtype. Furthermore, the conditions needed for identifiability of these effects are violated if there are unobserved common causes of both disease subtypes (Stensrud et al., 2020).
Such common causes are expected to exist in our study and in many disease subtype studies (Nevo et al., 2021).

The Subtype-Free Average Causal Effect
Our focus in this paper is on the common scenario of two disease subtypes (K = 2), as in our study. We propose a new estimand, inspired by the SACE (Rubin, 2006). We define the SF -ACE on the difference scale by

SF -ACE
(1) to be the SF -ACE on the RR scale. The SF -ACE (1) is the average causal effect among the individuals that would have been free of disease subtype 2 under either exposure level (i.e., regardless of smoking status). Similarly to SF -ACE (1) , let

SF -ACE
(2) be the average causal effect on the difference scale of A on Y (2) among the individuals that would have been free of disease subtype 1 under either exposure level. Let also SF -ACE

Interpretation of the SF -ACE
As explained in Section 4.2, competing events present a challenge in defining causal estimands that shed light on the scientific question of interest. The SF -ACE presents another source of information. A positive SF -ACE (1) D value means that in the sub-population who would have been free of MSI-high CRC regardless of their smoking status, we will expect more non-MSI-high CRC cases had the entire sub-population smoked compared to the scenario the entire sub-population did not smoke before age 60.
Therefore, the SF -ACE offers a way to circumvent non-zero causal effects due to exposure-induced subtype switching. Nevertheless, the above appealing advantage does not come without a price. One disadvantage is that the SF -ACE is defined in a latent subset of the population. This is a general problem with principal stratification approaches.
A second, related, issue is that when contrasting SF -ACE (1) and SF -ACE (2) , we contrast effects defined in two overlapping, but not identical, sub-populations. The SF -ACE (1) is defined within the sub-population {Y (2) (0) = Y (2) (1) = 0}, and the SF -ACE (2) is defined Thus, a non-zero θ D captures changes that might be due to the different population, and as a result, the exact source of heterogeneity in the causal effect is not revealed. Instead, a non-zero θ D presents evidence of heterogeneity and calls for further investigation. Lastly but importantly, in rare disease scenarios, such as CRC, this problem is less concerning because the two sub-populations {Y (2) (0) = Y (2) (1) = 0} and {Y (1) (0) = Y (1) (1) = 0} have almost perfect overlap, as those who would not have CRC of any subtype under any exposure value belong to both groups and comprise the vast majority of the population.

Identifiability
Similar to the SACE in truncation-by-death scenarios, the SF -ACE is not identifiable from the data under the standard assumptions of SUTVA, weak ignorability and positivity. One commonly-made assumption when targeting the SACE is Monotonicity (Zhang and Rubin, 2003), stating that treatment cannot hurt survival. In our study, the definition and plausibility of Monotonicity-like assumptions are more nuanced, and, as it turns out, resulting in new insights regarding principal causal effects and their identification in a competing risks setup. We start with the definition of Subtype Monotonicity (S-Monotonicity). It is clear from the above that the SF -ACE RR is identifiable from the data under SUTVA, weak ignorability, positivity, and S-Monotonicity for both subtypes. The following proposition establishes identification of SF -ACE D , and presents identification formulas under either standardization or IPTW (Hernán and Robins, 2020).

Assumption 1 Subtype Monotonicity (S-Monotonicity
Proposition 1 Under SUTVA, weak ignorability, positivity, and S-Monotonicity for sub- using standardization, or via IPTW by The proofs and analogous expressions for SF -ACE (2) , θ D , and for causal contrasts on the RR scale are all given in Section A of the Appendix. Estimation based on these identification formulas is discussed in Section 5.5.

Relaxing S-Monotonicity
Assuming S-Monotonicity for both subtypes might be too restrictive, even if we believe the exposure cannot protect from the outcome, as S-Monotonicity disallows subtype switching for any individual in the population. For example, consider the scenario the exposure is believed to increase non-MSI-high CRC risk but not MSI-high CRC risk. If individual i would have been diagnosed with MSI-high CRC when unexposed, then under S-Monotonicity for MSI-high CRC, it will be impossible for individual i to be diagnosed with non-MSI-high CRC when exposed. This is even though the exposure is a risk factor only for non-MSI-high CRC.
Therefore, we consider a weaker assumption, termed Disease Monotonicity (D-Monotonicity), that states that an individual who would have been diagnosed with disease subtype k under no exposure, would have also been diagnosed with the disease under exposure, but not necessarily with the same subtype.
D-Monotonicity keeps the premise that the exposure cannot protect from the disease while being less restrictive than S-Monotonicity. As with S-Monotonicity, D-Monotonicity can be modified for scenarios where the exposure is suspected to have a protective effect.
For each disease subtype, we may assume S-Monotonicity, D-Monotonicity, or neither.
With two disease subtypes, we have eight different combinations of possible assumptions.

Under each combination, certain potential outcome profiles of {Y
(1) i (1)} are assumed to be absent from the population. For example, under S-Monotonicity for subtype 1, there is no one in the population with {Y (1) i (1)} = {1, 0, 0, 1}. Note also that because the subtypes are mutually exclusive, it is impossible to i (1) = 1 regardless of the taken monotonicity assumptions. Table 1 presents the nine possible potential outcome profiles, and specifies whether each profile exists under four different combinations of assumptions. This table demonstrates that the exact implications of the different monotonicity assumptions on the potential outcomes in the population are quite explicit. The observed data, however, does not hold all the information required for classifying each individual into their potential outcome profile. Table A.1 in Section A of the Appendix specifies the possible profiles for each participant according to their observed (A, Y (1) , Y (2) ) under each of the four combinations of assumptions.

Potential outcomes profiles
Monotonicity assumptions

A sensitivity analysis approach
By Proposition 1, under S-Monotonicity for both disease subtypes, the SF -ACE (1) , SF -ACE (2) and θ are identifiable from the observed data distribution. If we replace S-Monotonicity with D-Monotonicity or with no monotonicity assumption for at least one subtype, then the causal effects are no longer identifiable, and SF -ACE Therefore, we take a sensitivity analysis approach, as is common when studying causal effects under untestable assumptions. Let λ 1 = Pr[Y (2) (1) = 1|Y (1) (0) = 1] be the proportion of individuals who would have been diagnosed with disease subtype 2 when exposed, out of the individuals who would have been diagnosed with disease subtype 1 when unexposed. That is, λ 1 is the probability of exposure-induced subtype switching from subtype 1 to subtype 2. Note that λ 1 = 0 under S-Monotonicity for subtype 1 and λ 1 ≥ 0 under D-Monotonicity for subtype 1. Let also λ 2 = Pr[Y (1) (1) = 1|Y (2) (0) = 1] be the subtype-switching probability from subtype 2 to subtype 1. While λ 1 and λ 2 are not identifiable from the observed data, they can be bounded. Specifically, because The following two propositions demonstrate that the subtype-switching probabilities play a key role in the identification of the SF -ACE whenever the S-Monotonicity assumptions are not met. Proposition 2 shows that (a) under D-Monotonicity for subtype 1 and S-Monotonicity for subtype 2, both the SF -ACE

Proposition 2 Under SUTVA, weak ignorability, positivity and D-Monotonicity for both subtypes, the SF -ACE
(1) D is identifiable from the data by using standardization, or via IPTW by .
The identification of SF -ACE (2) D is analogous.
The proof is given in Section A of the Appendix. Of note is that if one further assumes S-Monotonicity for one of the subtypes (say subtype 2), then an identification formula of

SF -ACE
(1) D is obtained by setting the appropriate λ k (say λ 2 ) to zero in (8) and (9). A second remark is that, technically, results (8) and (9) require D-Monotonicity for subtype 2 only, and similarly, the identification for subtype 2 requires D-Monotonicity for subtype 1 only. In Section A of the Appendix, we present further identification results for SF -ACE D under D-Monotonicity for one subtype only or without any monotonicity assumptions, and show that in these cases λ 1 and λ 2 do not suffice for identification of SF -ACE D .
Turning to the RR scale, Proposition 3 asserts that SF -ACE RR is identifiable as a function of λ 1 and λ 2 and the observed data without imposing any monotonicity assumptions.

The identification of SF -ACE
(2) RR is analogous.
The proof is given in Section A of the Appendix.
Propositions 2 and 3 serve as basis for our proposed sensitivity analyses. Researchers first decide what is the effect scale of interest, and then for which if any of the two subtypes the S-Monotonicity or D-Monotonicity assumptions are plausible. Then, depending on the taken assumptions, they may set λ 1 or λ 2 to zero, and estimate the SF -ACE (k) and θ as a function of λ 1 and/or λ 2 . Such an analysis can reveal how sensitive the results are to different levels of exposure-induced subtype switching. We implement these analyses in our study of smoking effect on CRC MSI subtypes in Section 7.

Estimation
Recall that for each individual i, the observed data are ( i ).
Estimation of SF -ACE (k) and θ from the data alongside the sensitivity analyses for the S-Monotonicity assumptions require estimation of π k (a, x), k = 1, 2, a = 0, 1 for the standardization-based formulas, or e(x) for the IPTW-based identification formulas. Expectations over X or (X, A, Y (k) ) are estimated by taking the sample analogues.
Non-parametric estimation of π k (a, x) or e(x) when the vector of measured confounders X is of high dimension or includes many continuous variables could be a difficult task, that will yield an estimator with a too large variance to be able to derive conclusions from the data. Therefore, parametric estimation might be a better option, and also aligns with the models used in practice and that were presented in Section 4.1. We first discuss parametric estimation, before implementing a DR estimator that will be consistent even if not all models are correctly-specified.
Starting with standardization, let π(a, x; η 1 ) = (π 0 (a, x; η 1 ), π 1 (a, x; η 1 ), π 2 (a, x; η 1 )) be the parametric model for π(a, X) with parameter vector η 1 . This could be, for example, the regression models (1) or (2). Upon estimatingη 1 , the estimated probabilities π(a, X i ;η 1 ) are calculated in the entire sample, for a = 0, 1, before plugged in the relevant identification formulas. For example, the proposed SF -ACE If there is an unmeasured covariate U affecting both disease subtypes (Y (1) , Y (2) ), model (1) might be misspecified, if one believes the model that is also conditioned on U is a multinomial regression model. However, this is not a main concern in studies involving rare disease, as in our study (Nevo et al., 2021).
Turning to IPTW, let e(x; η 2 ) be a parametric model for e(x), with parameter vector η 2 , for example a logistic regression. Then, after estimatingη 2 , the SF -ACE (1) D estimator is obtained by and similar estimators for SF -ACE (2) D and θ D , or for effects on other scales, can be obtained.
When using the maximum likelihood approach to fit models π(a, x; η 1 ) or e(x, η 2 ), the obtained estimators η 1 or η 2 are consistent and asymptotically normal. Thus, the continuous mapping theorem and the delta method can be used to show that the SF -ACE and θ estimators are also consistent and asymptotically normal. In practice, the bootstrap can be used to estimate the standard errors (SEs) and construct confidence intervals (CIs).
The standardization estimator or the IPTW estimator will be consistent only if their respective models are correctly specified. Therefore, we extend our estimation approach to DR estimation (Robins et al., 1994;Lunceford and Davidian, 2004). That is, we construct an estimator that will be consistent whenever at least one of the two models, π(a, x; η 1 ) or e(x; η 2 ), but not necessarily both, is correctly specified. DR estimators are frequently used in causal inference to reduce dependency on the chosen model. Furthermore, a DR estimator constructed as an augmented-IPTW estimator is also more efficient than IPTW estimators and have desirable theoretical properties (Robins et al., 1994).
Letη 1 andη 2 be the estimators obtained for η 1 and η 2 . As can be seen from the proof of Proposition 1 (Section A of the Appendix), under S-Monotonicity, the SF -ACEs and . Under weak ignorability and positivity, a DR augmented-IPTW (Bang and Robins, 2005;Lunceford and Davidian, 2004). The estimator (12) and analogous and similar equations to obtain estimators for the causal effects that are consistent by the continuous mapping theorem, if at least one of the models, π(a, x; η 1 ) or e(x; η 2 ), is correctly specified. Similar arguments can be used to construct DR estimators when conducting sensitivity analyses based on Proposition 2 and 3.

Simulations
We conducted simulation studies to evaluate the finite-sample performance of our proposed estimators and to compare the SF -ACE to existing estimands. Detailed information about the data generating mechanism and parameters used, as well as more detailed results are given in Section B of the Appendix. For each i = 1, ...n, potential outcomes Y i (a) were simulated using a multinomial regression model with two confounders, X 1i ∼ Ber(0.5), X 2i ∼ N (0, 1), and one unobserved covariate U i ∼ N (0, 1) affecting both subtypes. S-Monotonicity was imposed for both subtypes (without violating the multinomial model assumption, see Section B). The exposure A i was simulated using a logistic regression with To study different aspects of the proposed approach, we conducted three studies. In the first, we fixed the model parameters and varied the sample size between 5,000 and 50,000, to study the methods' performance in different sample sizes. In the second study, we fixed the sample size and varied the effect the unmeasured covariate U i had on Y (2) i to compare between the different causal effects. In the third study, we investigated the methods' performance when the fitted models were misspecified.
In each simulated dataset, we first fitted the multinomial regression (1), and a logistic regression for e(x). Then, we calculated the conditional estimators (3), the TE estimators (4), and the three SF -ACE estimators described in Section 5.5. We used the bootstrap with 200 repetitions to estimate the SEs and to construct Wald-type 95% CIs.
In the first study (Tables B.3 and B.4), the SF -ACE estimators had a small relative bias. The SE of all estimators were similar, and were generally well-estimated. Empirical coverage rates of the CIs were satisfactory. These results were also apparent for the second study ( Figures B.1 and B.2). In the second study, as the effect of U on Y (2) increased, absolute bias of the conditional estimand and the TE on the difference scale, compared to the true SF -ACEs, increased, while the SF -ACE estimators remained unbiased. As expected, because data were simulated under S-Monotonicity for both subtypes, the T E

Data analysis
We conducted several analyses to study the causal effects of ever smoking by age 60 on having CRC MSI subtypes by age 70, using the NHS and HPFS datasets described in Section 2. The main results are presented here, and additional results are given in Section C of the Appendix. In all of our analyses, we used the confounders detailed in Section 2 and Table C.7. Wald-type 95% CIs and p-values were calculated using the bootstrap with 200 repetitions.
As indicated in Section 2, for more than 50% of CRC cases the subtype status was unknown. To minimize the impact of potential selection bias, we used inverse probability weighting for missing subtype (Liu et al., 2018). We first fitted a model for the probability of observing the subtypes (Table C.8). The model included additional variables available for CRC cases (e.g., tumor stage and size). Then, for each CRC case with known subtype, we calculated the weight as the reciprocal of the probability of observing the subtype. The weights were truncated at the 99% percentile and then incorporated into our estimators.
We estimated the TE of smoking on each disease subtype by standardization (4), using a multinomial regression model (1) (Table C.9). The point estimates indicated smoking increases the risk of both CRC MSI subtypes (Tables 2 and C.10), although the results were not significant at the 5% level. The point estimates indicated that the effect on non-MSI-high subtype was larger than the effect on the MSI-high subtype on the difference scale and vice versa on the RR scale. These results are in agreement with the existing literature, because MSI-high is more rare, so while the RR is larger, the number of excess cases due to smoking can be smaller.   The S-Monotonicity assumption might be too restrictive, especially for the non-MSIhigh subtype, for which there is limited information about the relative effect strength (Amitay et al., 2020). From a subject-manner perspective, we could not overlook the possibility of exposure-induced subtype switching for non-MSI-high subtype, as the relative effect of smoking on MSI-high CRC is known to be larger. Therefore, we replaced the S-Monotonicity assumption for non-MSI-high (subtype 1) with D-Monotonicity by fixing λ 2 = 0 and varying the values of λ 1 between 0 and (Figure 1). As λ 1 increased, the estimated SF -ACE for non-MSI-high CRC increased, the estimated SF -ACE for MSI-high CRC decreased. For λ 1 ≥ 0.04, H 0 : θ D = 0 was rejected. This result serves as evidence for heterogeneity in the SF -ACE as long as we believe subtype switching is possible and not extremely rare. Given some a-priori evidence for smoking effects on both subtypes, we also considered an analysis under D-Monotonicity for both subtypes. Figure 2 presents the estimated effects on the difference scale as a function of λ 1 and λ 2 (for which the data did not restrict the range), along with 95% CIs. For nearly the entire range, the effect on non-MSI-high CRC was significant and the effect on MSI-high CRC was not. Sensitivity analysis results were similar for all three estimators (Figures 2, C.4 and C.5). As λ 2 increased, the larger was the minimal value of λ 1 for which heterogeneity can be concluded. For the RR (Figures C.6-C.8), under the combination of low λ 2 and large λ 1 the heterogeneity assumption was rejected.
As described in Section 2, there were 17,854 study participants for which the last available data were before age 70 (and they were not diagnosed with CRC and did not die).
We treated these people as CRC-free at age 70. To assess the impact of this decision, we repeated our analyses while removing these people. In an additional analysis, we again removed these people but used inverse probability of weighting to take them into account.
The results did not change materially, although the estimated effect sizes and the standard errors were generally larger (Tables C.13 and C.14).
To summarize our main findings, we have provided evidence for smoking effect on both CRC MSI subtypes and evidence for heterogeneity under reasonable causal assumptions.
Even though our case numbers were smaller and our follow-up was more limited than previous studies, the estimated SF -ACE lend support of a causal effect on non-MSI-high CRC, for which available evidence was limited.

Discussion
Questions of etiologic heterogeneity are central to our understanding of how diseases evolve.
Viewed through the lenses of competing events, existing approaches largely overlooked formal discussions of causal estimands and identifiability assumptions. On the other hand, the recent progress in causal inference methodology for competing events presented causal estimands that either provide little information about etiology heterogeneity or estimands requiring identifiability assumptions that are implausible in the study of CRC subtypes.
In this paper, to study the effect of smoking on CRC subtypes defined by MSI status, we proposed the SF-ACEs as causal effects that can complement existing methods and provide further insight towards the study of etiologic heterogeneity. One challenge with using multinomial regression with cohort data is the issue of right censoring due to loss to follow-up or administrative censoring. To assess its impact, we considered in this paper a number of analyses. The censoring problem would have been more pronounced if follow-up was longer than 10 years. Therefore, considering time-to-event outcomes is an important topic of future research. Another potential extension will consider continuous exposures instead of binary ones. Specifically, the generalization of the monotonicity assumptions, both in terms of definition and plausibility, is of interest. • Section A includes proofs and additional theory. Section A.1 presents the proofs of Propositions 1-3. Section A.2 illustrates the the connection between observed data and the potential outcomes. Section A.3 discuss identification results under weaker combination of assumptions. Section A.4 provides the complete expression for the DR estimator.
• Section B includes details on the simulations and presents the results. Section B.1 gives additional details about the data generating mechanism of the simulation. Section B.2 presents the details of three different simulation studies. Finally, Section B.3 presents tables and figures summarizing the different simulation results.
• Section C presents details on the data as well as additional analyses and results. Section C.1 includes descriptive information about the data used. Section C.2 gives additional information about the weights used for handling missing subtypes and about the models fitted for estimating effects. Finally, Section C.3 presents additional data analysis results.

A.1 Proof of identification results
Before providing the identification formulas for the SF -ACE, we first show that the SF -ACE equals to the T E when defined on the RR scale under S-Monotonicity for both disease subtypes.
When defined on the RR scale, the SF -ACE (1) equals to (A.14) Equality (A.13) holds by the S-Monotonicity assumption for subtype 2. Equality (A.14) holds because the subtypes are mutually exclusive (numerator) and by the S-Monotonicity assumption for subtype 1 (denominator).

A.1.2 Proof of Proposition 1
Starting from the definition of SF -ACE, we have

SF -ACE
(1) For the standardization-based identification formula, we continue by recalling that under positivity and weak ignorability, Pr[Y (k) (a) = j] is identified by for a, j, k = 0, 1 (Hernán and Robins, 2020). With Pr(Y (k) = j|A = a, x) previously denoted by π k (a, x) for j = 1 and hence equals to 1 − π k (a, x) for j = 0, we substitute these expressions in (A.16) to get

(A.18)
The identification formulas for subtype 2 are analogous. Using standardization,
The difference between the effects is given by substituting the expressions above in θ D =
Finally, define θ RR to be the difference between the effects in the RR scale. So θ RR is given by

A.1.3 Proof of Proposition 2
We start with the effect on subtype 1.

SF -ACE
In (A.19), the numerator follows from the fact the subtypes are mutually exclusive and the denominator is justified by the law of total probability. In (A.20), the numerator holds because of the law of total probability and the denominator follows by noting that the D-Monotonicity assumption for subtype 2 implies Pr[Y (2) (0) = 1, Now, similarly to the proof of Proposition 1, we can continue using standardization, by weak ignorability and positivity. Alternatively, as in the proof of Proposition 1, we can continue from (A.21) using IPTW to obtain

SF -ACE
(1) , again by weak ignorability and positivity, similarly to the standard IPTW-based identification.
The identification formulas for subtype 2 are analogous. Using standardization,
The difference between the effects is given by

SF -ACE
(1) Equality (A.22) holds because the subtypes are mutually exclusive. Equality (A.23) is by the law of total probability. As before, using standardization, we may write
The identification formula of subtype 2 is analogous. Using standardization,

SF -ACE
(2) The difference between the effects is given by θ RR = SF -ACE (1) RR .

An interesting distinction between the results in Propositions 2 and 3 compared to
Proposition 1, is that if one further assumes S-Monotonicity for one of the subtypes only, such that S-Monotonicity holds for one subtype only, say subtype 1 so λ 1 = 0, then the identification formulas for the two subtypes are no longer symmetric. Table 1 of the main text defines the nine possible potential outcome profiles. Table A Observed data Monotonicity assumptions (1, 0, 1) 1, 6 1,6,7 1,6 1,6,7 Table A.3: Connection between observed data and possible potential outcome profiles defined in Table 1 under four combinations of monotonicity assumptions.

A.3.1 S-Monotonicity for one subtype and no monotonicity assumption for the other subtype
We present here identification for the SF -ACE (k) D under S-Monotonicity for subtype 1 and no monotonicity assumption, of any version, for subtype 2. To this end, in addition to λ 2 , we define an additional sensitivity parameter λ 0 This sensitivity parameter is the probability to be CRC-free when exposed among those who would have been diagnosed with subtype 2 had they were unexposed. Under either S-Monotonicity or D-Monotonicity, λ 0 2 = 0.
We begin with the effect on subtype 1,

SF -ACE
(1) Because different assumptions were made for subtype 1 and subtype 2, the effects are not analogous. Turning to the effect on subtype 2,

SF -ACE
(2) Equality (A.27) holds because of the S-Monotonicity assumption for subtype 1. Equality (A.28) holds because the subtypes are mutually exclusive. Equality (A.29) holds because of the law of total probability. As before, standardization, IPTW or DR estimation can be used to identify and estimate the effects.

SF -ACE
(1) Equality (A.30) holds because the subtypes are mutually exclusive. Equality (A.31) holds because of the law of total probability. Equality (A.32) holds because of the law of total probability and because the subtypes are mutually exclusive.

SF -ACE
(2) As before, standardization, IPTW or DR estimation can be used to identify and estimate the effects.

A.4 Details about the doubly-robust estimator
The doubly-robust estimators that are plugged in Equations (A.16) (under S-Monotonicity for both subtypes), (A.21) (under D-Monotonicity for both subtypes), or their RR ana-logues are given by,

B Additional simulation details and results
Section B.1 gives additional details about the data generating mechanism of the simulation. Section B.2 presents the details of three different simulation studies. Finally, Section B.3 presents tables and figures summarizing the different simulation results.

B.1 Data generating mechanism
The data were simulated to resemble an observational study. For each individual i, we first simulated two independent measured confounders, X i = (X i1 , X i2 ), where X 1 ∼ N (0, 1) and X 2 ∼ Ber(0.5) and one normally distributed unmeasured covariate U i ∼ N (0, 1). To simulate the potential outcomes, we used the multinomial regression model for k = 1, 2, and We first simulated the potential outcome under no exposure Y i (0), for each individual i, using model (A.33). To make sure the data followed the S-Monotonicity assumption for both disease subtypes, the potential outcomes under exposure (a = 1) were created as follows: Individuals with Y i (0) > 0, namely that they would have been diagnosed with disease subtype 1 or 2 under a = 0, would automatically be diagnosed with the same disease subtype under a = 1. For individuals who were disease-free under a = 0 the disease subtype under a = 1 was simulated using Equation (A.33) with an adjusted probability using the Law of Total Probability, We simulated the actual exposure status for each individual A i using a logistic regression model, .
The observed outcome Y was obtained for each individual by setting For each individual i, the observed data were (A i , X i , Y i ).

B.2 Additional details about the simulation studies
As indicated in the main text, in each of the studies and scenarios described below, we first estimated π(a, x; η 1 ) using a multinomial regression model and e(x) using a logistic regression model. Then, we computed the estimators for the conditional statistical estimand, the total effect, and the three SF -ACE estimators.
Three studies were investigated. In Study I, we fixed the model parameters and varied the sample size between 5,000 and 50,000. In Study II, we varied δ 2 , the effect of the unmeasured covariate U on Y (2) , while keeping all the other parameters fixed. In Study III, we investigated the estimators' performance when the models fitted were misspecified. To achieve this, a transformation of X 2 was used to create the data, but the analysis assumed X 2 is (logit) linearly connected to the exposure or/and the outcome. Misspecification of the A ∼ X relationship was created by taking log |X 2 | when generating the data instead of X 2 . Misspecification of the Y ∼ A + X model was also created by replacing X 2 with log |X 2 | when generating the data. Table B.4 gives the model parameters and sample sizes used in the three studies. It also gives the obtained true SF -ACE values for each set of parameters yields. In all studies, the prevalence of subtype 1 ranged between 5%-6% and the prevalence of subtype 2 ranged between 2%-4%. The prevalence of the exposure was 49%.

B.3 Simulation Results
Tables B.5 and B.6 present the result for Study I. As can be seen from the Tables, the SF -ACE estimators had a smaller relative bias than the other estimators, even with a small sample size.
The absolute difference between non-SF -ACE estimators and the SF -ACEs on the difference T E and SF -ACE coincide on the RR scale. The relative bias ( Figure B.4) remained generally the same. As all three methods (standardization, IPTW, and DR estimation) were based on correctly-specified models, their finite sample bias was minimal.  Finally, Table B.7 presents results from Study III under four settings: (a) Both π(a, x; η 1 ) and e(x; η 2 ) are correctly specified, (b) The π(a, x; η 1 ) model is correctly-specified and model e(x; η 2 ) is misspecified, (c) The π(a, x; η 1 ) model is misspecified and model e(x; η 2 ) is correctlyspecified, and (d) both models are misspecified. As expected the standardization-based estimator had minimal bias when π(a, x; η 1 ) was correctly specified, and was biased when the model was misspecified. The IPTW-based model had minimal bias whenever model e(x; η 2 ) was correctly specified and produced considerable bias when the model was misspecified. As expected the DR estimator was only biased when both models were misspecified and had minimal bias under scenarios (a), (b), and (c) described above.

C Details on the data analysis and additional results
Section C.1 presents details on the data as well as additional analyses and results. Section C.1 includes descriptive information about the data used. Section C.2 gives additional information about the weights used for handling missing subtypes and about the models fitted for estimating effects. Finally, Section C.3 presents additional data analysis results. Table C.9 gives basic descriptive information about the different covariates X used in our analyses.  (1976-1995, 1996-2000, 2001-2012) and family history of colorectal cancer (yes or no). The estimates and 95% CI of the model coefficients are given in Table C.10. After fitting the model, we removed the individuals who were diagnosed with CRC but their MSI subtype was missing, and used weights based on the inverse of the logistic regression predictions to adjust our sample. We truncated the weights at the 99% percentile to reduce variance. The distribution of the final weights used are given in

C.3 Additional results
Turning to the analysis results. First, we estimated the effects under S-Monotonicity for both disease subtypes. The results for the SF -ACE D are given in Table 2 of the main text. Here, we give the results for the SF -ACE RR in Table C.13, and the results for θ in Table C.14.
We then continued to the analysis under D-Monotonicity for both disease subtypes. The results for the SF -ACE D and θ D using the DR estimator were given in Figures 2 and 3 of the main text. Figures C.6 and C.7 give the corresponding results for standardization and IPTW methods. (1) RR estimator was relatively stable while SF -ACE (2) RR increased as a function of λ 2 and decreased as a function of λ 1 .
Finally, we present the results of an analysis after removing participants for which the last available data were before age 70, with and without inverse probability weighting to take them into account (Tables C.15


s2fieldsofstudy: Biology Chemistry

text:
Bioconversion of vitamin D3 to bioactive calcifediol and calcitriol as high-value compounds

Biological catalysis is an important approach for the production of high-value-added compounds, especially for products with complex structures. Limited by the complex steps of chemical synthesis and low yields, the bioconversion of vitamin D3 (VD3) to calcifediol and calcitriol, which are natural steroid products with high added value and significantly higher biological activity compared to VD3, is probably the most promising strategy for calcifediol and calcitriol production, and can be used as an alternative method for chemical synthesis. The conversion efficiency of VD3 to calcifediol and calcitriol has continued to rise in the past few decades with the help of several different VD3 hydroxylases, mostly cytochrome P450s (CYPs), and newly isolated strains. The production of calcifediol and calcitriol can be systematically increased in different ways. Specific CYPs and steroid C25 dehydrogenase (S25DH), as VD3 hydroxylases, are capable of converting VD3 to calcifediol and calcitriol. Some isolated actinomycetes have also been exploited for fermentative production of calcifediol and calcitriol, although the VD3 hydroxylases of these strains have not been elucidated. With the rapid development of synthetic biology and enzyme engineering, quite a lot of advances in bioproduction of calcifediol and calcitriol has been achieved in recent years. Therefore, here we review the successful strategies of promoting VD3 hydroxylation and provide some perspective on how to further improve the bioconversion of VD3 to calcifediol and calcitriol.

vitamin D metabolism. It has been reported that supplementation with calcitriol is helpful for the treatment of liver [5] and kidney disorders [6], and calcifediol is used for the treatment of calcium disorders, such as rickets [7]. In addition, calcifediol or calcitriol supplementation is able to prevent cartilage diseases, improve egg production and egg quality of laying hen, and enhance reproductive performance of breeding swine in the livestock and poultry production [8][9][10].
At present, calcifediol and calcitriol are mainly produced by chemical synthesis, which requires complex protection and deprotection steps for the regioselective introduction of hydroxyl groups, especially for hydroxylation on C-1, limiting the industrial production [11,12]. For example, starting with hydrindenylpropanol, Lythgoe et al. [13] synthesized des-AB-cholestane-8β,25-diol, which could be further used as a substrate for the synthesis of calcifediol. The synthesis process of des-ABcholestane-8β,25-diol requires multiple steps and further conversion to calcifediol. Calcitriol can be obtained by condensation of 1α-hydroxylated phosphine oxide with 25-hydroxy ketone [14]. However, the substrates also require multiple steps for synthesis. Thus, compared with chemical synthesis, the bioconversion of VD 3 is a successive process that comprises two (25-and 1α) or one (only 25-) hydroxylation steps. In addition, the yield of chemical synthesis of calcifediol is low (the classical approach by photochemical ring opening of steroidal Δ 5,7 -dienes yields less than 1%) [15]. Therefore, increasing attention is being paid to the bioconversion of VD 3 to calcifediol and calcitriol. Though hydroxylated sterols could be synthesized by transgenic Arabidopsis plants with non-heme monooxygenase [16], VD 3 has not yet been reported as a substrate. Microorganisms have been a powerful tool for VD 3 hydroxylation research, and in some work, the VD 3 conversion rate was remarkable. Therefore, a bioconversion process with VD 3 as the substrate and microorganisms as the chassis cells is promising for its simplicity and low pollution.
In the past few decades, a large number of strains expressing VD 3 hydroxylases have been developed for VD 3 bioconversion. A common strategy is to screen organisms with VD 3 hydroxylation capabilities and then examine their calcifediol and calcitriol production abilities by fermentation. Furthermore, the characterization of VD 3 hydroxylases from different sources and their heterologous expression in different host cells are also topics of particular interest in this field. These VD 3 hydroxylases mainly refer to cytochrome P450s. In addition, steroid C25 dehydrogenase (S25DH) was also reported to have the ability to hydroxylate VD 3 at the C-25 position. However, S25DH is sensitive to oxygen, causing a high cost for calcifediol production [17]. CYPs are the primary enzymes currently used for bioconversion of VD 3 to its active forms. A variety of CYPs from mammals and bacteria (mainly from actinomycetes) that are capable of converting VD 3 to calcifediol and calcitriol have been identified, and highly active mutants have been constructed. Enhanced calcifediol and calcitriol production could further benefit from the rapid development of creating new enzymes and new biochemical reactions via computational redesign of enzymes.
In general, enzymes and whole cells are two categories of biocatalysts commonly used for bioconversion [18,19]. Here, we review advances in the characterization and engineering of VD 3 hydroxylase that plays a key role in VD 3 bioconversion, and progress in whole-cell bioconversion, which can use inexpensive and abundant raw materials and avoiding the addition of expensive cofactors.

VD 3 hydroxylases and their protein engineering
For the synthesis of complex compounds, including calcifediol and calcitriol, enzymatic conversion has high regio-and stereoselectivity when compared with chemical synthesis. VD 3 hydroxylases, most of which are CYPs, play a central role in the bioconversion of VD 3 into calcifediol and calcitriol. Cytochrome P450s (CYPs) are a class of heme-containing enzymes with similar structures. CYPs catalyze a variety of reactions, mostly monooxygenation reactions, and have a wide range of substrate specificities [20]. The main CYPs related to VD 3 hydroxylation are shown in Table 1 and Fig. 1.
Mammalian and bacterial CYPs have been widely utilized for the bioproduction of calcifediol and calcitriol, the active forms of VD 3 . Mammalian CYPs have attracted great attention as they are associated with VD 3 metabolism related diseases. Bacterial CYPs are also of great importance because of their enormous functional diversity. A unique property of bacterial CYPs is that carbon monoxide (CO), SKF-525-A and metyrapone inhibit the hydroxylase activities of bacterial CYPs. Therefore, reduced CO difference spectroscopy can be used to identify bacterial CYPs that function as VD 3 hydroxylases [21,22]. To simplify the purification of CYPs, especially mammalian CYPs, the heterologous expression of CYPs in model organisms can increase protein yield and reduce the cost [23]. Escherichia coli and yeast are ideal host cells for heterologous expression of CYPs, as they grow quickly and genetic tools for their modifications have been well developed. Given that the hydroxylation reaction of CYPs requires redox proteins to deliver electrons, it is of great importance to improve VD 3 conversion rates using appropriate redox partners and coenzyme regeneration systems. Because of the existence of abundant redox proteins, actinomycetes are also becoming popular for heterologous expression of VD 3 hydroxylases. With an increased understanding of the relationship between the protein structures of CYPs and their hydroxylation abilities, protein engineering has been widely applied to the evolution of CYPs in the laboratory, which aims to improve the activity of CYPs for enhanced VD 3 bioconversion.

CYPs from mammals
Mammalian VD 3 hydroxylases such as CYP27A1, CYP27B1, CYP2D25, CYP2J3 and CYP2R1 have been found to be able to catalyze C-25 or C-1α hydroxylation [24][25][26][27][28]. The CYP27A1 was the first mammalian CYP to be thoroughly studied. Saccharomyces cerevisiae expressing a modified rat CYP27A1 localized in microsomes could convert 1α(OH)VD 3 to calcitriol when adrenodoxin (ADX) and NADPH-adrenodoxin reductase (ADR) were coexpresseded [24]. The ADRs from yeast and rats were unable to transfer electrons for rat CYP27A1. However, the purified bovine ADR and ADX facilitated the reaction. Therefore, the mitochondria-targeting signals of bovine ADR and ADX were removed to construct an electron transfer chain between the cytoplasm and the modified CYP27A1 that was incorporated into microsomes. The modified rat CYP27A1 could convert 1α(OH) VD 3 to calcitriol with the help of the newly constructed electron transfer chain. In one study, human CYP27A1, which also belongs to the CYP27 family was expressed in E. coli to examine its enzyme activity with VD 3 , and only a low activity of human CYP27A1 toward VD 3 (with V max of 0.27 mol/min/mol P450) was observed in this reaction system [29]. Although the 25-hydroxylation activity of human CYP27A1 was notably lower than that of several other mammalian CYPs such as human CYP2R1 (k cat = 0. 61  In addition, ergosterol in the yeast membrane might inhibit the VD 3 bioconversion as a potential substrate for CYP27A1. Therefore, compared with yeast, E. coli may be a better host for heterologous expression of mammalian CYPs due to its lack of the P450 genes and steroids [30]. Mitochondrial CYP27A1 in human liver is considered to be one of the CYPs that performs C-25 hydroxylation, while CYP27B1, another member of the CYP27 family, is responsible for hydroxylation at the C-1α position. However, Uchida et al. found that, by adding purified ADR and ADX to the reaction system, mouse CYP27B1 could perform both 25-hydroxylation and 1α-hydroxylation, similar to CYP27A1, suggesting that 1α-hydroxylation and 25-hydroxylation of vitamin D 3 are closely related [23]. This observation led to the hypothesis that vitamin can be arranged in the opposite direction of the substrate-binding pocket of CYP27B1. In addition, molecular chaperone GroEL/ES greatly enhanced the expression of mouse CYP27B1 in E. coli, making great progress for heterologous expression of CYPs in the future. In addition to mouse CYP27B1, human CYP27B1 also has 1α-hydroxylation activity [25]. However, the variant human CYP27B1 with CYP27B1 T409I or CYP27B1 Q65H amino acid substitution caused vitamin D 1α-hydroxylase deficiency (Vitamin D-dependent Rickets Type 1). The identification of amino acids that interact with the substrate could provide structural clues for the protein engineering of CYPs. The Ser 408 residue in mouse CYP27B1 and the Thr 409 residue in human CYP27B1 were found to be responsible for the C-1α hydroxylation by forming a hydrogen bond with the 25-hydroxyl group of calcifediol [31]. However, the problem with the low expression of human CYP27B1 still exists. Different strategies were thus employed to optimize the expression of human CYP27B1. High expression of human CYP27B1, shown by conversion to calcitriol with a yield of above 6000 pmol/mg protein in 30 min, was finally achieved by coexpressing GroEL/ES and human ADX/ADR rather than bovine ADX/ADR [32]. The addition of cardiolipin to the reconstituted system markedly lowered (threefold) the K m of human CYP27B1 for calcifediol, suggesting that the activity of CYPs could also be influenced by phospholipid composition. Several other mammalian CYPs were also found to be relevant to VD 3 hydroxylation. Both human CYP2R1 and human CYP27A1 can catalyze the 25-hydroxylation of VD 3 . However, the 25-hydroxylation activity of CYP2R1 on VD 3 is five times higher than that of CYP27A1. In addition, a mutation in the human CYP2R1 gene results in vitamin D-dependent rickets [33]. These data indicate that CYP2R1 also plays an important role in calcifediol biosynthesis in human body [28]. The GroEL/ES chaperonin can also be used to increase the expression level of human CYP2R1 in E. coli [34]. To understand which amino acid residues are responsible for the substrate specificity of CYP2R1, the crystal structure of CYP2R1 in complex with VD 3 was solved. VD 3 binds in an elongated conformation with the aliphatic side-chain pointing to the heme. The substrate access channel is covered by the ordered B′-helix. The conserved active site in a closed conformation of CYP2R1 is mainly composed of hydrophobic amino acid residues. Porcine CYP2D25 was expressed in yeast by Hosseinpour [26] and Araya [35], and the substitution of five residues in substrate binding site 3 (SRS-3) inhibited its 25-hydroxylase activity. Besides 25-hydroxylase towards VD 3 , CYP2D25 is also able to convert calcifediol into calcitriol, 25,26(OH) 2 VD 3 , and 25,27(OH) 2 VD 3 . Just like CYP2R1 as VD 3 25-hydroxylase, the recombinant rat CYP2J3 shows a 25-hydroxylation activity when it was coexpressed with GroEL/ES in E. coli [27]. The turnover numbers were 3.3 and 22 toward VD 3 and 1α(OH)VD 3 , respectively. Human CYP2J2 shows 73% amino acid homology with CYP2J3. After the deletion of a putative membrane anchor region in the N-terminus and the coexpression of GroEL/ES, efficient expression of CYP2J2 was achieved. However, human CYP2J2 does not seem to be a principal VD 3 25-hydroxylase [36]. The CYP2C11, a male-specific hepatic recombinant microsomal vitamin D 25-hydroxylase, has a substrate preference toward VD 3 over 1α(OH)VD 3 [37]. The 25-hydroxylation activity of the CYP2C11 towards VD 2 was higher than that towards VD 3 (25-hydroxylation rate of 0.49 ± 0.03 nmol·nmol −1 ·1.5 h −1 for VD 2 compared with 0.07 ± 0.02 nmol·nmol −1 ·1.5 h −1 for VD 3 ). The examples of CYP2J2 and CYP2C11 imply that some CYPs from mammals take VD 2 as a prioritized substrate and thus might not be the best choice for the hydroxylation of VD 3 .

CYPs from microorganisms
CYPs existing in a variety of organisms are involved in many biosynthetic pathways [20,25]. Since the initial discovery of the VD 3 hydroxylation in Amycolata sp. and Streptomyces sp. strains, CYPs from actinomycetes have been widely used for the VD 3 bioconversion process [21,22,38,39]. The 25-hydroxylation activity of VD 3 by CYP105A2 (with the highest yield of 10%) was observed by heterologous expression of CYP105A2 from Amycolata autotrophica (Pseudonocardia autotrophica) in Streptomyces lividans [40]. Cytochrome P450SU-1 (CYP105A1) from Streptomyces griseolus was found to convert VD 3 to 25(OH)VD 3 when expressed in E. coli [41]. Although the amino acid sequences of CYP105A1 and CYP105A2 share 55% identity, there is a clear difference in catalytic activity between CYP105A1 and CYP105A2. Besides the 25-hydroxylation of VD 3 , CYP105A1 also shows 1α-hydroxylation of calcifediol. In addition, CYP105A1 can also catalyze the hydroxylation of VD 2 with a relatively high V max /K m value (0.142 L −1 ·min −1 /mol P450, compared to 30 L −1 ·min −1 /mol P450 for VD 3 ). The CYP107CB2 from Bacillus lehensis G1 was classified into the CYP107 subfamily and showed 25-hydroxylation activity of both VD 3 and 1α(OH)VD 3 [42].
CYPs that catalyze the hydroxylation of VD 3 share a common catalytic mechanism with most CYPs. The catalytic cycle occurs on the heme prosthetic group which is linked by a conserved cysteine [43], but the catalytic specificity of CYPs mainly depends on their substrate binding pockets [44] (Fig. 2). Based on the catalytic mechanism, the catalytic specificity of CYP for VD 3 can be improved by directed evolution of CYPs.
Structural analysis of CYP105A1 in complex with calcitriol found that Arg 73 and Arg 84 along with Arg 193 were located in the distal pocket of CYP105A1. The Arg 73 and Arg 84 may have an inhibitory effect on activity, while Arg 193 is required for the activity [39,45,46]. Therefore, a series of CYP105A1 variants focusing on residues Arg 73 and Arg 84 were made to achieve a higher activity (with relative kcat/K m compared with wild type ranging from 1 min −1 ·μM −1 to 43 min −1 ·μM −1 with 1α(OH)VD 3 and calcifediol as substrate). Compared with the wild-type CYP105A1, the activity of CYP105A1 R84A variant greatly Fig. 2 Schematic representation of the substrate binding pocket. The substrate calcifediol is located in the substrate binding pocket in a suitable position, and will be further catalyzed by the Cys-linked heme  [45]. The variant of CYP105A1 containing amino acid substitutions R73V and R84A exhibited much higher hydroxylation activities at positions C-1α and C-25 compared to the wild type; the CYP105A1 R84A variant also shows a kcat/K m value for 25-hydroxylation that is 435fold higher than that of the wild type; the CYP105A1 R84A mutant also shows a higher 1α-hydroxylation capability than any single mutant. Inspired by the engineering of CYP105A1 R73VR84A , arginine residues around the substrate entrance and active site of CYP105D7 were mutated to generate a double mutant CYP105D7 R70AR190A [47]. The CYP105D7 R70AR190A achieved an almost ninefold increase in the conversion rate of testosterone, suggesting that the arginine residues around the substrate entrance and active site play important roles in determining the functions of the CYPs. Interestingly, the CYP105A1 is not the only VD 3 hydroxylase found in P. autotrophica. Vdh (CYP107BR1) from P. autotrophica NBRC 12743 also has the capability of converting VD 3 into calcifediol and calcitriol [48]. Different from other VD 3 hydroxylases, the Vdh has no preferred regio-specific hydroxylation of VD 3 .
A Vdh variant designated as Vdh-K1 with four amino acid substitutions was obtained by means of random mutagenesis, with its activity an order of magnitude higher than that of wild-type Vdh. Yasutake et al. then analyzed the structures of Vdh and Vdh-K1 [49]. The wild-type Vdh showed an open conformation, and the distal heme pocket was exposed to solvent with or without substrate. In contrast, Vdh-K1 shows a closed conformation, increasing the substrate binding affinity and catalytic activity. In addition, the combination of VD 3 and calcifediol with Vdh-K1 is in an antiparallel orientation. It is worth noting that the four substituted residues of Vdh-K1 are dispersed throughout the protein but result in a closed conformation of the heme pocket. The Vdh T107A variant, which was made by modifying the putative ferredoxin-binding site of Vdh, also exhibited a closed conformation similar to Vdh-K1 bound to VD 3 [50]. The ferredoxin-binding surfaces of Vdh T107A seem to show a stronger positive potential that is conducive to the binding of ferredoxin when compared with the structures of Vdh-K1. Although the enzyme activity of Vdh-K1 was much higher than that of wild-type Vdh, its expression level in Rhodococcus erythropolis cells was low, in contrast with that of Vdh T107A , which was approximately 70% of Vdh WT . In addition, Vdh T107A has an improved activity that is comparable to that of Vdh-K1, and achieves a similar expression that is comparable to Vdh simultaneously.
The above examples are the modification of the substrate binding pocket. Nevertheless, the engineering technology of CYPs has been greatly developed. Drawing on these technologies, it may be possible to use allosteric regulation, decoy-based small molecules etc., for reference in the hydroxylation of VD 3 in the future. For example, the conformational transition of CYP3A4 is induced by the allosteric ligand progesterone (PGS), and CYP3A4-PGS conjugates thus acquire a twofold higher capacity to oxidize testosterone than CYP3A4 [51]. Inert decoy molecules that have a native substrate-like structure, can constrain substrates, such as benzene and benzene derivatives, into the suitable positions to facilitate the catalysis by targeting the enzyme substrate pocket to induce an active enzyme intermediate [52,53]. If the above technologies are combined into the catalytic reaction of VD 3 , it may be possible to expand the types of CYPs that can be used to catalyze VD 3 , inspiring future research that may have unforeseen benefits.
Compared with CYPs, S25DH has advantages as well as disadvantages. The catalytic activity of S25DH is higher than that of the CYPs reported thus far. In addition, the regioselectivity toward the C-25 position of VD 3 is higher in S25DH than that in CYPs, making the recovery of calcifediol much easier. However, S25DH is sensitive to oxygen, so the bioconversion process must be carried out under anaerobic conditions [54]. Jacoby et al. developed a platform for the overproduction of four steroid C-25 hydroxylases (S25DH 1 and three isoenzymes S25DH 2 , S25DH 3 , and S25DH 4 ) by coexpressing an essential chaperone in the betaproteobacterium Thauera aromatica K172 [55]. The crude extract from T. aromatica that overexpresseed S25DH 1 achieved a reaction rate that was 6.5-fold higher than that of the wild-type bacterial extract. However, the oxygen intolerance and high cost prevent S25DH from being applied for the industrial production of calcifediol.

Whole cell conversion of VD 3 into calcifediol and calcitriol
A whole-cell conversion system expressing the CYPs of interest has advantages over enzymatic conversion. It avoids adding the expensive cofactor NAD(P)H, and the hydrogen dioxide produced by the uncoupling reactions does not accumulat because of the presence of catalase in the cell. Therefore, increasing attention has been given to the whole cell conversion of VD 3 into calcifediol and calcitriol. initially, attempts were made to screen natural microbial strains that are able to convert VD 3 to calcifediol and calcitriol in nature. With the development of genome editing and synthetic biology, the production of calcifediol and calcitriol has been further improved by either overexpressing heterogenous VD 3 hydroxylase genes or altering the electron transfer chains.

Conversion of VD 3 to calcifediol and calcitriol by microbes isolated in nature
Sasaki et al. reported the bioconversion of VD 3 derivatives for the first time [21]. After screening approximately 300 Streptomyces strains, S. sclerotialus FERM BP-1370 and S. roseosporus FERM BP-1574 were found to be able to introduce hydroxy group at C-1α position of calcifediol and C-25 position of 1α(OH) VD 3 , respectively. Although VD 3 derivatives were used as substrates rather than VD 3 itself, this work demonstrated that microorganisms could be used for the production of hydroxylated derivatives of VD 3 . By further screening within bacterial and fungal strains, P. autotrophica FERM BP-1573 was chosen from twelve VD 3 -hydroxylating strains and achieved a calcifediol titer of 8.3 mg/L and a calcitriol titer of 0.17 mg/L in a 200 L fermenter [22].
With the ever-increasing number of newly discovered VD 3 -hydroxylating strains, the production of calcifediol and calcitriol has been further improved. The actinomycete Kutzneria albida produced 70.4 mg/L calcifediol and 2.0 mg/L calcitriol in an optimized conversion process [56]. When using Pseudonocardia sp. KCTC 1029BP as a whole-cell catalyst, calcifediol production further rose to 356 mg/L, and calcitriol production rose to 61.87 mg/L, with bioconversion yields of 59.4% and 30.94%, respectively [57,58]. More recently, Tang et al. obtained up to 830 mg/L calcifediol using Bacillus cereus zju 4-2, representing a significant increase in calcifediol production [59] (Table 2). In addition to calcifediol and calcitriol, some side products of VD 3 hydroxylation such as 2α,25(OH) 2 VD 3 could also be synthesized, increasing the difficulty of product recovery [60]. It is of great significance to reduce the formation of side products by isolating a microbial strain with high regioselectivity.
Though many attempts have been made to enhance the calcifediol and calcitriol production, the acquirable concentrations of calcifediol and calcitriol in the fermentation broth are not yet satisfactory. Therefore, a series of fermentation parameters have been optimized to further improve the production of calcifediol and calcitriol. The influence of bioconversion buffers, solubilizers, and metal salts on calcitriol production has also been examined, and optimized bioconversion medium was finally been obtained. The optimization of three key variables (aeration rate, resting cell concentration and temperature) in Pseudonocardia fermentation improves the VD 3 bioconversion [57,58].

Conversion of VD 3 to calcifediol and calcitriol by engineered strains
With the development of genomics and synthetic biology, engineering of microbial cells provides a new strategy to improve the production of calcifediol and calcitriol. Besides previously mentioned heterologous expression of CYPs for purpose of enzymatic activity characterization and structural analysis, design of engineered strains for whole-cell bioconversion could be applied for industrialscale production of calcifediol and calcitriol. Actinomycete chassis cells such as S. lividans and R. erythropolis have stable redox environments, which may explain why most bacterial VD 3 -hydroxylating CYPs have been discovered in actinomycetes; thus, these cells could be promising chassis cells for calcifediol and calcitriol production. When expressing CYP105A1 R73VR84A in S. lividans, approximately 7.7 mg/L calcifediol was observed in the culture [39]. R. erythropolis is another industrially important actinomycete that exhibits resistance to organic solvents, so it was used as a host cell to convert VD 3 to calcifediol and calcitriol [61] (Fig. 3). The permeability of VD 3 into the cytoplasm is a rate-limiting step for VD 3 bioconversion [62]. In addition, the electron carrier ferredoxin and the NADH-regeneration system also play essential roles in the conversion of VD 3 to its bioactive form. When Vdh was coexpressed with the redox partners AciBC and glucose dehydrogenase GlcDH-IV in nisin-treated R. erythropolis cells, 1176.5 μg calcifediol was obtained, representing an improved calcifediol production when compared with the wild-type strain [61]. Deletion of CYP-sb3 gene dramatically impaired the VD 3 bioconversion of by Sebekia benihana, suggesting that CYP-sb3a functions as a VD 3 hydroxylase. When CYP-sb3a was expressed in Streptomyces coelicolor that originally had no VD 3 hydroxylation ability, calcifediol and calcitriol were detected in the recombinant S. coelicolor cells [63]. B. megaterium is another popular host that has been used for VD 3 bioconversion (Fig. 4). When B. megaterium expressing human CYP27A1 was used as whole-cell catalyst, 80.81 mg/L VD 3 was converted to its bioactive form. In addition, B. megaterium overexpressing CYP109A2 T103A produced 282.7 mg/L calcifediol in 48 h under optimized conditions. Though there are still technical challenges for calcifediol and calcitriol production by engineered strain, designing better artificial strains is a promising strategy to further enhance the calcifediol and calcitriol production by designing better artificial strains.
With the help of the NADPH-regenerating system, approximately 40% VD 3 was converted to calcifediol in 3 hours. The CYP109A2 [64] and CYP109E1 [65] from Bacillus megaterium also exhibited high VD 3 bioconversion abilities. After 24 h of whole-cell catalysis, the conversion rates of CYP109A2 and CYP109E1 reached 76% and 95%, respectively. However, NMR analysis showed that CYP109E1 catalyzed the hydroxylation of VD 3 at both C-24 and C-25. Approximately twofold more calcifediol was produced in the CYP109E1 I85A variant obtained by site-directed mutation than in the wild-type   [66]. The highest yield by engineered strains was achieved by the protein engineering of Vdh (CYP107BR1). For the Vdh and its variants, when R. erythropolis cells expressing Vdh T107A were treated with nisin for VD 3 conversion, 573 mg/L calcifediol was produced in 2 hours, representing a significantly improved bioconversion rate. Compared with wild-type Vdh, a 1.7-fold of calcifediol production was observed in shorter time (2 h compared with 16 h).
Whether by isolated strains or engineered strains, the efficiency of the VD 3 bioconversion process is limited in different levels, mainly genes expression, enzyme activity, efficiency of CYP-redox partner interactions, the NAD(P) H regeneration system and the substrate transportation. Since various efforts have been made in optimizations, a combination of different methods could result in a higher yield. For example, with E. coli as the host, not only by expression of high enzymatic activity CYPs which are assisted with adequate Fdr and Fdx, the introduction of GroEL/ES system along with the mutation of genes include acrAB and tolC [67], whose expression products compose a efflux pump system which play positive roles simultaneously. The deletion of acrAB and tolC in E. coli BL21star(DE3) works similar to the nisin treatment for R. erythropolis, which can increase in the intracellular substrate concentration. In addition, in the fermentation process, the addition of Fe 2+ and 5-aminolevulinic acid (5-ALA) [68], which is precursor of heme can be adopted in most cases to ensure the functional expression of CYPs.
In addition to improving the transportation of substrates, there are other potential methods of improving VD 3 bioconversion. CYPs displayed on the cell surface of E. coli [69,70] will provide alternative solution for the transportation of VD 3 . In the meanwhile, this method could be also used for the screening of highly efficient VD 3 hydroxylases. The fusion of CYPs and redox partners could construct a self-sufficient P450 [71], which is promising for future methods of VD 3 bioconversion.

Conclusions
Calcifediol and calcitriol have bright market prospects in the poultry farming and pharmaceutical industries. Compared with the traditional chemical synthesis, bioconversion of VD 3 to calcifediol and calcitriol represents a promising and eco-friendly technology, and will contribute to the sustainable development goals of United Nations. Given that VD 3 hydroxylase plays a key role in the bioconversion of VD 3 to calcifediol and calcitriol, a series of CYPs have been modified by protein engineering to increase their catalysis of VD 3 hydroxylation. In addition, the VD 3 permeability and the electron transfer from NADH to VD 3 hydroxylase are also Fig. 4 B. megaterium is used as a host cell for VD 3 bioconversion. The electrons required for CYP activity are supplied by the endogenous cytochrome P450 reductase (CPR): NADPH-dependent diflavin reductase [76]. Quillaja saponin is used as a membrane-solubilizing agent of great importance in determining the rate of VD 3 bioconversion. Therefore, the lipid II-targeting antibiotic nisin, redox partners and NAD(P)H regeneration systems are employed to further improve the VD 3 bioconversion by enhancing the VD 3 uptake and optimizing the electron transfer pathway. High-throughput screening technology has shown its advantages in the field of strain screening [72], and high-throughput screening based on the calcifediol/calcitriol sensor can be used to screen microbial strains that are capable of converting VD 3 to calcifediol or calcitriol more efficiently.
With increasing pressure to conserve energy and protect the environment, increasing attention is being paid to the bioproduction of value-added chemicals in a more eco-friendly manner [73,74]. The "next generation industrial biotechnology" based on moderately halophilic bacteria, which can avoid sterilization and the use of fresh water, exhibits advantages over the current industrial biotechnology. Given that a large number of VD 3 -hydroxylating strains belong to actinomycetes, efficient screening of halophilic actinomycetes could be a promising strategy to obtain more ideal calcifediol/calcitriol-producing strains. In addition, engineering of the well-studied halophilic bacteria such as Halomonas bluephagenesis could be another promising strategy to further improve the calcifediol/calcitriol production in the future.