Differentiation of industrial hemp strains by their cannabinoid and phenolic compounds using LC × LC-HRMS

Cannabis is an ancient plant that has been used for therapeutic and recreational purposes. Nowadays, industrial hemp, a variety with low concentration of the psychoactive cannabinoid Δ9-tetrahydrocannabinol (THC) and high concentration of non-psychoactive cannabinoids, is getting more and more interest in the food, pharmaceutical, and cosmetic industry. However, cannabis not only contains cannabinoids as bioactive components but also other metabolites like terpenes and phenolic compounds, and the content of these interesting secondary metabolites greatly differs with the genetic variety of the plant. Due to the huge complexity of composition of the cannabis matrix, in this work, a comprehensive two-dimensional liquid chromatography (LC × LC) method has been developed as a very power separation technique coupling a pentafluorophenyl (PFP) and a C18 in the first and second dimensions. Two industrial hemp strains (cookie and gelato) were analyzed to determine the difference in their content of cannabinoids and phenolic compounds. To do this, a new demodulation process was applied for the first time to transform 2D raw data into 1D data which allowed carrying out the chemometric analysis needed to determine the statistical differences between the hemp strains. The cookie strain presented a total of 41 cannabinoid markers, while the gelato strain presented more representative phenolic compounds, in total 24 phenolic compounds were detected as potential markers of this sample. These differences in the chemical composition could determine the industrial destiny of the different hemp strains. Graphical abstract Supplementary Information The online version contains supplementary material available at 10.1007/s00216-022-03925-8.

According to the European Industrial Hemp Association, marihuana contains 1-20% of THC, while industrial hemp does not exceed 0.2% of the psychoactive cannabinoid [2,4]. Industrial hemp is currently considered as a balanced and complete food with health-promoting effects [2,15,16]. Regarding the content of cannabinoids in industrial hemp, these plants present a high concentration of non-psychoactive cannabinoids such as the neutral forms, and the acidic forms like cannabidiolic acid (CBDA) and cannabichromenic acid (CBCA), among many others [2]. Recently, the number of foods, pharmaceutical, and cosmetic products based on industrial hemp that are rich in non-psychoactive cannabinoids has drastically increased [17].
Besides of cannabinoids, industrial hemp contains other important families of secondary plant metabolites like phenolic compounds, terpenes, amides, amines, or phytosterols [4]. All these compounds play an important role together with cannabinoids not only in the organoleptic properties but also in the food and pharmaceutical applications due to their potential bioactive effects. Regarding the phenolic compound content, more than 20 compounds present in cannabis have been detected as flavones, flavonols, and prenylated flavones, which have been related to several health therapeutical activities [4]. Moreover, the presence of phenolic compounds in hemp produces a synergic action over some biological activities related to cannabinoids, producing an enhancement of these interesting bioactivities [18].
However, the chemical composition of hemp cultivars greatly varies depending on multiple factors such as the harvested year, environmental conditions, or the different parts of the plants [2,15,19,20]. One of the most relevant factor that affects the nutritional and second metabolite profile of hemp plants is the genetic diversity between the plants [5,21,22]. Thus, each C. sativa and indica strain presents different chemical compositions and the analysis of the phytochemical profile of the different genetic varieties is important to establish the relationship between composition, medicinal effects, and industrial use [7,22].
The analysis of hemp extracts is a challenge since there are more than 550 compounds identified in cannabis plants belonging to different chemical families [6] which implies the need for different analytical platforms such as gas chromatography and liquid chromatography coupled to mass spectrometry (GC-MS and LC-MS) for the analysis of all the different components [18,23]. Besides, NMR [24] and two-dimensional comprehensive gas chromatography (GC × GC) [25] have been used for the analysis of hemp extracts. Liquid chromatography coupled to diode array detection (LC-DAD) and LC-MS methods have been successfully used for the separation and identification of target cannabinoids [7,22,23,26,27]. Berman

Cannabis sativa Cannabis indica
LC-MS methods have also been optimized for the analysis of the phenolic fraction of hemp inflorescences [28]. However, there are still some analytical challenges to obtain a complete separation of the whole hemp profile due to the high number of compounds that build the phytochemical composition. For this reason, analytical techniques with high separation power are required to create footprints of the different commercial varieties to establish food authenticity and industrial aims. In this work, a comprehensive two-dimensional liquid chromatography (LC × LC) method together with a new demodulation process that transforms 2D data in 1D data is proposed for the first time for the differentiation of two varieties of industrial hemp according to their cannabinoid and phenolic profile.

