The study was carried out in E-R and Veneto, two neighboring regions of northern Italy with a population of approximately 4.4 and 4.9 million residents, respectively. Both these two regions oversee a regional healthcare system and have exclusive competence in regulating, financing, and organizing healthcare services and activities carried out within their territory. Data were extracted from the E-R and Veneto Regions healthcare administrative databases. They include comprehensive information about healthcare provision by the regional healthcare systems. Secure record-linkage procedures were carried out at the individual level to merge pseudonymized data related to: official notifications of SARS-CoV-2 infections; drug prescriptions; outpatient care; residence and vital status; acute hospital admissions; community hospital admissions; emergency room access; long-term care facilities; and integrated home care. In the E-R cohort, an individual risk of hospitalization and death score was also assigned using a previously developed standardized algorithm. This algorithm relies on a multivariable prediction model which estimates a punctual measure for the individual risk of hospitalization or death within the reference year. The risk of hospitalization and death is then scored according to a four-level scale: low (probability <6%), moderate (≥ 6% and < 15%), high (≥ 15% and < 25%), very high (≥ 25%). This score is assigned yearly based on demographic and residence characteristics, comorbidities, and access to a wide spectrum of healthcare resources in a multiyear period before the reference year, and is routinely available in E-R administrative databases (20). The extracted data include updates of databases up to November 2021 for E-R and up to December 2021 for Veneto.

The study is a multicenter retrospective population-based cohort study aiming to investigate the incidence and risk factors of access to healthcare services. It focused on a largely healthy population, at low risk of acute severe COVID-19, to ensure that a high fraction of outcomes is attributable to COVID-19. Eligible participants included all adult subjects aged ≥18 years at diagnosis and with continuous residence status in the two regions in the 365 days before diagnosis. Low risk of severe acute COVID-19 was characterized by the absence of all the following types of care, in the 365 days prior to the SARS-CoV-2 diagnosis: hospitalization; visits to the emergency room; integrated home care; residence in a long-term care facility; selected drug prescriptions (at least one drug within the selected list; see Supplementary Table 1); selected outpatient care (at least one visit/diagnostics within the selected list; see Supplementary Table 1). Additionally, in the E-R cohort, low risk of severe acute COVID-19 was assigned only to subjects with a low-to-moderate risk of hospitalization and death score (19). The population of interest included all consecutive adult individuals with confirmed SARS-CoV-2 infection (PCR or antigen tests) in the E-R and Veneto Regions between February 2020 and November 2020 (E-R) or December 2020 (Veneto) who, at the time of diagnosis, were at low risk of severe acute COVID-19 disease. None of these subjects was vaccinated at the time of SARS-CoV-2 diagnosis, as the vaccination campaign in Italy only started in late December 2020. Individuals entered the cohort on the day of SARS-CoV-2 diagnosis. Outcomes were assessed within 365 days before and after the diagnosis. In the follow-up period we distinguished an acute phase (AP) from diagnosis till day 30 and a post-acute phase (PAP) from day 31 till day 365. For patients who were hospitalized at day 30, PAP started on the first day after hospital discharge. Outcomes were assessed during the PAP. Individuals who were continuously hospitalized from the AP to more than 365 days after SARS-CoV-2 diagnosis, or who had moved their residence outside the region during the AP, were excluded. Duplicate or incomplete records were discarded. The study was carried out and reported according to the GATHER statement (21).

The outcome of interest was the access to specific healthcare services during the PAP, defined through a selected list of drug prescriptions and outpatient visits and diagnostics. This combined outcome was considered a proxy for the long-term effects of COVID-19 requiring medical attention (22). A broad range of drug prescriptions and outpatient care services was considered, to reflect the multifactorial nature of COVID-19, which can affect several organs and systems. Selected drug prescriptions included cardiovascular system (e.g., antithrombotics, antiarrhythmics, antihypertensives, beta blockers), antidiabetic, nervous system (e.g., antidepressants), and respiratory system (e.g., adrenergics and other drugs for obstructive airway diseases) drugs, oxygen and corticosteroids (22–24). Selected outpatient visits and diagnostics included ambulatory visits in cardiology, pneumology, angiology, neurology, psychiatry, rehabilitation-motor, nephrology, and diabetes, as well as other diagnostic and therapeutic procedures such as chest imaging, cardiac ultrasound imaging, pneumological diagnostics, electrocardiography, oxygen therapy, respiratory, and cardiological rehabilitation, peripheral vascular ultrasound imaging, training for cognitive disorders, hemodialysis, renal imaging, and glycated hemoglobin analysis (22–24). Only drugs and outpatient care provided by the regional healthcare system were included in the analysis. Drugs that are not reimbursed as well as private outpatient care not provided by the regional healthcare systems are not recorded. Details of drugs and outpatient care are reported in Supplementary Table 1, together with extraction criteria based on ATC classification and on E-R regional outpatient codes. Two analyses were performed, each with a different definition of the outcome variable. In the first analysis we considered a combined outcome including selected drug prescriptions and selected outpatient care, whichever came first. In the second analysis, each of the two outcomes of interest (drug prescriptions and outpatient care) was analyzed separately and the other was considered neither as a competing event nor as censoring. Hospitalization and death were always considered as competing events following a competing risk analysis framework. Hospitalization was considered a competing event due to drug prescriptions and clinical consulting not being tracked at the individual level during acute care stay: This prevented us from observing these events in administrative databases during the period a patient was hospitalized.

