Longitudinal relationship of liver injury with inflammation biomarkers in COVID-19 hospitalized patients using a joint modeling approach

The mechanisms underlying liver disease in patients with COVID-19 are not entirely known. The aim is to investigate, by means of novel statistical techniques, the changes over time in the relationship between inflammation markers and liver damage markers in relation to survival in COVID-19. The study included 221 consecutive patients admitted to the hospital during the first COVID-19 wave in Spain. Generalized additive mixed models were used to investigate the influence of time and inflammation markers on liver damage markers in relation to survival. Joint modeling regression was used to evaluate the temporal correlations between inflammation markers (serum C-reactive protein [CRP], interleukin-6, plasma D-dimer, and blood lymphocyte count) and liver damage markers, after adjusting for age, sex, and therapy. The patients who died showed a significant elevation in serum aspartate transaminase (AST) and alkaline phosphatase levels over time. Conversely, a decrease in serum AST levels was observed in the survivors, who showed a negative correlation between inflammation markers and liver damage markers (CRP with serum AST, alanine transaminase [ALT], and gamma-glutamyl transferase [GGT]; and D-dimer with AST and ALT) after a week of hospitalization. Conversely, most correlations were positive in the patients who died, except lymphocyte count, which was negatively correlated with AST, GGT, and alkaline phosphatase. These correlations were attenuated with age. The patients who died during COVID-19 infection displayed a significant elevation of liver damage markers, which is correlated with inflammation markers over time. These results are consistent with the role of systemic inflammation in liver damage during COVID-19.