Samples and reagents
Two commercial products of dried hemp inflorescences were purchased from a company that produces hemp products. The products (known as cookie and gelato strains) were certified as industrial hemp with a THC content of less than 0.2%. Both of them were indica dominant hybrid (60% indica, 40% sativa) strains, differentiated by their genetic precedence. While the cookie strain came from the cross of OG Kush with Durban poison strains, the gelato sample was obtained by crossing the Sunset Sherbet with the thin Mint Girl Scout cookie strains.
All solvents used were LC-MS grade. Acetone was purchased from Sigma-Aldrich (Darmstadt, Germany), acetonitrile and methanol were acquired from VWR (Darmstadt, Germany), and formic acid was bought from Fisher Scientific (Schwerte, Germany). Ultrapure water (resistivity 18.2 M Ω cm −1 ) was obtained from a Sartorius Ultrapure Water System (Goettingen, Germany).

Sample preparation
For the extraction of cannabinoids and phenolic compounds, the dried hemp inflorescences were freeze-dried to remove any remaining water with a vacuum drier (Alpha 1-2 LDplus, Martin Christ, Osterode am Harz, Germany) for 16 h and then mortared into powder. 250 mg of each sample were extracted using 37.5 mL acetone/water 70:30 (v/v) as extraction solvent. The mixtures were vortexed for 3 min, sonicated with ultrasonic bath for 30 min, and then centrifuged (Centrifuge 5804R, Eppendorf, Hamburg, Germany) for 10 min (3000 rpm, 5 °C). After the centrifugation, the extract was evaporated under nitrogen stream to remove acetone and then freeze-dried (64 h) to eliminate the water. The powder of both samples was weighed and stored at 4 °C in the darkness until its analysis. Prior to analysis, the extracts were dissolved in water/methanol (50:50, v/v).

µLC × LC-DAD-qTOF MS analysis
The chemical characterization of the hemp extracts was carried out using a two-dimensional liquid chromatography (2DLC) system (Agilent, Walbronn, Germany) coupled with an Agilent 6545 QTOF-MS system (Agilent, Santa Clara, USA). The first dimension ( 1 D) was built with a 1260 Infinity HiP micro ALS autosampler module (G1377A), a 1260 Infinity capillary pump (G1376A), a 1260 Infinity column compartment (G1316A), and a 1260 Infinity DAD detector module (G1315C). The second dimension ( 2 D) was equipped with a 1290 Infinity II high-speed pump (G7120A) and a 1290 Infinity II DAD detector (G7117B). Moreover, a 1290 Infinity binary pump (G4220A) was used to create a make-up flow rate (additional pump). The coupling of the 1 D and 2 D was carried out by an automated controlled 2 ports/4-position dual valve (G1170A) equipped with two 40 µL sampling loops. The operation and control of the system were done using the program OpenLAB ChemStation Edition (Version C.01.07 SR3, Agilent, Santa Clara, USA).
The QTOF-MS system worked with an Agilent Dual Jet Stream ion source. The source conditions were as follows: nitrogen gas temperature 325 °C, drying gas 9 L min −1 , nebulizer 30 psi, sheath gas temperature 300 °C, sheath gas flow 10 L min −1 , VCap 3500 V, and nozzle voltage 750 V. For the MS qTOF parameters, the fragmentor was set at 380 V while the skimmer and Oct 1 RF Vpp were fixed at 30 V and 750 V, respectively. The samples were analyzed in both ionization modes. The mass range was from m/z 100 to 1700. Data-dependent analysis of the top 10 ions was carried out using 20 eV as collision energy. For the operation of the QTOF system and the data acquisition, MassHunter Workstation LC/MS Data Acquisition (Version B.09.00, Agilent, Santa Clara, USA) was used. For the 2D data visualization, LC Image software (Version 2.7r3.1 LC × LC, GC Image, Lincoln, USA) was employed. MassHunter Qualitative Analysis Navigator (Version B.08.00, Agilent, Santa Clara, USA) was used for the MS data analysis.