Some clinically relevant variables were determined based on data recorded in administrative databases. COVID-19 severity (expressed on a four-level ordinal scale as low, mild, moderate or severe) was algorithmically assigned. The algorithm was based on respiratory system diagnoses (i.e., acute respiratory insufficiency, pneumonia, acute lower respiratory tract infections, other respiratory diagnoses), on ventilation procedures administered (i.e., oxygen therapy, non-invasive ventilation, invasive ventilation), and on intensive or sub-intensive care unit stay during hospitalizations in the AP. For subjects not hospitalized in the AP, the “low” level of severity was assigned. The algorithm is reported in Figure 1. Acute extrapulmonary complications (i.e., vascular, hemorrhagic, thrombotic, cardiac, neurological, septicemia and acute organ failure) were also identified, based on diagnoses during hospitalizations in the AP. These variables can occur only for hospitalized subjects. Criteria for the identification of respiratory diagnoses, ventilation procedures and acute extrapulmonary Complications, based on ICD-9-CM codes, are reported in Supplementary Table 2.

The statistical analysis was carried out separately in the two Regions, and results were pooled using meta-analysis. The frequency distributions of the categorical characteristics were described as absolute and percentage numbers. The numerical variables were described as the mean ± standard deviation and range. To evaluate potential risk factors, multivariable regression models were carried out. The incidence of outcomes in the PAP was analyzed with a competing risk approach: selected drug prescriptions and/or selected outpatient care were the outcomes of interest, whereas hospitalization and death were the competing events. Individuals who moved their residence outside of the E-R region during the PAP, as well as those who did not experience any outcome or competing event by the end of the PAP, were treated as having censored follow-up times. The incidence of outcomes over time was described using cumulative incidence functions curves (25). Uncertainty in curves was expressed with 95% confidence interval (CI) calculated with the asymptotic Aalen method (26). A multivariable Fine-Gray (FG) proportional subdistribution hazard regression model was used to assess the relationship between subjects’ characteristics and the hazard of outcomes over time, also accounting for the occurrence of competing events (27). The explanatory variables were: time period (1st wave: February–May 2020, intermediate period: June–September 2020, 2nd wave: October–December 2020); age class at diagnosis (18–39, 40–49, 50–59, 60–69, 70–79, ≥80 years); sex (male, female); Italian citizenship (yes, no); only in the E-R cohort, the risk of hospitalization and death score (low, moderate); COVID-19 severity during the AP (low, mild, moderate, severe); acute extrapulmonary complications occurring in hospital during the AP (yes, no). Moreover, the area of residence (8 Local Health Units in E-R and 9 in Veneto) was also used as an independent variable to account for potentially different health policies. Associations were measured using the subdistribution hazard ratio (HR) and the uncertainty in results was expressed with 95% CI. Cis for HRs were calculated with the Wald method based on normal approximation. Pooling of hazard ratios obtained in the E-R and Veneto Regions’ cohorts was carried out using random effects meta-analysis (MA). MA was performed with the inverse variance weights method and maximum likelihood estimator for between-study variance. Between-cohort heterogeneity was measured with the tau statistic (28) and its significance was assessed with the Cochran’s Q test. All statistical tests were two-sided. Analyses were carried out by E-R with SAS/STAT 15.1 (SAS Institute Inc., Cary, NC) and R 4.0.4 (The R Foundation for Statistical Computing, Wien) and by Veneto with SAS/STAT 13.1 statistics software.