Patients and methods
Data were collected from the medical reports of patients diagnosed with COVID-19 and admitted to the Complejo Hospitalario Universitario of Santiago de Compostela in Spain, from March 12, 2020 (date of first COVID-19 diagnosis) to April 11, 2020. A confirmed case of COVID-19 was defined as a positive result in the reverse transcription polymerase chain reaction (RT-PCR) test on samples obtained from nasal or throat swabs performed in accordance with WHO protocol 17 . The study time for each included patient was from admission until the patient died or was admitted to the intensive care unit.
Using Python 3.8, a specific program with a graphical user interface was designed through which each patient's clinical data were imported and automatically recorded from the electronic medical record. In all cases, the data were encrypted prior to filing. We collected all demographic data, admission history, and medical comments on the patient's progression, clinical constants (heart rate, blood pressure, and body temperature), pulse oximetry, radiological studies, and pharmacological prescriptions, including the date and time of their prescription. Among the laboratory determinations, we recorded all the parameters of the complete blood count, and blood biochemistry (glucose, urea, creatinine, uric acid, sodium, potassium, calcium, magnesium, albumin, cholesterol, triglycerides, AST, ALT, gamma glutamyl aminotransferase [GGT], Alkaline phosphatase [ALP], bilirubin, and coagulation), which included the international normalized ratio, prothrombin, blood gases (pO 2 , pCO 2 , O 2 saturation, bicarbonate, lactate), procalcitonin, lactate dehydrogenase, vitamin D, and hemoglobin A1c. We considered CRP, D-dimer, IL-6, and lymphocyte counts as a specific group of IMs and recorded the reference limits of normality for each parameter and laboratory measurement. An abnormal increase in transaminase levels has been considered any value that exceeds the normal range. The transaminase elevation profile has been grouped according to R-factor and similar to how drug-induced liver damage (DILI) is classified depending on whether the predominant elevation is due to AST and ALT (hepatocellular pattern), GGT and ALP(cholestatic pattern) or elevation similar to one another (mixed).The R-value is defined as serum alanine aminotransferase/ upper limit of normal (ULN) divided by serum alkaline phosphatase/ULN. By common convention, R ≥ 5 is labeled as hepatocellular, R < 2 is labeled as cholestatic, and 2 < R < 5 is labeled as "mixed" 18 .
For those cases in which specific treatment had been started before hospital admission, the start date of each drug was manually recorded. The program classified each drug into its corresponding pharmacological group, adding this as an additional variable in each registry. The records were processed and archived in a relational sqlite3 database, with independent tables for physical constants, pharmacological prescriptions, symptoms, exploratory signs, radiological signs, laboratory parameters, complications, and status (emergency room, hospital ward, critical unit, discharge or death), indexed by case number and start date and linked to a general table with demographic characteristics, antecedents, and form of infection.
Ethical issues. The study was conducted in accordance with the guidelines of the Declaration of Helsinki and the principles of good practice and was approved by the Institutional Review Board of the Galician health service on April 3, 2020April 3, (# 2020 Given that not all blood samples were analyzed with the same apparatus, the raw values of certain variables were not comparable. Therefore, prior to any statistical analysis, we divided the raw values for transaminases, IL-6, and CRP by their maximum normal values a form of weighting. For these variables, we thereby obtain the proportion with respect to the acceptable maximum (< 1 implies normal levels, and > 1 implies high levels). At this point, given that transaminases and IMs do not follow a Gaussian distribution, we estimated the exploratory curves for these variables over time for the deceased and living patients by the least absolute shrinkage and selection operator regression analysis 19 , given that this technique does not assume any parametric distribution for data.
We then normalized the transaminases, D-dimer, and IL-6 values using a Box-Cox transformation to facilitate statistical modeling calculations. We considered that transforming the CRP and lymphocytes values would not produce any improvement, so we decided to use the raw data for these covariates.
To assess the factors that could influence the transaminase values, we fitted a generalized additive mixed model (GAMM) 20 for each Box-Cox-transformed transaminase, adjusting for sex, exitus (binary for alive or deceased), antiviral treatment (binary for lopinavir/ritonavir and/or hydroxychloroquine administration), age, IM-exitus association, and time-exitus association. We also included random intercepts by patient to consider the variability within each individual (the values recorded for the same individual are expected to be more closely related to each other than to the values of another individual). The smooth effects of continuous covariates were considered by means of thin-plate spline smoothers. The GAMM regression model considered here was as follows (1): where "Transaminase" is AST, ALT, GGT, or ALP, and "IM" is CRP, D-dimer, IL-6, or lymphocytes. In the above equation, s(IM)*Exitus represents the association of "Exitus" according to the IM curve, i.e., a possible different non-linear effect of IM for deceased and living patients, and the same for s(Time)*Exitus. "ID" indicates the patient, and "s(ID, bs = "re")" indicates the random intercepts by the patient. In the graphs, the smoothed effects are shown centered and with a 95% confidence interval.
After these first models, a second GAMM for each transaminase was fitted, including only the statistically significant covariates from the previous models. Statistical significance was established for p-values less than 0.05 and, in the case of smoothed continuous covariates, for those ranges of the variable in which the estimated 95% confidence intervals do not include zero. For the categorical covariates, as well as for the continuous variables whose smoothed effect was linear, we reported the corresponding estimated β coefficient. The effective degrees of freedom (edf) were also reported for all covariates. The edf number represents the parameters needed to estimate the effect of a covariate, i.e., edf indicates the complexity of the functional form of the effect. When the effect is linear, the edf is approximately 1.
The temporal correlation between transaminase levels and IM was estimated using the bivariate copula generalized additive models for location, scale, and shape 21 . Such a correlation is presented by means of Kendall's τ. The model for the correlation of each transaminase-IM pair is determined by the following equation (2): in which, once again, "Transaminase" is AST, ALT, GGT, or ALP, and IM is CRP, D-dimer, IL-6, or lymphocytes. In the above equation, "s(Time)*Exitus" represents the association of "Exitus" according to the IM curve, i.e., a possible different non-linear effect of time for the deceased and living patients. All the Box-Cox-transformed covariates followed a Gaussian distribution, and the CRP and lymphocytes each followed a gamma distribution. In this type of model, a copula function is needed to model the correlation (see 21 for more details). For this study, we selected the best copula according to Akaike's Information Criterion and the Bayesian Information Criterion.
All the statistical analyses were performed using the free statistical software R 22 , version 4.1.0, using the MASS package for the Box-Cox transformation, the mgcv package for the GAMM regression, the GJRM package for the correlation analysis, and ggplot2 for all the plots presented in this paper.

Results
General description of mortality and associated factors. During the study period, 212 patients were admitted to the hospital with COVID-19. Table 1 shows the sample's demographic and clinical characteristics. The median age was 68 years, and 56% were men. In our series, 37 (17.4%) patients died. The patients who died were significantly older, and there was a significantly higher percentage of patients with a history of ischemic heart disease, diabetes, and a lower oxygen saturation rate during the hospitalization. There were no differences between the deceased and those who survived with respect to the treatment undergone, past HBV or HCV infection, and statin treatment. There were no differences in body mass index, although those who died tended to have a higher index. Fifty three (53%), 56%, 42% and 19% of the patients presented elevated AST, ALT, GGT and ALP respectively. Of those patients: 12%,50% and 38%,presented elevated transaminase levels with a hepatocellular, cholestatic, and mixed pattern, respectively. The surviving patients who presented a hepatocyte pattern had a significant increase in ALT. The median AST and ALP values were significantly higher in the patients who died. www.nature.com/scientificreports/ There was a significant difference in IL-6, D-dimer, CRP, and lymphocyte values between the patients who died and those who survived ( Table 1). The temporal progression of enzyme and IM levels in the survivors and the deceased is shown in Fig. 1, which represents the mean estimate with a 95% confidence interval of these values over time and as a function of the outcome. The left part of the figure shows the transaminase behaviors in these patient groups. The patients who died showed a greater increase in transaminase levels than those who survived, although this increase was only notable for AST and ALP. Similarly, IMs showed a greater increase in the patients who died (right side of Fig. 1), although only CRP and D-dimer were significantly elevated compared with the levels in the patients who survived. Lymphocyte counts showed an inverse behavior, rising significantly in the patients who survived and decreasing significantly in the patients who died. This separation occurred throughout the disease's progression.
The progression of transaminase levels at fixed intervals every 7 days during the first 3 weeks of hospitalization (Table 2) shows a characteristic elevation, from the first week, with respect to AST, ALT, GGT, and ALP, especially in those patients who died. The same progression for IMs also demonstrated a characteristic elevation for CRP and D-dimer, but from day zero, with a second elevation at the end of the observation period. This movement was less evident in the case of IL-6. Conversely, the lymphocytes count fell from the first week for the patients who died, with the differences being more intense between the second and third weeks compared with those of the survivors throughout the entire hospitalization period.

Multivariate analysis of transaminase values.
In a preliminary analysis, each transaminase figure was evaluated individually with each of the IMs and adjusted for age, sex, and death, as well as for having undergone Table 1. Baseline characteristics of the study population, stratified by outcome. Age, AST, ALT, GGT, ALP, Bilirubin, CRP, D-Dimer, IL-6, Lymphocytes, BMI, and oxygen saturation: Median (interquartile range). (*) The effective "n" was a total of 180: 23 exitus and 157 living; ** The effective "n" was a total of 181: 23 exitus and 158 living. (***) The effective "n" was a total of 121: 24 exitus and 97 living. (****) The effective "n" was a total of 93 patients. www.nature.com/scientificreports/ antiviral treatment. This analysis was repeated by performing a simultaneous adjustment between each of the transaminases and the 4 inflammation parameters considered, adjusting for the statistically significant variables in the first analysis. The results of this second analysis are shown in Table 3.
With regard to the time variable, ALT and ALP increased significantly, reaching a maximum from the second week in the patients who died. In addition, the AST levels decreased significantly and the GGT increased in those who survived. The latter drops significantly at the end of the observation period in both these patients who lived and in those who died.
ALT is associated with age, increasing in middle age and decreasing significantly after the eighth decade. Despite its lack of statistical significance, the trend was similar in the analysis of age and the remaining transaminases. The following variables were statistically significant in this analysis: AST increased with IL-6 (p < 0.001), ALT decreased with CRP, and GGT decreased with lymphocytes. These 3 variables (AST, ALT and GGT) increased in turn with D-dimer (Fig. 2).

Longitudinal correlations between transaminases and inflammation parameters.
To assess the potential effect of time, drugs, sex, and age on the correlation between transaminases and IMs, we used the regression model expressed in equation 2 (see "Statistical analysis" section). Table 4 shows the regression models, indicating only the statistically significant variables. Sex affected the correlation between D-dimer and AST and between D-dimer and ALT, increasing in women (β, 1.231 [p = 0.007] and β, 0.272 [p < 0.001], respectively). Sex also had an influence (but inverse), making the correlation lower in women between AST and lymphocytes, between ALT and CRP, between GGT and CRP, between GGT and lymphocytes, and between ALP and lymphocytes

Effect of time on the correlation between transaminases and inflammation parameters.
The temporal correlation between transaminases and IMs (Fig. 3) showed a different dynamic between the survivors and the deceased. In the former (Fig. 3, in blue), CRP and AST showed a positive correlation that weakened from the second week of hospitalization. In this patient group and from the second week on, CRP and ALT showed a negative relationship, as did CRP and GGT. In contrast, the correlation between CRP and ALP in the surviving patients was positive. D-dimer showed an initially positive correlation to AST and ALT that became negative from the second week in the patients who survived. In contrast, the relationship between D-dimer, GGT, and ALP in the surviving patients was positive. IL-6 and AST in the surviving patients showed an initially positive correlation that weakened after 15 days. The correlation between IL-6 and GGT was also positive in these patients.
In the patients who died (Fig. 3, in pink), the correlation was positive between CRP and AST and between GGT and ALP, again after the second week. Although the correlation between D-dimer and ALT was not www.nature.com/scientificreports/ statistically significant, it was positive between AST, GGT, and ALP in the patients who died. In this group of patients, the correlation was significant and positive for IL-6 and AST, GGT, ALP, and ALT, although at the end of the observation period the correlation became negative for ALT. Lastly, the lymphocyte counts in the survivors showed a negative correlation for AST and a positive one for GGT and ALP from the second week. In the patients who died, the correlation was statistically significant and negative for AST, GGT, and ALP and especially for ALT, although it became positive at the end of the observation period. Figure 4 shows the correlation between transaminases and IMs as a function of age; this correlation was lost in those older than 80 years. The AST and CRP correlation weakened in these older patients and became negative with ALT and GGT. With D-dimer and AST, ALT had a negative correlation in this age subgroup and again weakened in its relationship with GGT and FA, although the correlation with the latter continued to be positive. The correlation between IL-6, ALT, and GGT again became negative. In contrast to these markers in lymphocytes, the correlation expressed a positive trend as age advanced.

Discussion
There is currently an abundance of articles indicating the high frequency of liver disorders in COVID-19. Nevertheless, it is still unclear why the liver is so frequently affected and what mechanism is ultimately involved [2][3][4][5][6][7][8][23][24][25] . Understanding why the liver experiences these changes is key to understanding the pathophysiology of COVID-19.
As with other studies, our study found that a high proportion (60%) of patients with COVID-19 had elevated transaminase levels. Unlike other studies, the percentage of patients with a hepatocellular and those with a cholestatic pattern was similar 3,25 . It is well documented that patients with previous liver damage appear more likely to present alterations with infection 5,11,24,26-29 . Our cohort was characterized by a low incidence of underlying liver disease, as well as of HCV and HBV infection. In addition, they were "first wave" patients treated with antivirals Table 3. Results for the multivariate regression models for the weighted and Box-Cox-transformed values of transaminases. SE means standard error. The edf number represents the parameters needed to estimate the effect of a covariate, indicating the complexity of the functional form of the smooth effect. www.nature.com/scientificreports/ (lopinavir-ritonavir) and antimalarials (hydroxychloroquine) but not with immunosuppressive drugs such as corticosteroids. We believe that both characteristics facilitate the explanation of the potential mechanisms of damage, given that they are not masked by certain therapeutic interventions or previous pathology. Another common finding in other studies is that older patients with ischemic heart disease and diabetes died more frequently from COVID-19. There is also a trend toward higher mortality in patients with higher body mass index. The missing data in this study have probably prevented us from observing a significant effect for obesity on prognosis, as reported in other publications 2,4,5,30,31 . We also observed that mortality was higher in those patients with greater higher increase in IMs, as has been reported in other studies 32,33 .
We found that the patients who died showed a significant increase in AST and ALP levels, whereas the patients who survived showed a significant decrease in AST levels. The increase in ALP levels and their temporal dynamics is especially striking. The increases in this enzyme in COVID-19 have been described in the medical literature as moderate 3,11 . The peculiar dynamics of ALP in this study have been revealed by the statistical methodology employed, given that we performed a continuous temporal analysis, the effects of which can be obviated in an analysis at a specific moment.
The mechanisms of liver damage have been associated with a direct action of the virus, the inflammatory response, hepatotoxicity, and anoxia [9][10][11] . Although the present study was not designed to determine the direct effects of the virus on its target organ, the statistical methodology applied allows us to infer the indirect mechanisms of damage. A correlation analysis found an association between IMs and transaminases. From the second week of hospitalization, the patients who died showed positive and significant 2-to-2 correlations between CRP and AST and between GGT and ALP. The D-dimer correlations are equally significant for this with AST, GGT, and ALP, which also occurs with IL-6. In contrast, the surviving patients showed negative correlations between CRP, ALT, and GGT, an effect that was equally significant in those patients who survived, from the second week of hospitalization, between D-dimer, AST, and ALT. The relationship between Il-6 and transaminases in these surviving patients is irrelevant.
In the patients who died, an inverse correlation between lymphocytes and AST, GGT, ALP, and even ALT is especially striking, although this effect was lost with this enzyme at the end of the observation period. Also of particular interest is the markedly negative correlation between ALP and lymphocytes for the patients who died and the positive correlation for those who survived. In adjusting the correlations with covariates, antiviral drugs exert a modifying effect, increasing the correlation between GGT and D-dimer and decreasing the correlation between ALT and CRP; however, these effects did not affect the form of the relationship. On the other hand, these drugs were administered in a similar proportion to the patients who lived and to those who died, which makes it impossible to attribute a significant effect to their use, both in terms of non-mortality and hepatotoxicity. We therefore believe that the treatment administered did not play a role in the liver damage. In this work, a sequential alteration of inflammatory markers followed by enzyme elevation has been detected, which could suggest that the latter is a consequence of inflammation, although a direct cytopathic effect cannot be ruled out.
This would justify the performance of specially directed studies to validate this hypothesis. The observed cholestatic pattern of elevated ALP levels in the patients who died in our study is more biologically plausible given that the different distribution of specific receptors for the virus in hepatocytes and cholangiocytes is well documented and is predominant in the latter 10,11,24 . www.nature.com/scientificreports/ The effect of age is particularly interesting. In a regression analysis, age behaved as a variable associated with ALT; ALT levels increase early in life and decrease from the eighth decade. When we established the correlations between age and the correlation itself with the markers of inflammation and transaminases, we observed that they became negative after the eighth decade. There is no clear biological explanation for this phenomenon, although the immunological senescence observed in vaccine trials could point to the effect 34,35 . Patients with higher increases in inflammation parameters showed higher mortality 32,33 , which contrasts with the fact that older patients despite showing less inflammation die more frequently; however, it is also true that older patients are affected by higher rates of diabetes, obesity, and heart disease, which influence mortality.
Various studies with longitudinal measurements of transaminases have been published, and a number of them have included inflammation parameters [12][13][14][15][16] . In a cross-sectional study, Phipps et al. found a relationship between ferritin and IL-6 with severe liver damage in a large cohort of patients with COVID-19 3 . Another longitudinal study established a relationship between ferritin and AST elevation but not with other IMs 12 . Another study observed a relationship between leukocyte counts, decreased lymphocyte counts, and increased transaminase levels 35 . In the remaining publications, specific associations were found between inflammation and liver disorders on one hand and mortality on the other but not specifically between IMs and liver disorders [14][15][16]27 .
In the present study, we used innovative flexible regression joint modeling techniques that, applied in this area, allow us to assess the dynamics of the temporal correlation between IMs and elevated transaminase levels. Furthermore, the advantage of this statistical methodology is that it allows us to adjust the temporal correlation for covariates such as age, sex, and treatment administration. For the successful application of this methodology, an exhaustive data capture was crucial, in which the measurements have been taken continuously over time, employing a specifically designed database. The use of these techniques has allowed us to establish a dynamic www.nature.com/scientificreports/ temporal representation of transaminases different from that of other studies. As far as we know, ours is the first study to use this type of tool to analyze liver disorders in SARS-CoV-2.
The retrospective nature of this study limits our conclusions. Despite exhaustive data capture, certain data was lacking to establish more robust conclusions. We also lacked histological studies. A histology compatible with immune damage would reinforce our conclusions. In addition, we only analyzed a certain IM battery, which limits the conclusions to the analyzed markers, although the statistical techniques employed can be extended to a larger panel of markers.
In conclusion, the present study once again confirms the frequent alteration of transaminase levels in patients with COVID-19 infection. A joint-modeling statistical methodology established that the patients who died from COVID-19 had characteristically elevated AST and ALP levels and showed a significant and sequential correlation over time with IMs. The patients who survived were better protected from inflammation-mediated liver damage. This relationship between transaminases and IMs weakened with age.

Data availability
The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.