**Personalized Medicine in Epidemics**

Editor

**Rutger A. Middelburg**

MDPI Basel Beijing Wuhan Barcelona Belgrade Manchester Tokyo Cluj Tianjin

*Editor* Rutger A. Middelburg Public Health and Primary Care Leiden University Medical Center Leiden Netherlands

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Journal of Personalized Medicine* (ISSN 2075-4426) (available at: www.mdpi.com/journal/jpm/ special issues/Personalized Epidemics).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-4208-9 (Hbk) ISBN 978-3-0365-4207-2 (PDF)**

© 2022 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


Reprinted from: *J. Pers. Med.* **2021**, *11*, 678, doi:10.3390/jpm11070678 . . . . . . . . . . . . . . . . **99**

#### **Fu-Zong Wu, Yun-Ju Wu, Chi-Shen Chen and Shu-Ching Yang**


## **About the Editor**

#### **Rutger A. Middelburg**

Rutger Middelburg achieved a PhD in Epidemiology at the department of Clinical Epidemiology of the Leiden University Center. He has since spent over a decade researching clinical transfusion medicine and advocating the personalization of medicine. Since 2020, he has worked at the department of Public Health and Primary Care of the Leiden University Medical Center and as an editor at the Dutch drug bulletin (a member of the International Society of Drug Bulletins).

## **Preface to "Personalized Medicine in Epidemics"**

The only good medicine is personalized medicine. There is no "one size fits all"in the vastly divergent landscape formed by the endless numbers of unique patients that need to be treated. No two patients are the same; thus, neither is the optimal treatment for two different patients ever completely identical. Therefore, a good doctor needs to try to always maximize the personalization of the medicine they practice. At the same time, the medical world is continuously challengeed by new epidemics. The SARS-CoV-2 pandemic is an obvious example. Possibly less obvious, but no less important, the epidemic of obesity and its associated morbidities also continues to impact all aspects of medicine. Even lesser epidemics such as seasonal flu or lung cancer all impact our ability to personalize the treatments offered to patients.

Therefore, this reprint explores the impact of epidemics of all kinds on the personalization of medicine. With eleven high-quality scientific studies on this subject, this book represents an important compilation for medical doctors and researchers alike. The broad overview offered in this reprint is by no way meant to be exhaustive, but rather to provide examples from diverse fields of medicine, thereby emphasizing the importance of the subject for all.

> **Rutger A. Middelburg** *Editor*

### *Editorial* **Personalized Medicine in Epidemics**

**Rutger A. Middelburg**

Department of Public Health and Primary Care, Leiden University Medical Center, 2300 RC Leiden, The Netherlands; r.a.middelburg@lumc.nl

Before you lies the Special Issue "Personalized Medicine in Epidemics". As we stated in our call for papers, we were looking for papers which make a novel contribution towards optimizing the personalization of medicine during epidemics. In this context, we aimed to include papers covering all kinds of epidemics, whether big or small, and whether infectious or non-infectious in nature. Personalized medicine is important in all fields of medicine, and all epidemics influence our ability to practice personalized medicine—not only for those suffering from the epidemic disease itself, but also for those with other diseases, which tend to get less attention as a major epidemic unfolds. This broadly inclusive view of the role of personalized medicine during epidemics has led to the inclusion of eleven high quality papers, on a wide range of topics within this Special Issue.

Dopazo et al. [1] review opportunities and challenges of personalized medicine in the context of the COVID-19 pandemic. Boboc et al. [2] were more specific in their approach, by reporting on their experience with diabetes mellitus type 1 in children during the COVID-19 pandemic. Besides an increased incidence, they also report important differences in patient characteristics between a pre-pandemic control group and the cases occurring during the pandemic.

Sticking with diabetes, Leutner et al. [3] report on different risk profiles and individual risk factors that predispose diabetic patients (predominantly type 2) to a whole range of diabetic complications. They further differentiated these associations between men and women. Meng et al. [4] report on the role of IgG N-glycan profiles in the prediction of progression from either isolated diabetes mellitus type 2 or isolated hypertension to a combination of diabetes mellitus type 2, hypertension, and diabetic comorbidity. They further differentiated these associations between different Chinese ethnic groups.

Bawadi et al. [5] investigated the relationship between different measures for abdominal fat and hypertension in adolescent males. Chang et al. [6] studied the association of serum urate concentrations with progression of kidney disease. They show a clear association for women, but not for men, suggesting different clinical strategies could be warranted for men and women.

Returning to infectious diseases, Tsai et al. [7] show that hepatitis C infection is associated with worse outcomes for diffuse large B-cell lymphoma, suggesting that direct-acting antiviral agents might help improve prognosis for this group of patients. Yang et al. [8] found that tuberculosis survivors who experience lasting ventilatory function disorders are more likely to also experience more respiratory symptoms, more limitations in physical activity, and a worse decline in quality of life.

Wu et al. [9] associated lung cancer prognosis with mode of detection (i.e., screen detected or non-screen detected), smoking status, and several other potential risk factors for poor outcome. For non-smokers, the screened status was one of the predictors, while it was not for smokers. Further, the probability of lung cancer being screen detected was much higher in non-smokers. Since smokers are also less likely to engage in screening, the authors suggest smokers' prognosis might be improved by more effectively motivating them to participate in screening programs.

**Citation:** Middelburg, R.A. Personalized Medicine in Epidemics. *J. Pers. Med.* **2022**, *12*, 583. https:// doi.org/10.3390/jpm12040583

Received: 1 April 2022 Accepted: 2 April 2022 Published: 5 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Finally, Hsiao et al. [10] identified several modifiable predictors of dental problems in children with disabilities, and Liu et al. [11] demonstrate that tongue pressure decline can be used as an indicator for chewing and swallowing problems in older adults.

Together, these eleven papers demonstrate the wide variety of epidemics in which personalization of medicine is affected. As personalization of medicine is important in all fields of medicine, so are all fields of medicine affected by epidemics. We therefore, after the successful completion of this first Special Issue on "Personalized Medicine in Epidemics", now open our second Special Issue on this topic "Personalized Medicine in Epidemics 2.0" for submissions.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


## *Review* **Implementing Personalized Medicine in COVID-19 in Andalusia: An Opportunity to Transform the Healthcare System**

**Joaquín Dopazo 1,2,3 , Douglas Maya-Miles 4,5, Federico García 1,6,7 , Nicola Lorusso 1,8 , Miguel Ángel Calleja 1,9, María Jesús Pareja 1,10, José López-Miranda 1,11,12,13, Jesús Rodríguez-Baño 1,4,14,15 , Javier Padillo 1,4,16,17, Isaac Túnez 1,11,12,18,19,\* and Manuel Romero-Gómez 1,4,5,15,20,\***


**Abstract:** The COVID-19 pandemic represents an unprecedented opportunity to exploit the advantages of personalized medicine for the prevention, diagnosis, treatment, surveillance and management of a new challenge in public health. COVID-19 infection is highly variable, ranging from asymptomatic infections to severe, life-threatening manifestations. Personalized medicine can play a key role in elucidating individual susceptibility to the infection as well as inter-individual variability in clinical course, prognosis and response to treatment. Integrating personalized medicine into clinical practice can also transform health care by enabling the design of preventive and therapeutic strategies tailored to individual profiles, improving the detection of outbreaks or defining transmission patterns at an increasingly local level. SARS-CoV2 genome sequencing, together with the assessment of specific patient genetic variants, will support clinical decision-makers and ultimately better ways to fight this disease. Additionally, it would facilitate a better stratification and selection of patients for clinical trials, thus increasing the likelihood of obtaining positive results. Lastly, defining a national strategy to implement in clinical practice all available tools of personalized medicine in COVID-19 could be challenging but linked to a positive transformation of the health care system. In this review, we provide an update of the achievements, promises, and challenges of personalized medicine in the fight against COVID-19 from susceptibility to natural history and response to therapy, as well as from surveillance to control measures and vaccination. We also discuss strategies to

**Citation:** Dopazo, J.; Maya-Miles, D.; García, F.; Lorusso, N.; Calleja, M.Á.; Pareja, M.J.; López-Miranda, J.; Rodríguez-Baño, J.; Padillo, J.; Túnez, I.; et al. Implementing Personalized Medicine in COVID-19 in Andalusia: An Opportunity to Transform the Healthcare System. *J. Pers. Med.* **2021**, *11*, 475. https://doi.org/10.3390/ jpm11060475

Academic Editor: Rutger A. Middelburg

Received: 22 April 2021 Accepted: 21 May 2021 Published: 26 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

facilitate the adoption of this new paradigm for medical and public health measures during and after the pandemic in health care systems.

**Keywords:** personalized medicine; precision medicine; Covid-19; SARS CoV2; epidemiology; host genetics; viral genome

#### **1. Introduction**

Policymakers, health care leaders and physicians should improve our response to the SARS-CoV-2 pandemic by promoting interdisciplinary collaboration and moving the main research and innovation milestones to clinical practice. The management of this pandemic has been a great challenge addressing health care delivery in stressful conditions due to inadequate capacity, supply shortages, redesigning care, being more transversal and focused on the patient and the virus, thinking in a few days on open intensive care units distributed throughout the campus and managed by several specialties. Innovation and research on COVID-19 required an adaptive process to be translated to clinicians as fast as possible when usefulness was supported in evidence-based medicine. Complex adaptive systems that operate in unpredictable environments should be replaced. Personalized medicine implementation in clinical practice requires a multidisciplinary approach putting together people working on genome sequences, bioinformatics, geneticists or microbiologists and physicians in charge of the patients as well as the management of big-data (BD). A well-defined circuit is mandatory together with a multidisciplinary group able to meet frequently to solve complex clinical situations as severe COVID-19. Precision medicine and beyond, P4 medicine, including predictive, personalized, preventive and participatory medicine, could find the correct scenario in this pandemic. The Andalusian Regional Government allowed us to put together a large population-based health database together with human genomics and viral sequencing to be addressed by Artificial Intelligence (AI) methods to develop robust algorithms able to predict not just the natural history and progression of the disease but also antiviral therapy response [1] and immune response to vaccination [2].

#### **2. Impact of Human Genome on COVID-19**

The clinical course in patients with COVID-19 has been reported as highly heterogeneous. While most people will experience a mild or asymptomatic course, some others may develop progressive and life-threatening bilateral pneumonia and acute respiratory distress syndrome (ARDS). Identifying which patients are at risk to progress to a severe form could reduce the burden of COVID-19, which is currently overloading many health care systems around the world. Factors contributing to disease severity include anthropometric factors (e.g., age, gender, BMI), comorbidities (e.g., hypertension, diabetes) and unhealthy habits (e.g., smoking) [3]. Host genetics studies in COVID-19 have reported genomic variations associated with disease severity in chromosomes 1 (1q22.1), 2 (2p21.1), 3 (3p21.1–3), 6 (6p21.1), 8 (8q24.13), 9 (9q34.1–2), 12 (12q24.1–2), 17 (17q21.3), 19 (19p13.1–3) and 21 (21q21–q22) as well as in specific loci that have been manually selected [4–8]. The COVID Host Genetics Initiative has set up one of the largest communities that are currently generating, sharing and analyzing data to learn the genetic determinants of COVID-19 susceptibility, severity and outcomes. This initiative, which has currently released its fifth data freeze, including data from 46 studies across 19 countries worldwide and analysis looking for genetic determinants of severity (+6000 death or intubated patients vs +1.4 M controls), hospitalization (+13,000 hospitalized patients vs +2 M controls) and infection (≈ 50 K infected patients vs +1.7 M controls), is currently setting up a platform in which researchers will be able to explore the genetic variations that have a deeper impact in SARS CoV-2 infection and severity [9]. An overview of some of the strongest genetic associations described so far for COVID-19 infection and severity, along with some potential gene

candidates, is shown in Table 1. Chromosome 3 genetic variation in the 3p21 locus is the genetic variant that has shown the strongest association in terms of reproducibility to both COVID-19 infection and severity. This region, which is thought to be present in around 30% of people in South Asia and 8% in Europe, has been shown to increase between 1.5- and 2-times approximately an infected person's odds of developing severe COVID-19 [4,7,10]. Carriers of rs10490770, an SNP strongly linked to this chromosome region have increased the odds of several COVID-19 complications, including severe respiratory failure (odds ratio [OR] 2·0, 95%CI 1·6–2·6), venous thromboembolism (OR 1·7, 95%CI 1·2–2·4), and hepatic injury (OR 1·6, 95%CI 1·2–2·0) and higher odds of death or severe respiratory failure, which are especially relevant in patients ≤60 years (OR 2.6, 95%CI 1.8–3.9) compared to those >60 years (OR 1.5 (95%CI 1.3–1.9, interaction *p*-value = 0·04) [11]. rs11385942, another SNP strongly associated with this genetic variation, has shown no association to biomarkers of systemic inflammation, including the C-reactive protein, ferritin, IL-6 and circulating neutrophils and lymphocytes but has been recently associated to increased amounts of the enhanced complement activation, both with C5a and terminal complement complex [12,13]. The 3p21.31 locus contains 17 known protein-coding genes, including SLC6A20, LZTFL1, CCR9, FYCO1, CXCR6, XCR1, CCR1, CCR3, CCR2 and CCR5. None of them seems to have a straightforward connection to SARS-CoV-2 infection or severity to date, besides maybe SLC6A20, which is a transporter regulated by the main SARS-CoV-2 receptor ACE2. However, there are some indirect connections that might be worth highlighting. CCR2 encodes a C-C type chemokine receptor for a chemokine (CCL2) that mediates monocyte chemotaxis. Its expression has recently shown a strong association with the 3p21.31 severe COVID-19-risk variant in certain CD4+ memory T cell subsets and classical monocytes. Patients with severe COVID-19 illness have increased CCR2 expression in circulating monocytes as well as very high levels of CCR2 ligand (CCL2) in bronchoalveolar lavage fluid, leading to the hypothesis that excessive recruitment of CCR2-expressing monocytes may drive pathogenic lung inflammation in COVID-19 [14,15]. LZTFL1/BBS17 is a member of the Bardet-Biedl syndrome (BBS) and encodes a protein involved in protein trafficking to primary cilia. Mutations in LZTFL1 have been reported in human BBS patients, which develop a wide range of pathologies, including obesity, which is so far one of the comorbidities with a stronger link to both COVID-19 infection and severity. It has been recently shown that the deletion of LZTFL1 can cause pleiotropic defects in mice, including obesity. Interestingly, this work links obesity to the expression of LZTFL1 in the brain and demonstrates that the loss of this protein specifically in the brain can lead to leptin resistance [16]. SARS-CoV-2 has been shown to be able to infect the epithelial cells of the gastrointestinal glands of the stomach, duodenum and rectum of COVID-19 patients. The continuous positive detection of SARS-CoV-2 viral RNA in stools suggests that viral particles can be secreted by gastrointestinal cells infected with the virus, and two recent works have proven the ability of this virus to infect human and bat enterocytes in vitro [17–19]. CCR9 is a small intestinal chemokine homing receptor normally found on most mucosal T cells in the gastrointestinal tract, and that has been linked to celiac disease [20] It has been recently shown that the CCR9-CCL25 axis in mice promotes the development of a Th1 population with features of TRM cell that regulates the local immune environment and that CCR9 can exert a protective response against infections in the gastrointestinal tract [21]. CCR9 expression has also been observed in inflammatory cells that are recruited to the lungs and in peripheral blood eosinophils of asthmatic subjects and can be upregulated by stimulation with proinflammatory mediators in human eosinophils-derived cell lines [22]. FYCO1 plays a role in microtubule plus end-directed transport of autophagic vesicles through interactions with the small GTPase Rab7. Although the molecular mechanism of SARS-CoV2 virus infection and spread in the body is not yet disclosed, studies on other beta-coronaviruses show that, upon cell infection, these viruses inhibit macroautophagy/autophagy flux and cause the accumulation of autophagosomes. RNA viruses such as HBV and HCV also modify the autophagy machinery to favor viral replication, translation and propagation. [23] Experiments performed

in hepatocyte cell lines have shown that HCV infection causes inhibition of ras-related protein Rab-7 (Rab7)-dependent endosome–lysosome fusion and promotes the cleavage of the Rab7 adaptor protein RILP (Rab interacting lysosomal protein). This cleavage allows changing the location of this protein to the cell periphery, promoting the export of viral particles outside the cell [24]. XCR1 is a chemokine receptor for XCL1 (Lymphotactin or Lptn). This chemokine is produced predominantly by NK and CD8+ T cells and plays a key role in the tissue-specific recruitment of T lymphocytes [25]. Nasal co-administration of XCL1 and a protein antigen enhances antigen-specific antibody responses both in blood plasma and in mucosal secretions. Lptn as adjuvant-induced antigen-specific CD4+ Th1 and Th2- type cells and IgG1 > IgG2a = IgG2b = IgG3 antibody subclasses [26]. Viral macrophage inflammatory protein-II, a viral protein that inhibits C class chemokines, has been shown to be a potent antagonist able to block the signaling of the XCR1 receptor [27], and the TAT protein of HIV can strongly increase the expression of XCL1 in several cell types [28].

*J. Pers. Med.* **2021**, *11*, 475

**Table 1.** List of SNPs significantly associated with COVID-19 risk of infection, hospitalization, or critical illness that have been identified by three of the major genetic studies that have analyzed host genetics. Chr: Chromosome. LD: linkage disequilibrium.


7


*vs Population control (45.875) Severe Covid: Patients in critical care, being profound hypoxaemic respiratory failure the archetypal presentation.*

(1.18–1.42) 1.6 × 10−8

1.4 (1.25–1.48) 1.6 (1.35–1.87)

4.0 × 10−12 2.3 × 10−8

(1.17–1.41) 2.3 × 10−8

12 rs10735079 113380008 A/G OAS1/3 Severe Covid 1.3

19 rs2109069 rs74956615

4719443 10427721

A/G A/T 21 rs2236757 33252612 A/G IFNAR2 Severe Covid 1.3

DPP9 TYK2

Severe Covid Severe Covid

*J. Pers. Med.* **2021**, *11*, 475

#### **3. Role of Viral Genetic Variants in Covid-19**

Coronaviruses (CoVs) are positive ssRNA viruses, non-fragmented, 26–32 kb belonging to the coronaviridae family. At least four types have been described (α, β, γ and δ). Alpha and beta are pathogenic to mammalians, including humans. Until now, they have been associated with respiratory diseases—α: HCoV-229E and HCoV-NL63 and β: HCoV-OC43, HCoV-HKU1, SARS-CoV-1, MERS-CoV and SARS-CoV-2. The virus did infect and replicate in the cell-expressing ACE-2 receptors. Infection is diagnosed by RT-PCR detecting at least two of these four genes (envelope, spike, nucleocapsid or RNA-dependent RNA polymerase) [29]. The enormous international effort in SARS-CoV-2 sequencing and the subsequent data sharing in sequence databases, along with the popularization of the interactive epidemiological map viewer Nexstrain [30], has uncovered a huge variability spectrum in the viral sequences reported, which are estimated to accumulate nucleotide mutations at a rate of about 1–2 mutations per month [31]. As new sequences cumulated, an increasing number of studies started to describe mutations with an apparent impact on increased infectivity of the virus, such as the well-known D614G variant in the spike protein, or even on increased mortality [32]. Additionally, mutations in the spike have been associated with mABs evasion, which implies a potential risk for vaccine effectivity [33]. Beyond the spike protein, mutations in other proteins, such as the RNA-dependent-RNA polymerase, can reduce the copy fidelity that could result in resistance to antiviral treatments [34] On the other hand, mutations present at the receptor-binding site of the spike protein apparently cause reduced infectivity [33]. ORF8 deletion has also been associated with a milder clinical infection and less post-infection inflammation [35]. Actually, it has been suggested that superspreading events seem to be driving the SARS-CoV-2 pandemic evolution [36]. Moreover, sometimes a more drastic evolutionary event can occur, such as the recently described transmission between humans and mink, and back to humans [37], which triggered a preventive systematic mink slaughter in Denmark because of the presence of variants in the spike protein that might compromise the effectivity of a vaccine [38]. Another case of a new SARS-CoV-2 variant with an unexpectedly large number of mutations in several proteins is the new lineage, B.1.1.7, described in the UK, which is apparently associated with higher transmission rates and mortality [39], or the more recent strains B.1.351 from South Africa and B.1.1.28.1 from Brazil.

The potential effect of many of these mutations has been questioned as speculative and often based on a small number of cases with no much clinical information associated, which casts serious doubts on their actual significance [40]. However, there is an obvious scenario where viral mutations are known to have a potential effect: resistance to antiviral treatments. A noticeable example is remdesivir, an antiviral agent developed against the Ebola virus with demonstrated in vitro activity against SARS-CoV-2. In vitro studies have linked some drug-resistance variants, mainly amino acid changes on RNA-dependent RNA-polymerase (corresponding to residues F480; V557; F480 + V553; F480 + V557), with reduced susceptibility to remdesivir (between 2.4- and 6-fold changes). In clinical cases, some emergent variants promoting drug resistance, mainly in the RNA-dependent RNApolymerase region (D484Y), as a therapeutic target of the drug, have been reported in patients receiving remdesivir [41–43]. Indeed, new variants are continuously emerging every day and may affect drug binding sites and, consequently, may bear the potential to promote drug resistance or escape to current vaccines. Therefore, the benefit of combined genomic and epidemiological analysis for the investigation of health-care-associated COVID-19 seems apparent, as has recently been reported, enabling the detection of cryptic transmission events and identify opportunities for early interventions.

#### **4. Genetic Epidemiology of COVID-19**

In recent years, the European Center for Disease Prevention and Control (ECDC) has published several documents informing about both the roadmap for and state of the art on the current situation of the integration of microbiological data in Public Health surveillance and proposing the need to implement molecular typing and genomic sequencing data for

outbreak surveillance and control [44,45]. The use of genomic surveillance and molecular typing in public health surveillance involves the availability of complementary information to the epidemiological survey and identification of contacts (allowing traceability of the transmission chain). It allows knowing if the cases belong almost unequivocally to the same exposition or transmission line. Using this technique, exposures to a common source can quickly and easily be identified, as demonstrated in the recently reported UK SARS-CoV-2 variant [46]. The creation of a regional geographic database allows to know which pathogens are endemic and which are imported and the identification of new clades, strains or variants that are imported. The analysis of genetic data with phylodynamic methods allows making inferences about the characteristics of the individuals involved in the transmission of the infection and about how contact patterns and the dynamics of risk behaviors affect the flow of transmission through a population [39]. Thus, epidemiological surveillance has to monitor for abrupt changes in rates of transmission or disease severity as part of a systematic process of identifying, response and assessing the impact of variants. A recent example has been the emergence and rapid spread of the above-mentioned new SARS-CoV-2 B.1.1.7 variant with multiple spike protein changes and mutations in other genomic regions associated with higher transmissibility [46].

Although retrospective incorporation has made it possible at the local level to secure the data that had been identified by epidemiologists through surveys when the virus sequencing of confirmed cases has to be carried out faster, the pathways of infection have been prospectively traced, and chains of transmission have been precisely identified, at the hospital level [44,45,47–49]. Technological advances in the classification of pathogens according to the genomic sequence propel us into a new era of massive availability of genomic data at increasingly reduced costs. This information applied to current surveillance systems allows the pathogen causing an infection to be more accurately discriminated, improves the detection of outbreaks and circumscribes the scope and impact of these outbreaks as it allows defining transmission patterns at an increasingly local level. SARS-CoV-2 sequencing has also proven to be useful for studying reinfection. The first case of SARS-CoV-2 reinfection was reported on 24 August 2020 in a Hong Kong citizen re-infected while travelling in Europe. A few other COVID-19 reinfections have been published or deposited in repositories; however, as diagnosing reinfection requires whole-genome sequencing strategies to evaluate the differences between the first and the second strain and samples from both episodes may not be available, it is suspected that reinfection may be a more frequent event than reported.

Another important aspect of genomic surveillance of SARS-CoV-2 is related to vaccination. An effective vaccine should consider the natural variation of the pathogen in order to provide protection with coverage as extensive as possible. The evolution of viruses by mutating epitopes to escape from different pressures has been demonstrated in vitro in the presence of monoclonal antibodies [50] and also in clinical trials [51]. Additionally, some cases of viruses that evade the immune response elicited by vaccines have been described [52,53]. As expected, SARS-CoV-2 can escape in vitro from neutralizing antibodies against the spike protein [54]. However, the recent report of the escape in vitro from a neutralizing COVID-19 convalescent plasma with only three mutations [55], is a serious warning that cannot be ignored and points to the convenience of surveillance that considers immune aspects. The use of immunogenomics, bioinformatics and systems biology helps to understand the basis of interindividual responses to vaccines, both in terms of acquired immunity and adverse effects. The application of these concepts opens the door to the possibility of quantifying and predicting the protective immune response induced by vaccines according to the genomic signature, both of the microorganism and of the vaccine recipient itself [56]. In fact, a bioinformatic approach has recently been described that can predict candidate targets for immune responses to SARS-CoV-2 [57], providing crucial information for understanding human immune responses to this virus and for evaluating vaccine candidates [58]. In fact, these predictions, based on Artificial Intelligence (AI) methodologies [59], allows understanding the individual responses of

patients against the virus [60]. Thus, genomic surveillance and patient screening of risk variants need to be considered for personalized approaches to SARS-CoV-2 vaccination and to prevent possible future vaccine failures [61].

#### **5. Data Science in Health Data Sheet from Large Populations: An Opportunity for COVID-19**

Recent estimates suggest that, while more than 50 years were needed to duplicate all the medical knowledge by 1950, only 70 days were necessary for this increase by the end of 2020 [62], thus providing an idea of the real dimension of current clinical *BD* and their amazing growing pace. This increasing volume of data poses growing challenges for its management, but at the same time, offers an invaluable opportunity for discovery. In fact, the secondary use of stored clinical data is gaining importance progressively and provides a solid substrate for an increasing number of real-world evidence (RWE) studies [63]. Additionally, in parallel to the growth of clinical data repositories, the field of AI has recently experienced a remarkable development, particularly in the case of clinical applications [64]. The AI is starting to integrate into many aspects of medicine with the perspective of optimizing processes, diagnostic procedures and treatments, as well as helping to reduce medical errors [65].

It is worth noting that, despite the short time since the first COVID-19 outbreak, the activity in the development of applications for the retrospective analysis of electronic health records (EHR) by means of artificial intelligence (AI) applied to different aspects of the pandemic is remarkable [66]. With a rapidly growing amount of data available, predictors for different endpoints are being proposed based on different clinical data, that include a medical image [67], clinical text data [68] or, in general, clinical data contained in the EHRs [69]. One of the strengths of AI is its ability to "learn" how individual EHR variables (potential risk factors) can be used and combined among them to produce personalized risk predictions. While conventional approximations (e.g., Cox proportional hazards model) cannot properly combine heterogeneous data from different natures and often incomplete EHRs, modern AI techniques, based on supervised learning, can efficiently learn from such a complex variable dataset and generate risk predictors, as well as update their predictions as data evolve with time [70].

Beyond clinical data from EHRs, the abundance of genomic data in public databases, as well as international initiatives to rapidly increase the biological knowledge on the viral infection process, such as the COVID-19 disease map [71], is fostering innovative applications of AI in the field of drug repurposing [72]. It has recently been demonstrated that a combination of AI methods and mathematical models of the COVID-19 disease map<sup>2</sup> has been able of predicting all targeted drugs for COVID-19 treatment currently in clinical trials [73], opening the door to a new era in which AI-based *in silico* studies will become mainstream. AI may also help in the design of new randomized trials by selecting the most appropriate subpopulations for testing specific drugs according to the best fitted a-priory hypothesis based on the mechanisms of action of the drugs.

#### **6. Ethics, Data Science and Data Sharing in the Times of COVID-19**

Sharing data and results arising from public research projects promotes scientific progress. This concept, widely accepted among research communities, funding bodies, and regulatory agencies, has acquired an unprecedented dimension during the COVID-19 pandemic. The scientific and medical communities have both put in enormous effort to promote data sharing as fast as possible in order to advance our knowledge during the pandemic and design more effective ways to fight it back. Huge efforts have been performed to identify biological and non-biological factors behind the enormous heterogeneity of COVID-19 outcomes and to design strategies able to improve the standard of care given to patients. Sharing clinical and genomic data can improve research efficiency, especially in the genetics field, where the number of samples required to achieve enough statistical power is high. The analysis and re-use of large datasets also allow performing studies that integrate better genomic and phenomic variations, increasing research translationality and

reproducibility. Last but not least, it ensures transparency of previous studies while maximizing the utility of existing datasets. Several public resources have been set in place to access genomic data. 1000 Genomes Project [74], dbGaP [75], European Genome-phenome Archive [76] or the NHGRI AnVIL [77], a US cloud environment for the analysis of large genomic and related datasets, are some examples of databases that have provided services for the archiving and distribution of genetic and phenotypic data resulting from biomedical research projects during this pandemic.

Big data genomic and phenomic databases are necessary to promote personalized medicine but imply certain risks and challenges that need to be taken into consideration that can be roughly grouped into three main categories: issues associated with privacy, the occurrence of incidental findings and challenges associated with the safe management and sharing of genomic data. Researchers need to ensure that patients are well-informed about the benefits and potential risks of data sharing while educating participants about the importance of sample donation as the main pillar of scientific and medical progress. A great example is depicted in some sentences included in the 1000 Genomes Project consent template: "*there may be new ways of linking information back to you that we cannot foresee now. [...] We believe that the benefits of learning more about human genetic variation and how it relates to health and disease outweigh the current and potential future risks, but this is something that you must judge for yourself.*" Strategies to ensure patient privacy go from oversampling (recruiting more individuals than the final number to be included) to not collecting personal data besides sex. Aggregating data, such as allele frequency or allele-presence information, is another strategy that allows protecting participant privacy and also simplifies data sharing and storage. Genomic and phenomic databases usually incorporate controlled-access mechanisms to protect the privacy and confidentiality of research participants, limiting and/or restricting access to personal information. Implement technology safeguards to prevent unauthorized access, use, or disclosure of confidential and private data is a common feature of most of them. Data at EGA, for example, is collected from individuals that authorize data release only for specific research use, and strict protocols govern how information is managed, stored and distributed. EGA databases contain several measures to ensure the security of data, including a regular risk assessment and mitigation, audit logs, cryptography and communication security, among others [76]. IT security has become especially relevant with the transition from locked filing cabinets to digital databases, which bring enormous opportunities for big data analysis but also have an additional set of risks that need to be taken into account. NHGRI AnVIL, for example, uses NIST 800-53 Rev 4 security controls at the Moderate baseline and NIST 800-53 privacy controls documented, security protocols similar to those established in the industry [77].

The COVID-19 Host Genetics Initiative represents a good example of a huge worldwide effort made for COVID-19 genetic and phenotypic data sharing (+160 registered studies +19 countries). This initiative allows submitting individual-level data that includes genetic and clinical phenotype data and also study summary statistics. Individual-level data is shared via the European Genome-phenome Archive (EGA) or via NHGRI AnVIL (US studies). Researchers are allowed to download the metanalysis summary statistics as they are released by the consortia and also have a genome browser that allows them to explore all genetic variations found associated with infection and/or disease severity. Researchers can also apply for study-specific summary statistics and can also request access to the initiative's data deposited on EGA and AnVIL via their respective Data Access Committee, which is composed of the PIs of the studies that have deposited the data (EGA) or managed directly by the NIH (AnVIL). Several genetic studies from Spain, including some launched from Andalusia, are currently contributing with data and samples to this initiative, which is currently working together with the EGA and with the ELIXIR network to establish the EGA Federation network and ensure that data from all countries can be deposited within all national jurisdictions [78].

#### **7. Translating Personalized Medicine into Clinical Practice: The Andalusian Experience**

Andalusia, located in the south of Spain and with 8.5 million inhabitants, is the third most populated region in Europe, and it is larger than half of the countries of the European Union. Remarkably, Andalusia has the whole population under a unique universal electronic health record, thus forming the largest resource of this kind in the European Union. Under this scenario, all the decisions and strategies taken around this huge clinical database acquire enormous relevance. Since 2001, the data recorded by the Andalusian Public Health System (SSPA) are systematically uploaded to the Population Health Base (BPS), making it one of the largest repositories of highly detailed clinical data in the world (with over 13 million registries) [79]. BPS constitutes a unique and privileged environment to carry out large-scale RWE studies. Actually, one of the BPS missions is facilitating the discovery of new biomedical knowledge by means secondary use of clinical data [80], paying special attention to the evaluation of impact in personal data protection [81].

The Andalusian SARS-CoV-2 genomic surveillance project [82] set the ground for the implementation of a clinical circuit for controlling the COVID-19 pandemic as well as other potential future emergent viruses. This project engaged the 16 main tertiary hospitals in Andalusia, along with three research centers with genome sequencing facilities (IBIS, Genyo and CABIMER) and the Bioinformatics Area of the Progress and Health Foundation in a circuit of genomic data production. In parallel, the COVID-19 registries from Public Health and the BPS provide ongoing and retrospective clinical data, respectively, to the Bioinformatics Area, where the data are linked to the genomic data in a circuit of genomic data interpretation. Bioinformatics then provide (i) to the microbiologists at hospitals with information on the lineage, clade and relevant mutations in the virus; (ii) to Public Health with epidemiological data; (iii) to BPS with the viral sequences for further secondary clinical data studies; (iv) to the research community with the viral genome sequences through the European Nucleotide Archive (ENA).

The Andalusian Health Research and Innovation Strategy 2020–2023, presented on 2 September 2020, focused on the improvement of the wellbeing of citizens in the framework of Horizon Europe 2027, including a response to the impact of the challenges of the SARS-CoV2 pandemic. The Regional Health Ministry has been working and collaborating on different initiatives for some time: (a) To promote Digital Clinical records (Diraya®), integrating all the information on people into a Single Health Record and facilitating access to all the services and provisions of the health system, ensuring that all the relevant information is structured. As opposed to systems that merely assemble records, the design of the applications in Diraya shared tables, codes and catalogues; (b) To create the Bioinformatics Research Area [82], on 14 June 2016, to improve the technological support of personalized medicine, genomics and clinical genetics programs in the SSPA; (c) To create, by resolution of the management direction of the Andalusian Health Service 11 March 2018, the information system called BPS of the Andalusian Public Health System, which integrates clinical and epidemiological data from each patient; (d) The resolution of 9 March 2018, of the Public Business Entity Red.es, by which the Agreement with the Andalusian Health Service is published for the application of Information and Communication Technologies in the management of chronicity and continuity of care in the SSPA; (e) On 27 June 2019, at the meeting of the board of the Progress and Health Foundation (FPS), the creation of the B-D Area in Health of Andalusia was approved as an area integrated into the FPS in order to provide the SSPA of a platform of powerful, safe and data analysis tools oriented to health results and the optimization of healthcare processes based on personalized medicine; (f) During the last quarter of 2020, the coordination of the information and communication technologies (ICT) strategy of the SSPA was promoted; (g) The promotion of research on COVID19 in Andalusia, as of December 6, the number of research studies related to COVID-19, presented and/or evaluated in our Research Ethics Committees in Andalusia has been 277, with participation in 24 clinical trials addressing all the spectrum of COVID-19 from epidemiology, diagnosis, biomarkers, genetics (as a contributing study of the COVID-HGI), therapeutic interventions and vaccination. Some of them granted

under the specific support for financing research, development and innovation (R+D+i) in COVID-19; (h) The Ministry of Health and Families, through its General Secretariat for Research, Development and Innovation in Health, has also set up three working groups: (1) prospective studies on the evolution of the pandemic; (2) personalized medicine in Covid-19; and (3) supplement and nutritional intervention against the SARS-Cov2 virus.

Implementing personalized medicine in Covid19 included developing actions to define by means of BD and AI the interaction of genomics, epigenetics, metagenomics and viral sequencing in the development of events such as infection, severe disease, response to treatment and response to vaccination. A joint instruction was carried out on January 2020 from the General Secretariat for Research, Development and Innovation in Health and the Management Directorate of the Andalusian Health Service for the Management of samples in the approach to Personalized Medicine in COVID-19. Healthcare professionals will also have access to SARS-CoV-2 virus complete sequencing study by electronic biochemical request (MPA). The San Cecilio Clinical Hospital for Eastern Andalusia and Virgen del Rocio University Hospital for Western Andalusia were established as reference centers for receiving viral samples (Figure 1A) and sequencing them, respectively (Figure 1B), and the Bioinformatics Area process sequencing data (Figure 1C), joint with COVID registry metadata (Figure 1D), previously collected from the Hospitals (Figure 1E), and reporting back relevant epidemiological information to the COVID registry (Figure 1F) and information on lineages and variants to the Hospitals for supporting clinical decisions (Figure 1G). As previously stated, the clinical data of the Andalusian Health System is stored in the BPS (Figure 1H), but in this case, viral genomes are also stored in BPS (Figure 1I) linked to the rest of the patient's clinical data, offering an unprecedented opportunity for largescale secondary studies and implementation in clinical practice (Figure 1). Finally, the Bioinformatics Area submits the viral sequences to the ENA database, which is available to the scientific community. Since February, more than 2000 whole viral genomes have been sequenced, allowing the construction of a resource that depicts the evolution of SARS-CoV-2 along time and across the geography of Andalusia [82]. This systematic genomic surveillance system has allowed following the increase of the B.1.1.7 since February to become a majority or to detect new VOCs, such as the Brazilian lineage P1 or the South African B.1.351, and VOIs such as the Ugandan variant A.23.1 or others. *J. Pers. Med.* **2021**, *11*, x FOR PEER REVIEW 13 of 17

**Figure 1.** Circuit for COVID19 genomic surveillance. (**A**) Two reference Hospitals, San Cecilio and Virgen del Rocio, collect SARS-CoV-2 samples from Eastern and Western Andalusia, respectively. (**B**) Samples are sent to GENYO or IBIS sequencing facilities. (**C**) Sequencing data are sent to the Bioinformatics Area for processing and (**D**) linking to metadata from the COVID registry, (**E**) previously collected from the Hospitals. (**F**) Bioinformatics reports relevant epidemiological information to the COVID registry and (**G**) information on lineages and variants to the Hospitals for supporting clinical decisions. (**H**) Clinical data on COVID-19 patients recorded in the Hospitals is stored in the BPS. (**I**) Viral genomes are also stored in BPS linked to the rest of the patient's clinical data for further secondary use. **8. Concluding Remarks** The pandemic has pushed us to a new scenario promoting association and relationship between governments and the scientific community at the same time that emerged multidisciplinary teams to take care of this complex disease together with telemedicine to guarantee health care keeping at home. Deep sequencing, bioinformatic area and clinicians working on personalized medicine could help to better understand the interaction between the virus and the host. These tools should be available for physicians able to include in their everyday decision-making process. The increasing need for personalized **Figure 1.** Circuit for COVID19 genomic surveillance. (**A**) Two reference Hospitals, San Cecilio and Virgen del Rocio, collect SARS-CoV-2 samples from Eastern and Western Andalusia, respectively. (**B**) Samples are sent to GENYO or IBIS sequencing facilities. (**C**) Sequencing data are sent to the Bioinformatics Area for processing and (**D**) linking to metadata from the COVID registry, (**E**) previously collected from the Hospitals. (**F**) Bioinformatics reports relevant epidemiological information to the COVID registry and (**G**) information on lineages and variants to the Hospitals for supporting clinical decisions. (**H**) Clinical data on COVID-19 patients recorded in the Hospitals is stored in the BPS. (**I**) Viral genomes are also stored in BPS linked to the rest of the patient's clinical data for further secondary use.

medicine supported by scientific and objective data, big data and AI systems to create algorithms based on individual variables (genomic), the host and the guest (pathogen and patient subject). Public and private investment for the generation and transfer of

**Author Contributions:** Conception and design: J.D., D.M.-M., I.T., M.R.-G.; Analysis & interpretation; All authors contributed; Writing the article: J.D., D.M.-M., M.R.-G., I.T.; Critical revision & modifications: All authors contributed; Data collection: All authors contributed; Funding: All authors contributed. Literature search: All authors contributed; Figures, tables: D.M.-M., J.D. All au-

**Funding:** The authors included in this review have received funding for two COVID-19 projects (COVID GWAs, Premed COVID-19) from the Consejería de Salud y Familias of the Andalusian Government. DMM`s contract is supported by the Andalussian government (Proyectos Estraté-

gicos Fondos Feder PE-0451-2018).

thors have read and agreed to the published version of the manuscript.

#### **8. Concluding Remarks**

The pandemic has pushed us to a new scenario promoting association and relationship between governments and the scientific community at the same time that emerged multidisciplinary teams to take care of this complex disease together with telemedicine to guarantee health care keeping at home. Deep sequencing, bioinformatic area and clinicians working on personalized medicine could help to better understand the interaction between the virus and the host. These tools should be available for physicians able to include in their everyday decision-making process. The increasing need for personalized medicine supported by scientific and objective data, big data and AI systems to create algorithms based on individual variables (genomic), the host and the guest (pathogen and patient subject). Public and private investment for the generation and transfer of knowledge could support the development of high-quality translational and collaborative research to face a threatening situation similar to this terrible pandemic.

**Author Contributions:** Conception and design: J.D., D.M.-M., I.T., M.R.-G.; Analysis & interpretation; All authors contributed; Writing the article: J.D., D.M.-M., M.R.-G., I.T.; Critical revision & modifications: All authors contributed; Data collection: All authors contributed; Funding: All authors contributed. Literature search: All authors contributed; Figures, tables: D.M.-M., J.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors included in this review have received funding for two COVID-19 projects (COVID GWAs, Premed COVID-19) from the Consejería de Salud y Familias of the Andalusian Government. DMM's contract is supported by the Andalussian government (Proyectos Estratégicos Fondos Feder PE-0451-2018).

**Institutional Review Board Statement:** The circuit for COVID-19 genomic surveillance (Premed Covid) is being conducted according to the guidelines of the Declaration of Helsinki, and has been reviewed and approved by the Andalusian Ethics Biomedicine Committee (ethics id: 1954-N-20).

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We acknowledge all members of the "Grupo de Trabajo en Medicina Personalizada contra el COVID-19 de Andalucía."

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study, the writing of the manuscript or in the decision to publish it.

#### **References**


## *Article* **The Impact of SARS-CoV-2 Pandemic on the New Cases of T1DM in Children. A Single-Centre Cohort Study**

**Anca Andreea Boboc 1,†, Carmen Nicoleta Novac 1,†, Maria Teodora Ilie 1,† , Mara Ioana Ies, anu 1,† , Felicia Galos, 1,2,†, Mihaela Bălgrădean 1,2,† , Elena Camelia Berghea 1,2,\* ,† and Marcela Daniela Ionescu 1,2,†**


**Abstract:** Type 1 diabetes mellitus (T1DM) represents one of the most frequent chronic illnesses affecting children. The early diagnosis of this disease is crucial, as it plays a key role in preventing the development of a life-threatening acute complication: diabetic ketoacidosis. The etiopathogenetic role of viral infections has long been suggested and emerging data are pointing towards a complex bidirectional relationship between diabetes and COVID-19. The aim of this study is to assess the impact of the COVID-19 pandemic on the incidence and severity of new T1DM cases in children in Romania. We analyzed the differences between a group of 312 patients diagnosed with T1DM in the period 2003–2019 and a group of 147 children diagnosed during the pandemic. The data were investigated using statistical analysis of a series of relevant variables. The total number of newly diagnosed T1DM increased by 30.08% in the period March 2020–February 2021 compared to the previous years. The patients in the pandemic group had a higher mean age at the onset of T1DM, were less frequently living in an urban area, and presented a higher mean value of HbA1c. Diabetic ketoacidosis at the onset of T1DM was 67.40% more frequent, and a higher percentage of these patients presented with a severe form. The duration of T1DM symptoms did not differ significantly between the two groups. A number of 8 patients associated SARS-CoV-2 infection at the time of T1DM diagnosis.

**Keywords:** early diagnosis; type 1 diabetes mellitus; children; COVID-19; epidemics; diabetic ketoacidosis

#### **1. Introduction**

Type 1 diabetes mellitus (T1DM) is a metabolic disease characterized by the autoimmunemediated destruction of the pancreatic β-cells that leads to a deficit in the production of insulin with various repercussions on the intermediary metabolism. T1DM represents one of the most frequent chronic illnesses in pediatric population and has shown a continuously increasing incidence over the last decades [1]. The most frequent clinical manifestations of T1DM are weight loss, polyuria, and polydipsia; therefore, the rapid recognition of these symptoms plays a key role in making an early diagnosis and thus avoiding the progression towards diabetic ketoacidosis (DKA). DKA is a life-threatening acute complication that can be associated with the onset of T1DM. The incidence of DKA varies widely from one geographic region to another, ranging from 15% to 70% in Europe and North America [2].

The specific etiopathogenetic mechanisms involved in the development of T1DM are not completely understood so far, but they are known to entail a complex interaction between the genome, metabolic processes, immune system characteristics, environmental factors, and microbiome [3]. Among the environmental factors, viral infections have long

**Citation:** Boboc, A.A.; Novac, C.N.; Ilie, M.T.; Ies,anu, M.I.; Galos, , F.; B˘algr˘adean, M.; Berghea, E.C.; Ionescu, M.D. The Impact of SARS-CoV-2 Pandemic on the New Cases of T1DM in Children. A Single-Centre Cohort Study. *J. Pers. Med.* **2021**, *11*, 551. https://doi.org/ 10.3390/jpm11060551

Academic Editor: Rutger A. Middelburg

Received: 10 May 2021 Accepted: 11 June 2021 Published: 13 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

been considered triggers of the autoimmune process leading to the development of T1DM. Various mechanisms could explain the role that viral infections play in promoting the onset of T1DM in children. In patients with a predisposing background of genetic and immunological factors, the exposure to certain viruses might accelerate the final processes leading to the onset of the disease [4].

COVID-19 is a highly infectious respiratory disease determined by SARS-CoV-2 that spread rapidly across the world, having been declared a pandemic on 11 March 2020. In children, the disease has a milder course, being asymptomatic in many cases and having a mortality rate of 0.18% in patients hospitalized for COVID-19 [5]. In Romania, the first case of SARS-CoV-2 infection was reported on 26 February 2020. On 28 February 2021, the total number of confirmed cases exceeded 800,000 (41.77 per 1 million people, a prevalence that made the country rank 24th in the world by the number of infections). Over 40,000 (5%) of these infections were reported among children [6].

Although extensive evidence has been gathered since the beginning of the COVID-19 pandemic about the complex bidirectional relation between SARS-CoV-2 infection and diabetes, data regarding T1DM specifically remain sparse. A multicentre study reported an 80% increase in new T1DM cases in children, some of these patients having a confirmed history of SARS-CoV-2 infection [7]. Several studies reported a significant increase in the frequency of DKA present at the diagnosis of T1DM, with a higher percentage of the severe form of this metabolic complication [7–12]. It was hypothesized that an explanation for this phenomenon is the delay in the presentation of these patients caused by the parents' fear of accessing the healthcare system [9,10]. However, reports from the United Kingdom showed a short duration of symptoms in the majority of children newly diagnosed with T1DM [7]. More information about the complex interactions between T1DM and COVID-19 was gathered due to a study which showed that SARS-CoV-2 is able to infect cultured human pancreatic islets and replicate inside these cells, altering the subcellular morphology and the response to glycemic variations [13]. The pancreatic endocrine cells strongly express the angiotensin converting enzyme 2 (ACE2) receptor, a main binding site for SARS-CoV-2 [14].

Another possible mechanism that could play a role in the observed increased number of patients with T1DM and other autoimmune disorders during the pandemic could be the decreased vitamin D levels, a possible consequence of strict lockdown measures with people spending significantly less time outside. Vitamin D, due to its immunomodulatory properties, shows multiple beneficial effects in T1DM [15]. This idea is supported by a study reporting that 2000 UI cholecalciferol per day reduced T1DM incidence by 78% [16]. Moreover, the decision of the Finish authorities to fortify dietary milk products with vitamin D had led to a decrease in the number of T1DM cases [16]. From another point of view, the capacity of vitamin D to enhance the immune system might ensure a protective effect from developing severe forms of SARS-CoV-2 infection and reduce the risk of contracting the disease [17].

The aims of this study were to assess a possible increase in the number of T1DM new cases in children from Bucharest and the surrounding areas during the COVID-19 pandemic, compared to previous years, to analyze the differences in patient and disease characteristics between the pre-pandemic and pandemic periods and to evaluate predictors (including presentation time) of DKA associated with the onset of the disease and DKA severity.

#### **2. Materials and Methods**

#### *2.1. Study Population*

We performed an observational retrospective cohort study that included pediatric patients with T1DM from Marie Curie Emergency Children's Hospital, Bucharest. Our institution is a reference centre for pediatric T1DM, all the new cases from Bucharest and the south-east region of Romania (approximately one quarter of the 3,895,000 children in Romania are living here) being managed here. In order to analyze the temporary evolution of newly diagnosed pediatric T1DM cases, we extracted from the electronic hospital registry the numbers of all new T1DM cases diagnosed in patients younger than 18 years old, for each month of 2018, 2019, and 2020. Patients with other types of diabetes were excluded.

In order to further assess patients and disease differences between children diagnosed during (March 2020–February 2021; pandemic group) and before the COVID-19 pandemic (2003–2019; pre-pandemic group), we split the study sample into 2 groups. The pre-pandemic group included patients under 18 years old, diagnosed with T1DM, with extensive follow-up and complete and easily available medical records. The pandemic group included all the children diagnosed with T1DM in the period mentioned before; thus, the follow-up criteria being impossible to apply considering the recent onset.

#### *2.2. Patients and Disease Characteristics*

The data were collected from the digital and paper-based medical records of the hospital. For each patient, we gathered the following information: age at the onset of T1DM, gender, living area, glycosylated hemoglobin (HbA1c) levels, pH, and duration of T1DM specific symptomatology. The diagnosis of T1DM was made on the following criteria: fasting plasma glucose levels higher than 126 mg/dL or symptoms of hyperglycemia (polyuria, polydipsia, unexplained weight loss) with a random plasma glucose level ≥200 mg/dL or 2-h plasma glucose ≥200 mg/dL during an oral glucose tolerance test.

The pH was determined at the time of presentation, and a level lower than 7.3 was defined as DKA. The degree of severity of DKA was classified according to 2018 International Society for Paediatric and Adolescent Diabetes (ISPAD) guidelines: mild DKA (pH < 7.3), moderate DKA (pH < 7.2), severe DKA (pH < 7.1).

HbA1c was measured through an immunoturbidimetric method (DCCT standardized: Diabetes Control and Complications Trial and NGSP certificated: National Glycohemoglobin Standardization Program). The normal levels were 4.8–6.4%.

#### *2.3. Statistical Analyses*

Patients and disease characteristics were reported as absolute and relative frequencies for categorical variables, medians, and interquartile ranges (IQR) for non-normally distributed continuous variables and means and standard deviations for normally distributed continuous variables. Tests of normality were conducted with the Kolmogorov-Smirnov and Shapiro-Wilks tests. The variables were reported for the overall study sample, and separately, according to the time of diagnosis (pre-pandemic/pandemic group), according to DKA presence and DKA severity. Differences between the groups were tested using χ 2 statistics for categorical variables, Wilcoxon or Kruskal-Wallis rank sum tests for nonnormally distributed continuous variables, and *t*-tests or ANOVA for normally distributed continuous variables. The eight cases that presented SARS-CoV-2 infection at the time of T1DM diagnosis were described separately based on multiple variables.

Logistic regression was performed to determine the effects of patient characteristics and time of presentation (independent variables) on the likelihood of DKA being present at the time of T1DM diagnosis (dependent variable). Proportional odds logistic regression was performed to assess effects on the likelihood of increased DKA severity (3-level ordered variable) for the subset of patients presenting DKA. Univariable and multivariable analyses were conducted for both outcomes. Patients' characteristics were included in multivariable analysis, regardless of their significance in univariable analysis, in order to estimate the independent effect of time of presentation on the outcomes.

Furthermore, to address selection bias induced by the selection of patients in the pre-pandemic group, the analyses were repeated on a propensity score-matched dataset. Propensity scores for being diagnosed during or before the pandemic period were estimated with logistic regression, using all baseline measured characteristics. Patients were matched using 1:1 greedy nearest neighbor matching without replacement. The balance of individual characteristics between the pandemic and pre-pandemic groups, before and after matching, was assessed using standardized mean differences (SMD). In the matched dataset, containing patients with similar characteristics, the effects of time of presenta-

tion (pre-pandemic or pandemic) on presence of DKA and DKA severity were estimated. Residual differences in observed characteristics between the pre-pandemic and pandemic groups were accounted for by covariate adjustment. groups were accounted for by covariate adjustment. Restricted cubic splines were used to test for non-linear relationships between continuous independent variables (age, HbA1c) and the outcomes of interest in the model

Furthermore, to address selection bias induced by the selection of patients in the prepandemic group, the analyses were repeated on a propensity score-matched dataset. Propensity scores for being diagnosed during or before the pandemic period were estimated with logistic regression, using all baseline measured characteristics. Patients were matched using 1:1 greedy nearest neighbor matching without replacement. The balance of individual characteristics between the pandemic and pre-pandemic groups, before and after matching, was assessed using standardized mean differences (SMD). In the matched dataset, containing patients with similar characteristics, the effects of time of presentation (pre-pandemic or pandemic) on presence of DKA and DKA severity were estimated. Residual differences in observed characteristics between the pre-pandemic and pandemic

*J. Pers. Med.* **2021**, *11*, x FOR PEER REVIEW 4 of 13

Restricted cubic splines were used to test for non-linear relationships between continuous independent variables (age, HbA1c) and the outcomes of interest in the model development phase. The estimated associations were reported as odds ratios (OR) and 95% confidence intervals (95% CI). To assess if the predictive strength of multivariable models was increased by adding time of presentation as a covariate, the proportion of explained outcome variability was calculated using model Nagelkerke's pseudo-R<sup>2</sup> . development phase. The estimated associations were reported as odds ratios (OR) and 95% confidence intervals (95% CI). To assess if the predictive strength of multivariable models was increased by adding time of presentation as a covariate, the proportion of explained outcome variability was calculated using model Nagelkerke's pseudo-R2. The percentage of missing values for each variable was reported in the descriptive tables and regression analyses were conducted on complete cases.

The percentage of missing values for each variable was reported in the descriptive tables and regression analyses were conducted on complete cases. The data were statistically analyzed using Microsoft Excel ( Microsoft Corp. S.R.L, Michigan, USA), GraphPad Prism? (version 6.0; GraphPad Prism Inc., California, USA), (;

The data were statistically analyzed using Microsoft Excel (Microsoft Corp. S.R.L, MI, USA), GraphPad Prism? (version 6.0; GraphPad Prism Inc., CA, USA), R (R version 4.0.3; R Foundation for Statistical Computing, Vienna, Austria) and R Studio (version 1.1.463; R Studio, Inc., Boston, MA, USA) considering statistical significance at a two-sided *p*-value of 0.05. Nearest neighbor matching was performed using the MatchIt package in R. ), R (R version 4.0.3; R Foundation for Statistical Computing, Vienna, Austria) and R Studio (version 1.1.463; R Studio, Inc., Boston, Massachusetts, USA) considering statistical significance at a two-sided *p*-value of 0.05. Nearest neighbor matching was performed using the MatchIt package in R.

#### **3. Results 3. Results**

#### *3.1. Analysis of Monthly New Type 1 Diabetes Cases 3.1. Analysis of Monthly New Type 1 Diabetes Cases*

In March and April 2020 the number of new T1DM cases was lower than the same period of 2018 and 2019. Between May 2020 and February 2021 the monthly number of new T1DM cases was higher, with a mean of 13.2 cases/month compared to the same period of the last two years with a mean of 9.4 cases/month (*p*-value = 0.015). The total number of newly diagnosed T1DM increased by 30.08% from 113, between March 2019 and February 2020, to 147 cases reported between March 2020 and February 2021 (Figure 1). In March and April 2020 the number of new T1DM cases was lower than the same period of 2018 and 2019. Between May 2020 and February 2021 the monthly number of new T1DM cases was higher, with a mean of 13.2 cases/month compared to the same period of the last two years with a mean of 9.4 cases/month (*p*-value = 0.015). The total number of newly diagnosed T1DM increased by 30.08% from 113, between March 2019 and February 2020, to 147 cases reported between March 2020 and February 2021 (Figure 1).

**Figure 1.** The number of new T1DM cases per month for the periods March 2018–February 2019, March 2019–February 2020 and March 2020–February 2021.

#### *3.2. Demographic Analysis*

The study included a number of 459 pediatric patients with a mean age at diagnosis of 7.59 years, the age group with the highest percentage of patients being 8–11 years (32.68%). The male/female ratio was 53%/47% (245/214), the background was urban in 352 cases (77%), and rural in 107 (23%).

2020 and March 2020–February 2021.

#### *3.3. Differences between the Pre-Pandemic and the Pandemic Group of Patients* 10.00] years, *p* < 0.001). The sex ratio was lower in the pandemic group compared to the

*3.3. Differences between the Pre-Pandemic and the Pandemic Group of Patients* 

*J. Pers. Med.* **2021**, *11*, x FOR PEER REVIEW 5 of 13

**Figure 1.** The number of new T1DM cases per month for the periods March 2018–February 2019, March 2019–February

*3.2. Demographic Analysis* 

cases (77%), and rural in 107 (23%).

The patients in the pandemic group had a higher median age at the onset of T1DM than those diagnosed during the period 2003–2019 (7.20 [7.07, 7.30] years vs. 7.00 [3.00, 10.00] years, *p* < 0.001). The sex ratio was lower in the pandemic group compared to the prepandemic group (M/F = 1.19 and 1.04, respectively, *p* = 0.55). The percentage of patients living in an urban area decreased from 83.65% before 2020 to 61.90% during the pandemic (*p*-value < 0.001). There was no significant difference in the mean duration of T1DM symptoms before presentation between the two groups (25.39 ± 2.23 days vs. 24.52 ± 1.83 days, *p* = 0.51). In the pandemic group the values of HbA1c were higher than those of the patients from the pre-pandemic group (mean HbA1C 12.47 ± 2.19% vs. 11.32 ± 2.18%, *p* < 0.001) (Table 1, Figure 2). In the pandemic group, one patient was diagnosed with both autoimmune thyroiditis and celiac disease at the time of T1DM diagnosis. pre-pandemic group (M/F = 1.19 and 1.04, respectively, *p* = 0.55). The percentage of patients living in an urban area decreased from 83.65% before 2020 to 61.90% during the pandemic (*p*-value < 0.001). There was no significant difference in the mean duration of T1DM symptoms before presentation between the two groups (25.39 ± 2.23 days vs. 24.52 ± 1.83 days, *p* = 0.51). In the pandemic group the values of HbA1c were higher than those of the patients from the pre-pandemic group (mean HbA1C 12.47 ± 2.19% vs. 11.32 ± 2.18%, *p* < 0.001) (Table 1, Figure 2). In the pandemic group, one patient was diagnosed with both autoimmune thyroiditis and celiac disease at the time of T1DM diagnosis. **Table 1.** Characteristics of the study sample according to the time of diagnosis. **Group of Patients Pre-Pandemic Group Pandemic Group n (%)** *p* **Values a** 

The study included a number of 459 pediatric patients with a mean age at diagnosis of 7.59 years, the age group with the highest percentage of patients being 8–11 years (32.68%). The male/female ratio was 53%/47% (245/214), the background was urban in 352

The patients in the pandemic group had a higher median age at the onset of T1DM than those diagnosed during the period 2003–2019 (7.20 [7.07, 7.30] years vs. 7.00 [3.00,

**Table 1.** Characteristics of the study sample according to the time of diagnosis. Gender

**n (%)** 


<sup>a</sup> *p* values correspond to χ 2 statistics for categorical variables, Wilcoxon rank sum tests for non-normally distributed continuous variables (age, pH) and t-tests for normally distributed continuous variables (HbA1c), comparing the two subgroups. The *p* value assesses compatibility with the null hypothesis of no differences between subgroups. ous variables (HbA1c), comparing the two subgroups. The *p* value assesses compatibility with the null hypothesis of no differences between subgroups.

**Figure 2.** Comparison between the pre-pandemic (pre-PG) and pandemic group(PG) of patients based on several parameters. (**A**) Age of onset; (**B**) Symptoms duration; (**C**) HbA1c.

> The proportion of DKA at the onset of the disease increased during the pandemic with 67.40%, from 39.42% in the pre-pandemic group to 65.99% in the pandemic group (OR = 2.98, CI 95% = 1.97–4.49, *p* < 0.0001). In the pandemic group, a higher percent of DKA cases developed the severe form compared to the pre-pandemic group (42.27% vs. 26.83%, OR = 1.99, CI 95% = 1.13–3.51, *p* = 0.016). (Table 2, Figure 3).

eters. (**A**) Age of onset; (**B**) Symptoms duration; (**C**) HbA1c.


**Table 2.** DKA presence and severity at the time of diagnosis in the pre-pandemic and pandemic group. **Table 2.** DKA presence and severity at the time of diagnosis in the pre-pandemic and pandemic group.

The proportion of DKA at the onset of the disease increased during the pandemic with 67.40%, from 39.42% in the pre-pandemic group to 65.99% in the pandemic group (OR = 2.98, CI 95% = 1.97–4.49, *p* < 0.0001). In the pandemic group, a higher percent of DKA cases developed the severe form compared to the pre-pandemic group (42.27% vs.

*J. Pers. Med.* **2021**, *11*, x FOR PEER REVIEW 6 of 13

**Figure 2.** Comparison between the pre-pandemic (pre-PG) and pandemic group(PG) of patients based on several param-

26.83%, OR = 1.99, CI 95% = 1.13–3.51, *p* = 0.016). (Table 2, Figure 3).

**Figure 3.** DKA presence and severity in the pre-pandemic group of patients (pre-PG) and the pandemic group (PG). **Figure 3.** DKA presence and severity in the pre-pandemic group of patients (pre-PG) and the pandemic group (PG).

#### *3.4. Association of Patient Characteristics with DKA Diagnosis and Severity*

*3.4. Association of Patient Characteristics with DKA Diagnosis and Severity*  In order to further assess the development of DKA the sample was divided into two groups based on the presence and absence of DKA at the time of diagnosis. The gender distribution and the background distribution did not differ significantly between the two In order to further assess the development of DKA the sample was divided into two groups based on the presence and absence of DKA at the time of diagnosis. The gender distribution and the background distribution did not differ significantly between the two groups.

groups. The median age was lower in the group of patients that presented with DKA (7.00 [3.00, 10.00] vs. 8.00 [5.00, 11.00], *p* = 0.006). Also, in this group HbA1c had a higher mean value (11.96 ± 2.19 vs. 11.45 ± 2.27), and a significantly higher percent of patients were The median age was lower in the group of patients that presented with DKA (7.00 [3.00, 10.00] vs. 8.00 [5.00, 11.00], *p* = 0.006). Also, in this group HbA1c had a higher mean value (11.96 ± 2.19 vs. 11.45 ± 2.27), and a significantly higher percent of patients were diagnosed during the pandemic (44.1% vs. 20.9%) (Table 3).

diagnosed during the pandemic (44.1% vs. 20.9%) (Table 3).  **Table 3.** Characteristics of the study sample, overall and according to DKA diagnosis.


<sup>a</sup> *p* values correspond to χ 2 statistics for categorical variables, Wilcoxon rank sum tests for non-normally distributed continuous variables (age, pH) and t-tests for normally distributed continuous variables (HbA1c), comparing the two subgroups. The *p* value assesses compatibility with the null hypothesis of no differences between subgroups.

> The group of patients that presented DKA at the time of diagnosis was further divided into three subgroups based on DKA severity. Gender distribution, the median age at diagnosis, background distribution, mean HbA1c, and the time of diagnosis (i.e., prepandemic/pandemic) did not differ significantly between the two groups (Table 4).


**Table 4.** Characteristics of the DKA subset, according to DKA severity.

<sup>a</sup> *p* values correspond to χ 2 statistics for categorical variables, Kruskal-Wallis rank sum tests for non-normally distributed continuous variables (age, pH) and ANOVA for normally distributed continuous variables (HbA1c), comparing the three subgroups. The *p* value assesses compatibility with the null hypothesis of no differences between subgroups.

> In univariable analysis, a significant association with a higher likelihood of being diagnosed with DKA was found for younger age (OR 0.94 per year increase in age, 95% CI 0.90–0.98), higher levels of HbA1c (OR 1.11 per percent increase in HbA1c, 95% CI 1.02–1.21) and presentation during COVID-19 pandemic (OR 2.98, 95% CI 1.99–4.52 compared to the pre-pandemic period). All three variables remained significantly associated with the likelihood of DKA in multivariable analysis, with the same direction of effects and larger effect estimates: OR 0.87 per year increase in age (95% CI 0.82–0.92), OR 1.17 per percent increase in HbA1c (95% CI 1.06–1.30) and OR 3.23 for presentation during the pandemic, compared to the pre-pandemic period (95% CI 2.06–5.15). Gender and urban/rural background were not significantly associated with the likelihood of being diagnosed with DKA in our sample. Multivariable model Nagelkerke's pseudo-R2 for DKA diagnosis was 8% without time of presentation and 15% with time of presentation (Table 5). In sensitivity analysis (*n* = 439), excluding patients with active SARS-CoV 2 infection at the time of T1DM diagnosis, the effects of the variables remained similar in both magnitude and significance: OR 3.15 for presentation during the pandemic, compared to the pre-pandemic period (95% CI 2.00–5.05) in multivariable analysis.

**Table 5.** Logistic regression of variables associated with DKA diagnosis.


<sup>a</sup> Complete-case analysis: 446 patients included (13 excluded due to missing HbA1c data). <sup>b</sup> No significant non-linear relationships (modeled with splines) between the continuous covariates (Age, HbA1c) and outcome were detected. <sup>c</sup> Nagelkerke's pseudo-R<sup>2</sup> of the multivariable model increased from 0.08 to 0.15 when adding presentation time in the model.

> In the subgroup of patients with DKA, the only variable that demonstrated a significant association with increased DKA severity was presentation during the pandemic (OR 1.85 compared to the pre-corona period, 95% CI 1.07–3.23 in multivariable analysis). Multivariable model Nagelkerke's pseudo-R2 for DKA severity was 1% without time of presentation and 4% with time of presentation. (Table 6). In sensitivity analysis (*n* = 212), excluding patients with active SARS-CoV 2 infection at the time of T1DM diagnosis, presentation during the pandemic was similarly the only variable to display a significant effect on increased DKA severity: OR 1.99, compared to the pre-pandemic period (95% CI 1.13–3.51) in multivariable analysis.


**Table 6.** Proportional odds logistic regression of variables associated with increased DKA severity.

<sup>a</sup> Complete-case analysis: 218 observations included (two excluded due to missing HbA1c data). <sup>b</sup> No significant non-linear relationships (modeled with splines) between the continuous covariates (Age, HbA1c) and outcome were detected. <sup>c</sup> Nagelkerke's pseudo-R<sup>2</sup> of the multivariable model increased from 0.01 to 0.04 when adding presentation time in the model.

> The cohort of propensity score-matched patients included 294 patients and consisted of all the patients in the pandemic group and their matched "controls" from the prepandemic group. The balance of individual characteristics between the groups, before and after matching, is presented in Supplementary Table S1 and Supplementary Figure S1.

> The positive effect of presentation during the pandemic on presence of DKA at the time of T1DM diagnosis was significant in the matched cohort, in both univariable analysis (OR 3.2 compared to the pre-pandemic period, 95% CI 2–5.2, *p* < 0.001) and multivariable analysis with covariate adjustment (OR 3.4 compared to the pre-pandemic period, 95% CI 2.1–5.6, *p* < 0.001).

> In the subset of 152 patients with DKA from the matched cohort, presentation during the pandemic was significantly associated with increased DKA severity, in both univariable analysis (OR 2.03 compared to the pre-pandemic period, 95% CI 1.1–3.77, *p* = 0.024) and multivariable analysis with covariate adjustment (OR 2.03 compared to the pre-pandemic period, 95% CI 1.09–3.82, *p* = 0.027).

#### *3.5. SARS-CoV-2 Infection at the Onset of Type 1 Diabetes Mellitus*

Eight of the patients diagnosed with T1DM in 2020 presented SARS-CoV-2 infection at the same time. The diagnosis of COVID was made on the positive results of the RT-PCR for SARS-CoV-2 tests. These patients' ages ranged from 1 year to 17 years, with a mean value of 8.87 and the male/female ratio was 1/1. None of the patients presented any significant past medical history. Regarding the presence of type 1 diabetes or other autoimmune diseases in the family medical history, two patients' mothers were diagnosed with autoimmune thyroiditis and one patient's cousin was diagnosed with T1DM. In half of these cases, the duration of T1DM characteristic symptoms (polyuria, polyphagia, weight loss) was 14 days; in two cases, it was 7 days; one patient had a 60 days history of the specific clinical manifestations, and one other patient presented no symptoms before the diagnosis. Among the eight patients, seven presented DKA which was mild in three cases, moderate in three, and severe in one patient. The HbA1c had a mean value of 12.19%, with a minimum of 10.70% and a maximum of 13.80%. The COVID-19 infection was completely asymptomatic in four of the children, one patient aged 2 presented fever (maximum temperature = 40 ◦C), eating refusal, and somnolence, two patients experienced vomiting, and abdominal pain (that could have been determined by SARS-CoV-2 infection, as well as by DKA) and one patient was diagnosed with stomatitis and streptococcal pharyngitis. The presence of inflammation based on an increased level of C reactive protein was detected in five patients, with slightly elevated levels in three patients and a maximum value of 39.98 mg/dL in the patient that presented fever. Vitamin D deficiency was identified in two of the five measurements performed (Table 7).



#### **4. Discussion**

Given the paucity of available information on the subsiding complex relations between SARS-CoV-2 infection and other illnesses, it is important to explore the possible effects of the pandemic on the incidence, severity, pattern, and evolution of the existing conditions in order to adjust the healthcare response to the present situation and to anticipate what the future may hold. Taking that into consideration, this article assessed the evolution of the number and severity of new T1DM cases during the pandemic.

Between March 2020 and February 2021 there was an increase of 30.08% in T1DM cases compared to the same period of the previous year (147 vs. 113). This observation is in agreement to another study showing an increase in new T1DM cases in children in 2020 [7]. It is important to mention that the changes in the functioning and organization of the healthcare system as a response to the pandemic could not have contributed significantly to the increase in the number of new cases, because of the fact that the majority of new T1DM cases from the south-east region of Romania were referred to our hospital both before and during the pandemic.

However, in the first 2 months of the pandemic, the number of new T1DM cases was lower than that reported in the same period of the previous year, the same difference being reported in another study which analyzed the number of children diagnosed with T1DM in Italy in the first 2 months of the pandemic [10]. In those first weeks of the pandemic, in our country a strict lockdown was imposed, with the population being allowed to leave the house only for a few clearly defined reasons like buying food and essential products, going to work in cases where the job was an essential one, and for mandatory medical assistance. Therefore, the lifestyle of children went through significant changes, one consequence being a significantly reduced exposure to common infectious microorganisms. Considering the potential of viral infections to act as triggers for the development of T1DM, the decreased incidence of this category of diseases might be an important cause of the lowered T1DM new cases at the beginning of the pandemic.

The physiopathological mechanisms that could play a role in the increase of the T1DM cases are yet to be completely understood but important progress in this direction was made by gathering evidence that SARS-CoV-2 is able to infect and replicate in pancreatic β-cells, leading to impaired function [13].

The children included in this study and diagnosed between March 2020 and February 2021 had a higher risk of presenting with DKA, time of presentation being significantly associated with increased likelihood of DKA diagnosis and higher DKA severity when adjusting for known predictors. While younger age, higher levels of HbA1C and presentation during the pandemic were significantly associated with a higher likelihood of being diagnosed with DKA, presentation during the pandemic was the only variable significantly associated with increased DKA severity. Previous studies found that the risk of developing DKA is higher in females, children with migration background, lower socio-economic

status, children under 3 years old, early teenage years, patients without any first-degree relatives with T1DM and delayed diagnosis [18–20].

The addition of time of presentation to the multivariable explanatory model resulted in increased predictive performance, as reflected by a 2-fold increase in Nagelkerke's pseudo-R<sup>2</sup> for the DKA diagnosis model and a 4-fold increase for the DKA severity model. However, the absolute proportions of explained variability are relatively low even when including time of presentation, which indicates poor predictability of the outcomes given the variables we have considered. Although the purpose of this analysis is effect estimation rather than prediction, weak [18] predictive performance might indicate that other important predictors of outcome exist, which we did not consider in our analysis. The strong effect we have identified for time of presentation could be confounded by such an important unobserved predictor, and as such our results should be interpreted with caution.

An increase in DKA frequency and severity has been prior reported in other countries [7–12]. As DKA is a life-threatening condition that requires complex management, sometimes possible only in an intensive care unit, the increased frequency of this metabolic complication could lead to extensive pressure on an already overwhelmed healthcare system. On the other hand, pointing out that an increase in the incidence of a chronic pediatric disease might be one of the consequences of SARS-CoV-2 exposure might change the paradigm that children are not remarkably affected by the pandemic.

The duration of symptoms did not differ significantly during the pandemic compared to the previous years, this observation supporting the fact that an early diagnosis was made in these cases and being in contradiction with the hypotheses suggesting that the increase in the DKA frequency is caused by a delay in the diagnosis of T1DM [9,10]. A multicentre study performed in the United Kingdom also reported a short duration of symptoms in the majority of children [7].

Another interesting observation was that the mean age at diagnosis was higher in the pandemic group. Considering that it has long been recognized that DKA is more frequent in children of younger ages at the onset of T1DM [2], an increase in DKA frequency would be expected to associate a lower mean age of the patients.

One limitation of this study is represented by the fact that the pre-pandemic group included only the children that were under observation in our hospital during the evolution of T1DM and not all the cases diagnosed, while the pandemic group included all children diagnosed with T1DM in our hospital during the pandemic. This has been addressed in the effect estimation analysis for presence of DKA and DKA severity by using a propensity score matched cohort, in which every patient diagnosed in the pandemic period was matched to a patient diagnosed in the pre-pandemic period. In the matched cohort, the effect of presentation during the pandemic remained significantly associated with both higher likelihood of DKA at diagnosis and higher DKA severity. Another limitation is that the SARS-CoV-2 specific antibodies were not available for the new cases of T1DM. Additionally, the duration of symptoms might be insufficiently accurate as the information was gathered from the patients' parents or guardians. The main strength of this study, on the other hand, is represented by the large sample used: a relatively high number of patients were included in the analysis. Moreover, all patients included in the study were diagnosed with T1DM in the same hospital and all the blood tests referenced were performed in the same laboratory, which ensures consistency.

As explained at the beginning of this paper, it is highly necessary to further investigate the relation between the increase in the number of cases and severity of T1DM cases in children and SARS-CoV-2 infection. More data supporting this link could be gathered by analyzing the presence of SARS-CoV-2 antibodies as evidence of previous infection in children newly diagnosed with T1DM.

#### **5. Conclusions**

Considering the fact that T1DM is one of the most frequent chronic pediatric illnesses and that its main acute complication, DKA is a life-threatening condition, it is important to analyze the repercussions that the new global pandemic might determine. Between March and February 2020, in Bucharest (Romania) and the surrounding areas there was a 30.08% increase in the total number of new T1DM cases in children. The patients that presented during the pandemic had a lower mean age, higher values of HbA1c, and were living more frequently in a rural area compared to the previous years. The number of diabetic ketoacidosis presentations also showed an increase during the pandemic, patients diagnosed after March 2020 having a significantly higher risk of developing DKA, when adjusting for other predictors. Moreover, a higher percentage of children presented a severe form of this metabolic complication. In the subgroup of patients that presented with DKA, the time of diagnosis (i.e., during the pandemic) was the only variable that showed a significant association with increased severity. These findings are supporting the idea that an increase in the incidence and severity of T1DM in children could be one of the consequences of COVID-19 pandemic. The duration of T1DM symptoms did not differ significantly between the pandemic and the pre-pandemic group; therefore, a delay in the presentation is unlikely to have been the main cause of the increased number of DKA presentations. Among the patients with the T1DM onset during the pandemic, 8 of them also presented SARS-CoV2 infection at the time of diagnosis.

Taking all these results into consideration we believe that further investigations are needed in order to discover other possible particular characteristics of T1DM cases diagnosed during the pandemic. Also, we are taking into consideration that a prospective cohort study including these patients could bring relevant additional data.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/jpm11060551/s1, Figure S1: title, Table S1: title, Video S1: title.

**Author Contributions:** Conceptualization, A.A.B., C.N.N., M.T.I., F.G., M.B., E.C.B. and M.D.I.; Data curation, A.A.B., C.N.N., M.T.I., M.I.I., F.G. and E.C.B.; Formal analysis, A.A.B., M.T.I., M.I.I. and M.B.; Investigation, A.A.B., C.N.N., M.T.I. and F.G.; Methodology, A.A.B., C.N.N., M.T.I., M.B., E.C.B. and M.D.I.; Project administration, M.B., E.C.B. and M.D.I.; Resources, E.C.B. and M.D.I.; Software, M.T.I. and M.I.I.; Supervision, A.A.B., C.N.N., F.G., M.B., E.C.B. and M.D.I.; Validation, F.G., M.B., E.C.B. and M.D.I.; Visualization, A.A.B., C.N.N., M.T.I., F.G., M.B. and M.D.I.; Writing—original draft, A.A.B., C.N.N., M.T.I., M.I.I., E.C.B. and M.D.I.; Writing—review & editing, A.A.B., C.N.N., M.T.I., M.B., E.C.B. and M.D.I. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** The research was approved by the Ethics Committee of "Marie Curie" Emergency Children's Hospital.

**Informed Consent Statement:** Informed consent was obtained from all the patients included in the study.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Risk of Typical Diabetes-Associated Complications in Different Clusters of Diabetic Patients: Analysis of Nine Risk Factors**

**Michael Leutner 1,†, Nils Haug 2,3,† , Luise Bellach <sup>1</sup> , Elma Dervic 2,3 , Alexander Kautzky <sup>4</sup> , Peter Klimek 2,3 and Alexandra Kautzky-Willer 1,5,\***


**Abstract:** Objectives: Diabetic patients are often diagnosed with several comorbidities. The aim of the present study was to investigate the relationship between different combinations of risk factors and complications in diabetic patients. Research design and methods: We used a longitudinal, population-wide dataset of patients with hospital diagnoses and identified all patients (*n* = 195,575) receiving a diagnosis of diabetes in the observation period from 2003–2014. We defined nine ICD-10 codes as risk factors and 16 ICD-10 codes as complications. Using a computational algorithm, cohort patients were assigned to clusters based on the risk factors they were diagnosed with. The clusters were defined so that the patients assigned to them developed similar complications. Complication risk was quantified in terms of relative risk (RR) compared with healthy control patients. Results: We identified five clusters associated with an increased risk of complications. A combined diagnosis of arterial hypertension (aHTN) and dyslipidemia was shared by all clusters and expressed a baseline of increased risk. Additional diagnosis of (1) smoking, (2) depression, (3) liver disease, or (4) obesity made up the other four clusters and further increased the risk of complications. Cluster 9 (aHTN, dyslipidemia and depression) represented diabetic patients at high risk of angina pectoris "AP" (RR: 7.35, CI: 6.74–8.01), kidney disease (RR: 3.18, CI: 3.04–3.32), polyneuropathy (RR: 4.80, CI: 4.23–5.45), and stroke (RR: 4.32, CI: 3.95–4.71), whereas cluster 10 (aHTN, dyslipidemia and smoking) identified patients with the highest risk of AP (RR: 10.10, CI: 9.28–10.98), atherosclerosis (RR: 4.07, CI: 3.84–4.31), and loss of extremities (RR: 4.21, CI: 1.5–11.84) compared to the controls. Conclusions: A comorbidity of aHTN and dyslipidemia was shown to be associated with diabetic complications across all riskclusters. This effect was amplified by a combination with either depression, smoking, obesity, or non-specific liver disease.

**Keywords:** diabetes mellitus; cluster analyses; risk factors; micro- and macrovascular disease

#### **1. Background**

The prevalence of diabetes mellitus is constantly increasing, making this disease a global health problem. Efficient diabetes management is necessary in order to decrease the risk of both micro- and macrovascular complications such as cardiovascular diseases or diabetic nephropathy [1], especially when patients are also diagnosed with other comorbidities such as hypertension or depression [2,3]. Data from large trials conducted in recent

**Citation:** Leutner, M.; Haug, N.; Bellach, L.; Dervic, E.; Kautzky, A.; Klimek, P.; Kautzky-Willer, A. Risk of Typical Diabetes-Associated Complications in Different Clusters of Diabetic Patients: Analysis of Nine Risk Factors. *J. Pers. Med.* **2021**, *11*, 328. https://doi.org/10.3390/ jpm11050328

Academic Editors: George Patrinos and Rutger A. Middelburg

Received: 18 March 2021 Accepted: 15 April 2021 Published: 22 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

decades such as the Collaborative Atorvastatin Diabetes Study (CARDS) or the Action in Diabetes and Vascular disease: preterAx and diamicroN-MR Controlled Evaluation (ADVANCE) have fortified a more generalized approach to the management of diabetics by taking into account additional factors such as cholesterol levels or blood pressure [4,5], or risk factors that are known to increase cardiovascular related mortality rate in diabetics [6]. Accordingly, it has been shown that comorbidities such as dyslipidemia, depression, arterial hypertension, or fatty liver disease predispose diabetic patients to an increased risk of developing serious diabetes-specific complications later in life [7–10]. Given the wide range of common comorbidities in patients with diabetes that may lead to different trajectories of complications, clustering analyses may be especially useful for treatment stratification and personalization in this collective. Earlier studies that defined clusters of diabetic patients on the basis of parameters related to the etiology of diabetes such as insulin resistance or body mass index (BMI) were able to show that these clusters differentiated the risk of developing various diabetes-specific complications [11–15]. Interestingly, the diabetic cluster defined by insulin resistance as the predominant cause of diabetes showed the greatest risk of progressing to renal complications like chronic kidney disease, macroalbuminuria, or end-stage renal disease within the next decade compared to the other clusters defined by autoimmune diabetes, insulin deficiency diabetes, mild age-related diabetes, and mild obesity-related diabetes [11]. The importance of cluster analyses has also been highlighted by a study of Ryu et al., who developed a screening model including gender-specific characteristics for the estimation of an undiagnosed diabetes mellitus in high-risk patients for development of diabetes mellitus, namely patients with a positive family history of diabetes [16]. Additionally, more recently, Nedyalkova et al. showed that *k*-means clustering to detect clinical variables might be useful to stratify type 2 diabetics into distinct subgroups of risk factors [17].

Nevertheless, many aspects have yet to be explored in the novel field of clusterbased analyses for the risk assessment of diabetes-specific complications. Despite the growing number of studies investigating the effect of single concomitant diseases on diabetic complications such as depression, dyslipidemia, fatty liver disease or arterial hypertension [7–10], data on pooled comorbidities, and the respective impact on the disease trajectory are sparse. The aim of the present study was to investigate the relationship between different clusters including different combinations of selected diabetes-associated comorbidities such as overweight/obesity, dyslipidemia or arterial hypertension and the occurrence of typical micro- and macrovascular complications.

#### **2. Methods**

#### *2.1. Data*

We used a longitudinal, population-wide dataset on hospital stays in Austria spanning the period from 1997 to 2014. For each hospital stay, the dataset records a pseudonym for the patient treated, the date of admission and release, and the main and secondary diagnoses associated with the stay, in terms of level-3 ICD-10 codes. The age and sex of each patient are also available. Since each patient is assigned a consistent pseudonym, it is possible to track the development of the state of a patient's health throughout the years.

#### *2.2. Cohort*

The cohort consists of the 195,575 patients who did not receive a hospital diagnosis with ICD-10 code A00–N99 from 1997–2002 and were diagnosed with diabetes mellitus (ICD-10 E10–E14) in the observational period from 2003–2014. Moreover, each patient was at least 60 years old in 2014. The requirement that patients did not receive hospital diagnoses throughout the six years from 1997–2002 justifies our assumption that patients were in good health at the beginning of the observation period. The ratio of males in the cohort was 50,3%. The mean age of cohort patients in 2014 was 74 y (67–82) and 82 y (73–89) for males and females, respectively, the values in brackets denoting lower and upper quartiles. The descriptive characteristics of the study cohort are summarized in Table 1.

**Table 1.** Descriptive characteristics of the study cohort.


#### *2.3. Risk Factors and Complications*

We defined nine typical diabetes specific risk factors, and 16 micro and macrovascular complications in terms of level-3 ICD-10 codes. For example, E66 (overweight and obesity) was defined as a risk factor, whereas I64 (stroke) was defined as a complication. The full list of risk factors and complications is given by the labels on the horizontal axis in Figure 1.

**Figure 1.** Relationship of different clusters with typical diabetes-specific complications.

#### *2.4. Identification of Clusters*

The state of each patient's health is characterized by the set of risk factors and complications which they were diagnosed with during the observation period. Our approach is thus cross-sectional. By using a computational algorithm, we assigned each patient to one of 11 different clusters, each containing more than 5000 individuals. The clusters were constructed algorithmically so that patients in each cluster were diagnosed with combinations of complications during the observation period that were more similar to each other compared to the complications of patients in other clusters. We employed a divisive clustering algorithm that iteratively generates a binary clustering tree, where each node of the tree represents a subset of all patients of the cohort and the leaves of the tree represent the clusters. The result of the algorithm can be read like a decision tree, meaning that each cluster is defined by a set of one or more risk factors diagnosed to each patient in that cluster (inclusion criteria) during the observation period, and a set of risk factors each patient in that cluster was not diagnosed with (exclusion criteria) during that period. It is this property that distinguishes the algorithm from methods as *k*-means [18] or Ward's

method [19], and leads to highly interpretable results. By construction, there exists exactly one cluster without inclusion criteria (the 'healthy cluster' or cluster 0). The algorithm aims to maximize the mutual difference of the obtained clusters in terms of the risk of included patients for developing the different complications. A precise definition of the employed algorithm is given in the Supplementary Materials.

We refer the reader to [20,21] for two previous applications of a modified version of the same clustering approach, demonstrating its usefulness in describing health trajectories and its ability to distinguish patient clusters based on diagnoses. It was shown in [22] that in many cases, the underlying clustering algorithm performed similarly well as *k*-means, while offering greater interpretability of its results.

#### *2.5. Cluster-Based Relative Risk Analysis*

For each cluster *c* and each complication *d*, we computed the relative risk (RR) of a patient in that cluster for receiving a diagnosis of *d* during the observational period compared to patients in the 'healthy cluster'. This allows us to quantify the extent to which the combination of inclusion criteria of cluster *c* increases the risk of complication *d.* We refer to this fraction as the *enrichment* of complication *d* in cluster *c*. While the clustering procedure is performed with all patients, the cluster-based relative risk analysis is also performed separately for males and females.

A detailed description of the clustering algorithm is given in the Supplementary Materials.

#### **3. Results**

Figure 1 visualizes the results of the analysis. On the left of the figure, for each cluster *c* and each risk factor *r*, the probability of a patient in that cluster being diagnosed with risk factor *r* is color-coded from white (zero) via gray to black (one). Black signifies that *r* is an inclusion criterion for cluster *c*, meaning that all patients in that cluster have the corresponding diagnosis. For example, depression (F32) and hypertension (I10) are inclusion criteria for cluster 5. The gray color in the boxes corresponding to obesity (E66), recurrent depressive disorder (F33), and unspecific liver disease (K76) also indicate increased probability of these diagnoses compared to the baseline of cluster 0 (the 'healthy' cluster). Table 2 gives the demographic characteristics of the different clusters. The clusters with the most unbalanced sex distribution are cluster 10 with more males, and clusters 5 and 9 with more females. Cluster 10 has nicotine dependence as an inclusion criterion, whereas clusters 5 and 9 share major depression as inclusion criterion.

On the right of Figure 1, for each cluster *c* and each complication *d*, the enrichment of *d* in cluster *c* is color-coded, blue shading indicating an enrichment of less than one and red shading denoting an enrichment of more than one. Gray color stands for insignificant results. Significance was tested using a *p*-value of 0.05. Adjustment for multiple hypothesis testing via a Bonferroni correction was found to only bring minor changes to the results and therefore omitted. For example, in cluster 6, the red shading in the box corresponding to angina pectoris (I20) indicates that the risk of patients in cluster 6 of being diagnosed with this condition increases by a factor of around 7 compared to the baseline of the healthy cluster 0. On the other hand, in cluster 5 (aHTN and depression), the light blue shading in the box corresponding to retinal disorders (H36) means that patients in cluster 5 had a slightly lower risk of being diagnosed with H36 compared to the healthy baseline cluster 0.

According to Figure 1, six risk factors were defined as inclusion criteria in at least one cluster, namely obesity, dyslipidemia, smoking, depression, arterial hypertension, and non-specific liver disease. A diagnosis of one of these risk factors significantly increased the risk of complications for cohort patients. The risk of cardiovascular disease was particularly increased in those clusters where both dyslipidemia and arterial hypertension were inclusion criteria (clusters 6–10). It was highest in cluster 10 (RR for AP: 10.10 CI: 9.28–10.98), where patients were also diagnosed with nicotine dependence. Patients in that cluster were also at increased risk of loss of extremities (RR: 4.21 CI: 1.5–11.84). The risk of stroke (RR: 4.32 CI: 3.95–4.71) and polyneuropathy (RR: 4.80 CI: 4.23–5.45) was highest in cluster 9, where patients were diagnosed with dyslipidemia, arterial hypertension, and depression, and also had a high prevalence of obesity. Cluster 6 (aHTN and dyslipidemia), cluster 7 (aHTN, dyslipidemia, and overweight/obesity), cluster 8 (aHTN, dyslipidemia and non-specific liver disease) and cluster 9 (aHTN, dyslipidemia and depression) represented cohorts with exaggerated risk of being diagnosed with cardiovascular disease (CVD) and kidney disease.

**Table 2.** Numbers of individuals included in the different clusters. The reported ages refer to the end of the observation period (2014). The values in brackets give lower and upper quartiles, respectively.


#### *Sex-Specific Relationships*

As shown in Figure 2, females in cluster 6–10 were characterized by an exaggerated risk of being diagnosed with angina pectoris. In addition, the loss of extremities and angina pectoris (RR: 10.82 CI: 9.31–12.57) was extremely high for females in cluster 10 (aHTN, dyslipidemia and nicotine dependence). Similar results regarding angina pectoris could be observed for males in clusters 6–10 (Figure 3). Males in cluster 9 (aHTN, dyslipidemia and depression) additionally showed an exaggerated risk of CVD (RR for AP: 7.46 CI: 6.61–8.43), stroke (RR: 5.85 CI: 5.11–6.69), kidney disease (RR for chronic kidney disease: 3.28 CI: 3.06–3.53), and loss of extremities.

**Figure 2.** Relationship of different clusters with typical diabetes-specific complications in females.

**Figure 3.** Relationship of different clusters with typical diabetes-specific complications in males.

#### **4. Discussion**

Exploiting a cross-sectional clustering approach, in this study, we presented ten different risk phenotypes of diabetes mellitus (given our cohort age, these patients will be diagnosed predominately by diabetes mellitus type 2) based on commonly associated risk factors that allow for risk stratification for 16 typical diabetes complications. The combination of arterial hypertension and dyslipidemia was related to an increased risk of being diagnosed with diabetic complications. The risk of being diagnosed with microand macrovascular complications was particularly increased if an additional diagnosis of depression, nicotine dependence, obesity, or non-specific liver disease was present. In females, especially clusters 9 (aHTN, dyslipidemia and depression) and 10 (aHTN, dyslipidemia, and nicotine dependence) and in males, cluster 9 represents patients with a severely increased risk of being diagnosed with multiple diabetes complications.

In recent decades, evidence from large trials has greatly refined the therapeutic approach for diabetic patients. Data from the UK Prospective Diabetes Study (UKPDS) highlighted the central role of glycemic control in the prevention of diabetes-associated complications, showing, among other things, a risk reduction of 25% for microvascular disease upon intensified glucose-lowering therapy [23] and a 32% reduction for any diabetes-related endpoint for overweight diabetics receiving metformin [24]. Furthermore, additional factors such as cholesterol levels or blood pressure have been shown to influence macrovascular outcomes of patients with diabetes mellitus in the CARDS trial, a randomized controlled trial assessing the effect of an add-on low dose of atorvastatin therapy in type 2 diabetics [5], and the ADVANCE trial, which investigated the effect of a fixed angiotensin converting enzyme (ACE) inhibitor–diuretic combination in type 2 diabetics irrespective of baseline blood pressure [4]. Following up on this, more recent studies have identified several risk factors that increase the risk of micro- and macrovascular complications in diabetic patients [8,10,11]. Diagnoses such as arterial hypertension and dyslipidemia, in particular, are related to major macrovascular complications such as coronary artery disease or stroke [25,26] and microvascular complications like diabetic retinopathy or renal disease [27]. However, the risk of developing micro- and macrovascular diseases as a function of different combinations of risk factors is an important topic in type 2 diabetic patients, as diabetic patients often present with more than one comorbidity [28,29]. In the present study, we cross-sectionally identified five clusters of diabetic patients with different

combinations of common comorbidities of diabetes, each of them showing an exaggerated risk of being diagnosed with micro- and macrovascular diabetic complications. These five clusters share the diagnoses of dyslipidemia and arterial hypertension as inclusion criteria, alongside combinations of smoking, depression, obesity, or non-specific liver disease. An increased risk of being diagnosed with micro- and macrovascular complications could be observed in the subgroup of diabetic patients with arterial hypertension and dyslipidemia, who also had a positive history of smoking. In this specific cluster, there was an increased risk of being diagnosed with CVD, atherosclerosis, and amputation of the legs. In a sex-specific analysis, a close relationship between smoking and diagnosis of angina pectoris and amputation of the legs in women could be observed. Smoking has been associated with an increased risk of glycemic progression in diabetics [30] as well as with an increased risk of macrovascular complications, especially peripheral artery disease, coronary artery disease [31], and diabetic foot ulcerations findings [32], which are superimposable with the risk characteristics of subgroup "cluster 10". The underrepresentation of affections of the retina in nicotine-dependent diabetic patients is in accordance with a study by Cai et al., who noted a decreased risk of diabetic retinopathy and proliferative diabetic retinopathy [33]. Nevertheless, our results highlight that the combination of arterial hypertension, dyslipidemia, and smoking defines a specific cohort at a major risk of being diagnosed with diabetic complications. Next to that, cluster 9, characterized by type 2 diabetic patients diagnosed with arterial hypertension, dyslipidemia, and additional depression, also showed a close relationship to macro- and microvascular disease. Hence, men in this cluster were particularly at increased risk, especially of leg amputations. To date, co-occurrence of depression and type 2 diabetes has been identified as a major risk factor for chronic kidney disease, coronary heart disease, stroke, and increased all-cause mortality in a retrospective analysis investigating 933,211 patients. Although patients with depression were younger and had a higher estimated glomerular filtration rate (eGFR) than diabetics without depression in the study by Novak et al., they were characterized by a higher rate of comorbidities [8]. Accordingly, in the present study, we could show that the combination of depression with aHTN and dyslipidemia is especially related to an increased risk of being diagnosed with typical diabetes-specific complications. One possible explanation for the elevated risk of complications in depressed diabetic patients could be deficient glycemic control, as patients with depression are characterized by a worse glycemic profile compared to type 2 diabetic patients without depression [34]. The third high-risk cluster in the present study, cluster 8, comprised type 2 diabetic patients with a predominant co-occurrence of aHTN and dyslipidemia additional to the diagnosis of non-specific liver disease such as non-alcoholic fatty liver disease (NAFLD). NAFLD is in general closely related to metabolic disorders like diabetes mellitus and the respective complications. Evidence points to a crosstalk between inflamed abdominal fatty tissue and the liver, combined with a state of gut microbial dysbiosis accompanied by intestinal barrier dysfunction, leading to systemic inflammation, which consequently constitutes a driving factor in the development of vascular complication in diabetic patients [35]. Cluster 7 is defined by type 2 diabetic patients with arterial hypertension, dyslipidemia, and obesity. To date, several cluster analyses have reported an increased risk of diabetes-specific complications in type 2 diabetic patients suffering from obesity [13,36].

In this study, we clustered cross-sectionally, and hence one limitation is that causal relationships cannot be evaluated. A further limitation of the present study is that we did not have access to 4-point ICD codes that would have allowed analysis of the combinations of the underlying diseases with comorbidities. Nor did we have access to data on laboratory parameters, anthropometric parameters, or prescribed medications. A further limitation is that we defined our cohort with the ICD-codes E10–E14, which also includes diagnosis of diabetes mellitus type 1. However, given our cohort age, the diagnoses of the patients in the present study will be dominated by patients with diabetes mellitus type 2. Additionally, we cannot draw conclusions as to whether adequate diabetes therapy is related to changes in the risk of complications. Another limitation is that we had no access to data on lifestyle

habits such as physical activity and nutritional data, which are closely related to the development of insulin resistance and diabetes mellitus [37].

#### **5. Conclusions**

In conclusion, our results demonstrate that the combination of arterial hypertension and dyslipidemia is especially associated with an elevated risk of diabetes-associated complications. Furthermore, an additional diagnosis of one of the following diseases, namely depression, obesity, non-specific liver disease, or the habit of smoking dramatically aggravates the risk of being diagnosed with micro- or macrovascular diseases in a diabetic patient collective. Therefore, patients who present with the unfavorable combinations of risk factors defined by these five clusters should be screened regularly and an intensified treatment and observation regimen may be implemented in the management algorithm of type 2 diabetes mellitus. In particular, the increasing prevalence of smoking in women in recent years [38] calls for better smoking cessation programs. Additionally, the close relationship between depression and typical diabetic complications already demonstrated by earlier work and shown to be even more pronounced in males in the present study calls for regular screening programs for depression in all diabetic patients, allowing for establishment of appropriate treatment [39]. However, prospective clinical studies would be necessary to further investigate the relationship between different risk factor combinations and the development of diabetes-specific outcomes.

**Supplementary Materials:** Supplementary information is available online at https://www.mdpi. com/article/10.3390/jpm11050328/s1.

**Author Contributions:** Study design: M.L., N.H., P.K. and A.K.-W. Conceived the analytical method: N.H. Implemented the cluster analysis and produced the figures with assistance by E.D. Manuscript writing: M.L. and N.H. Contributed to the introduction and the discussion: L.B., E.D., A.K., P.K. and A.K.-W. All authors read, reviewed, and approved the final manuscript. A.K.-W. is the guarantor of this work. All authors have read and agreed to the published version of the manuscript.

**Funding:** The present study was funded by the Vienna Science and Technology Fund (MA16-045).

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Ethics Committee of the Medical University of Vienna (Nr.:1020/2020).

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data were provided by the Federal Ministry for Health of Austria.

**Conflicts of Interest:** There are no conflict of interest associated with the manuscript and the authors have nothing to disclose.

#### **References**


## *Article* **Glycosylation of IgG Associates with Hypertension and Type 2 Diabetes Mellitus Comorbidity in the Chinese Muslim Ethnic Minorities and the Han Chinese**

**Xiaoni Meng <sup>1</sup> , Manshu Song 1,2,\* , Marija Vilaj <sup>3</sup> , Jerko Štambuk <sup>3</sup> , Mamatyusupu Dolikun <sup>4</sup> , Jie Zhang <sup>1</sup> , Di Liu <sup>1</sup> , Hao Wang <sup>5</sup> , Xiaoyu Zhang <sup>1</sup> , Jinxia Zhang <sup>1</sup> , Weijie Cao <sup>1</sup> , Ana Momˇcilovi´c <sup>3</sup> , Irena Trbojevi´c-Akmaˇci´c <sup>3</sup> , Xingang Li 2,6 , Deqiang Zheng <sup>1</sup> , Lijuan Wu <sup>1</sup> , Xiuhua Guo <sup>1</sup> , Youxin Wang 1,2 , Gordan Lauc 3,7 and Wei Wang 1,2,6**


**Abstract:** Objectives: Hypertension and type 2 diabetes mellitus comorbidity (HDC) is common, which confers a higher risk of cardiovascular disease than the presence of either condition alone. Describing the underlying glycomic changes of immunoglobulin G (IgG) that predispose individuals to HDC may help develop novel protective immune-targeted and anti-inflammatory therapies. Therefore, we investigated glycosylation changes of IgG associated with HDC. Methods: The IgG N-glycan profiles of 883 plasma samples from the three northwestern Chinese Muslim ethnic minorities and the Han Chinese were analyzed by ultra-performance liquid chromatography instrument. Results: We found that 12 and six IgG N-glycan traits showed significant associations with HDC in the Chinese Muslim ethnic minorities and the Han Chinese, respectively, after adjustment for potential confounders and false discovery rate. Adding the IgG N-glycan traits to the baseline models, the area under the receiver operating characteristic curves (AUCs) of the combined models differentiating HDC from hypertension (HTN), type 2 diabetes mellitus (T2DM), and healthy individuals were 0.717, 0.747, and 0.786 in the pooled samples of Chinese Muslim ethnic minorities, and 0.828, 0.689, and 0.901 in the Han Chinese, respectively, showing improved discriminating performance than both the baseline models and the glycan-based models. Conclusion: Altered IgG N-glycan profiles were shown to associate with HDC, suggesting the involvement of inflammatory processes of IgG glycosylation. The alterations of IgG N-glycome, illustrated here for the first time in HDC, demonstrate a biomarker potential, which may shed light on future studies investigating their potential for monitoring or preventing the progression from HTN or T2DM towards HDC.

**Keywords:** hypertension and type 2 diabetes mellitus comorbidity; IgG; N-glycosylation; biomarkers

**Citation:** Meng, X.; Song, M.; Vilaj, M.; Štambuk, J.; Dolikun, M.; Zhang, J.; Liu, D.; Wang, H.; Zhang, X.; Zhang, J.; et al. Glycosylation of IgG Associates with Hypertension and Type 2 Diabetes Mellitus Comorbidity in the Chinese Muslim Ethnic Minorities and the Han Chinese. *J. Pers. Med.* **2021**, *11*, 614. https:// doi.org/10.3390/jpm11070614

Academic Editor: Rutger A. Middelburg

Received: 7 May 2021 Accepted: 21 June 2021 Published: 29 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Hypertension (HTN) is one of the leading modifiable risk factors for morbidity and mortality of cardiovascular disease (CVD) worldwide [1]. Likewise, type 2 diabetes mellitus (T2DM) has become one of the global disease burdens [2]. HTN and T2DM are closely interlinked due to similar risk factors including obesity, dyslipidemia, insulin resistance, endothelial dysfunction, vascular inflammation, atherosclerosis, and sequential inappropriate activation of renin-angiotensin system, etc. [3,4]. Previous studies have shown that hypertension is at least twice as frequent in patients with diabetes compared with those who do not have diabetes [5]. Accordingly, a considerable portion of the world population undergoes coexistent T2DM and HTN, which increases the risk of CVD and renal failure in individuals [6]. Compared with the normotensive non-diabetic individuals, patients with HTN and T2DM comorbidity (HDC) have a four-fold increased risk of CVD [7]. The number of patients diagnosed with HTN and T2DM comorbidity (HDC) is on the rise, resulting in increased medical expenses and resource utilization [8]. Thus, further studies of biomarkers for monitoring and preventing the progression from HTN or T2DM towards HDC are timely and called for.

The glycan attachment to individual proteins is driven by specific enzymes and plays an indispensable role in almost all physiological processes [9]. In general, N-glycome constructed in healthy individuals is fairly stable but can vary due to pathological conditions, especially during inflammation [10]. Altered glycosylation pathways can have significant effects on protein structure and function that occur during the transition from healthy to diseased states. Due to the absence of the genetic template and the lack of stringency within the biosynthetic pathway, the glycan structures of individuals are influenced by the intricate interplay of genetic polymorphisms and environmental factors [11]. Therefore, it is warranted to elucidate the role of altered glycome composition in the pathogenic immune response and to analyze the possibility of using N-glycosylation as a latent predictor for the development of multifactorial diseases, such as HDC.

As the most abundant antibody in human blood, immunoglobulin G (IgG) regulates systemic inflammation balance at multiple levels, and its function is largely influenced by glycosylation. It has been shown that the pro- and anti-inflammatory activities of IgG could no longer be demonstrated if N-glycans are completely removed [12]. Aberrant IgG N-glycosylation has also been demonstrated to be involved in the pathogenesis such as aging [13], HTN [14], and T2DM [15,16], which are relevant risk factors contributing to HDC. These observations make us speculate that changes of IgG N-glycosylation may be related to the pathogenesis of HDC by regulating inflammation. Understanding the underlying altered IgG N-glycosylation profiles that predispose individuals to HDC may help to develop novel protective immune-targeted and anti-inflammatory therapies and/or prevention approaches. However, the associations between IgG N-glycans and HDC have not been described in current literature.

There are 55 officially recognized ethnic minority groups within China in addition to Han majority (which constitute about 91.9% of the total population). The population size of Han Chinese and Chinese ethnic minorities are 1286.3 and 125.5 million, respectively (according to the seventh national census statistics of China conducted in 2020). The genetic background, culture, socioeconomic conditions, climate, geographical features, lifestyle, and dietary habit are substantially different between Han Chinese and other ethnic minorities. The Kirgiz (KGZ), Tajik (TJK), and Uygur (UIG) are three Chinese Muslim ethnic minorities, with a total population size of 11.94 million, inhabiting the Xinjiang Uygur Autonomous Region, the northwestern frontier of China, for generations, without intermarry with other ethnic groups [17]. These Chinese Muslim minorities share relatively similar lifestyles, diet-related habits, living environment, culture, and genetic backgrounds, thereby are well-defined isolated populations with substantial ancestry contribution from Western Eurasian people [18,19]. The prevalence of HTN and T2DM in the northwestern Chinese minority populations is on the rise and relatively higher than that in the Han Chinese, the majority ethnic group in China [20,21]. Investigating the IgG N-glycan profiles

of ethnic groups with different prevalence of HTN and T2DM could shed light on the potential population-specific and/or general IgG N-glycan patterns and their relevance to the pathogenesis of HDC.

In this study, we investigated the changes in IgG N-glycosylation profiles associated with HDC and their biomarker potential for distinguishing HDC from HTN, T2DM, and healthy individuals in the aforementioned populations. This work could further increase our understanding of HDC pathophysiology and may shed light on future studies investigating their potential for monitoring or preventing the progression from HTN or T2DM towards HDC.

#### **2. Materials and Methods**

#### *2.1. Subjects*

The participants from the three northwestern Chinese ethnic minorities (KGZ, TJK, and UIG) were recruited, and their demographic characteristics (gender, age, ethnicity, levels of education, occupation, and urban/rural residence) and anthropometric measurements (height, weight, body mass index [BMI], and blood pressure), as well as blood samples, were collected. All the Muslim participants were recruited from the isolated underdeveloped rural areas (Halajun town, Tagan town, and Minfeng town) of Xinjiang, China, where the respective Muslim minority "exclusively" inhabit for generations without co-resident Han or other ethnic groups, with either primary or secondary education. The occupations of all the UIG participants were farmers and of all the KGZ and TJK participants were herdsmen. A total of 484 samples from 165 KGZ, 146 TJK, and 173 UIG were included in this study, composed of 67 HDC, 183 HTN, 51 T2DM, and 183 healthy individuals (median ages for the four health condition groups: 63, 59, 54, and 51 years). The participants met the following inclusion criteria: (1) Over 18 years old; (2) self-reported the KGZ, TJK, or UIG ethnicity without intermarriage history with other ethnic groups within at least the past three generations. Individuals were excluded based on the following criteria: (1) Pregnant or lactating women; (2) history of severe hepatic, nephritic and autoimmune disease; and (3) diagnosis of diseases of mental illness or infectious disease.

In this study, a total of 399 Han Chinese participants from a community-based survey in Beijing were also included, composed of 72 HDC, 112 HTN, 50 T2DM, and 165 healthy controls (median ages for the four health condition groups: 66, 49, 54, and 46 years). The demographic characteristics, anthropometric measurements, and blood samples of Han Chinese population were collected in a cross-sectional study conducted in 2012, as described previously [13].

All 883 participants included in this study have signed written informed consent. The study was approved by the local community leaders and the Ethics Committee of Capital Medical University, Beijing, China (No. 2015SY21 and 2009SY16) and adhered to the principles of the Declaration of Helsinki.

#### *2.2. Diagnostic Criteria*

HTN was defined as an average systolic blood pressure (SBP) ≥ 140 mmHg and/or an average diastolic blood pressure (DBP) ≥ 90 mmHg, according to the Chinese guidelines for the prevention and treatment of HTN in adults [22]. Peripheral blood samples of the participants for analyses of biochemical indicators (fasting plasma glucose and blood lipids) and IgG N-glycan profiles were collected using EDTA anti-coagulated tube after an overnight fast. The diagnosis of T2DM was based on fasting blood glucose ≥ 7.0 mmol/L in this study, according to World Health Organization (WHO) Criteria [23]. HDC was defined as the co-existence of HTN and T2DM [24]. Dyslipidemia was defined as total cholesterol (TC) ≥ 6.2 mmol/L, triglycerides (TG) ≥ 2.3 mmol/L, low-density lipoprotein cholesterol (LDL) ≥ 4.1 mmol/L, or high-density lipoprotein cholesterol (HDL) < 1.0 mmol/L, according to the Guidelines for Prevention and Treatment of Dyslipidemia in Adults in China [25].

#### *2.3. Analysis of IgG N-glycosylation*

IgG was isolated from plasma samples from all participants using 96-well protein G monolithic plates, and N-linked glycans were released and labeled as described previously [26]. In brief, the isolated IgG was eluted from protein G with 1 mL of 0.1 M formic acid and instantly neutralized with 1 M ammonium bicarbonate. Then, IgG N-glycans were released using PNGase F, labeled with 2-aminobenzamide and purified with 96% acetonitrile and ultrapure water to remove proteins, salts and excess of reagents from deglycosylation and fluorescent labeling reaction. In the end, IgG N-glycans were analyzed by hydrophilic interaction liquid chromatography (HILIC) using ultra-performance liquid chromatography (UPLC) instrument as described previously [27].

In total, 24 chromatographic peaks, (i.e., IgG glycan peaks (GPs, GP1–GP24)) were measured and the amount of glycans in each peak was expressed as percentage of total integrated area [15], and the detailed structures of the N-glycan peaks were shown in Table S1. The 54 derived glycan traits (DGs, DG1–DG54) were calculated from the directly measured GPs [28]. These DGs were used to describe the relative abundance of galactosylation, sialylation, bisecting *N*-acetylglucosamine (GlcNAc), and core fucosylation. To remove the experimental variation from measurements, normalization and batch correction were conducted on the UPLC measured glycan data [29]. Therefore, the measurements were comparable across all samples.

#### *2.4. Statistical Analysis*

Based on our previous publications, in which the correlation coefficients of the relationship between the sialylated glycans (i.e., IgG 4 G1S (that is, Galactose and Sialic acid) and monosialylated glycan, respectively) and SBP in the Chinese Muslim ethnic minority (Kazakh) and the Han Chinese were of a similar magnitude: −0.373 and −0.372, respectively [30,31]. With 75% power and a 0.05 two-sided significance level, the required sample size is estimated to be 48 participants for each of the four health condition groups (i.e., HDC, HTN, T2DM, and healthy individuals) using StatsToDo (http: //www.statstodo.com/SSizCorr\_Pgm.php; accessed on 9 October 2020). The data from the three Chinese Muslim ethnic minorities with similar genetic and environmental backgrounds were pooled together for the statistical analysis, given that the sample size of the respective four health condition groups of the individual Chinese Muslim ethnic minority is insufficient.

Normality distributions of all glycan traits were double-checked by using the Kolmogorov– Smirnov test and Shapiro–Wilk test. When the continuous variable conformed to the normality distribution, means together with standard deviations were used in descriptive statistics. Otherwise, medians together with the interquartile ranges were used. The Chi-square test was applied to investigate the differences in categorical variables. Since the majority of the IgG N-glycan traits were not normally distributed, the Kruskal–Wallis test was conducted to determine the differences of directly measured and derived IgG N-glycan traits among the HDC, HTN, T2DM, and healthy controls.

The IgG glycosylation features of HDC may combine the IgG glycosylation features of both HTN and T2DM, thus a Venn diagram was constructed to illustrate shared IgG N-glycan features between the four health condition groups using online software (http://bioinformatics.psb.ugent.be/webtools/Venn/; accessed on 17 April 2021). The associations between HDC and IgG N-glycan traits were assessed through logistic regression analysis, with the IgG N-glycan traits as independent variables, adjusted for potential confounders (i.e., age, gender, ethnicity, BMI, and dyslipidemia). Given the potential internal associations among IgG N-glycan traits, which could induce multicollinearity in the statistical models, the least absolute shrinkage and selection operator (LASSO) method was conducted to reduce dimension before constructing the logistic regression models.

To evaluate differences in IgG N-glycan traits of different health conditions in the three Chinese Muslim minorities pooled samples and Han Chinese samples, the IgG N-glycan traits screened by the LASSO method were included to construct logistic regression models

for HDC versus HTN, HDC versus T2DM, and HDC versus healthy controls, respectively. Three types of logistic classification models were applied: (1) The baseline models included age, gender, ethnicity, BMI, and dyslipidemia as the covariates; (2) the glycan-based models included the screened IgG N-glycan traits as the covariates; and (3) the combined models incorporated the screened IgG N-glycan traits on top of the covariates included in the baseline models. The prediction error of the five-fold cross-validation procedure was used to assess the performance of the discriminant models (R package "*boot*") [32]. The discriminant ability of the models in delineating the HDC from the HTN, T2DM, and healthy individuals was evaluated by receiver operating characteristic (ROC) curve analysis (R package "*pROC*").

Data were analyzed and visualized using SAS (Version 9.4) and R packages (Version 4.0.4). For the logistic regression analysis, the Benjamini–Hochberg procedure was used to control the false discovery rate (FDR) [33]. All analyses were two-sided, and statistical significance was set at *p* < 0.05. Bonferroni adjustment was used for multiple comparisons among the four health condition groups to control familywise error rate (FWER) in a very stringent criterion and to compute the adjusted *p* values by directly adjust the significance level as *α* 0 = *α*/*m* (i.e., the number of simultaneously tested hypotheses). Thus, *p* < 0.05/6 (0.0083) was considered to be statistical significance in the comparative analysis of demographic and biochemical characteristics, and the levels of IgG N-glycans between the groups.

#### **3. Results**

#### *3.1. Demographic and Biochemical Characteristics*

A total of 883 participants (484 Chinese Muslim minorities and 399 Han Chinese) were analyzed in this study. Demographic and biochemical characteristics of the HDC, HTN, T2DM, and healthy individuals for the Chinese Muslim ethnic minorities and those for the Han Chinese are shown separately (Table S2) and pooled together in Table 1. In the pooling of the Chinese Muslim ethnic minorities, the age significantly differed in HDC compared to T2DM and healthy controls, and the SBP and DBP were significantly higher in the HDC group than those in both T2DM and healthy individuals, and the FBG was higher in the HDC than that in both the HTN and healthy individuals. In the Han Chinese individuals, the age, gender, BMI, SBP, DBP, TC, TG, LDL, FBG, and dyslipidemia significantly differed among the four health condition groups (i.e., HDC, HTN, T2DM, and healthy controls).

*J. Pers. Med.* **2021**, *11*, 614



Kruskal–Wallis test was used. χ2 test: Gender, dyslipidemia; BMI, body mass index; DBP, diastolic blood pressure; FBG, fasting blood glucose; HDC, hypertension and type 2 diabetes mellitus comorbidity; HDL, high-density lipoprotein cholesterol; HTN, hypertension; LDL, low-density lipoprotein cholesterol; SBP, systolic blood pressure; T2DM, type 2 diabetes mellitus; TC, total cholesterol; TG, total triglycerides; P25, the 25th percentile; P75, the 75th percentile; \* *p* < 0.05 was considered statistical significance; & *p* < 0.0083 was considered statistical significance between HDC and HTN; \$ *p* < 0.0083 was considered statistical significance between HDC and T2DM; # *p* < 0.0083 was considered statistical significance between HDC and Controls.

#### *3.2. The Association of IgG N-Glycans with HDC*

IgG N-glycome composition (i.e., GPs and DGs) was analyzed in all samples. Out of the total 78 glycan traits, 56 and 62 glycan traits were found to be abnormally distributed (*p*Kolmogorov–Smirnov < 0.05 and *p*Shapiro–Wilk < 0.05) in the pooled samples from the three Chinese Muslim ethnic minorities and those from Han Chinese, respectively (Table S3). The differences of directly measured and derived glycan compositions among HDC, HTN, T2DM, and healthy controls for the Chinese Muslim ethnic minorities and the Han Chinese, respectively, were shown in Figure 1 and Table S3. We found the increased relative abundance of GP4 (FA2), GP6 (FA2B), and GP11 (FA2[3]BG1) and the reduced relative abundance of GP14 (FA2G2), GP18 (FA2G2S1), and GP20 (FA2FG2S1) in the HDC group compared with the other three health condition groups (i.e., HTN, T2DM, and healthy individuals), with similar trends in the Chinese Muslim ethnic minorities and the Han Chinese (Figure 1).

Table S4 and its accompanying Venn Diagrams (Figure 2) displays and illustrates respectively, the results of the logistic regression analysis of IgG N-glycans after adjustment for potential confounders (i.e., age, gender, ethnicity, BMI, and dyslipidemia), and shared and unique IgG glycan features between the different health conditions. Noteworthily, 12 (i.e., 7 in HDC versus HTN, 3 in HDC versus T2DM, and 2 in HDC versus Controls) and 6 (i.e., 3 in HDC versus HTN, 2 in HDC versus T2DM, and 1 in HDC versus Controls) IgG N-glycan traits were found to be significantly different between four health condition groups within the three Chinese Muslim minorities pooled samples and Han Chinese samples, respectively, after controlling the aforementioned confounders and FDR, in the HDC pairwise comparison with the three health condition groups (i.e., HDC versus HTN, HDC versus T2DM, and HDC versus healthy individuals) (Table S4 and Figure 3). To elaborate further, it was found that within the pooled samples from the three Chinese Muslim minorities, seven glycan traits (increased relative abundance of GP5 [M5], GP6 [FA2B], and DG54 [BG2n/(FG2<sup>n</sup> + FBG2<sup>n</sup> )] as well as reduced relative abundance of GP16 [FA2G1S1], GP18 [FA2G2S1], DG3 [FGS/(F + FG + FGS)], and DG24 [GP8<sup>n</sup> ]) were significant differences between HDC and HTN (*p* < 0.05). Three glycan traits (increased relative abundance of GP5 and reduced relative abundance of GP16 and GP18) were significant differences between HDC and T2DM (*p* < 0.05). In addition, two glycan traits (GP5 and GP6), with increased relative abundance, were significantly different between HDC and healthy controls (*p* < 0.05).

For the pooling of three Chinese Muslim ethnic minorities, we found that one shared glycan trait—GP5, containing high mannose glycan (M5) in the HDC, was significantly higher than those in the other three healthy condition groups (i.e., HTN, T2DM, and healthy individuals). Additionally, two shared glycan traits (GP16 and GP18), containing corefucosylated galactosylated glycans with sialic acid (FA2G1S1 and FA2G2S1, respectively), were observed in lower levels in the HDC compared with HTN and T2DM, and a corefucosylated agalactosylated glycan containing bisecting GlcNAc eluting in GP6 (FA2B) was observed in higher levels in the HDC than those in the HTN and healthy individuals (Figure 3 and Table S4).

Similarly, for the Han Chinese, in comparisons between HDC and the other three healthy condition groups (i.e., HTN, T2DM, and healthy controls), one shared glycan trait (GP20 [FA2FG2S1]) containing core-fucosylated digalactosylated glycan with sialic acid of reduced relative abundance was of statistical significance. The relative abundance of DG2 (FBGS/(FBG + FBGS)) and DG53 (FG2n/(BG2<sup>n</sup> + FBG2<sup>n</sup> )) in the HDC were significantly higher than those in the HTN. Additionally, the increased relative abundance of GP24 (FA2BG2S2) containing sialylation of fucosylated galactosylated structures with bisecting GlcNAc was found in the HDC compared with the T2DM (Figure 3 and Table S4).

Chinese (Figure 1).

± 1.5 the IQR.

IgG N*-*glycome composition (i.e., GPs and DGs) was analyzed in all samples. Out of the total 78 glycan traits, 56 and 62 glycan traits were found to be abnormally distributed (*p*Kolmogorov–Smirnov < 0.05 and *p*Shapiro–Wilk < 0.05) in the pooled samples from the three Chinese Muslim ethnic minorities and those from Han Chinese, respectively (Table S3). The differences of directly measured and derived glycan compositions among HDC, HTN, T2DM, and healthy controls for the Chinese Muslim ethnic minorities and the Han Chinese, respectively, were shown in Figure 1 and Table S3. We found the increased relative abundance of GP4 (FA2), GP6 (FA2B), and GP11 (FA2[3]BG1) and the reduced relative abundance of GP14 (FA2G2), GP18 (FA2G2S1), and GP20 (FA2FG2S1) in the HDC group compared with the other three health condition groups (i.e., HTN, T2DM, and healthy individuals), with similar trends in the Chinese Muslim ethnic minorities and the Han

*3.2. The Association of IgG N-Glycans with HDC* 

**Figure 1.** Differences in IgG N-glycan traits between HDC, HTN, T2DM patients, and healthy controls for the pooled samples from three Chinese Muslim ethnic minorities (**A**) and the Han Chinese samples (**B**). Each box represents the 25th to 75th percentiles (interquartile range [IQR]). Lines inside the boxes represent the medians. The whiskers represent the lowest and highest values within the boxes **Figure 1.** Differences in IgG N-glycan traits between HDC, HTN, T2DM patients, and healthy controls for the pooled samples from three Chinese Muslim ethnic minorities (**A**) and the Han Chinese samples (**B**). Each box represents the 25th to 75th percentiles (interquartile range [IQR]). Lines inside the boxes represent the medians. The whiskers represent the lowest and highest values within the boxes ±1.5 the IQR. GP, glycan peak; HDC, hypertension and type 2 diabetes mellitus comorbidity; HTN, hypertension; T2DM, type 2 diabetes mellitus. significant differences between HDC and HTN (*p* < 0.05). Three glycan traits (increased relative abundance of GP5 and reduced relative abundance of GP16 and GP18) were significant differences between HDC and T2DM (*p* < 0.05). In addition, two glycan traits (GP5 and GP6), with increased relative abundance, were significantly different between HDC and healthy controls (*p* < 0.05).

**Figure 2.** Venn diagram of shared and unique IgG N-glycan features between the different health condition groups for the pooled samples from three northwestern Chinese Muslim ethnic minorities and those from the Han Chinese. The numbers of shared and unique IgG N-glycan traits in the three pairwise comparison groups are shown based on the result of logistic regression analysis adjusted for the covariates (for the pooling of three Chinese Muslim minorities, logistic regression adjusted for age, gender, ethnicity, BMI, and dyslipidemia [A]; For the Han Chinese, logistic regression adjusted for age, gender, **Figure 2.** Venn diagram of shared and unique IgG N-glycan features between the different health condition groups for the pooled samples from three northwestern Chinese Muslim ethnic minorities and those from the Han Chinese. The numbers of shared and unique IgG N-glycan traits in the three pairwise comparison groups are shown based on the result of logistic regression analysis adjusted for the covariates (for the pooling of three Chinese Muslim minorities, logistic regression adjusted for age, gender, ethnicity, BMI, and dyslipidemia [A]; For the Han Chinese, logistic regression adjusted for age, gender, BMI, and dyslipidemia [B]). HDC, hypertension and type 2 diabetes mellitus comorbidity; HTN, hypertension; T2DM, type 2 diabetes mellitus.

HDC, hypertension and type 2 diabetes mellitus comorbidity; HTN, hypertension; T2DM, type

BMI, and dyslipidemia [B]).

2 diabetes mellitus.

**Figure 3.** The associations between IgG N-glycan traits and HDC. AORs and 95% Cis for the association of IgG N-glycan traits with HDC versus HTN/T2DM/healthy controls adjusted for the covariates and false discovery rate (for the pooling of three Chinese Muslim minorities, logistic regression adjusted for age, gender, ethnicity, BMI, and dyslipidemia; For the Han Chinese, logistic regression adjusted for age, gender, BMI, and dyslipidemia). FDR was controlled by Benjamini–Hochberg method and *p* < 0.05 was considered statistical **Figure 3.** The associations between IgG N-glycan traits and HDC. AORs and 95% Cis for the association of IgG N-glycan traits with HDC versus HTN/T2DM/healthy controls adjusted for the covariates and false discovery rate (for the pooling of three Chinese Muslim minorities, logistic regression adjusted for age, gender, ethnicity, BMI, and dyslipidemia; For the Han Chinese, logistic regression adjusted for age, gender, BMI, and dyslipidemia). FDR was controlled by Benjamini–Hochberg method and *p* < 0.05 was considered statistical significance. AOR, adjusted odds ratio; BMI, body mass index; CI, confidence interval; DG, derived glycan; FDR, false discovery rate; GP, glycan peak; HDC, hypertension and type 2 diabetes mellitus comorbidity; HTN, hypertension; T2DM, type 2 diabetes mellitus.

AOR, adjusted odds ratio; BMI, body mass index; CI, confidence interval; DG, derived glycan; FDR, false discovery rate; GP, glycan peak; HDC, hypertension and type 2 diabetes mellitus comor-

For the pooling of three Chinese Muslim ethnic minorities, we found that one shared glycan trait—GP5, containing high mannose glycan (M5) in the HDC, was significantly higher than those in the other three healthy condition groups (i.e., HTN, T2DM, and

significance.

bidity; HTN, hypertension; T2DM, type 2 diabetes mellitus.

#### *3.3. The Comparison of IgG N-Glycans between the Chinese Muslim Ethnic Minorities Pooled Samples and the Han Chinese Samples*

We observed the increased relative abundance of sialylated glycans (i.e., GP24, DG2, and DG53) in the Han Chinese with HDC (HDC versus HTN or T2DM), while a reduced trend (i.e., GP16, GP18, and DG3) was presented in the Chinese Muslim ethnic minorities pooled samples with HDC (HDC versus HTN or T2DM). The HDC-associated changes in the sialylated glycans showed different trends between the Chinese Muslim ethnic minorities and the Han Chinese as a consequence of the presence or absence of bisecting GlcNAc in these sialylated glycans in the different ethnic groups.

Although the logistic regression analyses suggested that directly measured and derived IgG N-glycans between the three Chinese Muslim ethnic minorities and the Han Chinese in the three pairwise comparison groups (i.e., HDC versus HTN, HDC versus T2DM, and HDC versus healthy controls) did not share similarities in their results, we found that their respective results all presented low levels of core-fucosylation and galactosylation (i.e., GP16, GP18, and GP20), as well as high levels of bisecting GlcNAc (i.e., GP6, DG54, GP24, and DG53) in the three Chinese Muslim ethnic minorities with HDC and also in the Han Chinese with HDC (Figure 3 and Table S4).

#### *3.4. Development of the Classification Models and Discrimination of HDC*

Three types of logistic classification models (i.e., baseline models, glycan-based models, and combined models) were constructed for distinguishing HDC from HTN, T2DM, and healthy individuals, and were trained and evaluated by the five-fold cross-validation procedure and ROC curve analysis (Figure 4 and Table 2). The detailed information on the covariates included in the logistic classification models was shown in Table 2. Specifically, in the pooling of three Chinese Muslim ethnic minorities, the AUCs of the glycan-based models discriminating HDC from HTN and T2DM were 0.678 and 0.715, respectively, showing limited improved discriminative performance than the baseline models. However, the AUC of the glycan-based model discriminating HDC from healthy controls was lower than that of the baseline model (difference between AUCs, 0.093, *<sup>p</sup>* = 3.97 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ). A similar trend was observed in the Han Chinese for discriminating HDC from healthy controls (difference between AUCs, 0.115, *<sup>p</sup>* = 9.56 <sup>×</sup> <sup>10</sup>−<sup>3</sup> ) (Table S5).

Adding both IgG N-glycan and baseline traits to the combined models, the AUCs of the combined models differentiating HDC from HTN, T2DM, and healthy controls were 0.717, 0.747, and 0.786 in the pooling of three Chinese Muslim ethnic minorities, and 0.828, 0.689, and 0.901 in the Han Chinese, respectively (Table 2). The AUCs showed modest improvement in the combined models compared to the baseline models, with increments of 0.064 (*<sup>p</sup>* = 4.18 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ) for HDC versus HTN, 0.065 (*<sup>p</sup>* = 8.66 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ) for HDC versus T2DM, and 0.015 (*<sup>p</sup>* = 2.45 <sup>×</sup> <sup>10</sup>−<sup>1</sup> ) for HDC versus healthy controls, in the pooling of three Chinese Muslim ethnic minorities; and with increments of 0.064 (*<sup>p</sup>* = 2.22 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ) for HDC vs. HTN, 0.080 (*<sup>p</sup>* = 8.99 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ) for HDC vs. T2DM, and 0.030 (*<sup>p</sup>* = 2.73 <sup>×</sup> <sup>10</sup>−<sup>1</sup> ) for HDC vs. healthy controls, in the Han Chinese (Table S5). The discriminative performances of the combined models were shown to be superior to the baseline models. Although not all the AUCs of the combined models are statistically different from the baseline models, the addition of glycan variables into the baseline models can increase the discriminative power of the models. It suggests that IgG N-glycans could reinforce the discriminative potential of models. Additionally, all but the combined model discriminating HDC from HTN (Chinese Muslim ethnic minorities) showed lower prediction error than the baseline and glycan-based models in five-fold cross-validation (Table 2).

**Figure 4.** ROC curve illustrated the performance of the models in the classification of HDC from HTN, T2DM, and healthy individuals for the pooled samples from three Chinese northwestern Muslim ethnic minorities (**A**) and those from Han Chinese (**B**). When glycan traits are added into the combined models, the AUCs indicated that the IgG glycan profiles provide more information reinforcing the baseline models utilizing gender, age, ethnicity, BMI, and dyslipidemia as the covariates (with difference between 2 AUCs (95% CI) of 0.064 (0.002, 0.126), 0.065 (−0.009, 0.139), and 0.015 (−0.011, 0.042) in distinguishing HDC from HTN, T2DM, and healthy controls, respectively, in the pooling of three Chinese Muslim ethnic minorities, and of 0.064 (0.008, 0.118), 0.080 (−0.013, 0.172), **Figure 4.** ROC curve illustrated the performance of the models in the classification of HDC from HTN, T2DM, and healthy individuals for the pooled samples from three Chinese northwestern Muslim ethnic minorities (**A**) and those from Han Chinese (**B**). When glycan traits are added into the combined models, the AUCs indicated that the IgG glycan profiles provide more information reinforcing the baseline models utilizing gender, age, ethnicity, BMI, and dyslipidemia as the covariates (with difference between 2 AUCs (95% CI) of 0.064 (0.002, 0.126), 0.065 (−0.009, 0.139), and 0.015 (−0.011, 0.042) in distinguishing HDC from HTN, T2DM, and healthy controls, respectively, in the pooling of three Chinese Muslim ethnic minorities, and of 0.064 (0.008, 0.118), 0.080 (−0.013, 0.172), and 0.030 (0.003, 0.055), respectively, in the Han Chinese). AUC, area under the receiver operating characteristic curve; CI, confidence interval; HDC, hypertension and type 2 diabetes mellitus comorbidity; HTN, hypertension; T2DM, type 2 diabetes mellitus; ROC, receiver operating characteristic curve.

and 0.030 (0.003, 0.055), respectively, in the Han Chinese).

type 2 diabetes mellitus; ROC, receiver operating characteristic curve.

AUC, area under the receiver operating characteristic curve; CI, confidence interval; HDC, hypertension and type 2 diabetes mellitus comorbidity; HTN, hypertension; T2DM,

Adding both IgG N-glycan and baseline traits to the combined models, the AUCs of the combined models differentiating HDC from HTN, T2DM, and healthy controls were 0.717, 0.747, and 0.786 in the pooling of three Chinese Muslim ethnic minorities, and 0.828, 0.689, and 0.901 in the Han Chinese, respectively (Table 2). The AUCs showed modest improvement in the combined models compared to the baseline models, with increments of 0.064 (*p* = 4.18 × 10−2) for HDC versus HTN, 0.065 (*p* = 8.66 × 10−2) for HDC versus T2DM, and 0.015 (*p* = 2.45 × 10−1) for HDC versus healthy controls, in the pooling of three Chinese Muslim ethnic minorities; and with increments of 0.064 (*p* = 2.22 × 10−2) for HDC vs. HTN, 0.080 (*p* = 8.99 × 10−2) for HDC vs. T2DM, and 0.030 (*p* = 2.73 × 10−1) for HDC vs. healthy controls, in the Han Chinese (Table S5). The discriminative performances of the combined models were shown to be superior to the baseline models. Although not all the AUCs of



AUC, area under the curve; BMI, body mass index; CI, confidence interval; DG, derived glycan; GP, glycan peak; HDC, hypertension and type 2 diabetes mellitus comorbidity; HTN, hypertension; SD, standard

deviation; T2DM, type 2 diabetes mellitus; \* *p* < 0.05 was considered statistical significance; The prediction error of five-fold cross-validation procedure was presented in the form of mean ± SD.

52

#### **4. Discussion**

To the best of our knowledge, this is the first investigation on the association and discriminative potential of IgG N-glycan traits in HDC versus HTN, T2DM, and healthy individuals in the pooling of three Chinese Muslim ethnic minorities (KGZ, TJK, and UIG), inhabiting the northwestern frontier of China, as well as the Han Chinese population.

HTN and T2DM are common comorbidities that are multifactorial diseases driven by many common pathogenic factors, including modifiable risk factors (obesity, dyslipidemia, insulin resistance, endothelial dysfunction, vascular inflammation, atherosclerosis, and sequential inappropriate activation of renin-angiotensin system, etc.) [3,4], environmental factors, and genetic predisposition [4,34]. IgG contains complex-type biantennary N-glycan structures, which activate or inhibit IgG Fcγ receptors (FcγRs) by altering the binding affinity of IgG, thereby decisively affecting effector functions of IgG [35]. The pro- and anti-inflammatory actions of IgG are modulated by the specific N-glycan residues including galactose, fucose, sialic acid, and bisecting GlcNAc [9,36,37]. Cumulative evidence indicates that the alteration of terminal modification in IgG N-glycan plays a vital role in the anti- or pro-inflammatory process. Previous studies reported that aberrant glycosylation of IgG (low levels of sialylation and core fucosylation, as well as the high level of bisecting GlcNAc) have been demonstrated to be involved in a range of specific health conditions related to HDC, including aging, obesity, HTN, T2DM, dyslipidemia, and metabolic syndrome [13,14,31,32,38–40]. In the present study, the altered IgG N-glycosylation profile was also observed in the HDC patients, further suggesting that IgG N-glycome patterns are associated with the progression of inflammatory disorders, such as HDC, via inflammatory signaling. These findings present a situation where the role of IgG glycans in the pathogenesis of HDC extends beyond just a reflection of changes stemming from the onset of disease, but rather perhaps has some contribution towards influencing pathogenesis itself resulting from the altered inflammatory actions of IgG.

In this study, the higher level of GP5 (M5) and GP6 (FA2B), and the lower levels of GP16 (FA2G1S1) and GP18 (FA2G2S1) were found to be the population-specific IgG N-glycomic changes with statistical significance (*p* < 0.05) for the pooled samples from three Chinese Muslim ethnic minorities with HDC; and the lower level of GP20 (FA2FG2S1) were observed to be unique to the Han Chinese with HDC. In the pooled samples from the three Chinese Muslim ethnic minorities, HDC patients presented increased core-fucosylation with bisecting GlcNAc of IgG (GP6 and DG54), but decreased sialylation (GP16, GP18, and DG3) and monogalactosylated (GP16) and digalactosylated (GP18) core-fucosylated glycans without bisecting GlcNAc of IgG, when compared with HTN patients. In comparison with T2DM patients, HDC patients showed decreased galactosylation, sialylation, and core-fucosylation without bisecting GlcNAc of IgG (GP16 and GP18), but increased mannose structures of IgG (GP5). Additionally, compared to healthy individuals, HDC patients presented increased core-fucosylation with bisecting GlcNAc (GP6) and mannose structures of IgG (GP5). In the Han Chinese, the lower level of core-fucosylated galactosylated glycan containing sialic acid of IgG (GP20) was observed in the HDC compared with HTN/T2DM/healthy individuals, while sialylated glycans (reduced relative abundance of GP20 and increased relative abundance of GP24 [FA2BG2S2], DG2 [FBGS/(FBG + FBGS)], and DG53 [FG2n/(BG2<sup>n</sup> + FBG2<sup>n</sup> )]) presented inconsistent trends in HDC patients. Specifically, the common characteristic of these HDC-associated sialylated glycans is the presence of bisecting N-GlcNAc, and the increase of this structural feature in HDC outweighed the effects of a decrease in sialylation, which is in agreement with the alterations of sialyation observed in the study on systemic lupus erythematosus [29].

Fucosylation plays key roles in the function of IgG and highly fucosylated IgG shows a reduced magnitude of antibody-dependent cellular cytotoxicity (ADCC), which is a result of IgG-subclass IgG1 binding to each of the different Fc receptors [41]. Bisecting GlcNAc is the central branch of N-glycans, which is biosynthesized by the glycosyltransferase GnT-III encoded by the MGAT3 gene [42]. The presence of bisecting GlcNAc in N-glycans could suppress various terminal modifications (fucosylation, sialylation, and others) of N-

glycans [42], affecting the structure of glycans exerting noteworthy effects on inflammatory responses or immune function. IgG carrying glycans with bisecting GlcNAc was reported to have pronouncedly increased affinity to FcγRIIIa leading to an enhancement in ADCC, whereas fucosylation without bisecting GlcNAc of IgG can reduce ADCC [37]. Previous studies showed that such increased fucosylation with bisecting GlcNAc was associated with T2DM [15]. We also observed similarly increased fucosylation in structures with bisecting GlcNAc (GP6), and decreased fucosylation in structures without bisecting GlcNAc (GP16, GP18, and GP20) in the HDC patients, compared with HTN and T2DM patients. The increased fucosylation in structures with bisecting GlcNAc of IgG N-glycans could lead to the pro-inflammatory status of the body [37], thereby contributing to the progression of HTN or T2DM towards HDC. Thus, the present results indicate that the individuals with higher levels of fucosylated structures with bisecting GlcNAc of IgG may increase the risk for HDC, and that pathophysiological states of individuals with HDC are conditioned on a chronic inflammatory process.

The level of galactosylation has an essential impact on the inflammatory activities of IgG; lack of galactosylation was shown to associate with the pro-inflammatory autoimmune responses through activation of the complement cascade [43]. The down-regulation of galactosyltransferase activity in plasma cells and the activation of complement cascades were identified as potential mechanisms for the changes in IgG galactosylation [43]. As a modulator of inflammatory activity, high galactosylation of IgG Fc N-glycan promotes the association of the inhibitory receptor FcγRIIb and the lectin-like receptor dectin-1, which blocks the pro-inflammatory effector functions of C5aR and CXCR2, and enhances the anti-inflammatory activity of antigen-specific IgG1 immune complexes [44]. Previous studies have suggested that the decrease in galactosylation of IgG N-glycans play potential roles in the development and progression of chronic pro-inflammatory diseases, including HTN, T2DM, and dementia [14–16,45], which are in agreement with the alterations of IgG N-glycosylation patterns observed in the HDC patients in the present study. The decreased IgG galactosylation found in HDC presented a similar pattern found in the previous studies and thus our finding further supports the hypothesis that a decrease in galactosylation is a general feature of multiple chronic diseases associated with decreased anti-inflammatory/immunosuppressive potential of circulating IgG.

Besides fucosylation and galactosylation, sialylation mediates partly the inflammatory function of IgG. IgG shows different pro- and anti-inflammatory potentials on the basis of different Fcγ receptors that IgG preferentially binds to, and the affinity of the receptor depends on the composition of the glycans attached to the 297 asparagine [46]. The addition of sialic acid to IgG can decrease its affinity for FcγRIII, thereby reducing the degree of ADCC, converting its function from pro- to anti-inflammatory [47]. In the two comparison groups (HDC versus HTN, HDC versus T2DM), we observed that patients with HDC showed a significant decrease in sialylation (GP16, GP18, and GP20), indicating that sialylation is correlated with HDC through decreasing adaptive pro-inflammatory actions of IgG. However, sialylation with bisecting GlcNAc showed the reduced relative abundance in HDC patients, which was consistent with the suppression of core fucosylation expression due to the presence of bisecting GlcNAc. Bisecting GlcNAc inhibits the biosynthesis and modification of sialic acid of terminal epitopes in N-glycans, which results in the relative abundance of sialylation (GP16, GP18, GP20, GP24, DG2, DG3, and DG53) presenting inconsistent trends in HDC patients. Additionally, we found the decrease of GPs containing galactosylated sialylated glycan structures (GP16, GP18, and GP20) in HDC patients. Since terminal sialic acid residues are covalently attached to galactosylated glycans, the decrease of galactosylated sialylation could be a consequence of a decrease in galactosylation [48], which could provide a plausible explanation for the phenomenon of the co-directional changes of sialylation and galactosylation shown in the HDC patients. It is worthwhile to investigate the potential of the level of sialylation and galactosylation of IgG for monitoring the progression from HTN or T2DM towards HDC. Thus, future prospective cohort studies are called for to explore if there is a stepwise decrease trend in levels of sialylation and

galactosylation of IgG during progression from healthy individuals to HTN or T2DM and then to HDC.

Among the different Chinese ethnic minorities, the similarities and differences in the glycosylation patterns were found in both the HTN and T2DM patients, compared with the HDC patients, suggesting that those individuals who have developed T2DM or HTN do not show consistent alterations of IgG N-glycans. More work is required on understanding alterations of IgG N-glycans in HTN or T2DM individuals to see if their alterations over time show an inclination towards HDC. Such IgG N-glycans profiling will help to contribute towards developing potential glycan biosignatures for monitoring the progression from HTN or T2DM towards HDC.

The AUCs of the combined models, based on the glycan traits and the baseline characteristics, demonstrate higher discriminative performance to help discern between HDC, HTN, T2DM, and healthy individuals than the models based on either the glycan traits or the baseline characteristics alone. IgG N-glycan profiles reflect the integrative effects of genetic, metabolic, and environmental factors, contributing significantly towards their discriminative potential [11]. Although the addition of glycan traits to the baseline models did not substantially improve the discrimination in all the combined models, the glycan based models constructed based on IgG N-glycans alone displayed enough discrimination power and all the combined models showed improved discrimination performance. These results suggest that IgG N-glycan profiles may contribute significantly to much of the combined risk of these characteristics, where dynamic and populationspecific IgG N-glycomic changes can serve as functional effectors of the disease potentially increasing the diagnostic accuracy for complex diseases. The IgG N-glycosylation profiles associated with HDC described in this study may thus help to shed light on future studies investigating their biomarker potential for personalized monitoring and prevention of the progression from HTN or T2DM towards HDC.

The strength of this study is that we evaluated IgG N-glycomes in relation to HDC and explored similarities and differences in the IgG N-glycosylation profiles for the first time, particularly in different Chinese ethnic groups, including the Chinese Muslim minorities. The alterations of IgG N-glycosylation in HDC patients found in this study suggest a theoretical foundation and present practical direction for future studies to further elaborate the association of IgG glycosylation with HDC. These suggestions are, however, built upon the constraints of this study. Firstly, due to the small sample size of HDC in each ethnic minority in the present study, it is unfeasible to analyze the association of IgG N-glycans with HDC in each of the Muslim minorities separately. Consequently, the same health condition groups from the three Chinese Muslim ethnic minorities were pooled together for the analysis. Furthermore, the Han subjects included in the study were not co-resident in the Xinjiang Autonomous Region since the Muslim minorities were recruited from the isolated rural areas in Xinjiang, where the respective Muslim minority "exclusively" inhabit for generations without co-resident Han or other ethnic groups. Additionally, it is difficult to match the Han subjects with the Muslim minorities in terms of their profile of demographic variables due to the relatively small sample size of the respective health condition group of the Muslim minorities. To compensate for the potential impact of unmatched demographic variables on the results, we adjusted the potential confounding factors in the logistic regression analyses. Thirdly, the casecontrol nature of the present study was not able to provide the information about a causal relationship between the altered IgG glycosylation and the progression of HDC. A previous study reported that the polymorphisms at six loci (i.e., *FUT8*, *FUT6*/*FUT3*, *HNF1A*, *MGAT5*, *B3GAT1*, and *SLC9A9*) affect plasma levels of N-glycans [49], which could be used as the instrumental variables for future Mendelian randomization (MR) studies investigating the causal relationship between the altered IgG glycosylation and HDC. Fourthly, the information on medical treatment of the patients with HTN and/or T2DM was not completely collected. Consequently, the potential impact of medical treatment for HTN and/or T2DM on the associations between HDC and IgG glycans was unable to be

addressed in this current study. Lastly, only the five-fold cross-validation of the model is performed, external validation of model performance is still needed in further studies.

In the current study, the IgG N-glycans were analyzed by HILIC-UPLC method, which is a profiling method for total IgG N-glycosylation. The detailed differences between the glycosylation profiles of the individual IgG subclasses associated with HDC remain to be investigated, which would be performed in our coming study. Future prospective studies with larger samples in multiracial populations and/or MR studies examining causal inference are warranted to validate the potential of IgG N-glycosylation biomarkers for personalized monitoring and prevention of HDC. Clinical exploitation of IgG N-glycan biomarkers will be facilitated by the existing broad application in clinical laboratories of UPLC. This study shows the potential clinical application of N-glycans as dynamic biomarkers (e.g., of routine diagnostic and prophylactic test for HDC patients). Future longitudinal and prospective studies would provide a better understanding of the transition from HTN or T2DM to HDC and allow us to generate glycan-based markers for use in earlier detection of the risk of progression towards HDC. As such, healthy lifestyle intervention and therapeutic decision-making could be delivered timely for the HDC prophylaxis from the perspectives of predictive, preventive and personalized/precision medicine.

#### **5. Conclusions**

This study presents conclusive examples of population-specific N-glycosylation of IgG associated with HDC in three northwestern Chinese Muslim ethnic minorities the KGZ, TJK, and UIG, as well as a Han Chinese population, respectively. These findings offer an opportunity to understand how IgG N-glycosylation profiles may be developed as potential biomarkers indicative of pro- or anti-inflammatory states shared by HTN, T2DM, and HDC. Further prospective follow-up studies exploring the onset of HDC using glycomics are expected to provide a more comprehensive understanding of specific glycomic biomarkers that may serve as tools for personalized monitoring or to help prevent the progression from HTN or T2DM towards HDC.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/jpm11070614/s1, Supplementary Table S1: The structures of 24 directly measured IgG glycan peaks by UPLC, Supplementary Table S2: Demographic and biochemical characteristics of the study subjects in the four Chinese ethnic groups: Kirgiz, Tajik, Uygur, and Han, Supplementary Table S3: The description of IgG glycome composition and the levels of IgG glycan traits in the HDC, HTN, T2DM patients, and healthy controls for the pooled samples from Chinese Muslim ethnic minorities and the Han Chinese samples, Supplementary Table S4: The associations of IgG N-glycan traits with HDC vs. HTN/T2DM/healthy individuals for the pooled samples from Chinese Muslim ethnic minorities and the Han Chinese samples, Supplementary Table S5: The differences of AUC between three type of models.

**Author Contributions:** Conceptualization, M.S., Y.W., G.L. and W.W.; data curation, W.C., D.Z. and X.G.; investigation, M.D., X.Z. and J.Z. (Jinxia Zhang); methodology, X.M., M.V., J.Š., J.Z. (Jie Zhang), H.W., A.M. and I.T.-A.; software, D.L., X.L. and L.W.; supervision, M.S., G.L. and W.W.; visualization, X.M.; writing—original draft, X.M., M.S.; writing—review and editing, M.S., H.W., I.T.-A., X.L., Y.W. and W.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by the National Natural Science Foundation of China (81573215, 31460285, and 81673247), and the China-Australian Collaborative Grant (NSFC 81561128020-NHMRC APP1112767). Genos has received funding from European Commission FP7 grants MIMOmics (contract #305280), HTP-GlycoMet (contract #324400) and PainOmics (contract #602736), and H2020 grants GlySign (contract #722095), GlyCoCan (contract #676412), SYSCID (contract #733100), and IMforFuture (contract #721815), as well as funding from the European Structural and Investments funds for projects "New generation of high throughput glycoanalytical services (contract #KK.01.2.1.01.0003) and "Croatian National Centre of Research Excellence in Personalized Healthcare" (contract#KK.01.1.1.01.0010).

**Institutional Review Board Statement:** The study was approved by the local community leaders and the Ethics Committee of Capital Medical University, Beijing, China (No. 2015SY21 and 2009SY16) and adhered to the principles of the Declaration of Helsinki.

**Informed Consent Statement:** Informed consent was obtained from all participants included in this study.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy restrictions.

**Acknowledgments:** The authors thank the directors and relevant staff of relevant Health Bureaus for their support for recruiting the Kirgiz, Tajik, and Uygur participants, and these volunteers and community leaders for their participation and support.

**Conflicts of Interest:** Gordan Lauc is the founder and owner of Genos Ltd., Zagreb, Croatia, a private research organization that specializes in high-throughput glycomic analysis and has several patents in this field. Ana Momˇcilovi´c, Irena Trbojevi´c-Akmaˇci´c, Jerko Štambuk, and Marija Vilaj are employees of Genos Ltd., Zagreb, Croatia.

#### **Abbreviations**


#### **References**


## *Article* **Identification of Adolescents with Adiposities and Elevated Blood Pressure and Implementation of Preventive Measures Warrants the Use of Multiple Clinical Assessment Tools**

**Hiba Bawadi <sup>1</sup> , Manal Kassab <sup>2</sup> , Abdel Hadi Zanabili <sup>3</sup> and Reema Tayyem 1,\***


**Abstract:** Abstract: BackgroundThe burden of abdominal adiposity has increased globally, which is recognized as a key condition for the development of obesity-related disorders among youth, including type 2 diabetes, cardiovascular disease, and hypertension. High blood pressure (BP) and cardiovascular diseases increase the rates of premature mortality and morbidity substantially. **Aims**: to investigate the relation between abdominal adiposity and elevated BP among adolescent males in Jordan. **Methods**: Nationally representative sample of male adolescents was selected using multi-cluster sampling technique. Study sample included 1035 adolescent males aged 12 to 17 years. Multiple indicators were used to assess adiposity including waist circumference (WC) and total body fat (TF), truncal fat (TrF), and visceral fat (VF). Systolic blood pressure was measured to assess hypertension. **Results**: After adjusting for age, smoking status, and physical activity, the odds of having stage two hypertension increased 6, 7, and 8 times for adolescents who were on 90th percentile or above for Trf, VF, and WC, respectively. **Conclusion**: Elevated BP was significantly associated with total and abdominal adiposity among adolescent males in Jordan. Use of multiple clinical assessment tools is essential to assess abdominal obesity among adolescents.

**Keywords:** abdominal adiposity; adolescents; blood pressure; obesity; Jordan; boys

#### **1. Introduction**

Obesity is considered the most serious future-health threatening problem among children and adolescents. In Jordan, about 17.3% of children and adolescents were overweight, and 15.7% were obese, as documented by Zayed and colleagues (2016) [1]. Recent research indicates that abdominal or central obesity is a more robust predictor for several debilitating and deadly chronic diseases than overall adiposity [2,3]. Besides being considered an important predictor of adult obesity, abdominal obesity among children and adolescents is also linked to several life threatening diseases, including cardiovascular diseases, diabetes, orthopaedic problems, cancers, metabolic syndrome, and high blood pressure [4–7]. Obesity does this through a variety of pathways, some as straightforward as the mechanical stress of carrying extra pounds and some involving multipart alterations in both hormones and cellular metabolism [8].

In the last decade, elevated BP increased among children and adolescents [9]. Adolescents with elevated BP can develop several chronic diseases and body organ damage [10,11]. Children with high BMI are more prone to have elevated BP [12]. Premature mortality is attributed to elevated BP by increasing the incidence of cardiovascular disease [13]. Prevention of childhood obesity initiates reduction in BP which leads to substantial prevention of cardiovascular disease [14]. Several studies supported the association between

**Citation:** Bawadi, H.; Kassab, M.; Zanabili, A.H.; Tayyem, R. Identification of Adolescents with Adiposities and Elevated Blood Pressure and Implementation of Preventive Measures Warrants the Use of Multiple Clinical Assessment Tools. *J. Pers. Med.* **2021**, *11*, 873. https://doi.org/10.3390/jpm11090873

Academic Editor: Francesco Masedu

Received: 17 June 2021 Accepted: 23 August 2021 Published: 31 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

elevated BP and obesity among children and adolescents [15,16]. In Jordan, there are not enough data available about rates of abdominal obesity among children and adolescents, nor the association between abdominal adiposity and elevated blood pressure among this vulnerable group. The limited data available in Jordan used a single assessment tool for obesity with no consideration for identification and specification of risk and precision of preventive treatment. Therefore, this study aimed to investigate the relation between abdominal adiposity and elevated BP among adolescent males in Jordan by implementing a panel of assessment tools.

#### **2. Materials and Methods**

#### *2.1. Study Design*

A population-based, nationally representative cross-sectional study was conducted, and included 1035 male adolescents from 18 elementary and high schools in Jordan across the whole kingdom of Jordan. The study protocol was examined and approved by Jordan University of Science, Technology, and Ministry of Education-Jordan.

A list of all males' schools (341 schools) in Jordan was obtained from Ministry of Education (MOE). Three geographical zones were identified based on MOE classification; north, mid, and south Jordan. Distribution of Schools over the geographical zones was also obtained from MOE; 130 schools in north area, 159 schools in middle area, and 52 schools in south area. Based on this distribution, a sample of 5% of schools (n = 18 schools) were selected as follows: 7 schools in northern Jordan, 8 schools in mid-Jordan, and 3 schools in southern Jordan. Classes in each school and adolescents in each class were selected using simple random selection.

#### *2.2. Data Collection*

All adolescents were asked to complete a self-administered questionnaire including data about socio demographic status, dietary habits, and physical activity level. The questionnaire was pilot tested to check for any vague, poorly stated items.

#### *2.3. Measurements*

#### 2.3.1. Total Adiposity

Total adiposity was assessed using body mass index (BMI) and total body fat (TF) percentage. Anthropometric measurements were performed by a trained dietitian. Weight was measured to the nearest 0.1 kg (Tanita BF-350, Tanita Corporation, Tokyo, Japan). Height was measured to the nearest 0.1 cm with a portable Harpenden stadiometer. BMI was then calculated as the ratio of weight (kg) to height (m) squared (kg/m<sup>2</sup> ) [17]. BMI for age charts were used to assess overweight and obesity [18]. BMI and age were plotted to obtain adolescents' BMI-age percentiles. Data were interpreted as follows: <5th percentile were assessed as underweight, 5th–85th percentile were assessed as normal weight, 85th–95th percentile were assessed as overweight, and ≥95th percentile were assessed as obese [19].

TF% was estimated using bioelectrical impedance technique (Tanita BF-350, Tanita Corporation, Tokyo, Japan). Age and gender sensitive cut-off points for TF% were used to assess obesity [19].

#### 2.3.2. Abdominal Adiposity

Several indicators were used to assess abdominal adiposity; waist circumference (WC), truncal fat (TrF) %, and visceral fat (VF). WC was measured at the narrowest area of the torso, as seen from the anterior aspect, to the nearest 0.1 cm with a non-elastic tape measure snugly fitted to measure WC at the level of the natural waist. The WC value corresponding to ≥90th percentile of WC for gender and age were used as cut-off values to identify adolescents with abdominal obesity [20].

Truncal fat (TrF)% and VF were estimated using bioelectrical impedance technique (Tanita AB-140 ViScan; Tanita Corporation, Tokyo, Japan). Adolescents were asked to lie back and expose their abdominal region. Electrodes were attached to the impedance meter and placed across the navel of the subject with the electrodes facing down. The positioning line was aligned on the impedance meter with the laser to adjust the position. Values associated with ≥90th percentile were considered high.

#### *2.4. Blood Pressure Measurement*

Blood pressure (BP) was measured by a registered nurse using validated automated oscillometric device (Omron HEM-7051, Omron Corporation, Kyoto, Japan) [21]. Blood pressure was measured according to a standardized protocol. Blood pressure was measured first on the right arm and then at left arm. Both readings were taken after the participants were allowed to sit and rest for 5 min to relieve anxiety with their back supported, feet on the floor, arm supported, and cubital fossa at heart level. The same technique was followed for right and left arm measurement. Two readings were obtained at a 1-min interval, and the average of the two readings was recorded. Systolic BP was determined by the onset of the "tapping" KSs (K1) and the fifth KS (K5), or the disappearance of KS, as the definition of Diastolic BP. After measuring the mean values of BP, the BP percentiles were determined accordingly.

The following criteria were used to categorize adolescents according to their BP: normal BP for adolescents with systolic BP <90th percentile; prehypertension for adolescents with 90th to <95th percentile; stage 1 hypertension for adolescents with 95th to <99th percentile and stage 2 hypertension for adolescents with >99th percentile [22]. Only Systolic BP was used in this study as indicator for hypertension, as the literature suggests that diastolic hypertension rarely occurs without Systolic hypertension in children and adolescents [23,24].

#### *2.5. Statistical Analysis*

Statistical analysis was performed using SPSS for Windows version 11.5 (Statistical Package for the Social Sciences). Data were expressed as frequencies and percentages. Multivariate logistic regression was performed to investigate the association between adolescent's risk of elevated systolic blood pressure and anthropometric measurements. Odds ratios (OR) and 95% confidence intervals (95% CI) were calculated and *p*-value was set at *p* < 0.05.

#### **3. Results**

This study included 1035 adolescents aged between 12 and 17 years. Socio-demographic characteristics of Jordanian adolescents were presented in Table 1. Half of participants were from Mid Jordan (49%), followed by North (35%), and South (16%). Participant's school grades ranged from 7th grade up to 11th grade. Students' parenteral educational levels were high school diplomas and college education in most cases. Few students reported either an illiterate mother or father.

**Table 1.** Socio-demographic characteristic of Jordanian adolescents.


**Table 1.** *Cont.*


\* Jordanian Dinar (JD) = USD 1.408. † education level: ≤12 years: include elementary and high schools, ≥12 years: include after high school education.

Table 2 presents data related to systolic hypertension according to their dietary habits. Students' self-reported regular consumption of vegetables (*p* = 0.008), fruits (*p* = 0.033) and chocolate (*p* = 0.013) were related to their systolic blood pressure. On the other hand, milk and dairy products, meat, fish, legumes, canned fruit juice, soda, and nuts were not significantly associated with blood pressure. Table 3 shows adjusted odds ratios for adolescents' SHTN based on their general and abdominal obesity. After adjusting for confounding variables—age, smoking status, and physical activity—it was found that BMI, WC, TF, TrF, and VF are positively related to increased odds of SHTN.

**Table 2.** Prevalence of hypertension according to adolescents food group intake.


Different letters means that there is a statistically significant difference between the groups.


**Table 3.** Multivariate logistic regression between adolescent's systolic hypertension and anthropometric measurements.

\* BMI for age was obtained and categorized according WHO cutoff point underweight (for age BMI < 5% percentile), normal weight (5% percentile < BMI < 85% percentile), overweight (85% percentile < BMI < 95% percentile), obese (BMI ≥ 95% percentile) (Riley and Bluhm, 2012). † TOTAL FAT cut off point were determined. ‡ SHTN refer to systolic hypertension.

> Odds of pre-hypertension increased in a dose–response pattern among overweight (OR, 2.3; 95% CI, 1.4–3.8) and obese (OR, 2.7; 95% CI, 1.6–4.5) adolescents as compared to adolescents with normal body weight. Similarly, odds of stage 1 hypertension increased in a dose–response pattern across overweight (OR,1.8; 95% CI, 1.01–3.1) and obesity (OR, 2.7; 95% CI, 1.6–4.6) classes. A similar trend was observed for stage 2 SHTN.

> With regard to WC, and after adjusting for age, smoking status, and physical activity, adolescents with the highest 10th percentile of WC had 5 times higher odds of having stage 1 hypertension (OR, 95% CI, 2.3–10.6) and 8.6 higher odds of having stage 2 hypertension (OR, 95% CI, 4,5–16.4). A similar relation was found with total body fat. Obese adolescents had twice as high odds of having stage 1 hypertension (OR, 95% CI, 1.3–3.8) and 3.5 times higher odds of having stage 2 hypertension (OR, 95% CI, 2.3–5.5).

> With regard to TrF, and after adjusting for age, smoking status, and physical activity, adolescents with the highest 10th percentile of TF had 4 times higher odds of having stage 1 hypertension (OR, 95% CI, 2.2–7.1) and 5.5 higher odds of having stage 2 hypertension (OR, 95% CI, 3.4–9.1). Similar relation was found with VF. Adolescents with the highest 10th percentile of VF had 4.2 times higher odds of having stage 1 hypertension (OR, 95% CI, 1.8–9.4) and 7.2 times higher odds of having stage 2 hypertension (OR, 95% CI, 3.7–14.1).

#### **4. Discussion**

Health consequences of overweight and obesity in adolescence are strongly associated with risk factors for chronic health conditions such as cardiovascular disease, diabetes, and elevated blood pressure [4,6].

Fruit and vegetable consumption was significantly related to blood pressure (*p* = 0.008 and *p* = 0.033, respectively). This finding was consistent with the findings of several studies [25–28]. Apple et al. (2006) studied the contents of fruit and vegetables for vitamins, minerals, and fibers. Fruits and vegetables also contained potassium; and this increase consumption of potassium intake was associated with reduction in blood pressure. The increase in potassium intake had same lowering effect on blood pressure as a decrease in sodium intake. Potassium plays a major role in balancing out the negative effects of sodium. Whelton et al. (1997) recommended potassium for prevention and treatment of hypertension. Increasing serum levels of vitamins A, C, E [25], and D [28] were associated with lowering blood pressure. Meta-analysis suggested that increasing the dietary fiber intake had a lowering effect on blood pressure [27]. According to this study, chocolate was significantly related to blood pressure (*p* = 0.013). Studies explained chocolate's role to lower blood pressure [29–31]. Chocolate contained coca that include polyphenols, especially flavones. Strong effects of flavones on blood pressure as a vasodilator were applied by increasing the formation of endothelial nitric oxide.

Findings of the current study are in line with previously published studies suggesting that obesity is an independent risk factor for early hypertension [4,6,23,24]. Abdominal adiposity showed stronger association that between fatness and the risks of stage 1 and stage 2 systolic hypertension among study participants. Risks increased up to 8 times with abdominal adiposity indicators. Our findings were similar to those reported by Pausova et al. (2012) and Pazin et al. (2020), where an association was found between adolescents abdominal adiposity and high BP, especially among in boys [32,33].

The relationship between abdominal adiposity and elevated blood pressure is not fully explained yet. However, researchers linked abdominal adiposity to elevated blood pressure via insulin resistance and related inflammation [8,24,34–37]. Abdominal adiposity stimulates the state of insulin resistance and the inflammatory process [38]. Inflammation leads several changes vascular endothelial function, which in turn leads to the development of hypertension [39]. Visceral fat that accumulates around individual internal organs decreases body sensitivity to insulin, which contributes to oxidative stress, inflammation, vascular endothelial dysfunction, and hypertension [35,37].

Landsberg et al. (2013) explained that the pathophysiology of fat accumulation in the abdominal region led to an increase in BP through the stimulation of insulin resistance, increasing central nervous system activity, renin-angiotensin-aldosterone system activity, angiotensinogen from intra-abdominal adipocytes, aldosterone production, and the renal sodium reabsorption [40]. High adiposity can lead to adiposopathy in adolescence, with associated increases in inflammation and oxidative stress, changes in adipokines and endocrinopathy [39]. This manifests as cardio-metabolic risk factors in similar patterns to those in noted in obese adults [5]. Both obesity and cardio-metabolic risk factors have been shown to be associated with vascular changes indicative of early hypertension, atherosclerosis, as well as ventricular hypertrophy, dilation, and dysfunction [8,24].

Ke and colleagues (2018) clarified that high BP is being induced by VF due to high levels of free fatty acids excretion in liver via portal circulation after lipogenesis and lipolysis activity, gluconeogenesis, lipid synthesis, and insulin resistance [41]. This elevated amount of free fatty acids will induce hypertension and eventually atherosclerosis [42,43].

Many studies revealed that applying different anthropometric assessment methods is essential to evaluate body weight status accurately [44–47]. Body composition can be assessed at the molecular, cellular, and tissue levels using several different methods [44,46]. Specific limitations of anthropometry in the obese are obvious. These limitations include the inability to distinguish subcutaneous fat from visceral adipose tissue, which is helpful to assess disease risk [47] and accuracy may be lowered in the severely obese, due to difficulties finding the actual waistline or drooping abdominal fat can interfere with hip measurement [46]. Therefore, no single measurement is recommended to assess the body composition of obese people. Yet, each modality has benefits and drawbacks [46].

This study has several strengths, including the representation of the sample to the Jordanian adolescent males and multiple tools used to operationalize adiposity. However, the study is limited as it did not include intensive dietary data which were not to feasible in the school setting, especially with all the measurements being performed.

In conclusion, excess fat, especially abdominal fat, is associated with increased risk of systolic hypertension among adolescent males. The physiological relation between obesity and hypertension could not attribute to a single factor. It is critical to use multiple clinical assessment tools for obesity to optimizes identification of adolescents at risk; hence preventive measures can be timely considered. Findings of this study must be taken with high level of importance to prevent early hypertension and associated cardiovascular diseases risk later in adulthood.

**Author Contributions:** H.B.: proposal writing, data analysis, data interpretation, and writing manuscript. M.K.: data interpretation and writing manuscript. A.H.Z.: data collection and writing manuscript. R.T.: data interpretation and writing manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Deanship of Research in Jordan University of Science and Technology (project No. 2014/97).

**Institutional Review Board Statement:** This study was conducted according to the guidelines laid down in the Declaration of Helsinki and all procedures involving research study participants were approved by the IRB at Ministry of Education. IRB Approval number was 3/10/21299.

**Informed Consent Statement:** Written informed consent was obtained from all subjects/parents.

**Acknowledgments:** Authors would like to thank Fawzieh Hammad for her help in the data entry and analysis.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Sex-Specific Association of Uric Acid and Kidney Function Decline in Taiwan**

**Po-Ya Chang 1,\* , Yu-Wei Chang <sup>2</sup> , Yuh-Feng Lin <sup>3</sup> and Hueng-Chuen Fan <sup>4</sup>**


**Abstract:** An elevated serum urate concentration is associated with kidney damage. Men's uric acid levels are usually higher than women's. However, postmenopausal women have a higher risk of gout than men, and comorbidities are also higher than in men. This study examined the sex differences in the relationship between hyperuricemia and renal progression in early chronic kidney disease (CKD) and non-CKD, and further examined the incidence of CKD in non-CKD populations among patients over 50 years of age. We analyzed 1856 women and 1852 men participating in the epidemiology and risk factors surveillance of the CKD database. Women showed a significantly higher risk of renal progression and CKD than men within the hyperuricemia group. After adjusting covariates, women, but not men resulted in an hazard ratio (HR) for developing renal progression (HR = 1.12; 95% CI 1.01–1.24 in women and HR = 1.03; 95% CI 0.93–1.13 in men) and CKD (HR = 1.11; 95% CI 1.01–1.22 in women and HR = 0.95; 95% CI 0.85–1.05 in men) for each 1 mg/dL increase in serum urate levels. The association between serum urate levels and renal progression was stronger in women. Given the prevalence and impact of kidney disease, factors that impede optimal renal function management in women and men must be identified to provide tailored treatment recommendations.

**Keywords:** hyperuricemia; sex; renal progression; chronic kidney disease

#### **1. Introduction**

Uric acid is the final product of the breakdown of an exogenous pool of purines and endogenous purine metabolism [1], and most uric acid is excreted via the kidneys [2]. Elevated serum urate concentrations are associated with gout [3], diabetes mellitus [4,5], hypertension [6,7], metabolic syndrome [8] and cardiovascular mortality [9,10]. A metaanalysis of 1,134,073 participants showed a positive dose-response association between uric acid levels and cardiovascular disease mortality risk [11]. Hyperuricemia is common in chronic kidney disease (CKD) patients [12], and several studies have also shown a relationship between serum urate levels and kidney damage [13–16]. Urate-lowering therapy statistically significantly improved estimated glomerular filtration rate (eGFR) and serum creatinine and serum creatinine [17], and reduced the incidence of kidney failure events [18]. Uric acid crystallizations induce endothelial dysfunction [19,20], systemic inflammation [21], oxidative stress, and renin-angiotensin-aldosterone system activity [22], which may provide a pathogenic mechanism for chronic urate nephropathy.

Previous evidence has suggested that there are sex differences in uric acid levels; men have a higher uric acid level than women [23], and men have a higher prevalence of gout than women [24]. Sex differences decline with increasing age, but men still far exceed women in the prevalence of gout, even in older populations [24]. Age-related factors

**Citation:** Chang, P.-Y.; Chang, Y.-W.; Lin, Y.-F.; Fan, H.-C. Sex-Specific Association of Uric Acid and Kidney Function Decline in Taiwan. *J. Pers. Med.* **2021**, *11*, 415. https:// doi.org/10.3390/jpm11050415

Academic Editor: Rutger A. Middelburg

Received: 31 March 2021 Accepted: 12 May 2021 Published: 15 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and hormone are also associated with uric acid levels [25]. Premenopausal women have lower levels of uric acid than men or postmenopausal women because estrogens stimulate kidney clearance, secretion and reabsorption [25,26]. However, there is evidence that postmenopausal women are at higher risk of gout than men [27]. Moreover, women with hyperuricemia had more comorbidities than men with hyperuricemia. Women with high uric acid levels had a higher risk of developing metabolic syndrome [28], diabetes [29,30], hypertension [7], coronary artery disease [31,32], and cardiovascular events [9,23] than men. Sex differences have appeared in the physiological regulation of urate homeostasis. Estrogen may play an important role in the regulation of the expression or activity of urate transporters [33].

In addition, the effects of uric acid on renal function are sex-specific. Both high and low uric acid levels increase the risk of developing chronic kidney disease (CKD) in men, while only high levels are a risk factor in women [34]. A historic cohort study in Japan indicated that uric acid levels greater than 6.0 mg/dL are associated with kidney impairment in men with mild kidney dysfunction [13]. However, a screened study conducted in Japan showed that higher uric acid levels increase the risk of developing ESRD in women but not in men [35].

Sex-specificity exists in renal progression in patients with early CKD [36,37]. In fact, whether the incidence of CKD and renal progression differ between postmenopausal women and men with hyperuricemia is unclear. Early detection of loss of kidney function and early treatment can decrease the risk of complications [38], prevent disease progression and prevent high medical costs in later stages [39]. Therefore, we investigated sex differences in the relationship between hyperuricemia and renal progression in early CKD and non-CKD populations, and further examined the incidence of CKD in non-CKD populations among patients over 50 years of age.

#### **2. Materials and Methods**

#### *2.1. Study Population*

A cohort based on the Epidemiology and Risk Factors Surveillance of CKD database was recruited from 14 medical centers and community in Taiwan, namely Tri-Service General Hospital, Taipei Medical University Hospital, Shuang Ho Hospital, Taipei Chang Gung Memorial Hospital of the C.G.M.F., Cardinal Tien Hospital, E-DA Hospital, National Cheng Kung University Hospital, China Medical University Hospital, Show Chwan Memorial Hospital, Changhua Christian Hospital, Taipei Veterans General Hospital, Kaohsiung Medical University Chung-Ho Memorial Hospital, Kaohsiung Municipal Siaogang Hospital, and Kaohsiung Chang Gung Memorial Hospital from October 2008 to February 2016. A total of 10,823 participants aged ≥18 years who had at least two serum creatinine measurements were recruited in this study. For this study, we excluded patients according to the following criteria: renal function could not be assessed, the follow-up periods were less than 12 months, missing data on uric acid or younger than 50 years of age. In total, 6085 patients were enrolled, with 1728, 1980 and 2377 having non-CKD, early CKD and late CKD, respectively. In the present study, we examined the risk factors for renal progression in subjects with non-CKD and early CKD, as well as the incidence of CKD in subjects with non-CKD (Figure 1).

**Figure 1.** Flowchart of participants analyzed in this study. **Figure 1.** Flowchart of participants analyzed in this study.

The studied medical centers used the same medical laboratory criteria and protocols, and thus, the serum creatinine values from the different medical centers could be compared and standardized. We measured eGFR changes at the individual level and re-examined patients at the same medical center to control individual variation. After a thorough explanation of the study, all participants provided informed consent before data collection. This study was approved by the Ethics Committees of the Tri-Service General Hospital (TSGHIRB 100-05-197), Taipei Medical University (TMU-JIRB 20124036), Kaohsiung Medical University Chung-Ho Memorial Hospital (KMUHIRB 20120019), National Cheng Kung University Hospital (A-ER-101-117), Changhua Christian Hospital (CCHIRB 120405), Kaohsiung Chang Gung Memorial Hospital (101-1096B), Cardinal Tien Hospital (TMU-JIRB 201204035) and China Medical University Hospital (DMR101-IRB2- 273(CR-1)). All the research methods were implemented in accordance with guidelines approved by the Ethics Committees. The studied medical centers used the same medical laboratory criteria and protocols, and thus, the serum creatinine values from the different medical centers could be compared and standardized. We measured eGFR changes at the individual level and re-examined patients at the same medical center to control individual variation. After a thorough explanation of the study, all participants provided informed consent before data collection. This study was approved by the Ethics Committees of the Tri-Service General Hospital (TSGHIRB 100-05-197), Taipei Medical University (TMU-JIRB 20124036), Kaohsiung Medical University Chung-Ho Memorial Hospital (KMUHIRB 20120019), National Cheng Kung University Hospital (A-ER-101-117), Changhua Christian Hospital (CCHIRB 120405), Kaohsiung Chang Gung Memorial Hospital (101-1096B), Cardinal Tien Hospital (TMU-JIRB 201204035) and China Medical University Hospital (DMR101-IRB2-273(CR-1)). All the research methods were implemented in accordance with guidelines approved by the Ethics Committees.

#### *2.2. Study Variables and Definitions 2.2. Study Variables and Definitions*

The definition of CKD was based on the Kidney Disease Outcomes Quality Initiative guidelines [40]. eGFR was calculated using the chronic nephrology epidemiology collaboration (CKD-EPI) equation [eGFR (mL/min/1.73 m<sup>2</sup> ) = 141 × min (serum creatinine (SCr)/κ, 1)<sup>α</sup> × max SCr/κ, 1)−1.209 × 0.993age × 1.018 (if women) × 1.159 (if black)], κ = 0.7 (women) and 0.9 (men), α = −0.329 (women) and −0.411(men); min indicates the minimum of SCr/κ or 1, and max indicates the maximum of SCr/κ or 1 [41]. CKD stage 1 was defined The definition of CKD was based on the Kidney Disease Outcomes Quality Initiative guidelines [40]. eGFR was calculated using the chronic nephrology epidemiology collaboration (CKD-EPI) equation [eGFR (mL/min/1.73 m<sup>2</sup> ) = 141 × min (serum creatinine (SCr)/κ, 1)<sup>α</sup> <sup>×</sup> max SCr/κ, 1)−1.209 <sup>×</sup> 0.993age <sup>×</sup> 1.018 (if women) <sup>×</sup> 1.159 (if black)], <sup>κ</sup> = 0.7 (women) and 0.9 (men), α = −0.329 (women) and −0.411(men); min indicates the minimum of SCr/κ or 1, and max indicates the maximum of SCr/κ or 1 [41]. CKD stage 1 was defined as a eGFR value of <sup>≥</sup>90 mL/min/1.73 m<sup>2</sup> and the presence of kidney damage (i.e., proteinuria dipsticks ≥1+, urine protein-to-creatinine ratio (UPCR) ≥150, or urine albumin-tocreatinine ratio (UACR) <sup>≥</sup> 30); CKD stage 2 was defined as eGFR = 60–89 mL/min/1.73 m<sup>2</sup> and the presence of kidney damage (i.e., proteinuria dipsticks ≥1+, urine protein-tocreatinine ratio (UPCR) ≥150, or urine albumin-to-creatinine ratio (UACR) ≥30); CKD stage 3a was defined as eGFR = 45–59 mL/min/1.73 m<sup>2</sup> ; CKD stage 3b was defined as eGFR = 30–44 mL/min/1.73 m<sup>2</sup> ; CKD stage 4 was defined as eGFR = 15–29 mL/min/1.73 m<sup>2</sup> ; and CKD stage 5 was defined as eGFR <15 mL/min/1.73 m<sup>2</sup> [42]. CKD stages 1, 2, and 3a referred to patients with early CKD, whereas CKD stages 3b, 4, and 5 referred to patients with late CKD.

According to the baseline eGFR, renal progression was defined as an eGFR that was reduced by more than 25% [43]. Non-CKD patients were followed up from the baseline measurement of renal function to the development of renal progression and CKD, or the end of the study period. Early CKD patients were followed up from the baseline measurement of renal function to the development of renal progression, or the end of the study period. The patients were censored at the end of the follow-up period.

Uric acid status was classified as the absence of hyperuricemia (serum urate concentration values of <7.0 mg/dL for men and <6.0 mg/dL for women) or the presence of hyperuricemia (serum urate concentration values of ≥7.0 mg/dL for men and ≥6.0 mg/dL for women) [44,45].

Sociodemographic characteristics, comorbidities, and health-related behaviors, including age, sex, hypertension, diabetes mellitus, dyslipidemia, CKD, gout, stroke, cigarette smoking, and alcohol consumption, were obtained at baseline using a structured questionnaire. Serum urate levels and physical examination variables, including height, weight, body mass index (BMI), serum creatinine, baseline eGFR, systolic blood pressure (SBP), diastolic blood pressure (DBP), fasting glucose, and total cholesterol, were obtained via medical chart review. Cigarette smoking was divided into never smoking, current smoker (had smoked at least 100 cigarettes in the lifetime and is currently smoking cigarettes) or former smoker (had smoked at least 100 cigarettes in the lifetime but had quit smoking at the time of the interview) [46,47]. Alcohol consumption was divided into current alcohol consumers versus participants who never consumed alcohol.

#### *2.3. Statistical Analysis*

To determine the significant differences between women and men, participants were examined using the chi-square test and Student's t test for categorical and continuous variables, respectively. The incidence density method was used to estimate the uric acid status and sex-specific incidence rates of renal progression and CKD. Multivariate Cox proportional hazards analysis was used to estimate the hazard ratio (HR) and 95% confidence interval (CI) of the risk of renal progression and CKD associated with sex. Covariates included age, CKD, hypertension, diabetes mellitus, dyslipidemia, gout, stroke, BMI, serum creatinine, cigarette smoking, and alcohol consumption.

To induce a similar distribution of covariates between women and men, we used the inverse probability of treatment weighting (IPTW) analysis method, where the weight was the inverse (i.e., reciprocal) of the propensity scores [48]. To assess the balance of baseline covariates between women and men, we applied absolute standardized mean differences. In addition, the value of the absolute standardized mean difference was required to be lower than 10% between the sexes, indicating negligible differences among the covariates. In addition, the value of the absolute standardized mean difference was less than 10% between sexes, indicating negligible differences among the potential confounders [49]. The risks of renal progression and CKD among the groups of women and men were compared using Cox proportional hazards regression with IPTW. In order to compare the results of premenopausal and postmenopausal women, we further analyzed the sex differences in the relationship between hyperuricemia and renal progression and the incidence of CKD among patients before 50 years of age. Analyses were performed using SAS software (SAS Institute), version 9.4.

#### **3. Results**

Table 1 shows the demographic and clinical characteristics of the women and men. Of the 3708 participants, 1856 (50.05%) and 1852 (49.95%) were women and men, respectively. Before employing IPTW, the women exhibited significantly lower prevalences of hypertension (1124 (60.56%)), CKD (859 (46.28%)), gout (118 (6.36%)), stroke (87 (4.69%)), current smoking (36 (1.96%)), and alcohol consumption (63 (3.43%)) as well as low levels of serum urate (5.49 (SD 1.36)), serum creatinine (0.73 (SD 0.15)), and DBP (76.02 (SD 10.71)) than men (Table 1). After applying IPTW, the different sexes were well balanced with respect to key characteristics, such as fasting glucose, BMI, SBP, DBP, age, diabetes mellitus, total cholesterol, baseline eGFR, dyslipidemia, stroke, hypertension, CKD, and gout (all absolute standardized mean difference lower than 10%; Figure 2). The mean levels of uric acid in men was higher than that in women in each age group (18–29 group, *p* < 0.001; 30–39 group, *p* < 0.001; 40–49 group, *p* < 0.001; 50–59 group, *p* < 0.001; 60–69 group, *p* < 0.001; 70–79 group, *p* < 0.001; ≥80, *p* = 0.214), but levels of uric acid were increases with advancing age in women but not in men. Both in women and men, uric acid levels showed statistically significant differences in each age group (women, *p* < 0.001; men *p* = 0.004) (Figure 3).

**Table 1.** The baseline characteristics subjects by sex.


Abbreviations: CKD, chronic kidney disease; BMI, body mass index; SBP, systolic blood pressure; DBP, diastolic blood pressure. <sup>a</sup> Female <sup>≥</sup> 6.0 mg/dL; Male ≥ 7.0 mg/dL. \* *p* < 0.05, \*\* *p* < 0.01, \*\*\* *p* < 0.001.

Serum uric acid (mg/dL)

Abbreviations: CKD, chronic kidney disease; BMI, body mass index; SBP, systolic blood pressure; DBP, diastolic blood

pressure. <sup>a</sup>Female ≥ 6.0 mg/dL; Male ≥ 7.0 mg/dL. \* *p* < 0.05**,** \*\* *p* < 0.01, \*\*\* *p* < 0.001.

*J. Pers. Med.* **2021**, *11*, 415 6 of 13

Abbreviations: CKD, chronic kidney disease; BMI, body mass index; SBP, systolic blood pressure; DBP, diastolic blood

Alcohol consumption 63 (3.43) 353 (19.34) <0.001 \*\*\*

pressure. <sup>a</sup>Female ≥ 6.0 mg/dL; Male ≥ 7.0 mg/dL. \* *p* < 0.05**,** \*\* *p* < 0.01, \*\*\* *p* < 0.001.

**Figure 2.** Women and men by characteristics before and after propensity weighting.

**Figure 2.** Women and men by characteristics before and after propensity weighting. **Figure 2.** Women and men by characteristics before and after propensity weighting.

Men 7.6 7.06 6.87 6.86 6.89 6.93 6.81 **Figure 3.** Bar chart of uric acid concentrations in each age group.

Women 5.25 5.32 5.55 5.85 5.99 6.31 6.61

**Figure 3.** Bar chart of uric acid concentrations in each age group. Within the group without hyperuricemia, the incidence rates of renal progression in men and women were 4.86 and 3.79 per 100 patient-years, respectively; the incidence rates **Figure 3.** Bar chart of uric acid concentrations in each age group. Within the group without hyperuricemia, the incidence rates of renal progression in men and women were 4.86 and 3.79 per 100 patient-years, respectively; the incidence rates of CKD in men and women were 13.31 and 9.54 per 100 patient-years, respectively. Within the hyperuricemia group, the incidence rates for renal progression in men and women Within the group without hyperuricemia, the incidence rates of renal progression in men and women were 4.86 and 3.79 per 100 patient-years, respectively; the incidence rates of CKD in men and women were 13.31 and 9.54 per 100 patient-years, respectively. Within the hyperuricemia group, the incidence rates for renal progression in men and women were 5.44 and 6.15 per 100 patient-years, respectively; the incidence rates for CKD in men and women were 14.54 and 15.26 per 100 patient-years, respectively (Table 2).

of CKD in men and women were 13.31 and 9.54 per 100 patient-years, respectively. Within the hyperuricemia group, the incidence rates for renal progression in men and women After using IPTW, women showed a significantly higher risk of renal progression (HR 2.21; 95% CI 1.55–3.14) and CKD (HR 2.69; 95% CI 1.90–3.81) than men within the group without hyperuricemia. Women showed a significantly higher risk of renal progression (HR 2.14; 95% CI 1.35–3.40) and CKD (HR 3.56; 95% CI 2.04–6.25) than men within the hyperuricemia group (Table 3).

In women, the multivariate adjustment for age, CKD, hypertension, diabetes mellitus, dyslipidemia, gout, stroke, BMI, serum creatinine, cigarette smoking, and alcohol consumption resulted in an HR for developing renal progression of 1.12 (95% CI, 1.01–1.24) and for CKD of 1.11 (95% CI, 1.01–1.22) for each 1 mg/dL increase in serum urate levels. In men, there was no association between serum urate levels and developing renal progression and CKD (Figure 4).

**Table 2.** Incidence rate for risk of renal progression (25% decline in GFR) and CKD by uric acid in women and men.


Abbreviation: CKD, chronic kidney disease; CI, confidence interval; IR, incidence rate. <sup>a</sup> IR/100 person-years. <sup>b</sup> Female <sup>≥</sup> 6.0 mg/dL; Male ≥ 7.0 mg/dL.

**Table 3.** Risk of renal progression and CKD by hyperuricemia in women and men.


Abbreviation: CKD, chronic kidney disease; HR, hazard ratio; CI, confidence interval. <sup>a</sup> Adjusted for age, CKD, hypertension, diabetes mellitus, dyslipidemia, gout, stroke, BMI, serum creatinine, cigarette smoking, and alcohol consumption. <sup>b</sup> Adjusted for age, hypertension, diabetes mellitus, dyslipidemia, gout, stroke, BMI, serum creatinine, cigarette smoking, and alcohol consumption. <sup>c</sup> *p* < 0.05; <sup>d</sup> *p* < 0.01; <sup>e</sup> *p* < 0.001.

> We further analyzed the results of patients before the age of 50. Women did not have a higher risk of renal progression and CKD than men within the group without hyperuricemia (Table S1).

**Figure 4.** Risk of renal progression and CKD per 1 mg/dL increase in serum uric acid by sex. Adjusted for age, CKD, hypertension, diabetes mellitus, dyslipidemia, gout, stroke, BMI, serum creatinine, cigarette smoking, and alcohol consumption. **Figure 4.** Risk of renal progression and CKD per 1 mg/dL increase in serum uric acid by sex. Adjusted for age, CKD, hypertension, diabetes mellitus, dyslipidemia, gout, stroke, BMI, serum creatinine, cigarette smoking, and alcohol consumption.

#### **4. Discussion 4. Discussion**

The multicenter cohort study enrolled subjects with early CKD and without CKD to demonstrate that serum urate levels are an independent predictor of renal progression and incident CKD. We found significant differences between men and women. The effect of sex differences in the role of serum urate levels on the long-term prognosis of patients in terms of renal progression and incident CKD among early CKD and non-CKD has not previously been fully explored in middle-aged and elderly adults. To our knowledge, this is the first study to elucidate sex differences in serum urate levels for the prognosis of The multicenter cohort study enrolled subjects with early CKD and without CKD to demonstrate that serum urate levels are an independent predictor of renal progression and incident CKD. We found significant differences between men and women. The effect of sex differences in the role of serum urate levels on the long-term prognosis of patients in terms of renal progression and incident CKD among early CKD and non-CKD has not previously been fully explored in middle-aged and elderly adults. To our knowledge, this is the first study to elucidate sex differences in serum urate levels for the prognosis of patients with kidney damage among middle-aged and elderly adults with early CKD and without CKD.

patients with kidney damage among middle-aged and elderly adults with early CKD and without CKD. Serum urate levels are a well-known marker of hypertension, gout, arterial stiffness, cardiac hypertrophy, metabolic syndrome, and cardiovascular damage [3,8,9,50–53]. Our findings on the effect of serum urate levels on the risk of renal progression among individuals with early CKD are consistent with previous studies. In patients with renal disease, a higher uric acid concentration was a strong independent predictor of kidney failure in earlier stages of CKD [16]. Srivastava et al. described this relationship of uric acid concentration with all-cause mortality as J-shaped in individuals with CKD [16]. High uric acid was found to be a risk factor for CKD and a rapid decline in eGFR in patients with mild kidney dysfunction [13]. We observed that hyperuricemia is an independent risk Serum urate levels are a well-known marker of hypertension, gout, arterial stiffness, cardiac hypertrophy, metabolic syndrome, and cardiovascular damage [3,8,9,50–53]. Our findings on the effect of serum urate levels on the risk of renal progression among individuals with early CKD are consistent with previous studies. In patients with renal disease, a higher uric acid concentration was a strong independent predictor of kidney failure in earlier stages of CKD [16]. Srivastava et al. described this relationship of uric acid concentration with all-cause mortality as J-shaped in individuals with CKD [16]. High uric acid was found to be a risk factor for CKD and a rapid decline in eGFR in patients with mild kidney dysfunction [13]. We observed that hyperuricemia is an independent risk factor for incident CKD among individuals without CKD. In the general population, uric acid levels promote initial kidney damage and incident kidney disease [14,54,55].

factor for incident CKD among individuals without CKD. In the general population, uric acid levels promote initial kidney damage and incident kidney disease [14,54,55]. However, a recently meta-analysis of randomized controlled trials demonstrated that urate-lowering therapy was associated with slowing the decline rate of GFR, but there is insufficient evidence to support urate lowering in patients to improve kidney outcomes [12]. Previous Mendelian randomized study found that there is no evidence to support a causal relationship between serum uric acid levels and eGFR levels or incident CKD [56]. The Controlled Trial of Slowing of Kidney Disease Progression from the Inhibition of Xanthine Oxidase (CKD-FIX) demonstrated that allopurinol dose not slow down the progression of CKD in patients with severe CKD and high risk of progression [57]. The Preventing Early Renal Loss in Diabetes (PERL) trial did not show a beneficial effect of reduction of the serum urate level with allopurinol therapy on the kidney outcomes in persons with type 1 diabetes and early-to-moderate diabetic kidney disease [58]. The PERL and CKD-However, a recently meta-analysis of randomized controlled trials demonstrated that urate-lowering therapy was associated with slowing the decline rate of GFR, but there is insufficient evidence to support urate lowering in patients to improve kidney outcomes [12]. Previous Mendelian randomized study found that there is no evidence to support a causal relationship between serum uric acid levels and eGFR levels or incident CKD [56]. The Controlled Trial of Slowing of Kidney Disease Progression from the Inhibition of Xanthine Oxidase (CKD-FIX) demonstrated that allopurinol dose not slow down the progression of CKD in patients with severe CKD and high risk of progression [57]. The Preventing Early Renal Loss in Diabetes (PERL) trial did not show a beneficial effect of reduction of the serum urate level with allopurinol therapy on the kidney outcomes in persons with type 1 diabetes and early-to-moderate diabetic kidney disease [58]. The PERL and CKD-FIX trials do not support the clinical benefits of allopurinol in the treatment of asymptomatic hyperuricemia on kidney outcomes [57,58].

FIX trials do not support the clinical benefits of allopurinol in the treatment of asymptomatic hyperuricemia on kidney outcomes [57,58]. Although, the causal relationship between the effects of uric acid and CKD is still debatable. In a perspectives article argues that hyperuricemia is directly related to the Although, the causal relationship between the effects of uric acid and CKD is still debatable. In a perspectives article argues that hyperuricemia is directly related to the development and progression of CKD, and treatment with urate-lowering therapy may provide clinical benefits [59]. Hyperuricemia has a detrimental effect on kidney function [59].

development and progression of CKD, and treatment with urate-lowering therapy may Hyperuricemia accelerates renal progression by increasing renal renin and cyclooxygenase-2 expression and causing vascular smooth muscle cell proliferation [60]. Oxidative stress plays an important pathophysiological role in renal progression [22]. Furthermore, impaired kidney clearance is responsible for increased uric acid concentrations [16].

There are sex differences in hypertension and cardiovascular events in patients with high serum urate levels, and the rates in female patients are significantly higher than those in males [9,61]. The present study demonstrated that sex differences in the association between uric acid level and renal progression and incident CKD were stronger in women than in men. A previous study in Okinawa, Japan also found a significant association between uric acid and renal failure requiring renal replacement in women, but no significant association was observed in men [35].

It is well known that men have higher uric acid levels than women. Our findings also showed that men had higher levels of uric acid than women in each age group; however, we found that uric acid levels increased with advancing age in women but not in men. Women before the age of 50 did not have a higher risk of renal progression and CKD than men. Uric acid levels increase with age, and after menopause, the levels roughly equal the levels in men in later years [62]. A previous study in Japan showed that women uric acid levels increased until the age of 70 years or older as a result of changes in hormones accompanying menopause [63]. This finding provides a potential explanation for why women have a higher risk of kidney damage than men.

Young women have higher levels of estrogen, which can promote renal uric acid excretion and reduce uric acid levels [26,64]. The Third National Health and Nutrition Examination Survey reported that postmenopausal women receiving hormone replacement therapy may have reduced uric acid levels [25]. Further investigations are encouraged to investigate the mechanism by which hyperuricemic women have a higher risk of reduced renal function than men among elderly populations.

Our findings indicate that young men aged 18–29 and 30–39 have higher serum urate levels and decrease with age, while the opposite is the case in women. If young people aged 18–40 have normal kidney function and they develop hyperuricemia, the most likely clinical indication is genetic variants in urate transporters, such as SLC2A9 [65], rather than underlying chronic kidney disease. Additionally, there is extensive evidence that variants in ABCG2 cause hyperuricemia in young people [66–68]. Thus, predicting whether these patients develop CKD or whether hyperuricemia contributes to CKD progression later in life is impossible due to genetic components and environmental factors as previously described in a study in Taiwan [69]. Therefore, in this study, the sex-specific relationship between hyperuricemia and renal progression was studied, with a focus on patients over 50 years of age.

In addition, we found that younger women with hyperuricemia have a higher risk of kidney deterioration than men with hyperuricemia. A study of a military cohort of young individuals also found that the association between hyperuricemia and raised blood pressure was stronger in women than in men [70]. One possible explanation for the sex difference is that women with hyperuricemia may have higher levels of xanthine oxidase activity and reactive oxygen species productions than men, leading to more severe renal vascular inflammation in women than in men [71,72].

This study has several advantages. First, we had a large cohort that recruited patients from multiple centers throughout Taiwan. Second, we provided sex-specific data showing a difference in the relationship between hyperuricemia and renal progression. Third, this study used well-trained interviewers to collect data on demographic characteristics and health-related behaviors face to face. Finally, in order to reduce the effects of confounding in observational study, we used the IPTW analysis method to induce similar covariate distribution between women and men.

However, our study has some limitations. First, this research may have potential selection bias because our patients were voluntarily recruited. Second, certain test results may be underestimated due to the use of structured questionnaires to collect clinical disease data. Third, this study does not contain information on potential factors, including menopausal age and hormone replacement therapy. Fourth, this study uses structured questionnaires to collect alcohol consumption data, but does not include daily alcohol consumption. Therefore, it is difficult to accurately analyze the possible relationship between alcohol consumption with hyperuricemia and renal decline. Therefore, these factors should be considered in future research.

#### **5. Conclusions**

Given the prevalence and impact of kidney disease, factors that impede optimal renal function management in women and men must be identified to provide tailored treatment recommendations. Elevated serum urate levels are a strong risk factor for developing kidney disease. Therefore, our results suggest that individuals, especially postmenopausal women, with higher uric acid levels are at increased risk of developing renal progression and incident CKD.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/jpm11050415/s1, Table S1: Risk of renal progression and CKD by hyperuricemia in younger women and men (Less than 50 years old).

**Author Contributions:** Conceptualization, P.-Y.C.; methodology, P.-Y.C.; software, P.-Y.C.; validation, Y.-W.C. and Y.-F.L.; formal analysis, P.-Y.C.; investigation, P.-Y.C. and Y.-F.L.; writing—original draft preparation, P.-Y.C. and Y.-W.C.; writing—review and editing, P.-Y.C.; supervision, Y.-F.L. and H.-C.F.; funding acquisition, Y.-F.L. and H.-C.F. All authors have read and agree to the published version of the manuscript.

**Funding:** This study was funded by a grant from the Health Promotion Administration, Ministry of Health and Welfare, Taiwan, Republic of China (DOH101-HP-1103, DOH102-HP-1103, MOHW103- HPA-H-114-134101, MOHW104-HPA-H-114-134101); National Taipei University of Nursing and Health Sciences (109ntunhs-NT-02).

**Institutional Review Board Statement:** This study was approved by the Ethics Committees of Tri-Service General Hospital (TSGHIRB 100-05-197), Taipei Medical University (TMU-JIRB 20124036), Kaohsiung Medical University Chung-Ho Memorial Hospital (KMUHIRB 20120019), National Cheng Kung University Hospital (A-ER-101-117), Changhua Christian Hospital (CCHIRB 120405), Kaohsiung Chang Gung Memorial Hospital (101-1096B), Cardinal Tien Hospital (TMU-JIRB 201204035) and China Medical University Hospital (DMR101-IRB2-273(CR-1)). All research methods were implemented in accordance with guidelines approved by the Ethics Committees.

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** The data are not publicly available due to subject's privacy.

**Conflicts of Interest:** The authors declare that there is no duality of interest associated with this manuscript.

#### **References**


## *Article* **Poor Prognosis of Diffuse Large B-Cell Lymphoma with Hepatitis C Infection**

**Yu-Fen Tsai 1,2, Yi-Chang Liu 3,4, Ching-I Yang 3,5, Tzer-Ming Chuang <sup>3</sup> , Ya-Lun Ke <sup>3</sup> , Tsung-Jang Yeh <sup>3</sup> , Yuh-Ching Gau <sup>3</sup> , Jeng-Shiun Du <sup>3</sup> , Hui-Ching Wang 3,4 , Shih-Feng Cho 3,4, Chin-Mu Hsu <sup>3</sup> , Pey-Fang Wu <sup>6</sup> , Ching-I Huang 4,6, Chung-Feng Huang 4,6, Ming-Lung Yu 4,6 , Chia-Yen Dai 4,6 and Hui-Hua Hsiao 3,4,7,8,\***


**Abstract:** Background: Hepatitis C virus (HCV) in diffuse large B-cell lymphoma (DLBCL) is associated with a higher prevalence and distinctive clinical characteristics and outcomes. Methods: A retrospective analysis of adult DLBCL patients from 2011 to 2015 was studied. Results: A total of 206 adult DLBCL were enrolled with 22 (10.7%) HCV-positive patients. Compared to HCV-negative patients, the HCV-positive group had a poor performance status (*p* = 0.011), lower platelet count (*p* = 0.029), and higher spleen and liver involvement incidences (liver involvement, *p* = 0.027, spleen involvement, *p* = 0.026), and they received fewer cycles of chemotherapy significantly due to morbidity and mortality (*p* = 0.048). Overall survival was shorter in HCV-positive DLBCL (25.3 months in HCV-positive vs. not reached (NR), *p* = 0.049). With multivariate analysis, poor performance status (*p* < 0.001), advanced stage (*p* < 0.001), less chemotherapy cycles (*p* < 0.001), and the presence of liver toxicity (*p* = 0.001) contributed to poor OS in DLBCL. Among HCV-positive DLBCL, the severity of liver fibrosis was the main risk factor related to death. Conclusion: Inferior survival of HCV-positive DLBCL was observed and associated with poor performance status, higher numbers of complications, and intolerance of treatment, leading to fewer therapy. Therefore, anti-HCV therapy, such as direct-acting antiviral agents, might benefit these patients in the future.

**Keywords:** hepatitis C virus; diffuse large B-cell lymphoma; survival; fibrosis; performance status; liver toxicity

#### **1. Introduction**

Hepatitis C virus (HCV) infection is one of the major causes of chronic liver disease around the world with a prevalence rate of around 3% [1]. In addition to hepatic function, various extrahepatic manifestations, including vasculitis, glomerulonephrisis, and a wide range of lymphoproliferative disorders, have been associated with HCV infections,

**Citation:** Tsai, Y.-F.; Liu, Y.-C.; Yang, C.-I.; Chuang, T.-M.; Ke, Y.-L.; Yeh, T.-J.; Gau, Y.-C.; Du, J.-S.; Wang, H.-C.; Cho, S.-F.; et al. Poor Prognosis of Diffuse Large B-Cell Lymphoma with Hepatitis C Infection. *J. Pers. Med.* **2021**, *11*, 844. https://doi.org/ 10.3390/jpm11090844

Academic Editor: Stephen Opat

Received: 9 July 2021 Accepted: 25 August 2021 Published: 27 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

resulting in further morbidities and mortalities [2,3]. HCV infection are also associated with B-cell non-Hodgkin lymphoma (NHL) with distinct clinical presentation and outcome, especially in diffuse large B-cell lymphoma (DLBCL) and marginal zone lymphoma (MZL) in recent epidemiological studies [4–9]. Due to the fact that direct anti-HCV therapy induces hematologic response in indolent NHL patients with HCV infection, HCV infection has been linked to NHL, and front-line therapy of asymptomatic indolent NHLs is recommended [10,11]. DLBCL, the most common NHL, is an aggressive type of NHL with specific genetic features and clinical presentations [12]. Although HCV-positive DLB-CLs display specific clinical features and outcomes, the optimal treatment and timing of anti-HCV therapy in aggressive lymphoma is still uncertain [1,13–16]. Previous reports showed that HCV-positive patients with DLBCL exhibited worse overall survival (OS), and the incidence of severe hepatic toxicity in HCV-positive patients was significantly higher than that of HCV-negative patients [13,17]. However, all these studies only included small numbers of HCV-positive patients, and not all studies regarding HCV infection and DLBCL showed consistent results.

In addition, the seroprevalence of anti-HCV antigen in Taiwan is around 1 to 5.4% in the general population [18–20], and a retrospective study that investigated the association of HCV and lymphoma in southern Taiwan by Chuang et al. displayed that the incidence of HCV infection among lymphoma patients was significantly higher than that in healthy controls (11.0% vs. 1.8%, *p* < 0.001) [21]. The differences of HCV infection in NHL patients reflect geographic differences in HCV epidemiology [1,22,23]. The purpose of this study is to analyze clinicopathological characteristics, tolerance to chemotherapy, and clinical outcomes of DLBCL patients between HCV-positive and HCV-negative groups in order to determine the impact of HCV infection in this population.

#### **2. Materials and Methods**

#### *2.1. Study Population*

In this retrospective observational study, newly diagnosed DLBCL patients who were managed at Kaohsiung Medical University Hospital between January 2011 and December 2015 were enrolled. Clinical characteristics, including patients' demographics, baseline laboratory and biochemical data, hepatitis virus serology including hepatitis B virus (HBV) and HCV, and FIB (fibrosis)-4 index were collected and analyzed from medical records. Lymphoma was treated according to the treatment guideline in the hospital. The study was approved by the Institutional Review Board of Kaohsiung Medical University Hospital (KMUHIRB-E(I)-20210119).

#### *2.2. Risk Groups and Outcomes*

Patients were categorized as HCV-positive or HCV-negative groups based on the presence or absence of anti-HCV antibody. Since the interferon-free direct-acting antiviral agent (DAA) regimens were added to the reimbursement list in Taiwan for HCV treatment in 2017, patients were only followed until 2017 to eliminate the confounding factor of DAA treatment. All HCV-positive patients in this study did not receive interferon-free DAA prior to or during chemotherapy. All DLBCL patients with HBsAg-positive received antiviral agents for HBV prophylaxis. Performance status was evaluated by Eastern Cooperative Oncology Group (ECOG) scale [24]. The FIB-4 index was used to assess the severity of liver fibrosis as it is considered to be a useful fibrosis scoring system in HCV. A cutoff of >3.25 had a positive predictive value of 65% and a specificity of 97% to predict advanced fibrosis [25,26]. In Taiwan, liver cirrhosis was diagnosed by the image, either sonography or computed tomography or by FIB-4 index ≥ 6.5. The international prognostic index (IPI) score divided patients into four groups: low, low-intermediate, high-intermediate, and high for risk stratification [27], and the revised international prognostic index (R-IPI) classified patients into three groups: very good, good, and poor [28], which was calculated and used for analysis. The treatment response was evaluated by computed tomography (CT) or positron emission tomography–computed tomography (PET-CT) every 3–6 months

according to response criteria [29]. The first time of treatment response evaluation was at the time after 3–4 courses of chemoimmunotherapy. The liver toxicity was calculated by the common terminology criteria for adverse events (CTACE), version 4. Progressionfree survival (PFS) was defined as the length of time from the date of diagnosis to the date of progression or death. Overall survival (OS) was defined as the length of time from the date of diagnosis to the last date of follow-up or death. Liver function was monitored during every course of treatment during chemotherapy and around every 3 months after chemotherapy.

#### *2.3. Statistical Analysis*

Patient characteristics between HCV-positive and HCV-negative groups were analyzed using descriptive statistics and are presented as frequencies, percentages, and means with standard deviations. Student's t-test or nonparametric statistics were utilized to test for statistically significant differences in continuous variables, whereas the chi-square or Fisher's exact test was used for categorical variables. OS and PFS was estimated using the Kaplan–Meier method. Variables with a *p*-value of less than 0.05 in the univariate analysis of overall survival were subsequently subjected to multivariate analysis using a Cox regression model. A two-tailed *p*-value < 0.05 was considered statistically significant.

#### **3. Results**

#### *3.1. Baseline Characteristics*

A total of 206 DLBCL patients were enrolled in the study. Among them, twentytwo (10.7%) were HCV-positive. The characteristics at the time of diagnosis are listed in Table 1. Compared to HCV-negative patients, HCV-positive patients seemed to have a significantly lower platelet count (186.7 ± 68.8 × 103/µL vs. 236.2 ± 102.9 × 103/µL, *p* = 0.029) and a trend of older age (mean age 67.23 ± 17.13 vs. 61.24 ± 15.27, *p* = 0.088). HCV-negative patients had better performance status compared to HCV-positive patients (*p* = 0.011). Higher spleen and liver involvement rates were observed in HCV-positive patients compared to HCV-negative (liver involvement: 18.2% in HCV-positive vs. 4.3% in HCV-negative, *p* = 0.027, spleen involvement: 31.8% in HCV-positive vs. 13.6% in HCV-negative, *p* = 0.026). More HCV-positive patients had elevated aminotransferase compared to HCV-negative population, however, not to a statistical significant degree (31.8% vs. 16.3%, *p* = 0.073). There were no significant differences on sex, B symptoms, bone marrow involvement, IPI score, R-IPI score, HBsAg status, white blood cell (WBC) count, hemoglobin, albumin, lactate dehydrogenase (LDH) levels, alanine aminotransferase (GPT) levels, Beta-2 macroglobulin levels, and the percentage of receiving treatment between two groups.

#### *3.2. Response to Therapy and Outcomes*

Patients were managed by guidelines according to clinical status. Most (92.9%) of the HCV-negative patients and 86.4% of the HCV-positive patients received treatment respectively without statistical difference (*p* = 0.388) with most of them (74.3% in HCVnegative, 68.4% in HCV-positive) receiving an R-CHOP (rituximab, cyclophosphamide, daunorubicin, oncovin, and prednisolone) regimen. The details of the chemotherapy regimen were presented in Table 2, and there was no significant difference in treatment regimen between the HCV-positive and HCV-negative population. The HCV-positive population received fewer cycles of chemotherapy compared with HCV-negative patients (mean cycles of chemotherapy: 4.8 ± 3.5 in HCV-positive vs. 6.3 ± 3.4 in HCV-negative, *p* = 0.048) (Table 2). A lower evaluation of treatment response during the follow-up period was found in HCV-positive patients compared to the HCV-negative patient group (31.8% in HCV positive vs. 14.7% in HCV-negative, *p* = 0.041). Among these HCV-positive patients who were unable to be evaluated, four of seven patients received limited treatment with three of them receiving only one cycle of chemotherapy due to early mortality from sepsis

and/or other morbidities. Three of seven patients had supportive management due to old age (*n* = 2, over 75 years old) and severe liver cirrhosis (*n* = 1).


**Table 1.** Baseline characteristics of HCV-negative and HCV-positive DLBCL patients.

\*: *p* < 0.05, ECOG: Eastern Cooperative Oncology Group, R-IPI: revised international prognostic index, IPI: international prognostic index, SD: standard deviation, WBC: white blood cell, Hgb: hemoglobin, LDH: lactate dehydrogenase, GPT: alanine aminotransferase.


**Table 2.** Treatment detail, response, and liver toxicity in HCV-negative and HCV-positive DLBCL patients.

\* *p* < 0.05, R-CHOP: rituximab, cyclophosphamide, daunorubicin, oncovin and prednisolone; R-COP: rituximab, cyclophosphamide, oncovin and prednisolone; R-EPOCH: rituximab, etoposide, cyclophosphamide, daunorubicin, oncovin and prednisolone; R: rituximab.

The progression and recurrent events were similar in HCV-positive and HCV-negative groups in evaluable patients (27.3% in HCV-positive vs. 28.3% in HCV-negative, *p* = 0.922). Higher incidence of liver toxicity (86.4% in HCV-positive vs. 45.1% in HCV-negative, *p* < 0.001) and severe liver toxicity (≥grade 3 liver toxicity, 36.4% in HCV-positive vs. 9.8% in HCV-negative, *p* < 0.001) were observed in HCV-positive patients than HCV-negative patients. The incidence of liver toxicity was not associated with baseline GPT levels in HCV-positive patients (20.67 IU/mL, without liver toxicity vs. 37.16 IU/mL, with liver toxicity, *p* = 0.366). More patients in the HCV-positive group delayed or discontinued their treatment due to toxicities than in the HCV-negative group (40.9% vs. 12.5%, *p* = 0.001)

There was a trend of worse PFS in HCV-infected patients than non-infected, but no statistical significance (19.6 months in HCV-positive vs. NR in HCV-negative, *p* = 0.080), as shown in Figure 1. The 1-year, 2-year, and 3-year OS rates in HCV-positive and HCV-negative DLBCL patients were 67.0% vs. 76.4%, 50.2% vs. 71.6%, and 43.1% vs. 67.0%, respectively. The OS was worse in HCV-positive DLBCL patients (25.3 months in HCV-positive vs. not reached (NR) in HCV-negative, *p* = 0.049 by Kaplan–Meier method, Figure 2), which might result from early mortality, higher liver toxicity, and higher incidence of treatment delay/discontinuation with fewer treatment courses.

#### *3.3. Prognostic Factors for Overall Survival and Progression Free Survival in DLBCL Patients*

Based on univariate analysis via Cox regression model, we found that older age (HR: 1.03, 95% CI: 1.01–1.04, *p* = 0.003), poor performance status (ECOG ≥ 2, HR: 6.22, 95% CI: 3.66–10.56, *p* < 0.001), advanced stage (HR 4.36, 95% CI: 2.23–8.55, *p* <0.001), less chemotherapy cycles (HR: 0.77, 95% CI: 0.70–0.84, *p* < 0.001), and the presence of liver toxicity (HR: 2.46, 95% CI: 1.47–4.12, *p* = 0.001) contributed to poor OS. In multivariate analysis, four key factors: advanced stage, poor performance status, liver toxicity presentation, and less chemotherapy cycles were independent poor prognostic factors for overall survival in DLBCL patients (Table 3). Regarding the progression-free survival, we found that after multivariate analysis, poor performance status (ECOG ≥ 2, HR: 2.82, 95% CI: 1.56–5.09, *p* < 0.001), advanced stage (HR: 3.18, 95% CI: 1.81–5.58, *p* < 0.001), the presence of liver

toxicity (HR: 2.67, 95% CI: 1.69–4.22, *p* < 0.001), and less chemotherapy cycles (HR 0.90, 95% CI: 0.83–0.98, *p* = 0.015), contributed to poor PFS in DLBCL patients (Table 4). *J. Pers. Med.* **2021**, *11*, x FOR PEER REVIEW 6 of 13

*J. Pers. Med.* **2021**, *11*, x FOR PEER REVIEW 6 of 13

**Figure 1.** Progression-free survival of HCV-positive and HCV-negative DLBCL patients. NR: not reached, HR: hazard ratio, CI: confidence interval. **Figure 1.** Progression-free survival of HCV-positive and HCV-negative DLBCL patients. NR: not reached, HR: hazard ratio, CI: confidence interval. **Figure 1.** Progression-free survival of HCV-positive and HCV-negative DLBCL patients. NR: not reached, HR: hazard ratio, CI: confidence interval.

**Figure 2.** Overall survival of HCV-positive and HCV-negative DLBCL patients. NR: not reached, **Figure 2.** Overall survival of HCV-positive and HCV-negative DLBCL patients. NR: not reached, HR: hazard ration, CI: confidence interval. **Figure 2.** Overall survival of HCV-positive and HCV-negative DLBCL patients. NR: not reached, HR: hazard ration, CI: confidence interval.

*3.3. Prognostic Factors for Overall Survival and Progression Free Survival in DLBCL Patients* 

*3.3. Prognostic Factors for Overall Survival and Progression Free Survival in DLBCL Patients* 

Based on univariate analysis via Cox regression model, we found that older age (HR: 1.03, 95% CI: 1.01–1.04, *p* = 0.003), poor performance status (ECOG ≥ 2, HR: 6.22, 95% CI: 3.66–10.56, *p* < 0.001), advanced stage (HR 4.36, 95% CI: 2.23–8.55, *p* <0.001), less chemotherapy cycles (HR: 0.77, 95% CI: 0.70–0.84, *p* < 0.001), and the presence of liver toxicity

Based on univariate analysis via Cox regression model, we found that older age (HR: 1.03, 95% CI: 1.01–1.04, *p* = 0.003), poor performance status (ECOG ≥ 2, HR: 6.22, 95% CI: 3.66–10.56, *p* < 0.001), advanced stage (HR 4.36, 95% CI: 2.23–8.55, *p* <0.001), less chemotherapy cycles (HR: 0.77, 95% CI: 0.70–0.84, *p* < 0.001), and the presence of liver toxicity (HR: 2.46, 95% CI: 1.47–4.12, *p* = 0.001) contributed to poor OS. In multivariate analysis,

HR: hazard ration, CI: confidence interval.


**Table 3.** Univariate and multivariate analysis of risk factors for overall survival.

\*: *p* < 0.05, HR: hazard ratio, CI: confidence interval, ECOG: Eastern Cooperative Oncology Group. <sup>a</sup> : Multivariate analyses: variables with *p* < 0.1 in univariate analyses were included in multivariate analyses.

**Table 4.** Univariate and multivariate analysis of risk factors for progression-free survival.


\*: *p* < 0.05, HR: hazard ratio, CI: confidence interval, ECOG: Eastern Cooperative Oncology Group <sup>a</sup> : Multivariate analyses: variables with *p* < 0.1 in univariate analyses were included in multivariate analyses.

#### *3.4. Details of HCV-Positive DLBCL Patients during Follow-Up Period*

The details of HCV-positive DLBCL patients are listed in Table S1. Half of HCVpositive (11 of 22) DLBCL patients died during the follow-up period. Seven died of sepsis, two died due to the progression of hepatocellular carcinoma, one died of gastrointestinal bleeding related to liver cirrhosis, and only one patient died of lymphoma progression. Based on subgroup analyses of the death events in the HCV-positive population, no

difference of initial performance status, stage, IPI score, platelet count, coagulation profile (prothrombin time, PT), albumin, GPT, and LDH levels was found between the dead and alive groups. However, less chemotherapy cycles, higher percentage of severe liver toxicity, and treatment delay/discontinuation were observed in the death population, but not to a statistically significant degree. We found that only liver cirrhosis status contributed to more death events in HCV-positive patients. In HCV-positive patients, all were alive of those without liver cirrhosis (*p* = 0.035) or FIB-4 index less than 1.45 points (*p* = 0.032) (Table 5). In addition, only one of 22 HCV-positive patients started peg-interferon and ribavirin treatment for HCV at 4 months after the end of chemotherapy, and HCV viral load was 2758.2 IU/mL at that time. The HCV viral load became undetectable after peg-interferon and ribavirin status to now.


**Table 5.** Different variables in the dead and alive groups of HCV-positive DLBCL patients.

\*: *p* < 0.05, SD: standard deviation, ECOG: Eastern Cooperative Oncology Group, IPI: international prognostic index GPT: alanine aminotransferase, PLT: platelet, PT: prothrombin time, INR: international normalized ratio, LDH: lactate dehydrogenase, FIB-4: fibrosis-4.

#### **4. Discussion**

In our study, the prevalence of HCV infection was 10.7% (22/206 patients) in DLBCL patients. It was higher than the prevalence rate of the general population in Taiwan at 4.4% (6904/157,720) of anti-HCV positive [30]. There have been some hypothetical models to describe the possible pathologic role of HCV infection in aggressive B-cell lymphoma. The direct transformation mechanism is one kind of hypothesis accounting for HCV-associated lymphomagenesis. In vitro, the infection of B-cell lines with HCV leads to somatic mutations of several oncogenes and tumor-suppressor genes such as p53, betacatenin, and Bcl6 [31]. On the contrary, some studies support the role of HCV as an indirect transformation agent by chronically stimulating B-cell immunologic response and finally leading to lymphoma. Polyclonal or monoclonal B-cell proliferation can be detected in the blood or bone marrow of HCV patients and is associated with mixed cryoglobulinemia type II. The risk of lymphoma in patients with HCV-associated cryoglobulinemia is estimated to be 35 times higher than that in the general population [32]. All these theories are proposed to explain the high prevalence of HCV infection in lymphoma patients.

In previous retrospective studies, HCV-infected DLBCL patients shared distinctive clinical features. A higher percentage of HCV-positive DLBCL cases were associated with old age (> or =60) than HCV-negative DLBCL cases at diagnosis [14,33]. HCV-positive patients had more frequent extra-nodal involvement [13,33] and more frequent elevated lactate dehydrogenase (LDH) levels than other patients [13]. In our studies, we found that HCV-positive DLBCL patients had poorer performance status, lower platelet count, and higher probability of liver and spleen involvement than HCV-negative patients. In terms of OS and PFS, the influence of HCV infection on survival in DLBCL patients remains controversial. Studies conducted by Ennishi et al. reported a similar outcome in HCV-positive DLBCL patients compared to HCV-negative (3-year OS 75% in HCV-positive vs. 84% in HCV-negative, *p* = 0.07) [17]. On the contrary, two other studies presented an inferior outcome in HCV-positive patients [16,34]. Chen et al. found that three independent factors predicted a dismal OS: lower albumin level (<3 g/dL vs. ≥3 g/dL; HR = 13.21, 95% CI = 2.69–64.98, *p* = 0.001), presence of HCV infection (HCV-positive vs. HCV-negative; HR = 9.75, 95% CI = 1.97–48.34, *p* = 0.005), and poor IPI risk (high-intermediate or high vs. low-intermediate or low; HR = 5.56, 95% CI = 1.17–26.55, *p* = 0.031) [34]. The Fondazione Italiana Linfomi has conducted a multicenter retrospective study to explore a new prognostic system for HCV-associated DLBCL. Therefore, the HCV Prognostic Score, based on performance status, albumin level, and HCV-RNA viral load, was therefore introduced as a useful tool to predict the outcome of HCV-associated DLBCL [16].

In our study, we found that an inferior overall survival was observed in HCV-positive patients by the Kaplan–Meier method. However, after multivariate analysis, HCV infection status was not the key risk factor for overall survival. Instead, advanced stage, poor performance status, liver toxicity presentation, and less chemotherapy cycles were the main risk factors that predicted overall survival and progression-free survival. A higher percentage of HCV-positive patients' treatment response was unable to be evaluated during the follow-up period (31.8% in HCV-positive vs. 14.7% in HCV-negative, *p* = 0.041). Nearly half (42.9%, three of seven) of them did not receive any treatment due to elderly age with poor performance status and severe liver cirrhosis; the rest of the patients received treatment of no more than two cycles with no possible follow-up image due to sepsisrelated death. This highlighted the ability of HCV co-morbidity to hamper treatment completeness. The number of chemotherapy cycles was a major predictive factor for OS and HCV-positive patients; fewer chemotherapy cycles resulted in poor survival in the HCV-positive group. The phenomenon was also observed by Dlouhy et al., who presented higher number of complications and early deaths in HCV-positive DLBCL patients [33].

Regarding the HCV-positive DLBCL population, we found that liver cirrhosis was the main factor related to death events. Unlike the whole DLBCL population, age, stage, and performance status carried no vital role in predicting survival in HCV-positive DLBCL patients. Although fewer chemotherapy cycles, higher percentage of severe liver toxicity, and treatment delay/discontinuation were observed among those patients that died, the difference was not statistically significant compared to those who remained alive. While these findings highlighted the influence of the severity of liver cirrhosis for HCV-positive

DLBCL patients, the small sample size might limit the statistical findings. This result highlighted that the severity of cirrhosis was the major factor to determine survival in HCV-positive DLBCL patients and also contributed to high complications and early deaths in HCV-positive DLBCL patients.

We found higher liver toxicity and higher ≥ grade 3 liver toxicity in HCV-positive population, and the incidence of liver toxicity was not associated with initial GPT levels. One patient (4.5%, one of 22 HCV-positive) received peg-interferon and ribavirin treatment for high HCV RNA levels after chemoimmunotherapy. Since we did not have the comprehensive data of HCV viral load, it was difficult to be conclude as a HCV reactivation in our study. The increased risk of liver toxicity in the HCV-positive population was consistent with other studies. In a Japanese multicenter retrospective study of 553 DLBCL patients, HCV RNA levels significantly increased in the HCV-positive patients during immunochemotherapy, and more grade 3–4 hepatic toxicity (27%) was observed in patients with HCV-positive patients than HCV-negative patients (3%) [17]. In another Italian retrospective study of 156 HCV-positive DLBCL patients, they found that 85% of patients completed their therapeutic program without any interruption or dosage reduction and five patients (4%) had to discontinue chemotherapy due to severe hepatic function impairment (toxicity grade 3–4) [35]. These studies revealed that without initial liver dysfunction, HCV-positive patients experience a similar outcome compared to HCV-negative patients when treated with chemotherapy or immunotherapy [36]. Therefore, careful monitoring of hepatic function and viral load is suggested during the therapy, especially for those with initial liver dysfunction patients.

Some studies showed that DAA treatment could reverse the liver inflammation, fibrosis, and decrease the complications of compensated liver cirrhosis in chronic hepatitis C [37–39]. From the observation report of Occhipinti et al., seven DLBCL patients with HCV infection received different DAA regimens concurrently with immunochemotherapy [40]. All patients completed their scheduled treatment with no liver toxicity occurrence. Thus, the concomitant administration of DAA and immunochemotherapy for HCV-positive DLBCL seems safe and may prevent immunochemotherapy-induced liver toxicity. Further prospective studies are warranted to identify if concurrent DAA therapy with chemoimmunotherapy in HCV-positive DLBCL patients can lead to survival benefit, especially for those with advanced fibrosis.

There were some limitations in this study. First, our study was a retrospective study with limited numbers of HCV patients, which limits the scope of our conclusion. Secondly, HCV viral load was not routinely checked during 2011–2015. Therefore, we did not have the comprehensive data of HCV viral load in all HCV-positive patients. Thus, we only used the anti-HCV positivity as criteria for being HCV-positive. Thirdly, some patients had a shortened course of therapy due to mortalities, which hinders proper response evaluations, especially in HCV-positive patients. However, this study presents real-world data of clinical characteristics and outcomes of HCV-positive DLBCL patients before the introduction of DAA treatment and provided the association of HCV and DLBCL.

#### **5. Conclusions**

In our study, DLBCL patients have a higher prevalence (10.7%) of HCV infection than the general population in our area. Clinically, HCV-positive DLBCL patients had poor performance status, lower platelet count, higher probability of spleen and liver involvement, and received less chemotherapy cycles due to morbidity and mortality. With multivariate analysis, stage, performance status, liver toxicity presentation, and chemotherapy cycles were found to be the major factors for the OS and PFS in DLBCL patients. Inferior survival of HCV-positive DLBCL was contributed by poor performance status, higher numbers of complications, and intolerance of treatment, leading to fewer therapy. Severity of liver fibrosis was the main factor related to death in the HCV-positive population. As such, anti-HCV therapy such as direct-acting antiviral agents might benefit such patients in the future to improve outcomes.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/jpm11090844/s1, Table S1: the details of HCV-positive DLBCL patients.

**Author Contributions:** Conceptualization, Y.-F.T., H.-H.H., Y.-C.L. and C.-I.Y.; methodology, Y.-F.T., C.-I.Y. and P.-F.W.; validation, C.-I.H., C.-F.H., M.-L.Y. and C.-Y.D.; resources, T.-J.Y., Y.-C.G., J.-S.D., H.-C.W., S.-F.C., Y.-F.T., H.-H.H. and Y.-C.L.; data curation, T.-M.C. and Y.-L.K.; writing—original draft preparation, Y.-F.T. and H.-H.H.; writing—review and editing, Y.-F.T., H.-H.H. and C.-M.H.; supervision, P.-F.W., C.-I.H., C.-F.H., M.-L.Y. and C.-Y.D.; project administration, Y.-F.T. and H.-H.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by a grant from the Kaohsiung Medical University Hospital (KMUH108-8R24 and KMUH109-9R19) and Ministry of Health and Welfare (MOHW110-TDU-B-212- 124006).

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Institutional Review Board of Kaohsiung Medical University Hospital (KMUHIRB-E(I)-20210119).

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data used in the present study are available from the corresponding author upon reasonable request.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Association of Ventilatory Disorders with Respiratory Symptoms, Physical Activity, and Quality of Life in Subjects with Prior Tuberculosis: A National Database Study in Korea**

**Bumhee Yang 1,†, Hayoung Choi 2,† , Sun Hye Shin <sup>3</sup> , Youlim Kim 4,5, Ji-Yong Moon <sup>6</sup> , Hye Yun Park 3,\* ,‡ and Hyun Lee 6,\* ,‡**


**Abstract:** Tuberculosis (TB) survivors experience post-TB lung damage and ventilatory function disorders. However, the proportions of obstructive and restrictive ventilatory disorders as well as normal ventilation among subjects with prior TB are unknown. In addition, the impacts of ventilatory disorder and its severity on respiratory symptoms, physical activity limitations, and the quality of life in subjects with prior TB remain unclear. Subjects who participated in the Korean National Health and Nutritional Examination Survey 2007–2016 were enrolled in this study. We evaluated the impact of each ventilatory disorder and its severity on respiratory symptoms, physical activity limitations, and quality of life (measured by the EuroQoL five dimensions questionnaire [EQ-5D] index values) in subjects with prior TB. Among 1466 subjects with prior TB, 29% and 16% had obstructive ventilatory disorders and restrictive ventilatory disorders, respectively. Mild and moderate obstructive ventilatory disorders were not associated with respiratory symptoms, physical activity limitations, or EQ-5D index value compared with normal ventilation; however, severe obstructive ventilatory disorders were associated with more respiratory symptoms (adjusted odds ratio [aOR] = 13.62, 95% confidence interval [CI] = 4.64–39.99), more physical activity limitation (aOR = 218.58, 95% CI = 26.82–1781.12), and decreased EQ-5D index (adjusted coefficient = −0.06, 95% CI = (−0.12–−0.10) compared with normal ventilation. Mild restrictive ventilatory disorders were associated with more respiratory symptoms (aOR = 2.10, 95% CI = 1.07–4.14) compared with normal ventilation, while moderate (aOR = 5.71, 95% CI = 1.14–28.62) and severe restrictive ventilatory disorders (aOR = 9.17, 95% CI = 1.02–82.22) were associated with physical activity limitation compared with normal ventilation. In conclusion, among subjects with prior TB, 29% and 16% developed obstructive and restrictive ventilatory disorders, respectively. Severe obstructive ventilatory disorder was associated with more respiratory symptoms, more physical activity limitation, and poorer quality of life, while severe restrictive ventilatory disorder was associated with more physical activity limitations.

**Keywords:** tuberculosis; pulmonary function; respiratory symptoms; quality of life

**Citation:** Yang, B.; Choi, H.; Shin, S.H.; Kim, Y.; Moon, J.-Y.; Park, H.Y.; Lee, H. Association of Ventilatory Disorders with Respiratory Symptoms, Physical Activity, and Quality of Life in Subjects with Prior Tuberculosis: A National Database Study in Korea. *J. Pers. Med.* **2021**, *11*, 678. https:// doi.org/10.3390/jpm11070678

Academic Editors: Rutger A. Middelburg and José Carmelo Adsuar Sala

Received: 9 June 2021 Accepted: 16 July 2021 Published: 19 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Tuberculosis (TB) is the leading cause of death from a single epidemic disease worldwide [1]. Since recognition of TB as a deadly communicable disease, the mortality rate has decreased, and 85% of patients treated for their first episode of TB survive with treatment success (cure or completion) [1,2]. Despite the improved survival rate, approximately 50% of TB survivors experience post-TB lung damage and ventilatory disorders [3–5]. Depending on the severity of post-TB lung damage, they show various clinical courses, including no respiratory symptoms, dyspnea, or impaired quality of life (QoL) [6,7].

The ventilatory disorders of subjects with prior TB comprise obstructive and restrictive patterns [8]. Obstructive ventilatory disorder, which appears as airflow limitation on pulmonary function tests, is the most well-known form of lung damage after TB treatment [5,8]. Post-TB lung damage, including cavity, bronchiectasis, or distorted airway, can cause obstructive ventilatory disorder, which can lead to dyspnea, chronic obstructive pulmonary disease, and reduced exercise capacity [9,10]. Although obstructive disorder has been recognized as post-TB lung damage, restrictive ventilatory disorder is also found in subjects with prior TB and can lead to dyspnea or chest pain [11]. The restrictive pattern has been suggested to be a consequence of excessive fibrosis, fibrotic bands, or bronchovascular distortion in the process of post-TB lung repair [12]. The proportions of obstructive and restrictive ventilatory disorders as well as normal ventilation among subjects with prior TB remain unclear. Additionally, the impacts of ventilatory defects and severity of defects on respiratory symptoms, physical limitations, and QoL in subjects with prior TB have not been investigated.

This study sought to determine the composition of ventilatory defects among subjects with prior TB using a nationally representative database in South Korea. Furthermore, this study investigated the impact of each ventilatory disorder and its severity on respiratory symptoms, physical activity limitations, and QoL among subjects with prior TB.

#### **2. Material and Methods**

#### *2.1. Study Population*

We used the data from the 2007–2016 Korea National Health and Nutrition Examination Survey (NHANES), a nationally representative health survey collected by the Korean Ministry of Health and Welfare. Health-related questionnaires, health examinations, and spirometry results were used in this study. Previous pulmonary TB was defined based on formal reading of a chest X-ray or a history of physician diagnosed pulmonary TB. We classified subjects with prior TB into three groups according to spirometric pattern (Figure 1). The study protocol was approved by the Institutional Review Board of Chungbuk National University Hospital (application no. 2021-01-041).

#### *2.2. Measurements*

sis [16].

*2.3. Definitions of Ventilatory Disorder* 

*2.2. Measurements*  Data on age, sex, body mass index (BMI), smoking history, physical activity limitations, occupation, EuroQoL five dimensions questionnaire (EQ-5D) index value, and spirometric results were obtained from the Korea NHANES database. Respiratory symptoms, including cough, sputum, and dyspnea, were measured qualitatively (presence or Data on age, sex, body mass index (BMI), smoking history, physical activity limitations, occupation, EuroQoL five dimensions questionnaire (EQ-5D) index value, and spirometric results were obtained from the Korea NHANES database. Respiratory symptoms, including cough, sputum, and dyspnea, were measured qualitatively (presence or absence). To assess physical activity limitations, we used a questionnaire, "Do you experience physical activity limitations due to respiratory disease?", and this was also measured qualitatively (yes or no).

absence). To assess physical activity limitations, we used a questionnaire, "Do you experience physical activity limitations due to respiratory disease?", and this was also measured qualitatively (yes or no). The EQ-5D enables the respondent to classify his or her health according to five dimensions. These dimensions define health in terms of mobility, self-care, usual activity, pain/discomfort, and anxiety/depression. Each dimension is divided into three levels, i.e., The EQ-5D enables the respondent to classify his or her health according to five dimensions. These dimensions define health in terms of mobility, self-care, usual activity, pain/discomfort, and anxiety/depression. Each dimension is divided into three levels, i.e., no problem/some or moderate problems/extreme problems. The information derived from the EQ-5D self-classifier can be converted into a single summary index: the EQ-5D index [13]. The EQ-5D index ranges between 0 (worst imaginable health state) and 1 (best imaginable health state).

no problem/some or moderate problems/extreme problems. The information derived from the EQ-5D self-classifier can be converted into a single summary index: the EQ-5D index [13]. The EQ-5D index ranges between 0 (worst imaginable health state) and 1 (best imaginable health state). Spirometry was performed according to the recommendations of the American Tho-Spirometry was performed according to the recommendations of the American Thoracic Society and European Respiratory Society [14]. Since post-bronchodilator spirometry was not available in Korea NHANES database, pre-bronchodilator spirometry results were used in our study. Absolute values of forced expiratory volume in 1 s (FEV1) and forced vital capacity (FVC) were obtained, and the percentage of predicted values (% predicted) for FEV<sup>1</sup> and FVC were calculated using the reference equation obtained on analysis of a representative Korean sample [15]. Comorbidities of asthma, diabetes mellitus, hyper-

racic Society and European Respiratory Society [14]. Since post-bronchodilator spirometry was not available in Korea NHANES database, pre-bronchodilator spirometry results

forced vital capacity (FVC) were obtained, and the percentage of predicted values (% predicted) for FEV1 and FVC were calculated using the reference equation obtained on analysis of a representative Korean sample [15]. Comorbidities of asthma, diabetes mellitus, hypertension, dyslipidemia, cardiovascular disease, osteoporosis, osteoarthritis or rheumatoid arthritis, and depression were self-reported based on previous physician diagnotension, dyslipidemia, cardiovascular disease, osteoporosis, osteoarthritis or rheumatoid arthritis, and depression were self-reported based on previous physician diagnosis [16].

#### *2.3. Definitions of Ventilatory Disorder*

Normal ventilation was defined as pre-bronchodilator FEV1/FVC ≥ 0.70 and FVC ≥ 80% predicted. Obstructive ventilatory disorder was defined as pre-bronchodilator FEV1/FVC < 0.70 [17]. For cases with obstructive ventilatory disorder, FEV<sup>1</sup> ≥ 80% predicted, FEV<sup>1</sup> of 50–79% predicted, and FEV<sup>1</sup> < 50% predicted were classified as mild, moderate, and severe, respectively. Restrictive ventilatory disorder was defined as FEV1/FVC ≥ 0.7 and FVC < 80% predicted. For cases with restrictive ventilatory disorder, FVC ≥ 70% predicted, FVC of 60–69% predicted, and FVC < 60% predicted were classified as mild, moderate, and severe, respectively [18,19].

#### *2.4. Outcomes*

We compared respiratory symptoms, physical activity limitations due to respiratory diseases (hereafter physical activity limitations), and QoL (measured using the EQ-5D index) between subjects with prior TB with different ventilatory disorders. We also analyzed the impacts of ventilatory disorder severity on respiratory symptoms, physical activity limitations, and QoL in subjects with prior TB.

#### *2.5. Statistical Analysis*

All analysis was performed using survey commands in STATA 15.1 version (StataCorp LP, College Station, TX, USA) to account for the complex sampling design and survey weights. For each variable, we calculated prevalence and 95% confidence interval (CI) by group.

The associations between ventilatory disorders and respiratory symptoms (cough, sputum, or dyspnea) and physical activity limitations were analyzed using logistic regression analysis: Multivariable analysis was adjusted for age (categorized as ≥65 years or not), sex, BMI, smoking amount (pack-years), education level (categorized as high school or less vs. college or above), and family income (categorized as low vs. high). Linear regression analysis was performed to assess the association between ventilatory disorders and EQ-5D index scores: Multivariable analysis was adjusted for covariates as mentioned above.

A two-sided *p* value < 0.05 was considered significant. To account for multiple comparisons, post hoc Bonferroni correction was applied (normal ventilation vs. obstructive ventilatory disorder, normal ventilation vs. restrictive ventilatory disorder, and obstructive ventilatory disorder vs. restrictive ventilatory disorder), in which a *p* value of 0.05 corresponds to 0.017 (0.05/3).

#### **3. Results**

#### *3.1. Baseline Characteristics*

As shown in Table 1, the post-TB patient group included 1466 patients (54.9%) with normal ventilation, 783 patients (29.3%) with obstructive ventilatory disorders, and 420 patients (15.8%) with restrictive ventilatory disorders. The mean ages of subjects with normal ventilation, obstructive ventilatory disorders, and restrictive ventilatory disorders were 53.4, 64.4, and 59.6 years, respectively (*p* < 0.001). The proportion of males was highest in subjects with obstructive ventilatory disorders, followed by those with restrictive ventilatory disorders and those with normal ventilation (76.2%, 52.0%, and 49.5%, respectively, *p* < 0.001). The subjects with obstructive ventilatory disorders had a lower BMI than those with normal ventilation or restrictive ventilatory disorders (22.8 kg/m<sup>2</sup> , 23.7 kg/m<sup>2</sup> , and 23.8 kg/m<sup>2</sup> , respectively, *p* < 0.001). The subjects with obstructive ventilatory disorders had higher rates of current or ex-smoking than those with normal ventilation or restrictive ventilatory disorders (71.8%, 44.5%, and 39.4%, respectively, *p* < 0.001). The proportion of subjects with low family income was highest among obstructive ventilatory disorders, followed by restrictive ventilatory disorders, and normal ventilation (68.3%, 55.9%, and

43.0%, respectively, *p* < 0.001). The subjects with obstructive ventilatory disorders had the highest prevalence of asthma (10.9%), diabetes mellitus (19.1%), and hypertension (50.8%) among subjects with prior TB, while subjects with restrictive ventilatory disorders had the highest prevalence of dyslipidemia (48.4%) and osteoporosis (11.9%). Regarding spirometric results, FVC (L) and FVC % predicted were lowest in subjects with restrictive ventilatory disorders, while FEV<sup>1</sup> (L), FEV<sup>1</sup> % predicted, and FEV1/FVC were lowest in those with obstructive ventilatory disorders among subjects with prior TB (*p* < 0.001 for all).


**Table 1.** Baseline characteristics of subjects with prior TB by spirometric pattern.

Data are presented as weighted mean (95% confidence interval) or weighted percentage (95% confidence interval). *p* values are comparisons of three groups. TB, tuberculosis; BMI, body mass index; FVC, forced vital capacity; FEV1, forced expiratory volume in 1 s. The results of Bonferroni correction with three comparisons are provided as superscripts in Table (a *p* value of 0.05 corresponds to 0.017 [0.05/3]). a Indicates statistical significance for the comparison of normal ventilation and obstructive ventilatory disorder. <sup>b</sup> Indicates statistical significance for the comparison of normal ventilation and restrictive ventilatory disorder. <sup>c</sup> Indicates statistical significance for the comparison of obstructive ventilatory disorder and restrictive ventilatory disorder.

#### *3.2. Comparison of Symptoms, Physical Activity, and Quality of Life*

As shown in Table 2, respiratory symptoms including sputum (18.2%, 15.1%, and 9.5%, respectively, *p* = 0.004) and dyspnea (3.8%, 2.4%, and 0.9%, respectively, *p* < 0.020) and physical activity limitations (27.7%, 13.0%, and 5.1%, respectively, *p* < 0.001) were most frequently observed in subjects with obstructive ventilatory disorders, followed by those with restrictive ventilatory disorders, and those with normal ventilation. Cough was observed most frequently in subjects with obstructive ventilatory disorders followed by

those with normal ventilation, and those with restrictive ventilatory disorders (11.6%, 5.4%, and 5.1%, respectively, *p* < 0.001).


**Table 2.** Comparison of symptoms, physical activity, and quality of life according to spirometric pattern.

Data are presented as weighted mean (95% confidence interval) or weighted percentage (95% confidence interval). *p* values are comparisons of three groups. EQ-5D, EuroQoL five dimensions. The results of Bonferroni correction with three comparisons are provided as superscripts in Table (a *p* value of 0.05 corresponds to 0.017 [0.05/3]). <sup>a</sup> Indicates statistical significance for the comparison of normal ventilation and obstructive ventilatory disorder. <sup>b</sup> Indicates statistical significance for the comparison of normal ventilation and restrictive ventilatory disorder. <sup>c</sup> Indicates statistical significance for the comparison of obstructive ventilatory disorder and restrictive ventilatory disorder.

> The EQ-5D index values, denoting QoL, were lower among subjects with obstructive ventilatory disorders and restrictive ventilatory disorders than in those with normal ventilation (0.91, 0.91, and 0.94, respectively, *p* = 0.002). Regarding the individual EQ-5D component arm, the subjects with obstructive ventilatory disorders had the highest rates of mobility difficulty (24.1%, *p* < 0.001) and self-care limitations (7.7%, *p* = 0.005) among subjects with prior TB; however, those with restrictive ventilatory disorders had the highest rates of difficulty with usual activities (19.1%, *p* < 0.001) and pain/discomfort (36.3%, *p* = 0.015) among subjects with prior TB (Table 2).

#### *3.3. The Impact of Obstructive Ventilatory Disorder and Its Severity on Respiratory Symptoms, Physical Activity Limitations, and EQ-5D Index in Subjects with Prior TB*

As shown in Supplementary Table S1, the subjects with obstructive ventilatory disorders were 1.63 (95% CI = 1.05–2.82) times more likely to have any respiratory symptoms compared to those with normal ventilation in the fully adjusted model; additionally, obstructive ventilatory disorders were closely associated with sputum (aOR = 1.85, 95% CI = 1.03–3.37) and dyspnea (aOR = 4.19, 95% CI = 1.03–17.14). Likewise, the subjects with obstructive ventilatory disorder had more physical activity limitations (aOR = 6.59, 95% CI = 1.98–21.93) compared to those with normal ventilation. Although subjects with obstructive ventilatory disorders were more likely to have a lower EQ-5D index compared to those with normal ventilation in the unadjusted model (coefficient = −0.02, 95% CI −0.04–−0.07), this was not significant in the adjusted model.

In the analyses according to obstructive ventilatory disorder severity, mild and moderate obstructive ventilatory disorders were not associated with more respiratory symptoms, more physical activity limitations, or lower EQ-5D index values in the unadjusted and adjusted models. However, severe obstructive ventilatory disorder was significantly associated with more respiratory symptoms (aOR = 13.62, 95% CI = 4.64–39.99), more physical activity limitations (aOR = 218.58, 95% CI = 26.82–1781.12), and lower EQ-5D index values (adjusted coefficient = −0.06, 95% CI = −0.12–−0.01) (Table 3).


**Table 3.** The impact of obstructive ventilatory disorder severity on respiratory symptoms, physical activity limitations, and EQ-5D index value in subjects with prior TB.

Data are presented as a ratio (95% confidence interval) or a difference estimate (95% confidence interval). \* Age, sex, body mass index, smoking amount (pack-years), education (categorized as >high school or ≤high school), and family income (categorized as low or high) were adjusted. EQ-5D, EuroQoL five dimensions; TB, tuberculosis.

#### *3.4. The Impact of Restrictive Ventilatory Disorder and Its Severity on Respiratory Symptoms, Physical Activity Limitations, and EQ-5D Index Value in Subjects with Prior TB*

Restrictive ventilatory disorder was not significantly associated with increased respiratory symptoms (cough, sputum, or dyspnea), physical activity limitations, and EQ-5D index values compared with normal ventilation in adjusted models (Supplementary Table S1).

In the analyses according to restrictive ventilatory disorder severity, the subjects with mild restrictive ventilatory disorders were more likely to have respiratory symptoms (aOR = 2.10, 95% CI = 1.07–4.14) compared to those with normal ventilation, specifically sputum (aOR = 2.80, 95% CI = 1.34–5.83). In comparison, moderate (aOR = 5.71, 95% CI = 1.14–28.62) and severe (aOR = 9.17, 95% CI = 1.02–82.22) restrictive ventilatory disorders were associated with more physical activity limitations compared to those with normal ventilation (Table 4).

**Table 4.** The impact of restrictive ventilatory disorder severity on respiratory symptoms, physical activity limitations, and EQ-5D index value in subjects with prior TB.


Data are presented as a ratio (95% confidence interval) or a difference estimate (95% confidence interval). \* Adjusted for age, sex, body mass index, smoking amount (pack-years), education (categorized as >high school or ≤high school), and family income (categorized as low or high). EQ-5D, EuroQoL five dimensions; TB, tuberculosis.

#### **4. Discussion**

To the best of our knowledge, this is the first study evaluating respiratory symptoms, physical activity limitations, and QoL according to ventilatory disorder type and severity in subjects with prior TB. Among subjects with prior TB, approximately 29% and 16% developed obstructive and restrictive ventilatory disorders, respectively. Severe obstructive ventilatory disorders were associated with more respiratory symptoms, more physical activity limitations, and poorer quality of life. Severe restrictive ventilatory disorder was associated with more physical activity limitations.

TB survivors frequently experience structural and functional lung sequelae that vary in severity [7]. For example, approximately 24–35% of TB survivors experience obstructive ventilatory disorders [20–22]. In agreement with previous reports, 29% of the subjects with prior TB in this study had obstructive ventilatory disorders. Thus, development of obstructive ventilatory disorders can cause a concerning health-related burden in TB survivors. Subjects with prior TB also can experience restrictive ventilatory disorders. Post-TB survivors often show a fibrotic pattern on chest imaging due to the sequelae of pulmonary TB including destruction of lung parenchyma [23]. Restrictive ventilatory disorders occur in subjects with prior TB due to volume loss, lung scarring with reduction of pulmonary compliance, and an increase in elastic retraction pressure [24,25]. In contrast to the literature on obstructive ventilatory disorders, only a few studies have described restrictive ventilatory disorders, and the prevalence was reported to be 31–42% among TB survivors [26,27]. The small number of patients in these studies (*n* = 107 and *n* = 33) limited the study findings [26,27]. Thus, our study had the advantage of confirming these findings with a larger number of subjects using a nationwide database.

As shown in previous studies [20,28], clinical factors associated with poorer QoL (e.g., old age, male sex, smoking history, lower BMI, and lower education level) were more common in subjects with prior TB with obstructive ventilatory disorders than in those with normal ventilation. Approximately 72% of subjects with prior TB with obstructive ventilatory disorders in this study were current or ex-smokers. In line with this finding, smoking is a well-established factor associated with obstructive ventilatory disorders in subjects with prior TB [20]. However, little is known as to whether development of obstructive ventilatory disorders associated with higher symptomatic burdens in post-TB survivors, as most previous studies focused on the presence of obstructive ventilatory disorders and their severity after TB treatment [5,11,29]. From this perspective, our study has confirmed that the development of severe obstructive ventilatory disorders is associated with more respiratory symptoms, more physical activity limitations, and poorer QoL compared to patients with normal ventilation; however, these findings were not significant in patients with mild-to-moderate obstructive ventilatory disorders. These results indicate that regular health check-ups with pulmonary function measurement is necessary after completion of TB treatment to detect obstructive ventilatory disorders early and provide appropriate treatment to prevent further lung function impairment. Recent studies also support this suggestion in showing clinical improvement after bronchodilator treatment in TB-destroyed lung patients with obstructive ventilatory disorder [30,31].

Despite the prevalence of restrictive ventilatory disorders after TB, to the best of our knowledge, no studies have evaluated the association between restrictive ventilatory disorder and respiratory symptoms, physical activity limitations, and QoL in TB survivors. Our study revealed that respiratory symptoms and QoL were not significantly affected by restrictive ventilatory disorders in subjects with prior TB, while physical activity limitations were significant in subjects with prior TB with moderate-to-severe restrictive ventilatory disorders. Restrictive ventilatory disorders might be an underappreciated cause of functional impairments and respiratory symptoms [32,33]. One study showed that 35.4% of subjects with restrictive ventilatory disorders reported at least one chronic respiratory symptom [33]. One reason why our study results contrast with previous findings might be that most subjects with restrictive ventilatory disorders in our study had a mild degree of restrictive abnormality; thus, the number of subjects with moderate-to-severe restric-

tive ventilator disorders was relatively small. Accumulating evidence has shown that restrictive ventilatory disorders are related to physical activity limitations, which is in line with our study results [34,35]. The significant association of moderate-to-severe restrictive ventilatory disorders with physical activity limitations, but not with respiratory symptoms, suggests that restrictive ventilatory disorders influence physical activity limitations through mechanisms that are at least partly independent of respiratory symptoms. Decreased lung or chest wall compliance and increased elastic work of breathing might be one mechanism underlying physical activity limitations in patients with advanced restrictive ventilatory disorders [36].

This study has several limitations. First, this study was performed among a representative sample of Koreans. Thus, our data might not be generalizable to other ethnic groups or populations. Second, obstructive ventilatory disorders were defined by pre-bronchodilator spirometric results. This might lead to overestimation of the prevalence of obstructive ventilatory disorders; however, our estimates were similar to previous studies [20]. Third, the relatively small number of subjects with prior TB with moderate-to-severe restrictive ventilatory disorders might lead to statistical nonsignificance when analyzing the impact of ventilatory disorder severity on respiratory symptoms or QoL. Fourth, because this retrospective study used the Korea NHANES database that had been established before we designed this study, power calculations for sample size could not be performed during the design phase of the study. Fifth, there is a possibility that some subjects already had ventilatory disorders before *Mycobacterium tuberculosis* infection; however, in the Korea NHANES database, spirometry results before TB were not available. Sixth, as the Korea NHANES database does not provide some important clinical information that can affect respiratory symptoms and quality of life, such as medication use or exacerbation history, we could not adjust for these factors.

#### **5. Conclusions**

Among TB survivors, 29% had obstructive ventilatory disorders and 16% had restrictive ventilatory disorders. Severe obstructive ventilatory disorders were associated with an increased health-related burden, including more respiratory symptoms, more physical activity limitations, and poorer QoL, while severe restrictive ventilatory disorders were associated with more physical activity limitations. Further research is needed to establish strategies for early diagnosis and adequate treatment of ventilatory disorders in TB survivors.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/jpm11070678/s1, Table S1: The impact of obstructive and restrictive ventilatory disorders on respiratory symptoms, physical activity limitations, and EQ-5D index values in subjects with prior TB.

**Author Contributions:** Conceptualization, H.Y.P. and H.L.; data curation, B.Y., H.C., and Y.K.; formal analysis, B.Y. and H.C.; funding acquisition, B.Y.; H.C., H.Y.P. and H.L.; methodology, J.-Y.M., H.Y.P. and H.L.; supervision, H.Y.P. and H.L.; visualization, B.Y. and S.H.S.; H.C., H.Y.P. and H.L.; writing original draft, B.Y.; H.C., H.Y.P. and H.L.; writing—review and editing, All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by a National Research Foundation of Korea (NRF) grant funded by the Ministry of Science, Information, and Communications Technologies (MSIT) (No. 2020R1F1A1070468 and No. 2021M3E5D1A0101517621 to HL) and by the Korea Medical Device Development Fund grant funded by the Korean government (Ministry of Science and ICT, Ministry of Trade, Industry and Energy, the Ministry of Health and Welfare, the Ministry of Food and Drug Safety) (Project Number: KMDF\_PR\_20200901\_0214-2021-03 to HL).

**Institutional Review Board Statement:** The study protocol was approved by the Institutional Review Board of Chungbuk National University Hospital (application no. 2021-01-041).

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors have no conflict of interest to declare.

#### **References**


## *Article* **Impact of Smoking Status on Lung Cancer Characteristics and Mortality Rates between Screened and Non-Screened Lung Cancer Cohorts: Real-World Knowledge Translation and Education**

**Fu-Zong Wu 1,2,3,4,5 , Yun-Ju Wu <sup>2</sup> , Chi-Shen Chen <sup>6</sup> and Shu-Ching Yang 1,\***


**Abstract:** This was a retrospective hospital-based cohort study of participants diagnosed with lung cancer in the lung cancer register database, and our goal was to evaluate the impact of smoking and screening status on lung cancer characteristics and clinical outcomes. According to the hospital-based lung cancer register database, a total of 2883 lung cancers were diagnosed in 2883 patients between January 2007 and September 2017, which were divided into four groups according to smoking and screening status. A comparison was performed in terms of clinical characteristics and outcomes of lung cancer between the four groups. For non-smokers, age, gender, screened status, tumor size, targeted therapy, and curative surgery were independent prognostic factors of overall survival for lung cancer subjects. However, screened status and gender were not significant prognostic factors for lung cancer survival in smokers with lung cancer. For the non-smoker group, about 4.9% of lung cancer subjects (N = 81) were detected by screening. However, only 0.97% of lung cancer subjects (N = 12) were detected by screening in smokers. This could be attributed to smokers' negative attitudes and low socioeconomic status preventing LDCT lung cancer screening. In summary, our real-world data suggest that effectively encouraging smokers to be more willing to participate in lung cancer screening programs with screening allowance and educational training in the future is an important issue.

**Keywords:** lung cancer screening; smoking; willingness to screen; education

#### **1. Introduction**

Lung cancer is the most common cause of cancer-related death among both men and women worldwide in recent years, and about 70% of lung cancers are diagnosed at an advanced stage with poor survival outcomes [1–3]. Recent systematic reviews and meta-analyses across 31 randomized control trials (RCTs) have demonstrated the prolonged, substantial survival benefit of low-dose computed tomography (LDCT) for lung cancer screening in high-risk heavy smokers, with a 20% reduction in lung-cancer-specific mortality compared with the control group [4].

In recent years, several previous studies have emphasized the importance of a lung cancer screening program in the population with a high prevalence of non-smoking-related

**Citation:** Wu, F.-Z.; Wu, Y.-J.; Chen, C.-S.; Yang, S.-C. Impact of Smoking Status on Lung Cancer Characteristics and Mortality Rates between Screened and Non-Screened Lung Cancer Cohorts: Real-World Knowledge Translation and Education. *J. Pers. Med.* **2022**, *12*, 26. https://doi.org/10.3390/jpm 12010026

Academic Editor: Rutger A. Middelburg

Received: 19 October 2021 Accepted: 15 November 2021 Published: 2 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

lung cancer (reported prevalence of 1.55 to 2.34%) [5–10]. These studies have addressed that implementing a mass LDCT lung cancer screening program targeting the high-risk population could reduce the lung cancer mortality rate in the population with a high prevalence of non-smoking-related lung cancer.

Smoking causes genetic damage to epithelial cells in the lung and suppresses immune surveillance, which leads to an increased risk of lung cancer development [11,12]. Therefore, smoking is a well-known major risk for lung cancer development.

To the best of our knowledge, few studies have investigated the interaction of smoking and screening status on lung cancer characteristics and mortality rate. Therefore, this study was designed to provide an integrated approach to studying the relationship among smoking status, screening status, willingness to screening, and lung cancer characteristics with the prognostic outcome of subjects in a hospital-based lung cancer register cohort.

#### **2. Materials and Methods**

#### *2.1. Data Collection*

This was a retrospective cohort study based on the hospital-based lung cancer register database at Kaohsiung Veterans General Hospital between January 2007 and September 2017. The focus in this present study is on clinical characteristics and outcomes of lung cancer between screened and non-screening patient groups. The Institutional Review Board of Kaohsiung Veterans General Hospital approved this retrospective study, and thus informed consent was waived (VGHKS19-CT2-09). All research was performed in accordance with the relevant guidelines and regulations and all Institutional Review Board requirements. The study included 2883 patients aged 40–80 years with lung cancer diagnosed at Kaohsiung Veterans General Hospital, divided into four groups according to smoking and screening status. The study flowchart is illustrated in Figure 1. Patients were classified into four groups: group 1 included 1570 lung cancer subjects who had never been screened for lung cancer with LDCT and had no history of smoking; group 2 included 81 lung cancer subjects who had been screened for lung cancer with LDCT and had no history of smoking; group 3 included 1220 lung cancer subjects who had never been screened for lung cancer with LDCT and had a history of smoking, and group 4 included 12 lung cancer subjects who had been screened for lung cancer with LDCT and had a history of smoking. The subjects were classified according to the stringent definition of screening criteria. Positive screening status was defined as patients who had no clinical symptoms but had incidentally detected lung cancers by a self-paid LDCT exam between the ages of 40–80 years. Non-smoker was defined as an adult who had never smoked or who had smoked less than 100 cigarettes in his or her lifetime. Positive smoking history was defined as current smokers or a cessation of smoking within the previous 5 years, and who had smoked at least 100 cigarettes in his or her lifetime according to National Health Interview Survey. The information on survival and lung cancer deaths is based on the hospital-based center register data.

Clinical characteristics were recorded and compared for analysis in the four groups. The hospital-based lung cancer register database included the following patient demographics and clinical characteristics that are summarized in Table 1 (age, sex, tumor size, histopathologic type, adenocarcinoma spectrum classification, lung cancer clinical stage, lung cancer death, survival time, mortality rate, curative surgical treatment, targeted therapy, smoking habit, betel nut consumption, and alcohol drinking habit). The histological diagnosis was described according to the World Health Organization classification. All the patients were staged according to the seventh edition of the TNM staging system [13].

*J. Pers. Med.* **2022**, *12*, 26

**Table 1.** Demographic features of enrolled subjects for group (N = 2883) comparison based on smoking and screening status (means, standard deviations, and ANOVA results).


invasive adenocarcinoma; SD: standard deviation. Mean age at diagnosis, years

Median age at diagnosis,

**Figure 1.** Flowchart of patient recruitment. Patients were classified into four groups according to smoking and screening status. **Figure 1.** Flowchart of patient recruitment. Patients were classified into four groups according to smoking and screening status.

#### **Table 1.** Demographic features of enrolled subjects for group (N = 2883) comparison based on smoking and screening *2.2. Statistical Analysis*

status (means, standard deviations, and ANOVA results). **Group 1 (N = 1570) Group 2 (N = 81) Group 3 (N = 1220) Group 4 (N = 12)** *p***-Value 1 vs. 2 1 vs. 3 1 vs. 4 2 vs. 3 2 vs. 4 3 vs. 4**  (mean, SD) 66.36 ± 12.78 59.41 ± 7.41 69.15 ± 13.00 63.33 ± 11.83 <0.0001 <0.0001 <0.0001 1 <0.0001 1 0.695 years (range) 66 (40–99) 66 (42–77) 71 (41–99) 72 (42–83) Gender (n, %) <0.0001 0.437 <0.0001 <0.0001 <0.0001 <0.0001 1 Male 589 (37.5%) 24 (29.6%) 1175 (96.3%) 12 (100%) Female 981 (62.5%) 57 (70.4%) 45 (3.7%) 0 (0%) Smoking 0 (0%) 0 (0%) 1220 (100%) 12 (100%) Alcohol consumption 39 (2.5%) 4 (4.9%) 405 (33.2%) 6 (50%) <0.0001 1 <0.0001 <0.0001 <0.0001 <0.0001 0.481 Betel nut consumption 8 (0.5%) 0 (0%) 187 (15.3%) 3 (25%) <0.0001 1 <0.0001 0.003 <0.0001 0.005 1 Histology <0.0001 0.009 <0.0001 1 <0.0001 0.613 1 Adenocarcinoma 1269 (80.8%) 79 (97.5%) 744 (61%) 9 (75%) The patient demographics and clinical characteristics are expressed as mean ± standard deviation (SD) or median (interquartile range, IQR) and frequency (%) for group comparison. Multiple group comparisons were completed by analysis of variance (ANOVA) for normally distributed data and Kruskal–Wallis test for skewed data. We used the Kolmogorov test to assess the normality assumption in analysis of variance. We used the post hoc Bonferroni test to analyze the differences between each group. Cox regression analysis was performed to calculate and compare the survival rate between screened and non-screening groups in smokers and non-smokers. The proportional hazards assumption refers to the fact that the hazard functions are multiplicatively related. A multivariate analysis was used to estimate the hazard ratios (HRs) based on the Cox regression model, adjusted for age, gender, smoking, alcohol consumption, betel nut consumption, tumor size, curative surgery, targeted therapy, histology, and screening status. All statistical analyses were performed with SPSS 18.0 for Windows (SPSS Inc, Chicago, IL). A *p*-value of 0.05 was regarded as significant.

#### Small cell carcinoma 79 (5%) 1 (1.2%) 178 (14.6%) 0 (0%) **3. Results**

Mean survival days 676.03 ± 600.47 892.05 ±

Squamous cell carcinoma 170 (10.8%) 1 (1.2%) 273 (22.4%) 2 (16.7%)

Other 52 (3.3%) 0 (0%) 25 (2%) 1 (8.3%) Adenocarcinoma spectrum <0.0001 <0.0001 1 1 <0.0001 <0.0001 1 AAH 0 (0%) 6 (7.4%) 0 (0%) 0 (0%) AIS 0 (0%) 7 (8.6%) 0 (0%) 0 (0%) MIA 0 (0%) 9 (11.1%) 0 (0%) 0 (0%) IPA 1269 (100%) 57 (70.4%) 744 (60.9%) 9 (75%) The baseline characteristics of lung cancer subjects classified into four groups according to the screened and smoking status are shown in Table 1. There were statistically significant differences across the four groups for the clinical characteristics and outcomes demonstrated by one-way ANOVA. Bonferroni post hoc test results (inter-group comparison) are also shown in Table 1.

01

#### Stage <0.0001 <0.0001 <0.0001 <0.0001 <0.0001 0.95 <0.00 *3.1. Non-Smokers in Two-Group Comparison*

Carcinoma in situ 4 (0.3%) 14 (17.3%) 0 (0%) 0 (0%) I 326 (20.8%) 54 (66.7%) 129 (10.6%) 8 (66.7%) II 64 (4.1%) 4 (4.9%) 60 (4.9%) 1 (8.3%) III 334 (21.3%) 2 (2.5%) 302 (24.8%) 2 (16.7%) IV 842 (53.6%) 7 (8.6%) 729 (59.8%) 1 (8.3%) Curative surgery rate 459 (29.2%) 35 (83.3%) 211 (17.3%) 10 (83.3%) <0.0001 <0.0001 <0.0001 <0.0001 <0.0001 1 <0.00 01 Targeted therapy 319 (20.3%) 5 (11.9%) 197 (16.1%) 0 (0%) 0.008 0.982 0.028 0.418 1 1 0.898 Mean tumor size (mm) 41.25 ± 23.36 16.16 ± 13.74 51.12 ± 26.60 26.75 ± 19.99 <0.0001 <0.0001 <0.0001 0.249 <0.0001 0.976 0.004 Deaths 1159 (73.8%) 8 (9.9%) 1047 (85.8%) 2 (16.7%) <0.0001 <0.0001 <0.0001 <0.0001 <0.0001 1 <0.00 01 For non-smokers, the general characteristics and outcomes of the lung cancer subjects between the screened and non-screened groups are summarized in Supplemental Table S1. Table 2 summarizes the mortality and survival analysis according to screened and nonscreened status for non-smokers. For non-smokers, patients in the screened group had lower overall mortality than those in the non-screened group. One- and five-year mortality rates increased significantly from the screened group to the non-screened group (*p* < 0.001 for all). The mean survival time for the screened group was 892.05 ± 516.24 days (median 825 days), and the mean survival for the non-screened group was 676.03 ± 600.47 days (median 517.5 days).

337.21 <0.0001 0.003 <0.0001 1 <0.0001 0.869 1

516.24 444.85 ± 468.43 646.08 ±


**Table 2.** Clinical outcomes and mortality profiles of non-smokers with lung cancer according to screening status.

Cox regression was used to examine the association between screening status and survival and identify predictors of survival among non-smokers. For non-smokers, in the survival analysis for lung cancer subjects using the Cox regression model for multivariate effects, the hazard ratio for lung cancer mortality was determined adjusting for age, gender, smoking, alcohol consumption, betel nut consumption, tumor size, curative surgery, targeted therapy, histology, and screening status (Table 3). Multivariate survival analysis using Cox's regression model showed that age (HR = 1.011, *p* < 0.001) and tumor size (HR = 1.012, *p* < 0.001) were associated with unfavorable survival for lung cancer subjects, as shown in Table 3. However, Cox's regression model showed that a screening status (HR = 0.480, *p* = 0.040), gender (HR = 0.861, *p* = 0.034), targeted therapy (HR = 0.839, *p* = 0.031), and curative surgery (HR = 0.196, *p* < 0.001) were identified as independent, favorable prognostic factors of overall survival for lung cancer subjects, as shown in Table 3.



Abbreviations: CT: confidence interval; histology: adenocarcinoma versus other histology types (reference); gender (reference group: male); alcohol consumption (reference group: no alcoholic drinks); betel nut consumption (reference group: no betel nut consumption); screened (reference group: unscreened status); targeted therapy: (reference group: no targeted therapy); curative surgery (reference group: no curative surgery).

#### *3.2. Smokers in Two-Group Comparisons*

For smokers, the general characteristics and outcomes of the lung cancer subjects between the screened and non-screened groups are summarized in Supplemental Table S2. Table 4 summarizes the mortality and survival analysis according to screened and nonscreened status of smokers. For smokers, patients in the screened group had lower overall mortality than those in the non-screened group. One- and five-year mortality rates increased significantly from the screened group to the non-screened group (*p* < 0.001 for all). The mean survival time for the screened group was 646.08 ± 337.21 days (median 683 days), and the mean survival for the non-screened group was 444.85 ± 468.43 (median 304 days, *p* = 0.064). However, the *p*-value from the test may not reach statistical significance with a small sample size in the screened group (N = 12).

Cox regression was used to examine the association between screening status and survival and identify predictors of survival among smokers. For smokers, in the survival analysis for lung cancer subjects using the Cox regression model for multivariate effects, the hazard ratio for lung cancer mortality was determined adjusting for age, gender, smoking, alcohol consumption, betel nut consumption, tumor size, curative surgery, targeted

therapy, histology, and screening status (Table 5). Multivariate survival analysis using Cox's regression model showed that age (HR = 1.014, *p* < 0.001) and tumor size (HR = 1.011, *p* < 0.001) were identified as independent, unfavorable prognostic factors of overall survival for lung cancer subjects, as shown in Table 5. Multivariate survival analysis using Cox's regression model showed that targeted therapy (HR = 0.792, *p* = 0.023) and curative surgery (HR = 0.202, *p* < 0.001) were identified as independent, favorable prognostic factors of overall survival for lung cancer subjects, as shown in Table 5.

**Table 4.** Clinical outcomes and mortality profiles of smokers with lung cancer according to screening status.


**Table 5.** Multivariate Cox regression model of prognostic factors of smokers with lung cancer.


Abbreviations: CT: confidence interval; histology: adenocarcinoma versus other histology types (reference); gender (reference group: male); alcohol consumption (reference group: no alcoholic drinks); betel nut consumption (reference group: no betel nut consumption); screened (reference group: unscreened status); targeted therapy: (reference group: no targeted therapy); curative surgery (reference group: no curative surgery).

#### **4. Discussion**

The focus of the current study was to investigate the relationship between smoking status, screening status, willingness to screen, and lung cancer characteristics and the prognostic outcome of subjects in the hospital-based lung cancer register cohort. In this study, we demonstrated five major findings. The first one is that 57% of cancer subjects are non-smoking-related lung cancers, according to the hospital-based lung cancer register cohort. The second finding is that about 4.9% of lung cancer subjects (N = 81) are detected by screening in non-smokers. Third, only 0.97% of lung cancer subjects (N = 12) are detected by screening in smokers. Fourth, age, gender, screening status, tumor size, targeted therapy, and curative surgery are independently significant prognostic factors for lung cancer survival in non-smokers. Fifth, age, tumor size, targeted therapy, and curative surgery are independently significant prognostic factors for lung cancer survival in smokers.

In this present study, more than half (57%) were non-smoking-related lung cancers, according to the hospital-based lung cancer register cohort. Our study result is in line with previous studies in Asian countries. This finding suggests an increasing trend in the prevalence of non-smoking-related lung cancer in the Asian population [14–16]. In addition, if we only adopted the National Lung Screening Trial (NLST) criteria in the Asian population, approximately 70% of lung cancer subjects might have been missed, according to a previous study in Japan [17]. Therefore, not only heavy smokers, but also non-smokers at high risk should be enrolled in the target population for the LDCT

lung cancer screening program. For the non-smoker group, about 4.9% of lung cancer subjects (N = 81) are detected by screening in non-smokers. However, only 0.97% of lung cancer subjects (N = 12) are detected by screening in smokers. This could be contributed to smokers' negative attitudes and low socioeconomic status preventing LDCT lung cancer screening [18–20], as shown in Figure 2. Previous studies have shown that heavy smokers are less willing to be screened for lung cancer than former or non-smokers, because smokers do not believe in the survival benefits of early stage diagnosis [18,19]. In addition, results from our previous studies show that only about 28.75% of the study cohort had a smoking habit in our self-paid LDCT screening program [9]. Only 15% of the self-paid LDCT screening groups were eligible for NLST criteria [9,21]. Among all-screen detected lung cancer subjects, up to 87.1% (N = 81) of lung cancer subjects were non-smokers in the current study. As the number of screen-detected lung cancer subjects in smokers is still small (N = 12), larger studies will be needed to further explore the cause in the future. The percentage of screened lung cancer subjects in the non-smoker group was significantly higher than that of those in the smoker group (4.9% versus 0.97%; *p* < 0.001, shown in Figure 2). *J. Pers. Med.* **2021**, *11*, x FOR PEER REVIEW 7 of 10 smokers do not believe in the survival benefits of early stage diagnosis [18,19]. In addition, results from our previous studies show that only about 28.75% of the study cohort had a smoking habit in our self-paid LDCT screening program [9]. Only 15% of the self-paid LDCT screening groups were eligible for NLST criteria [9,21]. Among all-screen detected lung cancer subjects, up to 87.1% (N = 81) of lung cancer subjects were non-smokers in the current study. As the number of screen-detected lung cancer subjects in smokers is still small (N = 12), larger studies will be needed to further explore the cause in the future. The percentage of screened lung cancer subjects in the non-smoker group was significantly higher than that of those in the smoker group (4.9% versus 0.97%; *p* < 0.001, shown in Figure 2).

**Figure 2.** Screening status percentages according to lung cancer subjects based on smoking status. **Figure 2.** Screening status percentages according to lung cancer subjects based on smoking status.

These findings suggested that making heavy smokers more willing to participate in lung cancer screening programs is an important issue in the future [22]. A previous study demonstrated educational disparities in smokers, regardless of age and gender [23]. In line with the previous studies, our study suggests that we should pay more attention to shared-decision-making educational plans in heavy smokers. These findings suggested that making heavy smokers more willing to participate in lung cancer screening programs is an important issue in the future [22]. A previous study demonstrated educational disparities in smokers, regardless of age and gender [23]. In line with the previous studies, our study suggests that we should pay more attention to shared-decision-making educational plans in heavy smokers.

In this study, we aimed to investigate the difference of prognostic factors between smoker and non-smoker lung cancer subjects. In this study, a screened status is an important favorable prognostic factor for lung cancer in non-smokers. However, our study results did not find screened status to be an important prognostic factor for lung cancer In this study, we aimed to investigate the difference of prognostic factors between smoker and non-smoker lung cancer subjects. In this study, a screened status is an important favorable prognostic factor for lung cancer in non-smokers. However, our study results did not find screened status to be an important prognostic factor for lung cancer in

in smokers. Recent studies have demonstrated that low-dose computed tomography lung cancer screening trials have revealed a reduction in lung cancer mortality in sub-

12) were detected by screening in smokers. The rate is much lower than that of non-smoking lung cancer patients participating in lung cancer screening, as shown in Figure 2. These findings suggest that increasing smokers' participation in lung cancer

smokers. Recent studies have demonstrated that low-dose computed tomography lung cancer screening trials have revealed a reduction in lung cancer mortality in subjects with a smoking history [24,25]. In our study, only 0.97% of lung cancer subjects (N = 12) were detected by screening in smokers. The rate is much lower than that of non-smoking lung cancer patients participating in lung cancer screening, as shown in Figure 2. These findings suggest that increasing smokers' participation in lung cancer screening will improve lung cancer survival. However, several studies have shown that smokers are less likely to participate in lung cancer screening programs due to their low socioeconomic status or low willingness to screen [18–20]. As a result, the self-paid lung cancer screening program in this study has difficulty achieving significant effects in smokers. Through policy subsidies and shared-decision-making educational practices, barriers can be overcome, and the rate of smokers receiving lung cancer screening can be effectively increased. Therefore, screening can gradually reduce the lung cancer mortality of smokers in the real world.

In the non-smoking group, we found that female gender was a favorable prognostic factor predicting better survival. This could be due to the different biologic mechanisms of lung cancer that depend upon smoking status [26,27]. Previous studies have demonstrated a high prevalence of non-smoking-related lung cancer in the Asian population in the screening program, especially in women, usually manifesting with persistent subsolid nodules with indolent behavior [9,28,29]. In addition, about 30% of screen-detected lung cancer subjects met the diagnostic criteria for multiple primary lung cancers, according to previous studies [30,31]. However, these non-smoking lung cancer subjects have a better prognostic outcome of lung cancer. On the other hand, screening may also cause overdiagnosis and overtreatment during the screening process [32,33]. To maximize the benefit of LDCT screening and minimize the potential harm in over-diagnosis and overmanagement among non-smokers, a lung-cancer-risk-prediction algorithm with shared decision-making plans for indeterminate pulmonary nodule management is mandatory for a successful LDCT lung cancer screening program [34,35]. For smokers, effectively encouraging smokers to be more willing to participate in lung cancer screening programs with screening allowances and educational training programs (quitting smoking) in the future is an important issue. Promoting a personalized lung cancer screening program requires the support of government policies, primary care physicians, and public education for the incorporation of patient lung cancer risk, preferences, and socioeconomic status to resolve the complex trade-offs in the screening process.

This study benefits from real-world data through the prolonged implementation of the self-paid LDCT lung cancer screening program in this hospital-based cohort. This study describes why and how more non-smokers (especially women) could benefit from being diagnosed with lung cancer at an early stage. However, there are three major limitations in this study that could be addressed in future research. First, this study did not explore differences in attitudes and behavior towards screening for lung cancer between smokers and non-smokers. Real-world data from this study show that the smoking group had a very low proportion (0.97%) of lung cancer detected by screening, and only 12 people were diagnosed in this group (group 4). Our data also indirectly reflect smokers' negative attitudes and behavior towards participating in lung cancer screening. This finding is supported by previous studies that described that heavy smokers are less willing to be screened for lung cancer than former smokers or smokers [18,19]. Second, our self-paid, screened population may not represent the screening population eligible for NLST criteria. Therefore, the study result may not be generalized to the population at other institutions or with other demographics. Third, potential residual confounding in this study lowers the certainty of the evidence.

#### **5. Conclusions**

In summary, this study aimed at investigating the different impacts of smoking and screening status on lung cancer characteristics and mortality rates. For non-smokers, female gender and screened status are two independent, important, and favorable prognostic factors for lung cancer survival. However, these findings could not be demonstrated in smokers, which may be due to limited numbers of lung cancer subjects classified as smokers who have a negative attitude and low socioeconomic status preventing LDCT lung cancer screening. Effectively encouraging smokers to be more willing to participate in government-subsidized lung cancer screening programs in the future is an important issue.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/jpm12010026/s1, Table S1: Clinical characteristics of 1651 non-smokers with lung cancer according to screening status, Table S2: Clinical characteristics of 1232 smokers with lung cancer according to screening status.

**Author Contributions:** Conceptualization, F.-Z.W. and S.-C.Y.; methodology, F.-Z.W.; software, F.-Z.W. and Y.-J.W.; formal analysis, F.-Z.W. and Y.-J.W.; investigation, F.-Z.W. and Y.-J.W.; writing original draft preparation, F.-Z.W. and S.-C.Y.; writing—review and editing, F.-Z.W.; visualization, F.-Z.W., Y.-J.W. and C.-S.C.; supervision, F.-Z.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by grants from the Kaohsiung Veterans General Hospital, KSVGH110-120, KSVGH111-108, KSVGH111-D02-3, Taiwan, R.O.C, MOST 109-2314-B-075B-006-, MOST 110-2314-B-075B-008-.

**Informed Consent Statement:** Informed consent was waived due to the retrospective study design.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** We thank all doctors who responded to our screening and investigation. This study is based in part on data from the Cancer Registry Database provided by the Cancer Center of the Kaohsiung Veterans General Hospital.

**Conflicts of Interest:** All authors have declared no conflict of interest.

#### **References**


## *Article* **Dental Treatment Needs and Related Risk Factors among School Children with Special Needs in Taiwan**

**Szu-Yu Hsiao 1,2, Ping-Ho Chen <sup>2</sup> , Shan-Shan Huang <sup>1</sup> , Cheng-Wei Yen <sup>1</sup> , Shun-Te Huang 1,3, Shu-Yuan Yin <sup>4</sup> and Hsiu-Yueh Liu 3,5,\***


**Abstract:** The purpose of this study was to assess dental treatment needs (TNs) and related risk factors of children with disabilities (CD). This cross-sectional study recruited 484 CD, 6 to 12 years of age, from 10 special education schools in Taiwan. Dental status and TNs were examined and evaluated by well-trained dentists and based on the criteria set by the World Health Organization (1997). The results indicated that 61.78% required restorative dental treatment due to their dental caries. On average, each participant had 2.72 teeth that required treatment, and 6.38 surfaces required restoration. One-quarter of the participants (24.79%) required 1- or 2-surface restoration, and one out of three (36.98%) had more complex TNs (including 3 or more surfaces to be filled, pulp care, extraction, and more specialized care). The significant risk factors associated with restorative TNs among CD were those whose parents had lower socioeconomic status, frequent sweets intake, insufficient toothbrushing ability, and poor oral health. Most of the CD had extensive unmet TNs for their caries and required complex treatment to recover the function of their teeth. Encouraging parents/caregivers to take their children for dental treatment, promoting awareness of the importance of dental hygiene, giving assistance to brushing their teeth after eating, and controlling and/or modifying sweet diet habits are necessary to reduce CD's dental caries, especially those with lower socioeconomic status parents/caregivers.

**Keywords:** disability; children; oral health; caries; dental treatment

#### **1. Introduction**

People with disabilities usually suffer from a significantly higher prevalence of poor dental hygiene, plaque accumulation, dental caries, gingivitis, and periodontal disease than the ordinary population, and it gets worse with increasing age [1–6]. The top three medical care departments most frequently visited by people with intelligence disabilities in Taiwan are internal medicine (24.4%), psychiatry (16.7%), and dentistry (13.8%), as opposed to surgery, family medicine, and obstetrics and gynecology, which are utilized mainly by the general population [7]. People with disabilities have higher dental treatment needs (TNs) caused by secondary conditions as opposed to the general population, who seek and receive healthcare services regularly. According to Healthy People 2010 in the USA, secondary conditions refer to the problems related to medical, social, emotional, family, or community problems that a person with a primary disabling condition may encounter in his/her life [8]. It often aggravates and/or lessens their life quality in terms

**Citation:** Hsiao, S.-Y.; Chen, P.-H.; Huang, S.-S.; Yen, C.-W.; Huang, S.-T.; Yin, S.-Y.; Liu, H.-Y. Dental Treatment Needs and Related Risk Factors among School Children with Special Needs in Taiwan. *J. Pers. Med.* **2021**, *11*, 452. https://doi.org/10.3390/ jpm11060452

Academic Editor: Rutger A. Middelburg

Received: 17 April 2021 Accepted: 21 May 2021 Published: 23 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of physical, psychosocial, and social functions and increases the burden of health care on their parents, family members, or other caregivers with limited resources. Koritsas & Iacono reported that people with disabilities experienced an average of 11.3 secondary conditions (including medical complications, psychiatric disorders, environmental and quality-of-life issues, and difficulties with access to medical care/centers) during their development [9]. Secondary conditions often cause significant limitations, including reading difficulties, communication, physical fitness–conditioning, personal hygiene–appearance, weight, dental and oral hygiene, and memory problems [9]. Dental and oral hygiene is one of the significant limitations which is caused by secondary conditions.

A previous study reported that the greatest unmet health needs of people with disabilities are unmet dental treatment [10]. Studies have shown that children with disabilities (CD) usually receive less restorative care compared to the general population [11–15]. People with disabilities cannot express their physical discomfort properly because of their physical, intellectual, and psychological barriers. Long-term neglect of TNs, with increasing severity and complexity of oral diseases, leads to delay of time-sensitive, necessary treatment [16]. A previous study reported that only 32.37% of decayed teeth received restorative treatment among 6- to 12-year-old CD, and it was significantly lower than the ordinary population (47.72%) in Taiwan [17–19]. It reveals that the dental TNs for those children are unmet.

Dental treatment is a basic component of rehabilitation for CD, and it is difficult for dentists to perform [20]. The majority of children with special needs can be adequately treated using non-pharmacologic behavior management such as the tell-show-do technique [14,21]. Disabled people who have extensive and severe dental problems very often cannot cooperate well during the dental treatment process, and the treatment has to be performed by pharmacological behavior management techniques such as nitrous oxide/oxygen sedation, oral sedation. or general anesthesia (GA) to achieve higher quality treatment [14,22,23]. In order to avoid oral diseases' worsening and to provide CD with effective dental treatment services, it is helpful to provide appropriate baseline information regarding trends in unmet needs and related risk factors so as to make proper decisions to improve the oral health status of CD. Thus, the present study was carried out in an attempt to determine the unmet dental TNs and related factors of children with various disabilities in primary schools for CD in Taiwan.

#### **2. Method**

#### *2.1. Study Design and Sample Characteristics*

This cross-sectional study recruited 484 children, aged from 6 to 12 years old (mean age = 9.47 ± 2.06 years old), from 10 out of 18 special education primary schools in Taiwan. According to the definition of the Physically and Mentally Disabled Citizens Protection Act [24], people with disabilities refers to those who are limited or restricted from engagement in ordinary living activities and participation in society. All the participants have been evaluated by a committee composed of professionals from medicine, social work, special education, and employment counseling. After evaluation, they receive certificates that prove their classification and severity of disabilities after the processes of evaluation and assessment.

The CD were categorized according to their disability certificates in this study, which included vision disability (VD); hearing, voice or speech mechanism disability (HVD); intellectual disability (ID); and multiple disabilities (MD). MD means at least two or more types of impairment. The participants of this study included the certified MD children with ID and at least one or more co-occurring conditions. These conditions included vision disability, hearing mechanism disability, limb disability, etc.

#### *2.2. Ethical Approval*

This study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Institutional Review Board of Kaohsiung Medical University Hospital (Protocol number: KMUH-IRB-950125). The study purpose, procedures, and contents of the survey were explained to the parents or guardians of all CD, and written consent was obtained from parents or guardians who agreed to allow their children to participate.

#### *2.3. Data Collection*

#### 2.3.1. Oral Examination

Examination of the oral health status of study participants was carried out using a disposable dental mirror (Prosperity Island Medical Dressing Co, Ltd., Changhua, Taiwan), a standard dental explorer (CP-11.5B6, Hu-Freidy, Chicago, IL, USA), and a flashlight, and with the help of nursing staff in the classroom, lobby, auditorium, or other open space of the schools. An oral examination of the oral health status of study participants was performed in accordance with the methods and criteria of the World Health Organization [25], carried out by six well-trained dentists (S.-T.H., S.-Y.H., C.-Y.J., B.-M.C., C.-C.C.,R.-C.T.) who were evaluated prior to the survey by a senior dentist (gold standard: S.T.H.) with abundant clinical experience, and achieved acceptable reliability and inter-examiner agreement with a Kappa score between 0.80 and 0.84 (0.81, 0.84, 0.80, 0.80, and 0.82) between each two dentists.

The oral examination and records were determined on the basis of a maximum of 20–28 teeth for primary, mixed, or permanent dentition up to 12 years old. Unmet dental TNs were evaluated using dmt + DMT and dms + DMS indices to record numerically the number of unmet dental treatment teeth and surfaces. The dmt + DMT index is a sum of decayed teeth number (dt + DT) and missing teeth number (mt + MT) of mixed dentition. The dms + DMS index is a sum of decayed teeth surfaces (ds + DS) and missing surfaces (dm + DM) of mixed dentition. These two indices provide the consequences of untreated caries and present an aggregate value of current dental TNs.

Plaque and gingivitis indices were evaluated on 6 indicative surfaces (buccal side of 4 teeth and lingual side of 2 teeth) with the naked eye and a CPI explorer (CP-11.5B6, Hu-Freidy, Chicago, IL, USA). The 4 buccal surfaces selected for inspection were from 2 posterior teeth (tooth numbers 55 or 17/16 and 65 or 26/27) and 2 anterior teeth (tooth numbers 51 or 11 and 71 or 31). The 2 lingual surfaces selected for inspection were from 2 posterior teeth of the mandible (tooth numbers 85 or 47/46 and 75 or 36/37). The value of the poorest surface situation prevailed.

Plaque index (PI) was assessed according to criteria modified from Greene and Vermillion [26] and marked as: 0 = no plaque; 1 = plaque covering not more than one-third of the exposed tooth surface; 2 = plaque covering more than one-third, but not more than two-thirds, of the exposed tooth surface; and 3 = plaque covering more than two-thirds of the exposed tooth surface. The tooth surface with the worst plaque of six teeth was recorded as a whole-mouth oral hygiene score for each participant. Gingivitis index (GI) was assessed by inspection of the following and tendency toward spontaneously bleeding criteria: 0 = healthy; 1 = mild gingivitis—no bleeding on probing; 2 = moderate gingivitis bleeding on probing; 3 = severe gingivitis—ulceration [27]. The most-inflamed gingival surface of six teeth was identified as a whole gingivitis severity score for each participant. The PI and GI were classified into yes (score ≥ 1) and no (score = 0) for statistical analysis, respectively.

The dental TNs of each child were categorized as simple, moderate, and complex based on the following criteria [12,28]: simple = child had no caries and required no restorative treatment, but required assistance with oral hygiene and preventive treatment such as oral hygiene instruction (OHI), scaling, application of topical fluoride, and/or fissure sealants; moderate = one or more teeth caries and required one and/or two surface restorations; complex = one or more teeth caries and required three or four surface restorations/stainless steel crowns, endodontic therapy/crown, extraction, and/or prosthodontics.

#### 2.3.2. Questionnaire

The self-report questionnaire was completed by the parents or guardians of the participants. The questionnaire consisted of close-ended items and was constructed in three

parts: demographic data and family characteristics, sweet intake habits, and tooth-brushing habits of the participant. Sweet foods were considered as highly fermentable carbohydrates made with rich sugar to add a sweet taste such as chocolate, candy, cake, baked goods, ice cream, carbonated beverages, juice, milk with sugar, and chewing gum with sugar. The definition of the classification of parents' educational and occupation levels referred to the publications of Liu et al. [18]. The parents' educational levels were classified as follows: both low, one low and one high, or both high. The parents' occupation levels were divided into both unskilled, one unskilled and one skilled, or both skilled.

#### *2.4. Statistical Analysis*

All data were analyzed with JMP version 14 (SAS Institute, Cary, NC, USA). Categorical variables in a group were compared using Pearson's χ 2 test and Fisher's exact test and were presented as frequency and percentage, and the differences between numerical variables were analyzed using a *t*-test and one-way analysis of variance (ANOVA) and presented as the mean and standard deviation (SD). Both univariate and multivariate logistic regression models were estimated to assess the unadjusted and adjusted associations of a preset independent variable with TNs of dental caries. Only the independent variable, as the risk factor that was found to be significantly associated with dental TNs in the univariate logistic regression, was included in the multiple logistic regression models. A significance level of the *p*-value was set at 5%.

#### **3. Results**

The younger children had statistically significantly higher TNs of decayed teeth and surfaces than the older ones (all *p* < 0.001). We found that children with HVD had the highest TNs of teeth (3.48 ± 3.74) and surfaces (8.89 ± 11.57) compared to children with VD (1.87 ± 2.87 and 4.14 ± 8.59) (*p* = 0.023 and *p* = 0.030). The children with mild/moderate disabilities or with HVD had a statistically significantly higher decayed teeth and surfaces need for dental treatment than those with the other severities of disabilities or classifications of disabilities (*p* = 0.017 and *p* = 0.016). TN of surfaces statistically significantly decreased when the parents had higher educations or the parents' occupations tended to be skilled (*p* = 0.023 and *p* = 0.018). The results also showed the CD who asked for sweets, frequently had sweets, had independent tooth-brushing abilities, or infrequently brushed their teeth every day tended to have statistically significantly higher decayed teeth and surfaces for treatment (all *p* < 0.05) (Table 1).

**Table 1.** Dental treatment needs of teeth and surfaces by demographics and dietary and tooth-brushing habits among children with disabilities.



**Table 1.** *Cont*.

Overall, one-fourth of the children needed moderate treatment for one or two surface fillings, and more than one-third of the children needed complex treatment such as three or more surface fillings, pulp care, or crowns (Table 2). Dental treatment modalities had a statistically significantly positive association with age, severity of disability, parents' educational and occupation levels, asking for sweets, frequency of sweets intake, sweets as a reward for behavior control, and tooth-brushing ability of the participants (all *p* < 0.05). The dental TN modality of children with mild/moderate disability was statistically significantly more complex than those with severe/profound disability (*p* = 0.009). The higher the educational or occupational level of the parents, the less and simpler the dental TN modalities of their children (*p* = 0.015 and *p* = 0.045). Dental treatment modalities had a statistically significant positive association with the behaviors of asking for sweets and frequency of eating sweets (*p* = 0.021 and *p* = 0.001) among CD.

In order to better understand the participants, we further compared the oral status, demographics, sweet intake, and daily tooth-brushing habits among different disabilities, as shown in Table 3. We found the VD children had the lowest proportion of plaque and gingivitis (*p* = 0.036 and *p* < 0.001, respectively) than other disabled children. The MD children asked for sweets less, and their daily tooth-brushing is hard to perform by themselves. Less of the VD children's behavior was controlled by giving sweets as a reward than the others. Compared with other disabled children, the HVD children had a statistically significantly lower rate (8.24%) of brushing their teeth three or more-than-three times a day (*p* < 0.001).



Multiple logistic regression models showed the major risk factors for restorative TNs among CD were their intake of sweets at least once a week (AOR = 2.45, 95% CI = 1.48–4.11, *p* = 0.001) or more frequently than those who never or sometimes consumed sweets. The other risk factors—significantly related to children's poor oral health were having gingivitis (AOR = 1.94, 95% CI = 1.22–3.11, *p* = 0.006) and insufficient ability to brush their teeth without assistance from parents/caregivers (AOR = 1.95, 95% CI = 1.19–3.22, *p* = 0.008) (Table 4).


**Table 3.** Demographics, dietary and tooth-brushing habits, and oral hygiene status by classification of disability among children with disabilities.


**Table 4.** Logistic regression models of risk factors associated with restorative treatment needs among children with disabilities.

<sup>a</sup> COR: crude odds ratio. Data analysis by univariate logistic regression model. Dependent variable was with dental treatment needs. <sup>b</sup> AOR: adjusted odds ratio. Data analysis by multiple logistic regression model, adjusted participants' gender and age. Variables found with statistically significant associations in univariate logistic regression analysis were included in the multiple logistic regression models. Dependent variable was with treatment needs.

#### **4. Discussion**

The main findings of this study identified a higher caries prevalence and extensive unmet dental TNs among CD in special schools in Taiwan. High frequency of sweet intake and inadequate tooth-brushing ability were the critical risk factors attributed to CD' needs for extensive restorative care. Even though some CD have better independent abilities, they have the potential need for attention and assistance from their parents and/or caregivers to maintain and/or protect their oral hygiene with daily tooth-brushing behavior. This is beneficial to reduce the high dental caries treatment needed for their teeth and teeth surfaces.

The oral hygiene of CD is poor. This can be confirmed from our results that nearly 80% of CD had a mild-to-severe plaque index, and approximately 50% of them had gingivitis. Due to these CD being exposed to such a poor oral environment, 61.78% of them developed dental caries. If plaque is not completely removed from the teeth surfaces every day, a mild form of periodontal disease, which is called gingivitis, occurs. In our study, 47.93% of CD

had inflammation and bleeding of the gums. This results make us more convinced that poor oral hygiene will affect subsequent oral diseases such as dental decay and gingivitis, which is consistent with previous findings [12,18,28]. However, good oral hygiene can be achieved via regular removal of plaque by thorough tooth-brushing every day. If the plaque declines, tooth decay and gingivitis will simultaneously diminish.

Dental caries remains a primary health problem among CD. Unmet decayed and missing teeth in 6–12-year-old Taiwanese CD, with a mean age of 9.47 years old, have decreased, with the mean caries experience index declining from 3.36 teeth (72.77%) in 2005 [17] to 2.72 teeth (61.78%) in our study (2021). The same index in 9-year-old children also declined from 2.70 (73.55%) in 2007 to 2.66 (47.09%) in 2012 [17,29] among the ordinary population. The decreasing value of caries experience index among CD (10.99%) in our study was less than in those without disabilities (26.46%) [17,29]. In comparison with other studies, the caries prevalence in this study (61.78%) is higher than the result in Australia (56%) but lower than that in Iran (73.6%) [12,28]. The studies of 5–16-year-old CD in Australia and Iran reported less unmet dental treatment teeth (1.53 and 2.09), which was half to two-thirds of our findings (2.75) [12,28]. Furthermore, we compared the treatment modalities and found the CD in the present study had a higher need of complex treatment (36.64%) than the CD in the Australian (21%) and Iranian (25.1%) studies [12,28]. Fewer CD had caries lesions in the present study than in Iran, but they needed more complex dental treatment. It appeared that the CD in our study tended to have more severe tooth caries and required higher dental TNs than CD from other countries [28].

The moderate and complex dental TNs among CD were unmet in this study. Most of the patients receiving dental treatment under GA were people with special needs [14,21,30], especially when the patient needed to have extensive and complicated treatment [31]. The most common treatments (extraction, restoration, and pulp therapy) under GA over the past 10 years in Taiwan are related to a high proportion of multiple dental caries (86.4%) [32,33]. Our result is consistent with other studies, wherein children who need to be treated under GA usually have a higher unmet decayed teeth treatment of more than 10 teeth [32,33]. There were two-thirds (59.87%) out of 36.98% of our participants who might need to be treated under GA due to their uncooperative behavior in specialist clinics or regional hospitals that have critical care facilities. If the CD receive their first dental treatment as early as possible, they can achieve more effective dental rehabilitation [32].

There was a clear negative correlation between decayed teeth and surfaces for complex restorative TNs and severity of disability. This observation is in line with previous studies that children with profound disabilities acquire partial or complete assistance from others as it pertains to maintaining their oral health due to lack of adequate manual dexterity ability [12,22]. However, this is not always true, according to our findings. Our results showed the complex restorative TNs in HVD (71.26%) were higher than MD (58.02%) and VD (52.86%), which was also contrary to the results of a study in Saudi Arabia [34]. The studies of 5–16-year-old children with HVD in Saudi Arabia had restorative TN of 66.02%, which was lower than MD and VD (76.93% and 74.29%) in the present study. In this study, the HVD and VD children had the better capability to learn oral hygiene skills and take care of themselves in comparison with other groups. However, HVD children had ordinal vision and tended to be more susceptible to the allure of various sweets than VD children. The HVD children also had a better independent ability to access, buy, and obtain sweets by themselves than other disability groups. If HVD children do not perform the positive behavior of cleaning their teeth after sweets under supervision by parents/caregivers, it is likely to result in the repercussions of increasing the risk of more caries teeth and surfaces and higher restorative TNs.

Obstacles blocking these participants from achieving lower TNs in this study are improper sweet habits and inadequate tooth-brushing behavior. Manual dexterity is not a guarantee of better oral health and lower TNs among CD. The effectiveness of brushing is limited if CD receive less brushing assistance or supervision when they brush their teeth [18,34–36]. The higher the education level the parents/caregivers had, the more

aware they were about how to sufficiently complete the oral care needs of the CD [18,36,37]. If parents have higher occupation levels, they will have the sufficient financial capacity to take care of their children's oral hygiene to reach better oral health and decrease dental TNs [36,37]. However, the negative effect of caries TNs among CD from frequent sweets intake is superior to that from insufficient tooth-brushing ability [38]. The results of our study showed that frequent tooth-brushing could minimize the severity and prevalence of TN. The consumption of sweets in small amounts, along with other fermentable carbohydrates consumed frequently. will increase the caries risk and TNs, rather than large amounts eaten occasionally. Sweet limits, not only controlling the frequency, but also advising to give healthy or low-carcinogenic alternative foods such as fresh fruit and vegetables and/or low-carcinogenic foods by parents/caregivers contributed to a lower prevalence of untreated dental caries and TN rates among CD [39].

Before treatment, how to prevent and decrease the number and severity of TNs is important. Our study showed that tooth-brushing can, in part, diminish the association between giving sweets as a reward in behavior control, asking for sweets, and frequency of sweets intake on dental decay outcome in CDs [38]. Sticky foods with sugar and/or fermentable carbohydrates can stay in the oral environment for longer periods, thus increasing the potential and risk for tooth decay. Similar to ordinary children, having their children brush their teeth at least twice a day is deeply associated with parents'/caregivers' self-efficacy or confidence [40]. More frequent tooth-brushing might compensate for the inadequate tooth-brushing ability of children, as seen in previous studies [18]. To address the high TNs among CD in the present study, we need to provide promotion courses to encourage the parents/caregivers with lower educational levels, especially those in unskilled occupation levels, to implement effective preventive measures such as brushing after eating and correctly choosing healthy snacks for their children.

There are several limitations to the present study. One is that we have difficulty concluding regarding the causation between dental TNs and related risk factors by selfreported questionnaires, as is the case in most observational cross-sectional studies. Second, the observations were assessed outside the dental clinic, with limited accessibility and difficulties for patients' handling. Third, a small number of parents/caregivers might answer questions in such a way as to meet social expectations and thus cause answer bias. In addition, we collected the data from special schools, not including bedridden and homebound groups; therefore, our results may underestimate the real consideration of all disability groups and reflect the status of school CD. However, our findings may help to better understand the TNs of CD and propose some measures to improve their oral care.

#### **5. Conclusions**

The present findings illustrate TNs with an enormous need for restorative treatment, concentrated on 61.78% of CD. Improper sweet habits and inadequate tooth-brushing ability were found in the risk factors for having the highest odds of requiring restorative treatment. Study results reveal that oral problems among these children do not get much attention from their parents/caregivers, especially those with lower education or occupation levels. It is imperative to encourage low-educational-level parents/caregivers to take their children for dental treatment, teach them why and how to brush teeth after eating, and help them modify the sweet intake habits of their children to decrease their children's caries and improve their children's oral health.

**Author Contributions:** Conceptualization, S.-Y.H. and H.-Y.L.; methodology, P.-H.C. and S.-T.H.; software, P.-H.C.; validation, H.-Y.L.; formal analysis, S.-S.H. and C.-W.Y.; investigation, S.-Y.H. and S.-T.H.; data curation, S.-S.H. and C.-W.Y.; writing—original draft preparation, S.-Y.H. and S.-Y.Y.; writing—review and editing, H.-Y.L.; visualization, S.-Y.Y.; supervision, H.-Y.L.; project administration, H.-Y.L. The authors have confirmed that all authors meet the International Journal of Personalized Medicine criteria for authorship credit, as follows: (1) substantial contributions to conception and design, acquisition of data, or analysis and interpretation of data; (2) drafting the

article or revising it critically for important intellectual content; and (3) final approval of the version to be published. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Bureau of Health Promotion, Department of Health, Taiwan.

**Informed Consent Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Institutional Review Board of Kaohsiung Medical University Hospital (KMUH-IRB-950125). The participants' legal representatives provided written informed consent for the participants' involvement in the study.

**Acknowledgments:** The authors gratefully acknowledge the dentists from Tomlin children's dental clinic and Dada clinics, nurses of special schools, children, and their parents for conducting oral examinations, cordial help, participation, and cooperation, respectively. We would like to express our sincere gratitude to the administrators and staff of the facilities for people with disabilities for their support and assistance, and to all the wonderful children with disabilities for their generous participation.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Decreased Tongue Pressure Associated with Aging, Chewing and Swallowing Difficulties of Community-Dwelling Older Adults in Taiwan**

**Hsiu-Yueh Liu 1,2 , Jen-Hao Chen 3,4, Kun-Jung Hsu 3,5,6 , Ching-Teng Yao <sup>7</sup> , Ping-Ho Chen 3,8 , Szu-Yu Hsiao 3,9 and Chun-Li Lin 10,\***


**Abstract:** Personalized tongue pressure (TP) training focuses on improving swallowing. This study aims to establish the TP values of different age levels and compare changes between different swallowing status among community-dwelling elders. In this cross-sectional study, 1000 participants, aged 60 years old and above, were recruited from community care centers. All participants were classified into non chewing and/or swallowing difficulties (NCSD) and with chewing and/or swallowing difficulties (CSD) groups and their diseases and dieting status were recorded using a structured questionnaire. A disposable oral probe was used to measure TP by asking participants to compress it against the hard palate with maximum voluntary effort. Among 1000 elders, 63.10% had CSD and their TP (from 31.76 to 18.20 kPa) was lower than the NCSD group (from 33.56 to 24.51 kPa). Both groups showed the same tendency for TP decline with increasing age. Decline of TP makes CSD elderly have a poor appetite, eat a soft or liquid diet, and take longer to eat a meal (all *p* < 0.050). The secondary risk factor dominating TP decline for NCSD and CSD elders is having an education level less than primary school and an abnormal eating assessment, respectively. Our results demonstrated that TP decline has a significant relationship with age changes. Education level and an abnormal eating assessment score are closely associated with TP decline. A series of TP values can be used as a reference indicator of personalized medicine during the aging process among community-dwelling older adults.

**Keywords:** tongue pressure; aging; epidemics; community; elderly

#### **1. Introduction**

Taiwan's accelerated rate of aging issues have become an imperative topic in the country's long-term care policy and practice [1,2]. In response to the challenge of the rapidly growing older population, the new policy of long-term care, called LTC 2.0, was launched,

**Citation:** Liu, H.-Y.; Chen, J.-H.; Hsu, K.-J.; Yao, C.-T.; Chen, P.-H.; Hsiao, S.-Y.; Lin, C.-L. Decreased Tongue Pressure Associated with Aging, Chewing and Swallowing Difficulties of Community-Dwelling Older Adults in Taiwan. *J. Pers. Med.* **2021**, *11*, 653. https://doi.org/ 10.3390/jpm11070653

Academic Editor: Rutger A. Middelburg

Received: 17 June 2021 Accepted: 9 July 2021 Published: 11 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

planned, developed, and implemented in 2017 to prevent and delay the functional decline of the community-dwelling older population [3]. Aging has been confirmed to increase the decline of swallowing function by several studies [4–9]. Decreasing motor function in the lips, muscle strength and tongue pressure (TP), and loss of taste and smell with age increases and raises swallowing time, worsens swallowing performance, and impairs swallowing function [10]. A previous study in Taiwan found that 21.8% of the community's elderly in Taiwan aged 65–95 years old tended to choke at least three times a week during eating, and 12.8% of the study population were assessed as having swallowing disorders [8]. Based on this proportion, it is reasonable to estimate that approximately 480,000 elders in Taiwan have swallowing disorders. Swallowing disorders and/or difficulties reduce the willingness and ability to eat food or drink liquids they once loved, causes nutritional imbalance, weight loss, dehydration, and also increases the risk of and mortality from aspiration pneumonia in the elderly [9].

The aging process sees diminished muscle mass, a decline in muscle strength, and causes frailty and sarcopenia. As for most muscle groups, muscle reduction also occurs in the tongue muscles. Tongue muscle loss reflects tongue muscle strength [11,12]. The performance of tongue strength is considered the main driving force for moving the bolus during the chewing and swallowing of food. Swallowing efficacy and safety may worsen due to muscle fatigue and weakness of the swallowing-related tongue muscles as a result of repetitive swallowing during a long meal time over 30 min [10,13]. Due to the development and progress of science, the strength of the tongue muscle can be measured via a mechanical pressure device exerted by the tongue against the hard palate. These devices provide an objective method of assessing tongue strength for community screening and physical research investigation [14,15]. A meta-analysis conducted by Adams et al. reported average TP decreased by about 8–9 kPa from young to old adults overall [16]. Another analysis further reported that mean TP value ranges of healthy individuals ascertained by JMS (JMS TPM-01/02, JMS Co., Ltd., Hiroshima, Japan) [17] were significantly higher in healthy individuals less than 60 years (37.5 kPa to 41.1 kPa) relative to those 60 years or older (29.4 kPa to 31.1 kPa) [18].

Over the past three decades, tongue muscle strength has been used as a useful indicator to evaluate the physically frail, the eating and swallowing ability of the elderly, to screen patients with or without penetration, aspiration, or residue during eating, and assist patients (such as those with stroke, Parkinson's disease, and amyotrophic lateral sclerosis) with management targeted therapy goals or disease progression [13–15,19,20]. Personalized TP training is a new focus for improve swallowing. TP training or rehabilitation therapy at home is feasible, and it may become a new focus of personal precision medicine in the future. However, a large-scale TP database of older adults is currently unavailable in Taiwan. Moreover, the TP database may vary by country, subjects, and age groups. Before using TP to achieve personal swallowing precision medicine, it is necessary to establish a database of TP values among community-dwelling elderly. Therefore, this study aims to establish the TP values of different age levels among the community-dwelling elderly, and compare the TP changes between the elderly with different chewing and/or swallowing difficulties (CSD) status. These TP values of different age levels and CSD status may be used as an effective indicator for screening the changes of tongue function of the elderly and be an optimal treatment goal for patients with chewing and/or swallowing difficulties.

#### **2. Methods**

#### *2.1. Study Design and Population*

#### 2.1.1. Study Design

This observational cross-sectional study was carried out from October 2019 to November 2020. Community-dwelling older adults residing in Kaohsiung, Taiwan were enrolled in this study. The study adopted a multi-stage stratified cluster sampling design. The sampling probability is based on the probability proportional to size, with the community care centers as the sampling unit. First, Kaohsiung City has 38 districts stratified into

seven clusters based on the urbanization stratification of townships in Taiwan [21]. We merged the first and second clusters into urban areas (8 districts with 92 care centers), third and fourth clusters into town areas (9 districts with 89 care centers), and fifth to seventh clusters into rural areas (21 districts with 112 care centers). Second, 38 community care centers were randomly selected from 293 community care centers based on the population probability of older adults (60 years old or above), which were 47.76% (19 care centers), 31.39% (11 care centers) and 20.65% (7 care centers) of the 3 district clusters, respectively. Third, we recruited the eligible older adults from each selected community care center. A total of 1000 participants were recruited from urban areas (*n* = 463), new town areas (*n* = 339), and rural areas (*n* = 198).

#### 2.1.2. Sample Size

The sample size calculation considered the population of 616,280 older adults living in the city. We wanted to identify the sample size for a larger population and used the Slovin or Yamane formulas [22]. This formula, based on Krejcie and Morgan's recommendation [23], used 0.50 as an estimate of the population proportion to maximize variance, which will also produce the maximum sample size. With the population size of 616,280 older adults, at 95% confidence level and 5% margin of error, this resulted in our minimal sample size of 384 participants. We selected double the number of older adults necessary for research and considered 30% drop-out rate (as the rejection). Finally, the theoretical sample size in this designated study was at least 998 after original calculation. The actual sample size was adjusted to 1000 for further grouping and implementation. This study successfully recruited 1000 older adults (aged from 60 to 95 years old; 215 males, mean age of 74.75 ± 7.51 years; 785 females, mean age: 72.86 ± 7.58 years).

#### 2.1.3. Sample Selection and Randomization

The participants enrolled in this study (1) were equal to or older than 60 years; (2) had no craniofacial deformities or syndromes, no neuromuscular diseases, no history of head or neck cancer, had not undergone radiation therapy, and had no underlying neuromuscular diseases that were known as diseases that affect tongue strength. The exclusion criteria for the participants were (1) cognitively not fit for understanding and communication; (2) the coexistence of cognitive problems that affect language understanding, their inability to complete the tasks or the chance they may disrupt accurate TP measurement, and other problems that may impair the test.

#### 2.1.4. Ethical Approval

Ethical approval was obtained from the Human Experiment and Ethics Committees of Kaohsiung Medical University (Protocol number: KMUHIRB-F(I)-20190104). All procedures performed were in accordance with the ethical standards of the institutional and/or national research committee, and with the 1964 Helsinki Declaration and its later amendments or comparable ethical standards. Prior to oral examination and oral function evaluation, the purpose and content of the study protocol was thoroughly explained to the participants by the principal investigator. After agreeing to participate in the study, all participants signed a written informed consent.

#### *2.2. Data Collection Methods*

#### 2.2.1. Questionnaire Interview and Swallowing Assessment

A structured questionnaire was used to collect data on basic demographic characteristics (i.e., age, sex, education level), and diet status (i.e., type of diet and duration of meal time) referred to in the publication of Liu et al. [24]. Questionnaire data were collected by thoroughly trained interviewers using Mandarin or Taiwanese during face-to-face interviews in accordance with standard protocol. All interviewers participated in a 120-min training course on the standard process and data collection criteria to prevent information bias during interviews.

Swallowing difficulty screening used a self-administered questionnaire, called the Eating Assessment Tool (EAT-10) [25]. This self-administered tool is widely used to assess swallowing difficulty. According to the suggestion of EAT-10 normative data, a participant who scored 3 or more points on the EAT-10 was considered as having swallowing difficulty [25,26]. Each participant was also screened using the repetitive saliva swallowing test (RSST) within 30 s [27]. The RSST is a safe and simple screening test for functional swallowing difficulty. A participant who swallowed their saliva less than 3 times during 30 s of test was considered as having swallowing difficulty [27].

#### 2.2.2. Oral Examination

The oral examination of all study participants was performed by three well-trained and calibrated senior dentists according to the guidelines of the World Health Organization [28]. Intra-examiner reliability tests of tooth decay were carried out on these dentists during the oral data collection. The kappa coefficient was 0.81 for intra-examiner agreement. The participant's oral health examinations were performed in a group activity room in the care centers between 9.00 a.m. and 11.00 a.m. under natural light using disposable dental mirrors and CPI explorers, without radiographs.

Information was collected on the kind of denture and the condition of natural teeth, which included the caries status, location, and number. The dentist counted the number of functional natural teeth (FNT), fixed artificial teeth (FAT), removable artificial teeth (RAT), and functional teeth (FT) [29,30]. All examinations excluded third molars. Remaining natural teeth, teeth that were sound, decayed, filled, or filled but decayed were regarded as FT. Teeth with grade III mobility, residual roots, or extensive crown destruction (i.e., more than three-fourths of the clinical crown destroyed) were excluded. FAT was defined as fixed artificial teeth, including abutment teeth, pontics, and implant-supported prostheses. FAT with grade III mobility were excluded. RAT was regarded as artificial teeth on removable dentures being worn during the dental examination. FT included FNT, FAT, and RAT. The total number of FT ranged from 0 to 28. The participants with more teeth was helpful for assessment of TP [11,26,31]. According to the number of FT, participants were classified into 0–9, 10–19, and ≥20 teeth groups to analyze the relationship between FT and TP [11].

#### 2.2.3. Tongue Pressure Measurement

A wireless TP measurement device (OCPT-168, Shi-Heng Technology Co., Ltd., Taipei, Taiwan) was used to assess TP [32]. The device is a personal use device which consists of a disposable position mouthpiece and a pressure transducer (Figure 1a). During the test, the participant sits with foot support, keeps head upright, and eyes on a horizontal target (Figure 1b). A mouthpiece is used for each participant. The TP mouthpiece is placed behind the central incisor along the central groove of the tongue blade. The participant is asked to press their tongue against the hard palate as hard as possible for 5 s and a maximum TP is obtained on the digital screen of the device. Each participant was tested 3 times with 30-s rest interval in between. The obtained maximum value was recorded as maximum TP of each participant. For the data analysis, participants' TP less than the median value (27 kPa) was considered as decreased TP.

**Figure 1.** (**a**) A tongue pressure measurement device with wireless mobile application control function and a disposable oral positioning mouthpiece. (**b**) The tongue pressure device and its usage in the measurement of tongue pressure. **Figure 1.** (**a**) A tongue pressure measurement device with wireless mobile application control function and a disposable oral positioning mouthpiece. (**b**) The tongue pressure device and its usage in the measurement of tongue pressure.

#### *2.3. Statistical Analysis 2.3. Statistical Analysis*

In order to clarify the association between TP and CSD, the participants were divided into non CSD (NCSD) and CSD groups. If the participants had fewer than 20 FNT an RSST equal or less than three times, or an EAT-10 score equal or more than 3 or selfreported as having swallowing difficulty [25,26], they were classified into the CSD group. The participants of the NCSD group were the healthy older adults, who had no chewing or swallowing difficulties. Further to this, the participants in both groups (NCSD and CSD) were categorized into 6 age levels (60–64 years, 65–69 years, 70–74 years, 75–79 years, 80–84 years, and 85 years or above) to identify age-related changes in TP. In order to clarify the association between TP and CSD, the participants were divided into non CSD (NCSD) and CSD groups. If the participants had fewer than 20 FNT an RSST equal or less than three times, or an EAT-10 score equal or more than 3 or self-reported as having swallowing difficulty [25,26], they were classified into the CSD group. The participants of the NCSD group were the healthy older adults, who had no chewing or swallowing difficulties. Further to this, the participants in both groups (NCSD and CSD) were categorized into 6 age levels (60–64 years, 65–69 years, 70–74 years, 75–79 years, 80–84 years, and 85 years or above) to identify age-related changes in TP.

The participants' demographic characteristics, questionnaire data, and TP were entered into Microsoft Excel (Microsoft, Redmond, WA, USA). IBM SPSS Statistics 20.0 (IBM, Armonk, NY, USA), which was then used to perform descriptive and inferential statistical analysis. Demographic characteristics were expressed as numbers and percentages, and TP was presented as means and standard deviations. Chi-square test was used to compare the demographic characteristics between NCSD and CSD groups. The *t*-test and ANOVA were used to compare the TP value between or among different variable groups. Pearson correlation analysis was used to measure the correlation between age and TP in both groups. To determine factors involved in TP, multivariate linear regression analysis was used. In addition, in order to obtain the best explanatory models of age changes on decreased TP, univariate and multivariate logistic regression analysis was used. In all analyses, the significance level was set at 5% and the confidence interval was The participants' demographic characteristics, questionnaire data, and TP were entered into Microsoft Excel (Microsoft, Redmond, WA, USA). IBM SPSS Statistics 20.0 (IBM, Armonk, NY, USA), which was then used to perform descriptive and inferential statistical analysis. Demographic characteristics were expressed as numbers and percentages, and TP was presented as means and standard deviations. Chi-square test was used to compare the demographic characteristics between NCSD and CSD groups. The *t*-test and ANOVA were used to compare the TP value between or among different variable groups. Pearson correlation analysis was used to measure the correlation between age and TP in both groups. To determine factors involved in TP, multivariate linear regression analysis was used. In addition, in order to obtain the best explanatory models of age changes on decreased TP, univariate and multivariate logistic regression analysis was used. In all analyses, the significance level was set at 5% and the confidence interval was 95%.

#### 95%. **3. Results**

**3. Results**  There were 63.10% of older adults who had CSD. The distribution of CSD had a statistically significant increase with age, lower education level, more chronic diseases, fewer There were 63.10% of older adults who had CSD. The distribution of CSD had a statistically significant increase with age, lower education level, more chronic diseases, fewer teeth, lower RSST times, and a higher score of EAT-10 (all *p* < 0.001) (Table 1).

teeth, lower RSST times, and a higher score of EAT-10 (all *p* < 0.001) (Table 1).


**Table 1.** General characteristic of participants.

yrs: years old.

Both males and females aged 60–64 years old had the highest average TP compared with other age groups, as shown in Figure 2. The average TP was 30.63 ± 11.77 kPa in the NCSD group, and this was statistically significantly higher than the CSD group (25.17 ± 13.61 kPa) (*p* < 0.001) (Table 2). The NSCD and CSD participants whose education level was university or above, had the highest TP (34.67 kPa and 29.37 kPa) compared with those who had other education levels (all *p* = 0.001) (Figure 3). The TP value of CSD had a positive tendency with the number of functional teeth, but did not show any statistically significant difference. In the NCSD and CSD groups, TP statistically significantly decreased with age increase (all *p* < 0.001). The scatter plots and bar charts revealed there was a negative significant relationship between TP and age in both the NCSD and CSD groups (all *p* < 0.001) (Figure 4).


**Table 2.** Tongue pressure of participants.

yrs: years old; CSD: chewing or swallowing difficulties; <sup>a</sup> : Data are expressed as interquartile range.

**Figure 2.** The mean maximum tongue pressure of (**a**) male and (**b**) females among age groups. Each bar showed the mean value and standard error. \* *p* < 0.05, \*\* *p* < 0.01. N.S.: non-significance.

≥3 151 22.93 13.93 22.93 13.93

≥3 800 22.93 13.93 30.63 11.77 25.62 13.02 yrs: years old; CSD: chewing or swallowing difficulties; a: Data are expressed as interquartile range.

**Figure 3.** Comparison of maximum tongue pressure of participants without/with chewing and/or swallowing difficulties by different education levels. Each bar shows the mean value and standard error. \* *p* < 0.05, \*\* *p* < 0.01. **Figure 3.** Comparison of maximum tongue pressure of participants without/with chewing and/or swallowing difficulties by different education levels. Each bar shows the mean value and standard error. \* *p* < 0.05, \*\* *p* < 0.01. *J. Pers. Med.* **2021**, *11*, x FOR PEER REVIEW 11 of 18

**Figure 4.** Correlation between maximum tongue pressure and age of both participants (**a**) without and (**b**) with chewing and/or swallowing difficulties. Tongue pressure of participants (**c**) without and (**d**) with chewing and/or swallowing difficulties by age groups. Each bar shows the mean value and standard error. \* *p* < 0.05, \*\* *p* < 0.01. N.S.: non-significance. **Figure 4.** Correlation between maximum tongue pressure and age of both participants (**a**) without and (**b**) with chewing and/or swallowing difficulties. Tongue pressure of participants (**c**) without and (**d**) with chewing and/or swallowing difficulties by age groups. Each bar shows the mean value and standard error. \* *p* < 0.05, \*\* *p* < 0.01. N.S.: non-significance.

In order to better understand the participants, we further compared the dietary habits between the NCSD and CSD groups, as shown in Table 3. Higher TP values of CSD

TP decreased, CSD participants obviously felt that they needed to take longer to eat a meal

than before (*p* = 0.006).

In order to better understand the participants, we further compared the dietary habits between the NCSD and CSD groups, as shown in Table 3. Higher TP values of CSD participants had a statistically significant benefit with their independent ability to eat, various types of diet, and shorter duration of meal time (all *p* < 0.050) while eating a meal. As TP decreased, CSD participants obviously felt that they needed to take longer to eat a meal than before (*p* = 0.006).



CSD: chewing and/or swallowing difficulty.

Multiple linear regression models were performed to clarify the risk factors of decreased TP among different target groups (Table 4). In the three regression models, we found the only common risk factor in the total population and for NCSD and CSD groups was age level. Other risk factors which affected TP were lower education level, fewer FT and a higher EAT-10 score.

Finally, we reconfirmed the most important risk factor for decreased TP among older adults was age level through multiple logistic regression models (Table 5). Whether the older adults had CSD or not, the odds ratio of decreased TP which was 1.54 (AOR = 1.54, 95% CI: 0.98–2.43, *p* = 0.062) in the 65–69 years old group increased to 6.61 (AOR = 6.61, 95% CI: 3.58–12.57, *p* < 0.001) in the 85 years old or above group compared with the 60–65 years old group.


#### **Table 4.** Factors affecting tongue pressure of participants.

yrs: years old; CSD: chewing and/or swallowing difficulties; CI: confidence interval. Adjusted gender.



yrs: years old; CSD: chewing and/or swallowing difficulties; CI: confidence interval. a COR = Crude Odds Ratio. Data analysis by simple logistic regression model. b AOR = Adjusted Odds Ratio. Data analysis by multiple logistic regression models. Adjusted gender.

#### **4. Discussion**

Personalized TP training or rehabilitation therapy is a new measure to improve swallowing with a basis of customized design. In the present study, we successfully established a large amount of TP data in NCSD (healthy) and CSD adults over the age of 60 by a quantitative method. The results of this study show evidence that TP declines are in parallel for NSCD and CSD adults at every 5-year increase in age. But the TP of CSD adults saw an accelerated decline from 60 years old (31.76 kPa) until 85 years or older (18.20 kPa) than NCSD adults (33.56 kPa and 24.51 kPa). Furthermore, this study established the standard and normal reference ranges of TP values at 60 years (31.47 kPa), 70 years (26.22 kPa), and 80 years and above (20.83 kPa) of overall participants. If an elder's TP value is less than 30 kPa, 25 kPa, and 20 kPa, respectively, he or she may have CSD. These three TP values can also be used as reference indicators for community screening, clinical diagnosis, and treatment effectiveness of CSD. In the future, the elderly can refer to the TP values currently found under parameters such as age, gender, education level, CSD status, etc., to set personalized training goals for rehabilitation.

Besides the 60–64 and 75–79 years old groups, there was no difference in TP values between genders of most age groups in the present study, as seen in several previous studies [6,17,33,34]. The TP values in males was rather gradual with a slightly higher value (within 2.5 kPa) than in females in most age groups. A research reviewed numerous studies and concluded that the impact of TP in healthy older adults between different genders only exists in the participants who are less than 60 years old [18]. From the present study results, we can clearly see that the TP of male older adults at 60–64 years and 75–79 years are significantly higher than females. This may because the height of a male's anterior tongue have already decreases by their sixties [34]. Then, a man may experience greater loss of total muscle mass at this time, accompanied with a rapid decline in muscle strength at age 75 or older, until they finally reach the same level of pressure as females [5,17].

This study clearly demonstrates that TP decreased with natural aging in older adults. Past studies showed that mean TP values in healthy individuals at an age less than 60 years was 39.3 kPa, which was higher than those at age 60 years or above (30.3 kPa) [18]. Furthermore, we compared our TP values with several studies of similar age composition. These studies included the participants' age between 60 years to 95 years and the mean age of participants was between 70–75 years. We found that our mean TP conducted by TP wireless application (30.63 kPa) was similar with previous studies which were conducted using JMS (29.5 kPa to 31.8 kPa) [15,17,18,35,36]. This result confirmed that our TP data have high reference value. Hara et al. reported that the TP values of healthy people in their forties, fifties, sixties, seventies, and eighties were 36.48 kPa, 34.59 kPa, 33.16 kPa, 29.92 kPa, and 23.60 kPa, respectively [37]. From the above data, we observed that the TP values in 60–64 year olds was lower than those less than 60 years old. We also observed that our TP value for 60–64 year olds (33.56 kPa) was very close to the relevant study of Hara et al. (33.16 kPa). We confirmed that this study's results provide quantitative evidence of the age-related changes in TP during the aging process.

CSD is another risk factor of rapid TP decline under natural aging. From 75 years of age, TP continued to rapidly decrease with larger differences with increasing age between the NCSD and CSD groups until the later stage of aging. We observed that the TP decline rate from 60–64 years old to 80 years old or older of CSD (42.70%) participants were 1.5 times that of NCSD participants (26.97%). Wang et al. reported that swallowing function decreases significantly above the age of 75 years old in Taiwan [8], supporting the results of the present study. However, the former TP studies showed no evidence of independent differences in CSD attributable to age due to only a small sample size of the elderly group evaluated or a lack of clear health conditions of the participants [6,18,38]. We collected chronic disease, teeth number, swallowing, and eating status of each participant through a large sample size and excluded the relationship between diseases and TP during data analysis. Unquestionably, CSD accelerates the decline rate of TP. The study provides a supplement to explore the interaction between age and swallowing difficulty in TP.

This study observed that decreased TP strongly affects the eating status of CSD. The elderly who have swallowing difficulties commonly encounter swallowing muscle fatigue, longer mealtimes, low food consumption, and malnutrition, supporting the results of the present study [36,39]. A series of negative effects during the eating process may originate from weak TP. A low TP makes normal swallowing hard to achieve. As noted in previous studies, the percentage of maximum TP used during swallowing in older adults was 53.8%, which was higher than young adults at 38.8%. This could be regarded as an indispensable compensatory mechanism to complete safe swallowing [7] because posterior TP is needed to generate propulsion to transfer solid bolus into the pharynx [12]. Based on the current findings, order adults with CSD had lower TP than those without CSD and a decreasing TP is linked to CSD. This finding reveals that a lack of adequate tongue muscle strength is one of the eating disorders in older adults with CSD.

There are some limitations that need to be noted for the current study. First, we cannot well observe the impact of disease on TP because those who were actively willing to participate in this study may tend to be healthier than those who did not participate. Due to an inability to reach the complete population of elderly, the results of the current study may be overestimated. Second, the presence of fewer males than females may elicit selection bias. However, the number of each age group in this study were larger than in previous studies, and we believe that the results of the study can compensate for part of the sample bias.

#### **5. Conclusions**

TP decline is prevalent among community-dwelling older adults, and aging is the major factor in TP decrease in parallel with ordinary and CSD adults. Older adults with CSD have an accelerated decrease in their TP, especially when they are 75 years old or above. This study successfully established a series of TP values in NCSD (from 33.56 to 24.51 kPa) and CSD (from 31.76 to 18.20 kPa) in Taiwan and suggests that TP can be used as a reference indicator of aging progression among community-dwelling older adults. The authors have confirmed that all authors meet the Journal of Personalized Medicine criteria for authorship credit, as follows: (1) substantial contributions to conception and de-sign, acquisition of data, or analysis and interpretation of data; (2) drafting the article or revising it critically for important intellectual content; and (3) final approval of the version to be published.

**Author Contributions:** Conceptualization, H.-Y.L., C.-T.Y. and C.-L.L.; methodology, H.-Y.L.; software, H.-Y.L.; validation, P.-H.C.; formal analysis, K.-J.H.; investigation, H.-Y.L., J.-H.C., K.-J.H. and C.-T.Y.; resources, H.-Y.L.; data curation, J.-H.C. and C.-L.L.; writing—original draft preparation, J.-H.C., K.-J.H. and S.-Y.H.; writing—review and editing, H.-Y.L. and C.-L.L.; visualization, C.-T.Y. and P.-H.C.; supervision, C.-L.L.; project administration, H.-Y.L. and C.-L.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ministry of Science and Technology, Taiwan (grant number: MOST108-2218-E- 010-005/MOST109-2224-E- 010-001).

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no potential conflict of interest with respect to the research, authorship, and/or publication of this article.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Journal of Personalized Medicine* Editorial Office E-mail: jpm@mdpi.com www.mdpi.com/journal/jpm

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-0365-4207-2