Descriptive characteristics of the cohort are reported in Table 1. Most subjects were diagnosed for SARS-CoV-2 in Emilia-Romagna in October–November 2020 (80.3%) and in Veneto in October–December 2020 (91.8%), during the second epidemic wave in Italy. The average age at diagnosis was 41.1 ± 14.1 years (range: 18–91) in Emilia-Romagna and 41.8 ± 14.1 years (range: 18–97). Those aged more than 60 years were a small part (9.1 and 10.0%, respectively). About half of the included individuals were female (50.3 and 49.9%) and the most part had Italian citizenship (88.7 and 89.6%). At baseline, the risk of hospitalization and death score estimated in E-R was low for 98.7% and moderate for 1.3% of subjects. During the AP, 96.1% of low-risk individuals in Emilia-Romagna and 98.7% in Veneto experienced low severity COVID-19. Conversely, 0.2% and less than 0.1% of subjects had mild severity, 1.4 and 0.6% had moderate severity, and 2.2 and 0.7% experienced severe COVID-19 disease. Those who were hospitalized during the AP (in acute care hospitals or community hospitals) were 4.3 and 1.6%. Acute extrapulmonary complications have occurred in 0.3 and 0.2% of individuals (Table 1). The higher frequency of individuals who were not hospitalized and who were assigned the lowest severity level in the Veneto cohort likely reflects the wider use that was made of diagnostic testing in 2020 in that Region to detect positive cases even among asymptomatic people (29).

The total follow-up time in the PAP was equal to 26,826.7 person-years in E-R and 68,333.1 in Veneto. During this time, 9,208 (26.2%) of low-risk individuals in E-R and 19,769 (22.2%) in Veneto experienced at least one outcome of interest during the PAP. 4,633 (13.2%) and 9,409 (10.6%) subjects, respectively, had at least one selected drug prescription. On the other hand, 6,523 (18.6%) and 14,108 (15.9%) individuals had at least one selected outpatient care visit or diagnostic procedure, respectively. In the competing risks analysis, as reported in Figure 3, the cumulative incidence of the combined outcome at 11 months in the PAP was equal to 24.9% (95% CI = 24.5–25.4%) in E-R and to 21.2% (95% CI = 20.9–21.5%) in Veneto. Drug prescriptions were less frequent than outpatient care: at 11 months in the PAP, the incidence of the former was 11.9% (95% CI = 11.5–12.2%) in E-R and 9.5% (95% CI = 9.3–9.7%) in Veneto, whereas the incidence of the latter was 17.9% (95% CI = 17.5–18.3%) and 15.4% (95% CI = 15.1–15.6%), respectively. Monthly outcomes incidence data are reported in Supplementary Table 3. Competing events such as hospitalization and mortality occurred during the PAP in 1,292 (3.7%) and 26 (0.1%) individuals in E-R, and in 2,837 (3.2%) and 80 (0.1%) individuals in Veneto. The frequency of individuals who accessed healthcare services during the PAP is reported in detail in Table 2. Among selected drug prescriptions, cardiovascular system drugs (6.0% in E-R and 4.7% in Veneto), followed by corticosteroids for systemic use (4.8 and 3.9%) and respiratory system drugs (2.3 and 1.8%) were administered more frequently during the PAP. Among selected outpatient ambulatory services, the most frequent during the PAP were cardio-respiratory visits and procedures (11.6% in E-R and 9.1% in Veneto), followed by diabetic ambulatory visits or procedures (3.9 and 4.3%), and rehabilitation-motor visits (2.0 and 2.3%).