Data treatment
For the 2DLC data treatment, firstly, a feature list was created with the software MS-Dial 4.7 (http:// prime. psc. riken. jp/ compms/ msdial/ main. html). After that, in order to convert the 2D data in 1D data, the feature list was "demodulated" using a home-made program. This "demodulation" program applies an algorithm that recognized the retention times, m/z values, and intensities of the features. Giving the modulation time of the LC × LC analysis, the algorithm combines all the areas of a given m/z value that follow a Gaussian distribution along the modulations. Briefly, the program identifies the initial modulation of a peak when it recognizes a m/z value at certain intensity and combines all the areas of the following m/z values with higher intensity that eluted at exactly the given modulation time. Lastly, it recognizes the last modulation point of the modulated peak when the intensity of that m/z value decreases. At the end, the program provides the sum of the areas of the different modulated points for a given feature at the retention time where the maximum intensity was found for each m/z value. The program together with the source code and an in-depth explanation will be published separately in the future.
For the statistical analysis, the software Simca 16.0.2 (Sartorius Stedim Data Analytics AB, Umeå, Schweden) was used. Both non-supervised and supervised methods consisting of principle component analysis (PCA) and partial least squares discriminant analysis (PLS-DA) were performed. After that, a suspected targeted analysis of the data was carried out in order to identify cannabinoid-like and phenolic compounds. Finally, in order to distinguish both samples in terms of cannabinoids and phenolic compounds, a cluster was performed.

µLC × LC separation
2DLC provides the separation power that cannot be achieved by conventional 1DLC analysis. The reason for this high separation power is the possibility to analyze a sample by two separation mechanisms that present different selectivity for the analytes. Accordingly, the analytes that are not possible to be separated and coelute in the first column (or first dimension, 1 D) can be separated in the second column (second dimension, 2 D). Therefore, 2DLC is the analytical tool of choice for the analysis of very complex samples. In particular, LC × LC is the 2DLC mode preferred to do non-targeted analysis, since the complete sample is separated by both dimensions. In the last years, the application of LC × LC methods for the analysis of complex food and plant samples has greatly increased, showing the expansion on the use of this technique [29]. In this work, a µLC × LC method has been developed for the analysis of the extract of industrial hemp inflorescences [6]. The µLC × LC method was optimized for the separation of cannabinoids and phenolic compounds present in the industrial hemp.
For the 1 D, a µLC system was used to achieve high reproducible and robust gradients at low flow rates, typically used in 1 D [30]. Different column combinations were tested for the 2DLC separation of the hemp extract. The combination of HILIC × RP has been successfully used for the separation of very complex phenolic compound and other secondary metabolite mixtures [31][32][33]; therefore, firstly, hydrophilic interaction liquid chromatography (HILIC) column was tested as 1 D separation, although it was not possible to achieve an efficient separation to be modulated in the 2D system probably due to the non-polarity of the cannabinoids (data not shown). Then, reversed phase (RP) columns combined with 2 D HILIC stationary phases (RP × HILIC) like C18 × NH 2 , C18 × HILIC, or C18 × Cys were tested. Using C18 in the 1 D, a good separation was achieved for the hemp; however, the coupling of this C18 separation in the 1 D with HILIC produced a high breakthrough and a poor separation of the fractions transferred from the 1 D into the 2 D (Figure S1a-c). This effect is due to the high strength mismatch between the mobile phases of RP and HILIC (i.e., the weak solvent in RP is water, which is the strong solvent for HILIC). This fact together with the high sensitivity of HILIC to the injection solvent and the fast analysis carried out in the 2 D made that the analytes transferred to the 2 D diluted in the 1 D solvent were not focused and retained in the 2 D HILIC column. Finally, a RP × RP combination was checked. The main advantage of the RP × RP coupling is the good mobile phase compatibility between the two separation modes used in both dimensions. For this approach, the 1 D separation was carried out in a PFP column that presents good properties for the separation of aromatic substances like cannabinoids and phenolic compounds. On the other hand, a short C18 column performed the 2 D separation. C18 is the most common stationary phase used in the 2 D [34,35] due to its beneficial properties to perform fast analysis while keeping a high resolution, and parameters that are required in the 2 D to carry out the separation of each 1 D fraction before the next fraction are injected in the 2 D column. The flow rate used in the 1 D was optimized at 50 µL/min. Smaller flow rates were not optimal for the 1 D separation providing a high dead times and therefore long total analysis time as can be observed in Figure S2a-c. Higher flow rates than 50 µL/min were not tested to avoid the collection of large 1 D fractions that would disrupt the 2 D separation. Although 50 µL/min is not the optimal flow rate for columns with internal diameter of 2.1 mm, the reproducibility of the 1 D separation at this flow rate was very precise ( Figure S2d). On the other hand, a modulation time of 0.5 min was established in order to achieve a compromise between the undersampling effect and the minimum analysis time required to carry out the 2 D separation. The combination of two theoretically correlated separation mechanisms could drive to a non-orthogonal 2DLC separation and therefore, the separation of the compounds in the 2D space is limited to the diagonal. In fact, this effect occurred in the present PFP × C18 analysis before the optimization (Figure S1d), where the same 2 D gradient was used during the whole 2DLC analysis ( 2 D full gradient). To improve the orthogonality, different mobile phases were tested in both dimensions (methanol was used as organic solvent in PFP, and acetonitrile was used for the separation in the C18 2 D). Besides, the 2 D gradient was tailored according to the 1 D eluted fractions along the whole 2DLC analysis developing a 2 D multi-segment shifting gradient. The highly improved orthogonality obtained for the separation of the cookie and gelato hemps after the optimization of the tailored 2 D gradient can be observed in Fig. 2. To qualify the gain in the separation after the optimization and to quantify the effect of the 2 D multi-segment shifting gradient, the peak capacity and the orthogonality of both 2DLC methods, the 2 D full gradient (Figure S1d), and the 2 D shifting gradient (Fig. 2) were calculated. The practical peak capacity ( 2D n c,practical ) was calculated according to Li et al. [36], the orthogonality was estimated following the asterisk equations (A O ) [37], and the corrected peak capacity ( 2D n c,corr. ) was calculated multiplying the 2D n c,practical and the A O . For the 1 D, a peak width media of 1 and 0.8 min was calculated for the full and shifting gradient, respectively, while the peak width media of the 2 D was 3.3 and 0.4 s for the two respective settings. The peak capacities of the 1 D ( 1 n c ) for the full gradient and the shifting gradient methods were 59 and 75, respectively. On the other hand, the peak capacities of the 2 D ( 2 n c ) in the full gradient method was 10 and for the shifting gradient method was 62. Therefore, the 2D n c,practical values were 442 and 3080 for the 2 D full gradient and the 2 D shifting gradient, respectively. Regarding the orthogonality, a gain of 23% was obtained after the optimization of the 2 D gradient (A O = 40% for 2 D full gradient and A O = 63% for the 2 D shifting gradient). To have a more realistic value about the peak capacity, the 2D n c,practical was corrected by the real 2D space occupied by the separated compounds (A O ). The 2D n c,corr. values for the full gradient and shifting gradient were 170 and 1940, respectively. Although, as mentioned above, the peak capacity value should not be considered as a real number of separated peaks, this value is very useful not only to compare the improvement between different 2DLC methods carried out with the same setup during the optimization process but also to confirm the separation gain that 2DLC offers in comparison to conventional 1D separation for this kind of complex samples. In this work, the peak capacity achieved by 2DLC is much higher than the individual peak capacity obtained by the corresponding 1D alone. However, it is also important to remark that the 2DLC method is always accompanied by a big optimization and development effort and it involves specific instrumentation as well as a difficult data treatment as will be discussed in the following section. For this reason, the application of a 2DLC method should be always justified by a great gain in peak capacity.

Data treatment in µLC × LC
Although 2DLC has been established as a suitable alternative to conventional one-dimensional systems for the separation of very complex matrices in the last years, there is still a limitation on its use related to the data treatment. 2D raw data presents a matrix where each compound or feature is modulated several times due to the fractionation of each 1 D peak. Therefore, the signal of a single compound is divided with a time frame corresponding to the modulation time in the raw data. The modulation of the signal makes it not possible to apply automatic data treatment to create a feature list that allows the typical workflow for the identification and statistical analysis. Different efforts have been done to solve this limitation, for example, by the compression of the data applying a selection of the region of interest (ROI) in the MS dimension, a one-dimensional wavelet analysis of the previous m/z ROI values, and finally a time windowing strategy of the compressed data [38], or by the development of a new 2DLC concept called LC + LC, where the modulation time is increased to avoid the modulation of the peaks, that way, each peak is analyzed only once by the 1 D and the 2 D, achieving only one signal per compound [39]. However, the reported solutions present some limitations like the complexity of the data treatment workflow or the loss of 1 D resolution due to undersampling effects. In this work, a new data treatment has been used for the first time to convert 2D data feature analysis into a 1D data file that only shows one intensity or area per detected feature. This socalled demodulation process allows the use of common data analysis strategies (e.g., PCA) for the comparison between the chemical compositions of the cookie and gelato hemp samples acquired by a comprehensive 2D µLC × LC-HRMS analysis.
The analyses of the gelato and cookie samples were done in triplicate. Then, MS-Dial software was used for the identification and the alignment of all the features detected in the µLC × LC-HRMS analysis of the cookie and gelato hemp strains. The result of this alignment provided 75,387 features which included the divided signals of the modulated analytes. This feature list was introduced in the new demodulation tool. Although the signal of one compound is divided in modulation time frames, it follows a typical 1D Gaussian trend peak. The demodulation tool was programed to recognize the Gaussian trend of one m/z value that appears with time intervals equal to the modulation time and to combine the areas of all the modulated signals in the feature that corresponded to the maximum peak height. That way, the 2D modulated data are now reconstructed into 1D peaks that combine the total area of the 2D modulated peaks, that is, the demodulation tool transforms 2D data in 1D data. A graphical example of the demodulation process can be observed in Figure S3, where the ion m/z 865.2004 was modulated four times in the total retention time of 20.05, 20.55, 21.15, and 21.61 min and four times more at 29.99, 30.47, 30.95, and 31.44 min. After the demodulation, it can be seen how the areas of the corresponding modulations were summed up in the retention time of the maximum peak height. This tool provides a huge advantage in the global 2DLC methodology, since up to now, some omic applications have limited the use of this high resolution and separation technique due to the lack of powerful program able to deal with the complex 2DLC data treatment.
After applying the demodulation process, the feature list was reduced to 37,961 features. However, this number was still too high and therefore, filters were applied to obtain high quality features. Firstly, features with an intensity lower than 1 × 10 5 counts were discarded. Then, features that presented a relative standard deviation (RSD) higher than 50% across the triplicates were eliminated. This threshold was chosen as a compromise between analytical error and 2DLC repeatability (affected by both dimensions), since slight differences in the 1 D retention time can have high effect in both the 2 D retention time and area. Finally, all the features with a signal-to-noise ratio lower than 3 were also discarded of the data matrix. After this filtering, 5,296 features were considered for the statistical analysis.
Chemometric analysis was applied to evaluate the statistical differentiation of the cookie and gelato samples. First, a non-supervised method such as PCA was applied ( Figure S4a) and after that, a PLS-DA was performed as a supervised method to describe the model ( Figure S4b). In the loading plots of both statistical analyses, it is possible to observe a substantial number of features accumulated in the edge of the plot, which could be responsible for the statistical differentiation of the samples. Moreover, in the PCA and PLS-DA score plots, the PC1 was able to explain 65.2% of the samples' variance. These results indicate that the two hemp samples presented a different chemical composition that could be used for authentication purposes as well as for the application in specific therapeutical, pharmaceutical, or cosmetic motivations depending on the composition of those hemp strains.

Differentiation of the cannabinoids and phenolic compounds in the cookie and gelato strains
The composition of industrial hemp is highly affected by the genetic variation of the plant, and consequently, the composition of bioactive compounds is different in each of the strains [4,22]. Thus, the variety in the chemical composition and, particularly, in bioactive compounds can increase the applications of each hemp strain to very specific functions in the industry. The interest of this work was to evaluate the composition of cannabinoids and phenolic compounds of two commercially available industrial hemps and to assess possible differences between them. To do that, a flagging approach was carried out in the whole demodulated and filtered data.
For the suspected detection of cannabinoids, the database provided in the work of Berman et al. [7] was used to search for the most studied cannabinoids using the corresponding m/z values of [M-H] − ions as well as their typical MS/MS product ions. Besides, the most common product ions were selected to find other possible cannabinoids in the hemp strains.
Most of the cannabinoids found could be grouped into two families. On one hand, cannabinoids presented the typical fragment ion at m/z 179.1078 corresponding to the pentylresorcinol structure ( Figure S5a). The second family of cannabinoid compounds was found to have a product ion at m/z 195.1031 corresponding to a hydroxylated pentylresorcinol ( Figure S4b). For this last type of cannabinoids, only few information was available. In fact, only Berman et al. have reported them although they only classified and identified them as additional phytocannabinoids and provided the molecular formula. In particular, they found 27 compounds that were tentatively identified as potential phytocannabinoids by accurate mass and fragmentation patterns related to the already identified phytocannabinoids. Four of these 27 compounds were phytocannabinoid isomers at m/z 329.2122 and 373.2021 that presented the m/z 195.1027 as one of the main fragments [7]. Therefore, this fragment ion was also considered for the flagging approach to tentatively identify other potential phytocannabinoids. Besides these two families of cannabinoids, other cannabinoids described in the literature were considered.
After the flagging approach was done, the selected features were submitted to statistical analysis to determine the pattern of cannabinoids between the cookie and gelato sample. The heat map of the cannabinoids present in both samples is shown in Fig. 3a. As can be observed, it was possible to distinguish the two samples by their content on cannabinoids. Besides, this heat map revealed a cluster group of closely related compounds (marked in dash line) that were characteristic of the cookie sample. Moreover, it was also possible to appreciate that the cookie sample was richer in cannabinoids, being considered as potential markers since they presented a higher trend in this strain. The identification of these potential markers is shown in Table 1. Among them, there were some cannabinoids previously reported. For instance, the peaks observed at 46.08 and 47.07 min with an m/z value of 357.2091 and a fragmentation pattern consisted of m/z 313.2189, 245.1545, and 191.1080 could be assigned to the acidic forms CBDA, CBCA, CBLA, or CBRA. All these compounds as well as the acidic form of THC, the THCA, are isomers. However, THCA was not considered as possible identification due to the low concentration of THC in the industrial hemp samples. CBDA has been shown as the major compound in industrial hemp, and it is usually more abundant than the corresponding neutral form CBD. This fact is explained by the biosynthesis of these compounds since the phytocannabinoids are synthesized in the plant as acids [6,22,40]. The extracted ion chromatogram (EIC) of the ion m/z 357.2091 presented the maximum intensity at 46.08 min, so the compound eluted at that retention time could be tentatively identified as CBDA. Other  [7]. Moreover, additional phytocannabinoids described in the database without the complete characterization or name were also identified by means of the deprotonated molecule as important compounds for the characterization of the cookie sample. Following the nomenclature provided in the database, these compounds were named as the isomers 11a, 11b, 11c, or 11d (m/z 329.2136); 12a, 12b, 12c, or 12d (m/z 373.2020); 13c (m/z 327.1975); and 14a or 14b (m/z 371.1835).
The rest of the selected compounds for the cookie sample were not previously identified. As mentioned above, these potential unknown cannabinoids could be divided in two groups according to the main product ion shown in their fragmentation pathway. In total, eight additional phytocannabinoids presented as main in-source fragment ion m/z 179.1078 and fourteen phytocannabinoids showed the m/z 195.1031 as major fragment in the outstanding compounds in cookie hemp. The retention time, accurate mass, molecular formula, and the MS/HRMS fragments of all of them are summarized in Table 2. A possible molecular structure of one of each phytocannabinoid group is proposed in Figure S5c-d.
On the other hand, a similar strategy was followed for the suspected analysis of phenolic compounds. In this case, Phenol-Explorer was used as database for the search of phenolic compounds in both samples [41]. This search was done by monitoring the aglycone ions which usually are the main fragments of phenolic compounds. The heat map of the phenolic compounds found in the samples is depicted in Fig. 3b. As it happened with the cannabinoid content, it was also possible to differentiate the samples according to their phenolic compound content. Interestingly, in this case, gelato sample showed a group of closely related compounds highlighted in comparison to the gelato strain (dashed line). Thereby, these compounds were considered as potential authentication markers for the gelato sample.
The potential phenolic markers of the gelato sample are listed in Table 2. These compounds belonged to different classes of phenolic compounds. One of the main phenolic compound classes was the procyanidins, which are polymeric phenolic compounds formed by the link of catechin or epicatechin units and constitute the second most abundant group of phenolic compounds in nature. They are responsible not only for some organoleptic properties like astringency but also several therapeutical properties have been related to them, like antioxidant, anticancer, cardioprotective, antimicrobial, antiviral, neuro-protective, and antiinflammatory activities, among others [42,43]. The presence of this phenolic class in C. sativa has been recently reported for the first time [23]. In that work, two procyanidin dimers and two trimers were identified. In this work, one procyanidin dimer and five different trimers were shown as characteristic compounds from the gelato strain. For example, the procyanidin trimers were tentatively identified thanks to the accurate mass of the precursor [M-H] − ion at m/z 865.2004 and to the typical fragmentation pattern of procyanidins consisting of the neutral losses of one or two units of (epi) cathechin leading to the precursor ions at m/z 577.1281 or 575.1200 and 289.0728 or 287.0557, respectively, depending on where these losses were produced (terminal or the intermediate units). Besides, typical fragment ions from the Retro-Dial-Alder reaction were observed like m/z 739.1757, 713.1755, 451.1026, 425.0867, and 125.0254 ( Figure S6). Another representative family of phenolic compounds described as markers of the gelato strain was the flavone  The peak ID refers to the peak labeling in Figure S7 [4,23]. Therefore, these methylflavone derivates could be tentatively identified as diosmetin-or chrysoeriol-related compounds. However, more investigation should be done to confirm the identity of these compounds although it can be concluded that all these compounds are derived from the same aglycone (m/z 299.0556) and they have a defined glycosidic pattern. The rest of the detected compounds belonged to the flavonol, flavone, isoflavone, and phenolic acid classes. Among the flavonols, quercetin-dihexoside and quercetin-rutinoside were identified by their [M-H] − ions (m/z 625.1387 and 609.1461, respectively) and the common fragment ion m/z 301.0352 ([C 15 H 9 O 7 ] − , − 0.6 ppm) corresponding to the quercetin aglycone. Three flavones were also highlighted in the gelato strain. In this case, the found aglycone fragment ion presented an m/z value of 285.0416 ([C 15 H 9 O 6 ] − , 4.0 ppm), which corresponds to the isomers luteolin or kaempferol. The three derivatives of these aglycones were luteolin or kaempferol dihexoside (m/z 623.1258), luteolin or kaempferol hexuronide (m/z 461.0750), and luteolin or kaempferol hexuronide-penstoside (m/z 593.1179). Quercetin, luteolin, and kaempferol derivates have been extensively described in C. sativa varieties [4,28]. The isoflavone detected in the group of highlighted compounds in gelato was tentatively identified as genistein hexuronide with a Finally, phenolic acids also form part of the phenolic profile of cannabis. Caffeic acid, cinnamic acid, benzoic acid, coumaric acid, and ferulic acid are some of the phenolic acids that have been identified in different cannabis plants [4,28,44]. Here, several phenolic acids were identified in this group like caffeic acid (m/z 179.0362), salicylic acid-hexoside (m/z 299.0771), and dihydroxybenzoic acid-hexoside or protocatechuic acid-hexoside (m/z 315.0735). Therefore, the chemical composition of the cookie strain presented a richer cannabinoid profile in comparison with the gelato strain, showing several cannabinoids that distinguished that sample. In contrast, the gelato strain was more related to the content on phenolic compounds, being some of them characteristic analytes of this sample.
The peak labeling of all the tentatively identified compounds in cookie and gelato samples can be observed in Figure S7.

Conclusions
In this work, a µLC × LC-HRMS method is used for the analysis of two commercial industrial hemp strains (i.e., cookie and gelato). Two reversed phase modes were coupled in the 1 D and 2 D. To increase the orthogonality of these correlated modes, a 2 D multi-segment shifting gradient was optimized that enhanced the orthogonality in 23% and produced a peak capacity eleven times higher. After the analysis, a demodulation process, which is able to transform 2D data into 1D data, was applied for the first time. With this process, the 2D data treatment of a very complex sample (75,387 features) was easily handled for the statistical and identification process. Due to the high variability in bioactive compounds between hemp varieties and strains, a suspected analysis of the phytocannabinoids and phenolic compounds present in both samples was carried out. Cookie sample presented a higher content and a higher number of characteristic cannabinoids that could be considered as potential markers of this strain. Among them, the acidic form of CBD (CBDA) or its isomer, varinic acid-C3 type, and neutral-C4 form type phytocannabinoids were identified together with several cannabinoids that presented a common MS/MS pathway not previously reported. On the other hand, the gelato sample was richer in phenolic compounds, among which an important number of potential markers were highlighted. Procyanidins and diosmetin or chrysoeriol glycosidic derivates were the major compounds that were characteristic of this sample. Therefore, this study reveals the different bioactive compound profiles between two industrial hemp varieties. The chemical characterization of them would be of great interest for pharmaceutical, food, or cosmetic applications that could be targeted to specific interests according to the properties offered by the compounds present in each strain. Another conclusion of this work is that the number of secondary metabolites in cannabis is exponentially growing more and more, and further studies are needed for achieving a high confident identification level of all the bioactive compounds responsible for all the interesting applications of cannabis.