According to the confounder-adjusted pooled analysis reported in Figure 4, the major risk factor was the level of COVID-19 severity during the AP. Those with mild severity had +174% hazard of combined outcome compared to those with low severity (HR = 2.74, 95% CI = 1.49–5.02), whereas those with intermediate severity had +260% hazard (HR = 3.60, 95%CI = 2.41–5.36). Finally, those with severe COVID-19 had +321% hazard (HR = 4.21, 95%CI = 3.22–5.48). Such differences among severity groups were not homogeneous in the two regional cohorts, being HRs for Veneto higher than those for E-R. The reason for such a heterogeneity is shown in Figure 3, which describes cumulative incidence curves by severity groups in the two regional cohorts. In the E-R cohort, the risk of the combined outcome at the end of the PAP was 23.5% (95% CI = 23.0–23.9%) for low severity patients, 37.2% (95% CI = 26.6–47.4%) for mild-severity ones, 56.9% (95% CI = 52.2–61.1%) for moderate-severity ones and 66.5% (95% CI = 63.0–69.7%) for severe patients. In the Veneto cohort, the same figures were equal to 20.5% (95% CI = 20.2–20.7%), 59.5% (95% CI = 41.6–73.5%), 67.8% (95% CI = 63.6–71.7%) and 74.3% (95% CI = 70.6–77.6%), respectively. The occurrence of an acute extrapulmonary complication during the AP was associated with a higher risk (HR = 1.84, 95% CI = 1.44–2.35). Age at diagnosis was also a risk factor for the combined outcome, as the risk in higher age groups was always greater than or equal to the risk in the group of individuals aged 18–39. In particular, those aged 60–69 or 70–79 had more than a two-fold hazard (HR = 2.10, 95%CI = 2.01–2.19 and HR = 2.31, 95%CI = 1.96–2.73, respectively) and those aged ≥80 had less than a two-fold hazard (HR = 1.73, 95%CI = 1.15–2.60). Male individuals had −18% hazard of outcome compared to females (HR = 0.82, 95%CI = 0.80–0.84). Italians had a slightly higher level of risk compared to non-Italian citizens (HR = 1.09, 95%CI = 1.05–1.14). Finally, there were only minor differences between individuals diagnosed in the second epidemic wave (HR = 0.92, 95%CI = 0.88–0.97) or between individuals diagnosed in the intermediate period (HR = 0.81, 95%CI = 0.73–0.90) and those diagnosed in the first epidemic wave. The analysis of risk factors for each separate outcome (drug prescriptions or outpatient care), is reported in Figures 5, 6. Overall, the risk factors were similar to those already described for the combined outcome, with the following three exceptions. Firstly, COVID-19 severity in the AP was associated with a higher increase in the risk of outpatient care, compared to the increase in the risk of drug prescriptions. Secondly, individuals with Italian citizenship had a lower hazard of drug prescriptions (HR = 0.90, 95%CI = 0.85–0.96), whereas the opposite was observed in relation to outpatient care (HR = 1.14, 95%CI = 1.09–1.20). Thirdly, the presence of acute extrapulmonary complications during the AP was strongly associated with drug prescriptions (HR = 3.83, 95%CI = 2.96–4.95), whereas its relationship with outpatient care was of minor relevance (HR = 1.36, 95%CI = 1.10–1.66). Other minor differences in risk factor intensity were present, although they did not alter the overall interpretation of the results (Figures 5, 6).

The present study focused on the incidence and determinants of access to selected healthcare services in a largely healthy population of subjects who had a diagnosis of SARS-CoV-2 in the Emilia-Romagna and Veneto regions, Italy. The analysis covers people diagnosed in the period from February to December 2020, when vaccination against COVID-19 was not yet available. The outcome (combination of outpatient care and drug prescriptions) occurred in a high percentage of subjects (25% cumulative incidence in E-R and 21% in Veneto) during the PAP, i.e., from day 31 to day 365 after diagnosis of SARS-CoV-2. It should be however noted that the incidence in a health population in the absence of COVID-19 was not measured, leaving uncertainty about the proportion of outcomes directly attributable to the disease. Considering the two outcomes separately, outpatient care was observed more frequently than drug prescriptions (18% vs. 12% in E-R and 15% vs. 10% in Veneto). The most frequently administered drugs were cardiovascular system drugs and corticosteroids, whereas the most common category of outpatient care was the cardio-respiratory one. Hospitalization occurred more rarely (about 3%). The cumulative incidence curve grows steadily for the combined outcome and for the outpatient care outcome. For drug prescriptions, growth appears steeper in the first month of follow-up in the PAP. The curves show significant differences if the analysis is stratified by severity of acute COVID-19. In particular, subjects with severe forms showed a much steeper increase in cumulative incidence in the first part of the follow-up, especially for drug prescriptions. Risk factors associated with the combined outcome were female sex, age over 40 years, moderate to severe acute COVID-19, and acute extrapulmonary complications occurring during the AP. Although there are differences, the risk factors associated with outpatient care and drug prescriptions (considered separately) are the same as those associated with the combined outcome. Having Italian citizenship is an exception, as it was a risk factor for the prescription of drugs and a protection factor for outpatient visits in the follow-up. This finding may be related to a different use of private outpatient care. Acute COVID-19 extrapulmonary complications is another exception, as it was strongly associated with drug prescriptions, but not with outpatient care.

The results of the two cohorts, although providing quite consistent results, show a greater frequency of outcomes in E-R than in Veneto Region. The meta-analysis of the factors associated with the outcome shows heterogeneity which is high for some variables such as the severity of the acute COVID-19. In particular, the negative effect of increasing COVID-19 severity on the need for outpatient care was more intense in Veneto than in E-R. These findings can be explained by the differences in diagnostic testing policies in place in the two Regions in early 2020, when Veneto Region was the first to implement intensive diagnostic testing and contact tracing procedures, leading to a higher share of asymptomatic individuals as opposed to paucisymptomatic or weakly symptomatic ones (29). The lower incidence of outcomes observed in the Veneto cohort for subjects classified at low severity is consistent with this hypothesis. Furthermore, the cumulative incidence for mild, moderate, and severe patients was higher in Veneto than in E-R. This may depend on different long-term case management policies for hospitalized subjects in the two Regions.

The study results showed a remarkable frequency of outpatient care and drug prescriptions in the post-acute follow-up period. These data suggest that long-term effects of COVID-19 needing clinical attention occur in a substantial percentage of cases, even among a previously healthy population with low or mild severity of acute COVID-19. The latter finding may be very important for clinicians and for healthcare policy makers as, thus far, an active follow-up has often been restricted to COVID-19 patients with moderate or severe infection (3). Our data may indeed suggest the need for a greater clinical attention in the follow-up after SARS-CoV-2 infection, also in a non-hospitalized low-risk population, as other recent reports have suggested (4, 18, 23). Furthermore, in subjects with severe acute COVID-19, a high frequency of outpatient visits and prescriptions is observed in the first part of the follow-up, suggesting that there may be a continuation of the AP (defined by NICE as ongoing COVID). Our results also contribute to the evidence on risk factors for the post-COVID syndrome and are in line with previous studies carried out in other populations of COVID-19 patients (3–6, 18, 23).

Based on the results of our study, the access to outpatient care and drug prescriptions was very frequent even in a low-risk population. The significant healthcare resource consumption related to such outpatient care and administration of drugs should therefore be considered by researchers when evaluating the healthcare burden after a SARS-CoV-2 diagnosis. Moreover, our data can be useful for setting benchmarks in the levels of access to selected healthcare services, as they are referred to a large unvaccinated population at low risk of severe acute COVID-19.

This study used healthcare administrative databases as data sources; therefore, only drug prescriptions and outpatient care provided by the regional healthcare systems were included in the analysis. In addition, there are no data on diagnoses of outpatient care and indications for drug prescriptions. The occurrence of selected drug prescriptions and selected outpatient care was indeed used as the outcome variable and considered as a proxy for long COVID. Therefore, the present analysis was only able to assess a variation in prescriptions and outpatient services in general terms, whereas no information on the incidence and risk factors of single drugs and outpatient services can be derived. Secondly, our results refer to unvaccinated subjects, but the epidemiological features of long COVID have changed significantly since the introduction of vaccines. Thirdly, unmeasured characteristics (e.g., environmental, lifestyle and genetic characteristics, SARS-CoV-2 variants) may have affected our results. Furthermore, outpatient care could result to some extent from the implementation of algorithms for the patient’s follow-up, regardless of persistence or recurrence of symptoms or the emergence of new ones. This limitation could have caused an overestimation of the incidence of study outcomes, particularly those related to outpatient diagnostics. Another potential limitation is in the criteria for the selection of a healthy population at low risk of severe acute COVID-19 from administrative databases. It is indeed possible that some non-healthy subjects, especially those with chronic diseases of minor clinical relevance (e.g., asthma), were included in the cohort due to the absence of relevant access to healthcare services in the one-year period before SARS-CoV-2 diagnosis. Similarly, the criteria for acute COVID-19 severity, being based on information available only for hospitalized subjects, may have been influenced by different health policies and hospital capacity during different pandemic periods. Finally, no control group of subjects who were not diagnosed with SARS-CoV-2 infection was considered, leaving uncertainty on the incidence and risk factors of outcomes directly attributable to COVID-19.

The study data sources allowed for a multicenter population-based design (all cases of COVID-19 that met the inclusion criteria in two Italian Regions) and for a long follow-up period (one year from the diagnosis of SARS-CoV-2 infection). These are strengths compared to most published studies which present one or more limitations reducing the validity and generalizability of the results (3.17). The most frequent flaws of published studies refer to small or selected samples (e.g., focus on hospitalized or more symptomatic patients) and short follow-up periods (e.g., less than 12–24 weeks). Other limitations of the available studies relate to design (e.g., in surveys, people who are still unwell or who have had a long-lasting illness are more likely to participate and recall symptoms) and case definition (e.g., the use of serology as an inclusion test makes the dating of the infection inaccurate) (3.17). Moreover, the selection of low-risk subjects (i.e., no hospitalization, visit to the emergency room, prescription of specific drugs or specific outpatient visits in the year preceding the diagnosis of SARS-CoV-2 infection) makes the observed outcomes in the follow-up period largely attributable to the post-COVID syndrome.