**Preface**

This Special Issue, entitled "Feeding Strategies and Nutritional Quality of Animal Products," delves into the complex interplay between the diets of livestock, health and the nutritional value of the products they yield. In our ever-evolving world, the dynamics of food production and consumption have become increasingly vital. As the global population continues to grow, so does the demand for high-quality, nutritious food, especially animal products.

The quality of animal-derived food products is linked to the animals' diets. What animals eat significantly influences the composition of their products, impacting not only their health, but also the nutritional quality of their products. Researchers worldwide have been exploring innovative feeding strategies, ranging from natural foraging to precision nutrition, to enhance the nutritional quality of these products, while ensuring sustainable and ethical production practices.

This Special Issue collates a diverse array of studies, encompassing both fundamental research and practical applications. Authors' contributors investigate the impact of various feeding strategies on the levels of essential nutrients, bioactive compounds, and fatty acid profiles in animal products. They also explore the potential health benefits and drawbacks associated with the various dietary strategies employed for livestock.

The articles collected in this Special Issue demonstrate the relationship between feeding strategies and the nutritional quality of animal products, and aim to contribute to the application of sustainable feeding strategies and nutritious food production, ultimately benefiting both human health and the environment.

> **Arabela Elena Untea, Mihaela Saracila, and Petru Alexandru Vlaicu** *Editors*

### *Editorial* **Feeding Strategies and Nutritional Quality of Animal Products**

**Arabela Elena Untea \*, Mihaela Saracila \* and Petru Alexandru Vlaicu \***

Feed and Food Quality Department, National Research and Development Institute for Animal Biology and Nutrition, 077015 Balotesti, Romania

**\*** Correspondence<EMAIL_ADDRESS>(A.E.U<EMAIL_ADDRESS>(M.S.);

<EMAIL_ADDRESS>(P.A.V.)

Feeding strategies play a crucial role in determining the nutritional quality of animal products. The type and composition of an animal's diet directly influence the nutritional content of the products derived from them, such as meat, milk, and eggs. Various feeding strategies, including traditional and modern approaches, impact factors like the balance of macronutrients, essential vitamins, and minerals, as well as overall product quality and safety. Alternative feed sources containing bioactive compounds are being employed as agents that improve the health of animals and enhance the quality of animal products [1]. Incorporating high-quality and well-balanced feeds, optimizing grazing practices, considering supplementary nutrients, and even implementing sustainable and ethical feeding practices can lead to animal products with enhanced nutritional profiles and improved consumer satisfaction. It is evident that the careful selection and management of feeding strategies are essential to ensure the production of animal products that meet both human dietary needs and sustainability goals. Ongoing research and developments in animal nutrition continue to refine feeding strategies. Scientists are working to identify novel feed ingredients, optimize nutrient delivery, and develop precision feeding approaches that cater to the specific needs of animals at different stages of growth.

Among the papers included in this Special Issue, ten discuss different feeding strategies for ruminants and non-ruminant farm animals (nine research studies and one review paper). The research papers regarding the nutrition of ruminants focus on enhancing the diets of sheep and lambs, bioactive feed additives for dairy cattle, beef cattle, and newly weaned Nellore cattle, and nutritional strategies for the prevention of clostridial disease in cattle. The studies on non-ruminant feeding strategies not only center around the effects of microencapsulated probiotics on weaned piglets but also the effect of age at slaughter on the meat quality of male layer-type chickens. The exposure of piglets to lower concentrations of Zearalenone mycotoxin and a natural nitrite, as well as antioxidant sources for pork mince, are also considered. The scientific subject of the lone review paper mentioned above relates to obtaining valuable pigment additives for poultry nutrition.

The papers featured in this Special Issue explore three main topics, namely, animal health and well-being, growth performance, and product quality.

Animal health is a comprehensive concept that encompasses various factors that contribute to the overall well-being of animals, with a critical focus on gut health and the immune system. These aspects are interrelated and play a pivotal role in determining animal's resilience to diseases, growth performance, and the quality of animal products. Implementing appropriate feeding strategies is integral to maintaining optimal animal health and enhanced product quality.

For this topic, Barbosa et al. [2] explored the feasibility of substituting 0.018 g of lasalocid sodium with 1.500 g of dried and milled barbatimão bark or 0.300 g of barbatimão bark extracts (*Stryphnodendron* sp.) in the diets of feedlot lambs. The findings of this study highlighted the advantageous impact of the investigated feed additives on health status by reducing total cholesterol levels when compared with synthetic feed additives. These results require further investigation into the bioavailability of bioactive compounds derived

1

**Citation:** Untea, A.E.; Saracila, M.; Vlaicu, P.A. Feeding Strategies and Nutritional Quality of Animal Products. *Agriculture* **2023**, *13*, 1788. https://doi.org/10.3390/ agriculture13091788

Received: 25 August 2023 Revised: 1 September 2023 Accepted: 8 September 2023 Published: 9 September 2023

**Copyright:** © 2023 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/).

from such extracts. In the context of dairy cows, the persistent challenge of clostridial disease remains troubling despite the global utilization of various vaccine methodologies. Probiotics and holistic supplements are gaining prominence as preventive measures against clostridiosis. An innovative approach that has proved effective in bolstering animal resilience, particularly against *Clostridium perfringens* β-toxin, involves the amalgamation of microorganisms and bioactive substances. Nekrasov et al.'s study [3] discusses improved metabolism and heightened antioxidant status in cows. This comprehensive strategy offers valuable insights into the evolution of bovine clostridia prophylaxis. The use of biofunctional additives has emerged as a practical facet within the scope of animal biosecurity measures against toxins, paving the way for their integration into animal diets. This study marks a significant step in unraveling clostridia prevention strategies, holding potential applications in enhancing biosecurity through animal diets. Batista et al. [4] aimed to assess diverse strategies for supplementing immunomodulatory feed additives at (10 g/100 kg BW/day of NUTRA) before and after transportation in weaned Nellore calves raised on pasture and the subsequent impact on physiological parameters and performance as they progressed through the growth phase. Their findings suggest that the specific conditions of the experiment, such as the short duration of transportation stress and the relatively stress-free grazing environment, could have mitigated the detrimental effects of stress on the animals. Consequently, the anticipated advantages of employing the immunomodulatory feed additive might have been addressed. The impact of microencapsulated *Lactobacillus acidophilus* and *Lactobacillus plantarum* supplementation, individually and in combination, was studied by Lefter et al. [5] on weaning piglets. Their findings revealed a remarkable reduction in diarrhea incidence that was particularly evident in the *Lactobacillus plantarum* and combined *Lactobacillus acidophilus* + *Lactobacillus plantarum* groups. Furthermore, the use of microencapsulated probiotics exhibited distinct advantages in terms of intestinal health. These results collectively showcase the potential to mitigate instances of diarrhea while concurrently ameliorating both intestinal structure and microflora composition. Bulgaru et al. [6] concentrated on the effects of two ZENcontaminated diets—one below (75 ppb) and one above (290 ppb) the EU's recommended levels—on immune and oxidative responses in weaned piglets. Strikingly, exposure to ZEN at both concentrations exerted minimal effects on immune markers, oxidative stress, and colon inflammation. However, for a more comprehensive understanding of ZEN's influence on immune defense, broader studies on various organs are required. These findings contribute to the assessment of suitable ZEN concentrations in piglet feed and offer insight into the underlying mechanisms driving its effects.

However, achieving this intricate equilibrium is not a solitary endeavor. Instead, it hinges on the strategic application of prudent feeding approaches. The importance of these approaches cannot be exaggerated, as they are pivotal in nurturing animal welfare and enhancing the quality of resulting products. By tailoring diets to encompass the specific nutritional needs of animals, these strategies underpin the establishment and preservation of optimal health. Consequently, the animals are not only better equipped to endure potential health challenges but also thrive, achieving their growth potential and manifesting their genetic traits to the fullest, showing the profound correlation between the growth performance and nutritional quality of animal diets. This topic is also covered in the papers published in this Special Issue.

Arcos-Álvarez et al. [7] investigated the effect of olive oil treatment on various aspects related to lactating hair sheep, including their productive traits, and their treatment strategy displayed a positive outcome, as it did not exert any detrimental influence on the productive behavior, milk production, or chemical composition of lactating sheep. In feedlot lambs, Barbosa et al. [2] showed that substituting lasalocid sodium with dried and milled barbatimão bark or barbatimão bark extracts (*Stryphnodendron* sp.) demonstrated no detrimental effects on production performances. Similarly, Batista et al. [4], after supplementing immunomodulatory feed additives in weaned Nellore calves, found that their treatment

had no effect on animal performance. Likewise, in Lefter et al. study [5] on weaning piglets, probiotics supplements showed potential in enhancing production performance.

Feeding strategies have a direct impact on quality attributes in animal products. From influencing nutritional quality and tastes to determining fatty acid profiles, the nutritional composition of diets significantly shapes the sensory and nutritional characteristics of the final meat products, thus influencing consumer preferences and choices. The effect of animals' diets on product quality is a topic of continuous interest for researchers, farmers and producers, as shown in some of the papers published in this Special Issue.

According to Arcos-Álvarez et al. [7], including extra virgin olive oil in the diets of lactating sheep appears to offer a viable approach to mitigate saturated fatty acids while simultaneously elevating the levels of monounsaturated and polyunsaturated fatty acids. The utilization of extra virgin olive oil contributed to significant enhancements in specific fatty acid contents, including linoleic, linolenic, and eicosapentaenoic acids with potential positive impact on heart health. Soares et al. [8] presumed that the incorporation of 1.6 g/kg supplement dry matter of a blend of live yeast (*Saccharomyces cerevisiae* strains) and organic trace minerals, particularly chromium, could effectively replace the conventional monensin (30 mg/kg) and inorganic trace minerals during the finishing phase of beef cattle, and their study reports favorable effects on feed efficiency and carcass traits. In fact, the inclusion of a blend of live yeast and organic trace minerals exhibited notable benefits in terms of animal growth and meat color, all of which were attributed to elevated concentrate intake and heightened diet digestibility facilitated by the utilization of live yeast and organic minerals. The insights derived from the research of Popova et al. [9] contributed to understanding how the meat quality of slow-growing chicken male layertype chickens differs according to meat portions, particularly regarding chicken breasts and thighs. In terms of nutritional constituents, advancing age was linked to a noteworthy reduction in intramuscular fat content in thighs, with a similar trend observed in breast meat. Collectively, the meat from nine-week-old male layer-type chickens displayed certain drawbacks regarding fatty acid profiles. This aspect unveils potential avenues for future studies to explore diverse feeding strategies or housing systems to enhance this particular trait. The impact of introducing varying levels of ethanolic hawthorn berry extract in conjunction with a consistent concentration of fermented parsnip juice on factors like lipid stability, heme pigment conversion, residual nitrite content, and the growth of spoilage bacteria in refrigerated pork mince was investigated by Predescu et al. [10]. Coupling fermented parsnip juice with hawthorn extract exhibits a notable inhibitory impact on spoilage bacteria that break down nitrogen compounds, leading to enhanced stability in unsaturated fatty acids that are prone to oxidation in minced pork. In this study, the interplay between the hawthorn extract and the fermented parsnip juice appeared to be pivotal in governing color and lipid stability.

In response to the growing emphasis on natural products within the food and feed markets, the demand for organic ingredients is on the rise. Current trends in livestock nutrition research are leaning toward feed formulation diets that incorporate natural and organic feed additives [11]. These additives encompass active elements that impact plants, animals, and microorganisms, yielding valuable biocompounds. Notably, yeast pigments offer a wealth of natural colors, coupled with a broad spectrum of nutritional and medicinal attributes, as noted by Grigore et al. in their review [12]. Achieving these optimizations hinges on strategic approaches involving strain genetic engineering and process development, often utilizing cost-effective organic substrates. According to Grigore et al. [12], investigations into the effects of value-added yeast pigment additives on livestock well-being, productivity, and product quality are essential in substantiating their nutritional and medicinal potential. It is equally important to consider the perspectives and preferences of consumers regarding purchasing animal products derived from microbial pigment additives, as well as the need for an enhanced understanding of these additives and their implications.

In summary, the presented works contribute towards solving issues regarding the development of "Feeding Strategies and Nutritional Quality of Animal Products" in different livestock systems. The reported results addressing sheep, lambs, dairy cattle, weaned cattle, and piglets, as well as product quality, are of interest to specialists and scientists involved not only in research but also in daily farm support and management. The feeding strategies presented in the works published in this Special Issue address some important aspects in terms of animal health, production performance, and product quality.

**Author Contributions:** Conceptualization, A.E.U., M.S. and P.A.V.; writing—original draft preparation, A.E.U., M.S. and P.A.V.; writing—review and editing, A.E.U., M.S. and P.A.V.; visualization, A.E.U., M.S. and P.A.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** The GE would like to thank to the Romanian Ministry of Research, Innovation and Digitalization, project PED 631/2022 and the support provided by the program National Research Development Project to Finance Excellence (PFE)—8/2021, which facilitated the studies published in this Special Issue.

**Acknowledgments:** The Guest Editors would also like to thank all of the authors who contributed to this Special Issue, as well as the reviewers who reviewed these papers and the Managing Editors.

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

#### **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.


**2. Method of Study**

First, the morphological characteristics of Lingayen Gulf were briefly described using a satellite image and bathymetric survey data as well as the investigation of wave characteristics using the Global Wave Statistics [8]. Because of the shape of Lingayen Gulf with a large aspect ratio, oblique wave incidence with a large angle to the direction normal to the shoreline prevails in this gulf. Taking these characteristics into account, satellite images of Lingayen Gulf were collected and the shoreline changes along the east coast were determined, selecting the deltas of two rivers (Balili and Aringay Rivers) and a sand spit at

Santo Tomas with a supplementary field observation on 23 February 2018. Then, the wave field in Lingayen Gulf was calculated to explain the mechanism of occurrence of shoreline undulations along the east coast of the gulf due to high-angle wave instability [1] using the energy balance equation [9], based on the seabed contours of the gulf. Along the east coast of this bay, waves are obliquely incident to the direction normal to the shoreline over 45◦, causing the shoreline instability. In an area south of the Balili River delta, a sand spit has formed at a location where the shoreline configuration abruptly changes due to high-angle wave instability. Therefore, the results of the numerical simulation regarding the formation of a sand spit on a coast with a sudden change in shoreline configuration were shown and the comparison of the results with the measured data was made. Furthermore, the rate of longshore sand transport along the sand spit at Santo Tomas was estimated from the temporal change in sand volume by integrating the foreshore area of the beach. Finally, the morphological similarity of Lingayen Gulf to that of Suruga Bay in Japan, and the rate of longshore sand transport in both bays were compared to consider the future measures against beach erosion along the coasts in these bays.


*3.3.1. Clustering in MAC layer*

without the hidden/exposed terminal or deafness problem to ensure reliable message delivery. Although the infrastructure is an extra, it will be furnished on the highways extensively and applied in VANET in the near future. Therefore, compared with great and lifelong benefit, the infrastructures expense is of trifling importance at all. The efficient cluster based MAC and routing protocols can provide a more stable communication than a solution using V2V clustering. The optimum protocol should that take the advantages of fixed infrastructure and

V2V based clustering is a decentralized clustering where clusters are formed based on communication between vehicles. Additionally, the CH election will be based on V2V communication. There are several advantages of using V2V-centric clustering as compared with the infrastructure-centric VANETs. V2V-centric clustering can avoid the short commu‐ nication link period, high frequent hand-offs, fast channel fading, etc., that are caused by the high relative-speed difference between the fixed infrastructure and the fast-moving vehicles. Finally, the V2V-centric clustering performs better in active safety applications, which only requires exchanging messages among one hop vehicles within their transmission range.

V2V communications are expected to significantly improve transportation safety and mobility on the road. Several applications of V2V communications have been identified, from safety and warning applications, up to traffic control and driver assistance applications. In infra‐ structure centric clustering, all the communications is done via the infrastructure which causes a lot of control overhead and additional delay. Furthermore, it would be very cost intensive to build an infrastructure based communication all along the road structure. The V2V based clustering technique avoids the use of stationary base stations by building up VANETs, where all vehicles in a common transmission radius can exchange messages. However, CH selection carried out through V2V communications has some shortcomings, e.g., highly complex protocols, large computation and communication cost, need of additional devices and so on. Another important issue is that the connection between two adjacent CHs may be lost due to vehicles high speed, which drastically reduces the link quality. The hidden terminal problem where two vehicles are outside of each other's transmission radius, but both attempt to transmit to a vehicle that is within the radius of both. This issue is likely in pure V2V scenarios where there is no centralized communication system. The result of the hidden terminal problem is data collisions. By enabling vehicles to transmit/receive messages with each other via V2V as well as with infrastructure communications, VANETs could contribute to more safer and congestion free roads by providing correct and timely message to neighboring

Clustering can simplify essential functions like bandwidth utilization, routing, and channel access. In MAC layer, it can provide a fairer and reliable channel access to all vehicles in network. This can lead to increase in the reliability of packets and scalability of the network. In network layer, clustering for routing can find the closest vehicles to intended destination.

optimize the problem.

**3.2. V2V based clustering**

228 Contemporary Issues in Wireless Communications

vehicles and other related departments.

**3.3. Clustering in layers**

Introducing a cluster scheme already on the MAC layer additionally provides the possibility of a fairer medium access. When clustering applied in VANETs, it brings interesting research point such as broadcast storm that occurs when several vehicles are passing at a specific region at the same time, causing network congestion, packet collisions and delays in the medium access layer. A cluster-based MAC scheme is needed in V2V communication to overcome the lack of specialized hardware for infrastructure and the mobility to support network stability and channel utilization. In this case the CH can take over the responsibility to assign bandwidth to the CMs and therefore even QoS support can be improved. As the bandwidth can be assigned centrally fewer collisions have to be expected which consequently increases the reliability.

Many researchers have proposed cluster based multi-channel medium access control protocols to improve the performance and reliability of VANETs. In these protocols, clustering is used to limit channel contention and provide fair channel access within the cluster. On the other hand, multi-channel is used to increase the network capacity by the spatial reuse of the network resources and reduce the effect of the hidden terminal problem. Moreover, to optimize the communication range and the cluster size is very difficult especially in a highly dynamic environment such as VANETs. However, in order to overcome this situation some approaches divide the service area into a set of region units, and limit the number of vehicles in each region unit for the contentions of radio channels. Each region unit is then associated with a nonoverlapping radio channel pool. Since the number of vehicles in each region unit is limited, the contention period is reduced and the throughput is increased. However, these types of approach have low bandwidth utilization in case of sparse traffic. Some of clustering algo‐ rithms try to minimize the total number of clusters by creating hierarchical clusters with a diameter of at most four hops. In this section of the chapter, we compare well known cluster based MAC protocols in Table 3.

The MAC layer is divided into different cycles and each cycle is divided into contention based or contention free. In the current literature, several MAC protocols have been proposed to VANETs. Some of the well-known MAC protocols are ADHOC MAC [32], SDMA [33], VeMAC [34], DMMAC [35], STDMA [36], VeSOMAC [37] etc. These protocols are proposed for various scenarios and have many drawbacks such as hidden terminal problem, time unbounded, unreliability etc. There is a need for new MAC protocol in VANETs that can provide mobility (i.e., the MAC protocol should support vehicles to leave and join inter-vehicle communications at high speed), delay bounded (i.e., the communication must be delay bounded and real-time), scalability (i.e., VANET should scale itself according to the number of vehicles present), bandwidth efficiency (i.e., the radio resource should be utilized in an efficient and fair manner), cost (i.e., for cost-efficient and reliable communications, VANET should be fully decentral‐ ized), and fairness (i.e., every vehicle should get a fair chance to get the radio channel).The

**Table 3.** Comparison between various cluster based MAC protocols.

challenge of successfully deploying VANET services is to ensure timely and reliable data delivery for mobile vehicles.


*2.2. BEV Powertrain*

In the BEV powertrain layout illustrated in Figure 1b, the fuel tank, ICE and AMT of the CV powertrain are replaced with a high-voltage battery, an electric motor (EM) and a direct drive, respectively. Finally, power electronics enable the proper operation of the electric powertrain components. In this framework, speed *ωEM* and torque *TEM* of the EM can be evaluated at each time step following the same procedure illustrated for the CV case in Equations (1)–(3). A gear shift logic does not need implementation in this case, given the embedment of a direct drive. The EM electrical losses *lossEM*, including inverter losses as well, can be evaluated in this case by interpolating in a two-dimensional lookup table with *ωEM* and *TEM* as independent variables. Finally, the battery state-of-charge (SOC) over time can be evaluated as a function of the requested battery power *Pbatt* by adopting an equivalent circuit approach and following Equations (4)–(6) [34].

$$SOC(t) = \int\_{t\_0}^{t\_{end}} S \dot{O} \mathcal{C} \left[ P\_{batt}(t), SOC(t) \right] dt \tag{4}$$

with:

$$P\_{\rm hatt}(t) = \omega\_{EM}(t) \cdot T\_{EM}(t) + \log \text{s}\_{EM}[\omega\_{EM}(t), T\_{EM}(t)] + \log \text{s}\_{\rm aux}(t) \tag{5}$$

$$\begin{array}{l} \text{S\dot{O}C} \text{[}P\_{\text{batt}}(t)\text{,SOC}(t)\text{]}\\ = \frac{V\_{\text{OC}}[\text{SOC}(t)] - \sqrt{\{V\_{\text{OC}}[\text{SOC}(t)]\}^2 - 4 \cdot R\_{IN}[\text{SOC}(t)] \cdot P\_{\text{batt}}(t)}}{2 \cdot R\_{IN}[\text{SOC}(t)]} \cdot \frac{n\_P}{A h\_{\text{batt}} \cdot 3600} \end{array} \tag{6}$$

where . *SOC*, *t*0 and *tend* are the instantaneous rate of SOC, the initial time instant and the final time instant of the drive cycle, respectively. *lossaux* is the power requested by the accessories (e.g., air conditioning, lubrication), and it is modeled as having a constant value in this work. *RIN* and *VOC* represent the internal resistance and the open-circuit voltage of the battery pack, as obtained by interpolating in 1D lookup tables with SOC ad independent variables. *nP* is the number of cells in parallel according to the battery pack layout, while *Ahbatt* represents the battery pack energy capacity in ampere-hours. The factor of 3600 is considered here to convert energy units in ampere-seconds.


*Article* **Multi-Media Exposure to Polycyclic Aromatic Hydrocarbons at Lake Chaohu, the Fifth Largest Fresh Water Lake in China: Residual Levels, Sources and Carcinogenic Risk**

**Ning Qin 1,2, Wei He 2,3, Qishuang He 2,4, Xiangzhen Kong 2,5, Wenxiu Liu 2,6, Qin Wang 7 and Fuliu Xu 2,\***

**Abstract:** The residual levels of 16 priority polycyclic aromatic hydrocarbons (PAHs) in environment media and freshwater fish were collected and measured from Lake Chaohu by using Gas chromatography-mass spectrometry. Potential atmospheric sources were identified by molecular diagnostic ratios and the positive matrix factorization (PMF) method. PAH exposure doses through inhalation, intake of water and freshwater fish ingestion were estimated by the assessment model recommended by US EPA. The carcinogenic risks of PAH exposure were evaluated by probabilistic risk assessment and Monte Carlo simulation. The following results were obtained: (1) The PAH16 levels in gaseous, particulate phase, water and fish muscles were 59.4 ng·m−3, 14.2 ng·m−3, 170 ng·L−1 and 114 ng·g−1, respectively. No significant urban-rural difference was found between two sampling sites except gaseous *BaPeq*. The relationship between gaseous PAHs and PAH in water was detected by the application of Spearman correlation analysis. (2) Three potential sources were identified by the PMF model. The sources from biomass combustions, coal combustion and vehicle emission accounted for 43.6%, 30.6% and 25.8% of the total PAHs, respectively. (3) Fish intake has the highest lifetime average daily dose (LADD) of 3.01 × 10−6 mg·kg−1·d−1, followed by the particle inhalation with LADD of 2.94 × 10−6 mg·kg−1·d−1. (4) As a result of probabilistic cancer risk assessment, the median ILCRs were 3.1 × 10−5 to 3.3 × 10−5 in urban and rural residents, which were lower than the suggested serious level but higher than the acceptable level. In summary, the result suggests that potential carcinogenic risk exists among residents around Lake Chaohu. Fish ingestion and inhalation are two major PAH exposure pathways.

**Keywords:** PAHs; multi-media exposure; health risk; probabilistic risk assessment; Lake Chaohu

#### **1. Introduction**

Polycyclic aromatic hydrocarbons (PAHs) are globally concerned pollutants because of their widespread occurrence, strong persistence and long-range transportation potential [1]. Furthermore, they possess potential toxicity, mutagenicity and carcinogenicity [2–4]. Studies have shown that human cancer causes of skin, lungs and bladder have always been associated with PAHs [5–8], and 16 PAHs are included on the priority pollutants list of the

**Citation:** Qin, N.; He, W.; He, Q.; Kong, X.; Liu, W.; Wang, Q.; Xu, F. Multi-Media Exposure to Polycyclic Aromatic Hydrocarbons at Lake Chaohu, the Fifth Largest Fresh Water Lake in China: Residual Levels, Sources and Carcinogenic Risk. *Atmosphere* **2021**, *12*, 1241. https:// doi.org/10.3390/atmos12101241

Academic Editor: Alina Barbulescu

Received: 20 July 2021 Accepted: 19 September 2021 Published: 23 September 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/).

US EPA. PAHs have a wide variety of sources, including coal combustion, vehicle emission, coking industry and biomass burning [9–11]. After being emitted into the environment, PAHs may redistribute in environmental media and result in people being exposed to these pollutants through multiple pathways, including breathing in polluted air and particles, drinking water, dietary intake and dermal contact with contaminated soil [12,13]. Multi-media distribution and multi-pathway exposure render the assessment of PAH exposure complicated. Therefore, accurately evaluating the contribution of each exposure pathway, characterizing the carcinogenic risk and identifying the sensitive parameters in the exposure process are crucial to the management of PAH emission.

Water bodies act either as a sink [14] or as a source [15] for PAHs in the environment. The atmospheric PAHs can enter water system through wet deposition, dry deposition and gas exchange across the air–water interface [16–18]. Meanwhile, PAHs in water may accumulate in aquatic organisms by direct uptake from water through gills or skin or by the ingestion of suspended particles and contaminated food [19]. Residents living near the lake can be exposed to PAHs by inhaling polluted air and ingesting water and aquatic products.

Lake Chaohu is located in the Anhui Province, which belongs to one of the most developed areas in China, Yangtze River Delta Economic Zone. During the last decades, the PAHs' emission in China, especially in the above-mentioned areas, increased greatly due to the increasing energy demand associated with rapid population growth and economic development and to the low efficiency of energy utilization [20]. Chaohu is famous for its fresh water fish. It is also the drinking water source of large cities such as Hefei and Maanshan. PAH pollution in water system of Chaohu may increase the risk of residents' exposure through fish ingestion and drinking water. Therefore, in recent years, PAH exposure in Lake Chaohu has become a topic of concern. Some studies have been conducted on the PAH residual levels in environment media, source apportionment or environment behaviors [21–25]. Despite the progress in these directions, studies on three issues remain scarce. First, comparison among exposure contribution from different pathways are seldom reported. Second, carcinogenic risk due to total exposure remains unclear. Third, the factors influencing risk assessment are seldom studied. Hence, further studies on PAHs should be performed in order to obtain a comprehensive understanding of the risk profile of PAHs exposure among Chaohu residents.

In this research, EPA priority control PAHs were selected as target chemicals due to their extensive residence in the environment and their threat to public health. The contents of 16 l PAHs in air, particles, lake water and aquatic organisms were measured; physiological and behavior parameters influencing PAH exposure were collected; potential sources were identified; and cancer risks were calculated by US EPA model. The aims of this study were to elucidate the characteristics of multi-media PAH exposure of residents and to provide information for PAH management near Lake Chaohu.

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

#### *2.1. Sample Collection and Pretreatment*

Two sampling sites were selected near Chaohu City and Zhongmiao Town as urban and rural sites (Figure 1). Atmosphere samples were collected once a month from May 2010 to April 2011 by high-volume samplers. Polyurethane foam (PUF) disk and glass fiber filter (GFF) were used to collect gaseous phase and particulate phase PAHs, respectively. Water was collected from two sampling sites selected near the atmospheric sites. After shaking and mixing, a one liter aliquot of each collected water sample was filtered through a 0.45 μm glass fiber filter using a filtration device consisting of a peristaltic pump (80EL005, Millipore Co., Billerica, MA, USA).

**Figure 1.** The location of Lake Chaohu and sampling sites: (**a**) People's Republic of China, (**b**) Anhui Province and (**c**) Lake Chaohu.

In the laboratory, PUF and GFF samples were added with surrogate standards of 2-fluoro-1, 10-biphenyl and p-terphenyl-d14 (J&K Scientific, Beijing, China, 2.0 mg mL 1) before measurement. The PUF was Soxhlet extracted with 150 mL 1:1 mixture of n-hexane and acetone for 8 h. GFF was extracted by 25 mL hexane/acetone mixture (1:1) using a microwave accelerated reaction system (CEM Corporation, Matthews, NC, USA). Microwave power was set at 1200 W, and the temperature program was set to the following: ramp up to 100 ◦C in 10 min and held at 100 ◦C for another 10 min. Both PUF and GFF extracts were concentrated to 1 mL by rotary evaporation at a temperature below 38 ◦C and then transferred to a silica/alumina chromatography for cleanup. The elution solution was collected, concentrated, conversed to hexane solution and then added with internal standards (Nap-d8, Ace-d10, Ant-d10, Chr-d12 and Perylene-d12, J&K Scientific Ltd., Beijing, China).

The water samples were extracted by using a solid phase extraction (SPE) system (Supelco, Bellefonte, PA, USA). C18 cartridges (500 mg, 6 mL, Supelco, Bellefonte, PA, USA) were prewashed with dichloromethane (DCM) and conditioned with methanol and deionized water. A 1 L water sample was added with surrogate standards, passed through the SPE system and was extracted. The cartridges were eluted with 10 mL of dichloromethane. The volume of the extracts was reduced by a vacuum rotary evaporator in a water bath and was adjusted to a volume of 1 mL with hexane. Internal standards were added for analysis.

The fish samples were pretreated on the same day after being delivered back to the temporary laboratory. The muscles on both sides of the dorsal and chest were mixed. After obtaining the wet weight, the samples were freeze dried (FDU-830, Tokyo Rikakikai Co., Tokyo, Japan) and grounded into a granular powder with a ball mill (MM400, Retsch GmbH, Haan, Germany). Two gram powder samples were weighed into an extraction tube, and the surrogate standards were added to the samples to indicate recovery. After microwave extraction, the extracts were pressure filtered and concentrated to approximately 1 mL and cleansed by GPC instrument (GPC800+, Lab Tech Ltd., Hongkong, China) with a Bio Beads SX-3 column (300 mm × 20 mm, Bio-Rad Laboratories, Inc., Hercules, CA, USA). Subsequently the concentrate was loaded in a silica gel SPE cartridge (6 mL, 500 mg, Supelco Co., Bellefonte, PA, USA). The cartridge was eluted by hexane and mixed solution of dichloromethane and hexane. The extracts were concentrated to 1 mL, transferred to vials, added with internal standards and sealed for analysis. The details of experiment have been reported in previous research [19,26].

#### *2.2. Instrument Analysis*

The samples were analyzed by using Agilent 6890 gas chromatography and a 5976C mass spectrometer detector with a HP-5MS fused silica capillary column (30 m × 0.25 mm × 0.25 μm).

Helium was used as the carrier gas at a flow of 1 mL/min. The samples (1 μL) were injected by the autosampler under a splitless mode at a temperature of 220 ◦C. The column temperature program was as follows: 50 ◦C for 2 min, 10 ◦C/min to 150 ◦C, 3 ◦C/min for 240 ◦C, 240 ◦C for 5 min, 10 ◦C/min for 300 ◦C and 300 ◦C for 5 min. The ion source temperature of the mass spectrometer was 200 ◦C, the temperature of the transfer line was 250 ◦C and the temperature of the quadrupole was 150 ◦C. The compounds were quantified in the selected ion mode, and the calibration curve was quantified with the internal standard. There were three parallel samples in each species. The method blanks and procedure blanks were prepared following the same procedure.

The quantification was performed by the internal standard method. All of the solvents used were HPLC-grade pure (J&K Chemical, Beijing, China). All of the glassware was cleaned by using an ultrasonic cleaner (KQ-500B, Kun Shan Ultrasonic Instruments Co., Ltd., Kunshan, China) and heated to 400 ◦C for 6 h. In the sampling process, three parallel samples were been collected from each sample site. The laboratory blanks and sample blanks were analyzed with the true samples. A total of 16 priority control PAHs were measured. The PAH individuals, abbreviations as well as method recoveries in different environment media and aquatic organisms are shown in Table 1.

#### *2.3. Positive Matrix Factorization (PMF)*

In this study, positive matrix factorization method was applied in order to quantitatively identify the major sources. PMF is a useful factorization methodology that can determine source profile and contribution [27,28]. The PMF model can be expressed as follows:

$$X = GF + E \tag{1}$$

where *X* is the concentration matrix, consisting of n samples and m concentrations of the compounds (n × m); *G* is the factor contribution matrix; *F* is the factor profile matrix; and *E* (n × m) is the residual matrix. The elements of residual matrix are denoted as the following:

$$e\_{ij} = x\_{ij} - \sum\_{k=1}^{p} g\_{ik} f\_{ki} \tag{2}$$

where *xij*, *fki* and *gik* are the corresponding elements of *X*, *F* and *G*, respectively. Nonnegativity constraints are imposed on the contribution and profile matrices, and PMF

simultaneously weights individual data points based on uncertainty. *Q*(*E*) is an object function and a criterion for the model, defined as the following:

$$Q(E) = \sum\_{i=1}^{n} \sum\_{j=1}^{m} (e\_{ij} / s\_{ij})^2 \tag{3}$$

where *sij* is the uncertainty of the *j*th compound in the *i*th sample.

#### *2.4. Multi-Pathway Exposure and Risk Assessment*

The BaP equivalent concentration (*BaPeq*) and toxicity equivalency factors (TEFs) were used to express the effects of exposure to mixtures of PAHs on health [29]. *BaPeq* is directly derived from the mass concentrations of different PAHs using TEFs. Therefore, they can be directly compared and contrasted [30]. In order to evaluate the total exposure to dietary PAHs, *BaPeq* based on BaP toxicity was determined using the following equation:

$$BaP\_{eq} = \sum \mathbf{C}\_{\mathbf{i}} \times TEF \tag{4}$$

where *Ci* is the concentration of the PAH species in food, and *TEFi* is the toxic equivalence factor of the PAH's congener *i*. (Table 1).

In accordance with the *Exposure Factors Handbook* [31], the lifetime average daily dose (LADD) of PAH exposure through inhalation (air and particle), aquatic product ingestion, and water intake was calculated as follows:

$$LADD = \frac{\mathcal{C} \times IR \times EF \times ED}{BW \times AT} \tag{5}$$

where *C* is the concentration of PAHs in the environment media, and *IR* is the intake rate of PAHs through inhalation (IRinh, m3·day−1), water intake rate (IRwater, mL·day−1) and aquatic product intake rate (IRinh, g·day−1). *EF* is the exposure frequency (day·year−1); *ED* is the exposure duration (year); and *BW* is body weight (kg). *AT* is the average lifespan for carcinogens.

High uncertainty exists in risk assessment. Sample measurement errors were inevitable. There are also uncertainties in the parameters and estimates. In probabilistic risk assessment, exposure parameters are considered as random variables. In order to quantify experiment uncertainty and its impact on the estimation of expected risk, a 10,000 times Monte Carlo (MC) technique was used. The Crystal Ball software was employed to implement MC simulation.

#### **3. Results and Discussions**

#### *3.1. PAH Residual Levels in Environment Media*

The levels of PAHs in the environmental media from Lake Chaohu are presented in Table 2. Sixteen priority PAHs were all detected during both the gaseous phase and particulate phase. In comparison, the detection rates of PAHs with higher than four rings in the water phase were very low due to their poor hydrophilic. The total concentrations of 16 priority PAHs (PAH16) in gaseous and particulate phases were 59.4 ± 51.4 ng·m−3 and 14.2 ± 23.5 ng·m−3, respectively. The average gaseous PAH16 concentrations in urban and rural sites were 3.59 times and 4.95 times higher than in particles. The atmospheric PAH16 residual level in this study was lower than the values reported in Guangzhou (337 ng·m−3) [32] and in Tianjin (752 ng·m−3) [33], but it was greater than the highest level reported in mountain Taishan (9.07 ng·m−3) [34]. Compared with data reported abroad, the PAH16 level was higher than data reported in Chesapeake Bay (5.31~71.6 ng·m−3) [35], Athens (4.8~76 ng·m−3) [36] and in southwest Europe (0.32 ng·m−3) [37]. PAH16 in the muscle of fish from Chaohu Lake (also including snail and shrimp) was also comparable with data reported from other freshwater fish in Hebei (4.76–144 ng/g) [20] and less than data reported in Shanxi (160 ng/g) [38]. Generally speaking, the PAH16 contents

in environment media and fish indicated a low PAH pollution level in Lake Chaohu. In order to compare the toxicity of difference environment media, the concentrations were converted to *BaPeq* concentrations. Although the gaseous phase had much higher PAH16 content, the *BaPeq* was much higher in particles. The particulate *BaPeq* in urban and rural were 11.2 and 5.51 times higher than those in the gaseous phase. The difference between PAH16 and *BaPeq* can be attributed to PAH composition in gas and particles.

**Table 2.** Residual levels of PAH16 and *BaPeq* in environmental media.

PAH16: Tthe sum of 16 PAH components; GM: geometric mean; SD: standard deviation.

Spatial difference between urban and rural sites were compared by using the Kruskal– Wallis test. No significant difference was found between two sampling sites except gaseous *BaPeq*. A *p* < 0.05 significant difference was detected between urban and rural gaseous *BaPeq* concentrations. The results showed that the concentrations in most environment media were similar in urban and rural area. There may be two reasons accounting for the small spatial difference. First, both Chaohu City and Zhongmiao Town had small populations. No obvious different lifestyle was found between people in urban and rural areas. In particular, there is no obvious heating season in the area around Chaohu. Thus, local emission sources in urban and rural area were not obvious. Second, the city and the town were far away from the local thermal power plant and other industrial pollution sources, resulting in low local pollution levels.

The PAH compositions in environment media were illustrated in Figure 2. It can be observed that the water phase was dominated by the low molecular PAHs. PAHs with less than three rings accounted for 95.0% and 93.3% of total PAHs in water. The same ratios in gaseous phase were 87.1% and 82.9% in urban and rural samples. In the particulate phase, however, the PAHs with more than or equal to four rings contributed to 80.2% and 81.6% of the total PAHs. The high proportion of high molecular weight PAHs results in the increase in toxicity due to the high TEF of high molecular weight PAHs.

Spearman's correlation analysis was used to detect the relationship between content of PAHs in gaseous, particulate phase, water and in aquatic animals. As a result, a significant positive correlation was found between gaseous phase and water phase (Figure 3a). It is reported that wet deposition, dry deposition and gas exchange across the air–water interface are the three major ways that PAHs can enter the water system. The result suggested that gas exchange across the air–water interface is probably an important way for atmospheric PAHs to affect the aquatic system [39]. In contrast, the relationship between particulate and dissolved PAHs was not significant (*p* = 0.116). This can be partly explained by the solubility of different composition in particles and gas phase. The solubility of high molecular weight PAHs was lower than the solubility of those with low molecular weight. On the other hand, high molecular weight PAHs dominated the particle component. Thus, concentration levels of PAHs in particles had a weak correlation with their concentration levels in water. Significant positive correlation was also found between PAH in water and in fish tissues, which indicates the effect of environmental concentration on organisms.

**Figure 2.** PAHs composition in different environment media.

**Figure 3.** Spearman correlation between PAH levels in environment media (**a**) water and gaseous phase (**b**) and fish tissues and gaseous phase (**c**) water and fish tissues.

#### *3.2. Source Apportionment*

#### 3.2.1. Molecular Diagnostic Ratios

PAHs can be formed by multiple anthropogenic activities such as combustion of fossil fuels, or they can be formed naturally in the environment by oil seeps and plant debris and forest and prairie fires. Some methods have been established in order to identify PAH sources: for example, molecular diagnostic ratios (MDRs), the principal component analysis (PCA) method [37,40], the chemical material balance (CMB) model [41], the positive matrix factorization (PMF) method and stable carbon isotopic ratios analysis [42]. In this study, MDRs and PMF were used to identify the major sources and to obtain a reliable conclusion.

The MDR theory [17,43] is based on the hypothesis that some PAH ratios remain constant between the source and the receptor. The MDR method has been widely used in the identification of preliminary sources [44]. Based on the monitoring data, four normally used ratios were employed in this research. Common ratios used include Ant/(Ant + Phe) (mass 178), Fla/(Fla + Pyr) (mass 202), BaA/(BaA +Chr) (mass 228) and IcdP/(IcdP + BghiP) (mass 276).

The ratios of mass 178 increased from January and reached the peak in October (Figure 4). Generally, the ratios varied near the value of 0.1. The lowest value was obtained during the winter, and the highest value was achieved in the summer. For mass 178, a ratio < 0.10 usually is taken as an indication of petroleum, whereas a ratio > 0.10 indicates a dominance of combustion. From November to January, petroleum was the dominant source, and emission from combustion became the greatest contributor during the rest of the year. For mass 202, a ratio of 0.50 is usually defined as the petroleum/combustion transition, and point ratios between 0.40 and 0.50 are more characteristic of liquid fossil fuel (vehicle and crude oil) combustion, whereas ratios > 0.50 are characteristic of grass, wood or coal combustion. The ratios in our study indicated a strong influence of biomass and coal combustion before November, which was consistent with the result of mass 178. For mass 228, BaA/(BaA +Chr), ratios < 0.20 indicate petroleum sources, ratios from 0.20 to 0.35 indicate either petroleum or combustion and those > 0.35 imply combustion. The ratios from May to October were higher than 0.5, which indicated grass, wood or coal combustion. From November to April, the ratios were between 0.40 and 0.50, which was the characteristic of liquid fossil fuel combustion. For IcdP/(IcdP + BghiP), ratios < 0.20 likely indicates petroleum, those between 0.20 and 0.50 imply liquid fossil fuel combustion and ratios > 0.50 imply grass, wood and coal combustion. The ratios of mass 228 in our research indicated a mixture of petroleum or combustion. Combustion was the potential source from July to December. The ratios in our research suggested a potential combustion source [45].

The MDRs indicated that PAHs in Lake Chaohu were mainly from combustion and vehicle emission. We can also found that differences existed when we were using different PAH ratios. The use of PAH MDRs has been criticized in the past due to low accuracy. Overlap areas were reported between commonly applied ratios associated with different types of PAH emissions [44,46]. According to a study of MDRs based on the inventory and monitoring data over 20 years, it was found that the use of MDRs does not respond to known differences in atmospheric emission sources unless the source is strong [47]. due to the limitations of MDRs, a PMF model was also applied to detect the potential sources.

**Figure 4.** Seasonal and spatial variation of four ratios in atmospheric samples from Lake Chaohu.

#### 3.2.2. PMF Results

When the number of factors for PMF is three, good simulation has been achieved for most PAHs. Therefore, three components were extracted, and the source profiles of three factors are illustrated in Figure 5. It was found that factor one is predominately weighted by low molecular weight PAHs. Factor two is heavily weighted by the middle molecular weight PAHs, and factor three has a higher contribution in the PAHs with more than four rings. The first factor is predominately weighted by Acy, which has been proved to be a tracer of combustion of straws. In addition, factor one has high load on Fla, Pyr and Chr, which are the combustion products of firewood. Therefore, factor one appears to be biomass combustion. Factor two is predominately weighted by Phe and BbF. According to the literature, Flo, Phe and Ant are predominantly considered as coal combustion profiles. The high BbkF and Chr loads are also a typical sign of Chinese domestic coal emissions [2,17,48,49]. Factor two is supposed to be coal combustion emission. For factor three, BghiP has been identified as tracers of auto emissions, and IcdP is considered as a marker of diesel emission. It can be concluded that factor three represents emission from vehicle.

The percentages of the sources from the factors were estimated by PMF. These results showed that the sources from biomass combustions, coal combustion and vehicle emission accounted for 43.6%, 30.6% and 25.8% of the total PAHs, respectively. Compared with the result of MDRs, a similar conclusion was obtained by PMF and MRDs. It can also be observed that coal combustion plays an important part in local PAH emission. Our conclusion is different from some research studies conducted in China. Most domestic studies show that coal is the main source of pollution in China, which is related to the fact that coal is the main energy material in China. China is the largest coal producer and consumer in the world [50]. According to the National Bureau of Statistics, the coal production of China in 2018 was 3.5 billion tons, nearly half of the world's production. Thus, the high coal consumption produces high PAH contribution to the atmosphere. In this study, the atmospheric samples were collected from Mushan island and Chaohu City. The former is a rural area far away from the city, while the latter is a small city with a small population and is a suburb surrounded by rural areas. Compared with coal combustion, firewood combustion is a more important method supplying energy in the research area.

**Figure 5.** Factors of PMF analysis in Chaohu dustfall.

#### *3.3. Exposure through Different Pathways*

#### 3.3.1. Derivation of Exposure Parameters

Three physiological and exposure behavior parameters including bodyweight, water intake rate and inhalation rate were collected from *Exposure Factors Handbook of Chinese Population* [51]. The distribution modes of parameters were fitted by regression models. Specifically, first, we considered a normal distribution for BW and log-normal distributions for the water intake rate and inhalation rate because the normal and log-normal distribution models are the most widely applied in studies on the exposure parameters [52,53]. Second, the 5th, 25th, 50th, 75th and 95th percentiles were collected from the *Exposure Factors Handbook*. The distribution of the parameters was fitted using the Gaussian function. The fit curves and parameters are shown in Figure 6 and Table 3, respectively. It can be observed that relatively good regression results were obtained for all three parameters. In addition, the intake rates of fish consumption in urban and rural were considered as constants. According to a survey conducted between 2010 and 2013, the average daily intakes of freshwater products per capita were 19.0 g/day and 11.1 g/day for the urban and rural populations, respectively [54].

**Figure 6.** Fitting of (**a**) bodyweight, (**b**) water intake rate and (**c**) inhalation rate of urban and rural population in Anhui Province.

**Table 3.** Exposure parameters of adults for Monte Carlo simulation.

#### 3.3.2. Estimation of Exposure Doses

The lifetime average daily doses (LADD) of *BaPeq* exposure through inhalation (air and particle), aquatic product ingestion and water intake were calculated by 10,000 iterations of Monte Carlo simulation. The results were illustrated in Figure 7. *BaPeq* exposure through fish ingestion had the greatest contribution, with the average value of 3.01 × 10−6 (mg·kg−1·d−1). The LADD of residents in urban and rural were 2.97 × 10−6 (mg·kg−1·d−1) and 3.07 × 10−6 (mg·kg−1·d−1), respectively. Although the fish intake rate in urban residents was 1.71 times higher than that in rural, no obvious difference was observed between two population groups. Exposure through particle inhalation has a comparable contribution of 2.54 × 10−6 (mg·kg−1·d−1), which is one and two order of magnitude higher than exposure through gas inhalation and water ingestion, respectively. Our results are consistent with previous research. It has been reported that, for most non-occupationally exposed individuals, diet is the main route of exposure [55,56]. Inhalation of gaseous and particulate also had an ignorable contribution relative to the total *BaPeq* exposure.

**Figure 7.** *BaPeq* exposure through different pathways: (**a**) gaseous phase; (**b**) particulate phase; (**c**) water intake; and (**d**) fish intake.

#### *3.4. Carcinogenic Risk Assessment*

The probabilistic cancer risk for urban and rural residents were obtained based on the results of Monte Carlo simulation and presented in Figure 8a,b. The median total PAHs ILCR was 3.1 × 10−5 to 3.3 × 10−5 in urban and rural residents, which showed a slightly higher carcinogenic risk for rural residents than that for urban residents. The difference can be attributed to the variation of exposure behaviors. A one in a million chance of additional human cancer over a 70-year lifetime (ILCR = 10−6) is considered acceptable or inconsequential, and one in ten thousand or greater (ILCR = 10−4) is considered serious (US EPA). Most of the surrogate samples generated by Monte Carlo distributed between the range of 10−6 and 10−4, which indicates potential ILCR risk. Due to the different end points, the carcinogenic slopes of inhalation and ingestion are quite different. Therefore, fish ingestion had much higher contribution than particle inhalation. Inhalation (including gas and particles) accounted for 26.5% and 29.1% of the total risk in urban and rural areas, respectively.

Parameter sensitivity was quantitatively assessed by the Spearman's rank correlation coefficient of parameter and risk. The result was shown in Figure 8c,d. As the dominant exposure pathway, fish consumption had the highest risk contribution. PAH concentration in fish had the highest influence on the results of risk assessment. As the most important protection factor, the bodyweight was the second most sensitive parameter. PAH concentration associated with particles was the third most sensitive parameter in this research. It can also be observed that the sensitivity of exposure behavior parameters were relative low compared with media concentrations. This can be attributed to the small difference of these parameters among adults. As the most important behavior parameter, fish consumption rate was not treated as variables in the study due to the limited data support. However, the uncertainty of intake rate of fish should not be ignored due to the high contribution of fish intake on the total risk, which requires further research and more data support in the future.

**Figure 8.** Distributions of incremental lifetime cancer risk (**a**,**b**) and parameter sensitivity for urban and rural populations (**c**,**d**) derived using the Monte Carlo simulation.

#### **4. Conclusions**

PAHs in major environment samples were collected in Lake Chao for the entire year. Exposures through four pathways were estimated. ILCRs were characterized. The results found that the gaseous concentration had influence on the PAHs in the water according to the Spearman correlation analysis and may further affect the PAH content in fish tissues. Atmospheric transport is the source of the entire water system. The results of sources apportionment on atmospheric samples indicated the high contribution of biomass combustion. Probabilistic risk assessment suggested that inhalation and fish ingestion are two major pathways of PAH exposure, which are also the key processes in PAHs' risk control for people near Lake Chaohu.

**Author Contributions:** Investigation, W.H., Q.H., X.K. and W.L.; writing—original draft, N.Q.; writing—review and editing, F.X. and Q.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by China CDC Key Laboratory of Environment and Population Health, National Institute of Environmental Health, Chinese Center for Disease Control and Prevention (2021-CKL-01), by the Fundamental Research Funds for the Central Universities (FRF-TP-18-071A1) and by the National Natural Science Foundation of China (NSFC) (41977312, 41503104).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** We also thank the help from National Environmental and Energy Base for International Science and Technology Cooperation.

**Conflicts of Interest:** All authors declare that they have no conflicts of interest or financial conflicts to disclose.

#### **References**


**1.1 Genetic resources of** *Fragaria*

Genus *Fragaria* includes 20 species (**Table 1**) distributed in the Northern temperate and Holarctic zones [7, 8]. In *Fragaria*, four specific fertility classes are primarily associated with their ploidy level or chromosome number. Most widely distributed species *F. vesca* is diploid in nature (2n =14) and the most cultivated

#### **Table 1.**

*Species of Fragaria, ploidy and global distribution.*

species is *F. × ananassa* (2n = 56) is octoploid in nature. Diploid strawberries can be differentiated by their foliage, flower and fruits, inflorescence structure and plant habit. The *F. vesca* genome size was estimated at 197 Mb [9] and 206 Mb [10] and of *F*. *ananassa* 698 Mb [11]. Out of 20 species listed in **Table 1**, three species of *Fragaria* occur wild in the Himalayas in India [12–14]. These are: *1. F. nilgerrensis* Schecht. (Nilgiri strawberry): Habitat: A creeping herb found in Nilgiris and higher hills of north-eastern India. Uses: The fruit is pinkish, sub-acidic and juicy; *2. F*. *nubicola* Lindl. (Alpine Strawberry): Habitat: A herb distributed in the temperate Himalayas. Uses: Fruits are edible; and *3. F*. *vesca* L. (perpetual strawberry): Habitat: A herb found in the higher altitudes of temperate Himalayas. Uses: The red delicious fruits are edible.


**2. Theoretical and experimental approach to the molecular arrangement in organic liquids**

In this section the fundamentals of molecular structure and its formation mechanism in liquid organic phase are presented, including the experimental data and theoretical study (dynamic simulation theory) in combination with the authors analysis of this approach.

Although there are numerous papers and books devoted to this problem, none of them fully resolved the problem of organic liquids molecular arrangement, and the character of long-range molecular interaction in liquid bulk is still ambiguous. Therefore, in the first part of this chapter we tried to present the current state of considered problem in terms of theoretical and experimental approach to evaluate the role of DFT calculation in the solution of the presented problem.

However, at the beginning of this consideration we would like to review the dipole–dipole interaction concept briefly because it is the oldest and most common theory of liquid formation. In fact, this theory is based on the model of Coulomb interactions between polar systems or temporary dipoles caused by fluctuation in electronic distribution.

#### **2.1 Classical dipoles interaction model**

Cooled and pressed real gases become liquid due to the intermolecular interactions called van der Waals. This intermolecular interaction energy can be defined as the heat of a liquid evaporation or rather the difference between the vaporization temperature and the work of expanding one mole of gas at the atmospheric pressure. Value of this interaction at the boiling point of gases is 1–3 kJ/mol. As it was shown by the quantum mechanical calculations, the energy of the van der Waals interactions consists of electrostatic, inductive, and dispersive components [1].

#### *The Formation Mechanism and Structure of Organic Liquids in the DFT Challenges DOI: http://dx.doi.org/10.5772/intechopen.100429*

The so-called orientational interaction of polar molecules is the most important part of the electrostatic interaction. Its essence, known as the Keesome effect, consists in the orientation of two interacting polar molecules with identical dipole moments that leads to the minimization of the system energy. In this case, the head-to-tail orientation becomes the most advantageous configuration. The energy of the orientational binding can be calculated as the sum of the Coulomb attractions and repulsion of the pole charges in the dipoles, expressed in terms of the dipoles [2].

Inductive interaction, or the Debye effect [3], is the interaction of the constant dipole moment and the induced dipole moment, arising due to an additional charge separation. Inductive interactions occur with the formation of noble gases hydrates with the dissolution of polar substances in the non-polar liquids and they are valuable only for the molecules with significant polarizability (ex. molecules with conjugated bonds) [3].

The dispersion interaction, or the London effect [4], arises between the electrons in the interacting molecules. There are molecules that have no dipole moment – for example, homoatomic molecules of noble gases. The additivity of dispersion forces is manifested in adsorption, in the processes involving gas condensation, etc. Dispersion forces play the important role not only in the individual molecules, but also in the macroscopic particles (ex. colloids). However, these forces are relatively weak.

The intermolecular interaction is the summarizing action of attraction and repulsion forces. At large distances, the attraction prevails, and at very short distances the repulsion is the main contributor. The attractive forces of Van der Waals are long ranged [5, 6], and the attraction energy rather decreases with distance. Van der Waals interaction forces at the equilibrium distance is small: ~ 1–5 kJ/mol, which is considerably less than the chemical bond energy.

The main conclusion is that Coulomb and Van der Waals interactions can relate to the formation of liquid from the gaseous substances at the cooling or pressing. However, this theory is very limited for organic liquids. There are at least three fundamental reasons for that: i) the intermolecular forces in organic liquids are sufficiently stronger, accounting their vaporization enthalpy [7]; ii) there is no direct dependence of liquid stability on dipole and molecular volume of consisting molecules in many cases [7]; iii) the interaction of dipoles in organic fluids should have a certain orientation, while the molecules in liquid do not have any anisotropy [5, 6]; iv) this conception ignores the specific non-valence interaction in liquid systems [7].

## **2.2 Thermodynamic data source**

The concept of the specific interaction in liquid systems based on the thermodynamic data is developed in [7]. The regularities in the homological series exhibit no correlation between the molecular mass and the vaporization enthalpy in many cases. The numerous thermodynamic data indicate the existence of specific nonvalence interaction in organic liquids. In case of substances with the different variant of hydrogen bond (ex. with OH–, NHi–, C=O– or COOH–group) the specific interaction providing molecular structure can be explained by classical H–donor– H–acceptor interaction [7].

However, for the liquid saturated hydrocarbons, alkenes, and alkynes this approach is not applicable. The author of [7] has suggested the scheme of hydrocarbon interactions. The main idea is that the carbon atoms in chains have a different negative charge owing the shift of 2 *s*–electronic pair between the carbon atoms. The electron density of carbon with a larger negative charge transfers partially to the carbon atom with the reduced negative charge.

In other words, the alkane molecules can be bound in the liquid state so that the carbon atoms with an enhanced negative charge would be closer to the carbon atom with a reduced charge. The carbon atoms of the methyl groups in hydrocarbons that have the enhanced negative charge can participate in the specific interaction between molecules under the donor–acceptor mechanism resembling the hydrogen bond.

In author's opinion, the non-equivalence charges distribution in hydrocarbon chains leads to the appearance of charge shift not only for the carbon atom, but for the hydrogen atom as well. It means that the hydrogen atom of the neighboring molecules can form the donor–acceptor C‧‧‧H bond and even the dihydrogen bond. The comparison of thermodynamic data indicates the relatively high contribution of similar pair interactions in vaporization energy of organic liquids. These specific interactions are the nature of liquid saturated hydrocarbons formation mechanism and their stability without involving of the specific hydrogen bon. The evaluation of specific interaction energy is given in [7] and relates to the range 3–10 kJ/mol. These values correspond to the weak hydrogen bond [1, 2, 6].

The thermodynamic analysis of the vaporization enthalpy of unsaturated hydrocarbons leads to the conclusion that the presence of a double or a triple bond in the molecule results in a stronger electron density shift from one or two hydrogen atoms to the carbon atoms of neighboring molecule [7]. For example, the hydrogen atoms of the CH2–group in propene and CH–group in propyne possess enhanced positive charges, and the carbon atoms acquire higher negative charges. Therefore, the acceptor and donor properties of these substances are more pronounced than in the propane molecule (this statement are cited fully in accordance with [7]).

#### **2.3 Molecular light scattering (MLS) and X-ray data source**

In the 70s – 80s the professor at Moscow University M.I. Shachparonov has published a work, which was, unfortunately, limitedly known in the scholar community [8, 9]. In this paper he has predicted the existence of the specific interaction in the non-polar liquids, in which the hydrogen bond lacks. For benzene, the existence

**Figure 1.** *The stack of aromatic rings in liquid benzene.*

#### *The Formation Mechanism and Structure of Organic Liquids in the DFT Challenges DOI: http://dx.doi.org/10.5772/intechopen.100429*

of a new specific interaction in liquids – the formation of molecular π-complexes (molecules stacks) between aromatic rings – was suggested (**Figure 1**).

At present, this concept was confirmed and developed in the works of I.A. Abramovich and the coauthors [10–17], devoted to the structure of the liquid system in benzene and its substituted species. The specific interaction with different spatial geometry can form in a liquid bulk of molecules without hydrogen bonds. For example, in solution of *di*-chlorobenzene, there are strong interactions in neighboring molecules between chlorine atom and carbon in aromatic ring and between carbon atoms.

Since the arrangements of molecules in a crystal always correspond to the potential energy minima of the system, it can be used for the verification of the obtained data in LS experiments combined with the modeling procedure. This approach allowed to reveal the structure and characteristics of the mutual arrangement in molecules as well as the types of intermolecular contacts in the liquid phase [18, 19].

For the determination of intermolecular arrangement in liquid *di*-chlorobenzene and other compounds with chlorine atoms, the Cambridge crystallographic database for these organics in solid state were used. The intermolecular distances between molecules, obtained by MLS study, are 4–5 Å, which is considerably longer than the internal distances in molecules. However, the liquid system has a quite certain spatial geometry realized by this long-ranged molecular binding.

X-ray studies allow to define the internal parameters of molecules, but the intermolecular structure of organic liquids stays in an ambiguous conclusion: on one hand, in the liquid phase the distances between nearest molecules are longer than 4 Å, and on the other hand, the peak of distribution function was observed in dichloromethane at distance about 2 Å. However, in these experiments it was confirmed that the organic liquids have certain intermolecular composition and labile binding between identical molecules. However, this does not exclude the spatial transformation in liquid bulk [20–22].

#### **2.4 IR-data source**

The IR spectra of organic liquids in the middle- and high-frequency regions have "additional" bands [23] that cannot be expected by the normal coordinates' analysis [24–26]. The observed spectral phenomenon was interpreted in literature as a manifestation of the vibration anharmonicity or Fermi resonance (ex. [27, 28]). Such a version has a few contradictories; firstly, it is non-logic to mix the approach using the vibrational theory based on the classical mechanics of vibration (normal coordinates' analysis) and the quantum mechanical interpretation; secondly, many systems that are considered in this conception include heavy molecular weight atoms and their vibration anharmonicity is negligible; thirdly, the overtones and combination bands should have considerably lower intensities than the basic bands. At last, these bands remain in the solid state (in a low-temperature spectra), and they are observed for homological analogs [20]. We have suggested another variant of assignment based on the conception observing spectral features relating to the manifestation of specific interaction in liquids.

In the spectra of liquid CCl4, in accordance with the selection rules for the Td symmetry in IR spectra, only one stretching band should be active, but for the C3V symmetry, two bands of E– and A1–symmetry species are permitted. In the IR spectrum of liquid carbon tetrachloride two overlapping bands at 786 and 761 cm−1 having an approximately equal intensity were observed. These bands are well-resolved in IR spectrum of the low-temperature film recorded at 20 K [29].

In the gas phase two bands were also observed (at 795 and 776 cm−1), but they have a different counter structure: the first of them consists of overlapped components and another one is a single band. Therefore, we can conclude that carbon tetrachloride exists in gas phase not only in the single molecular shape that has Td symmetry, but also in the transformed shape having C3V symmetry. The first one is manifested by the band at 776 cm−1 and another one – by the band at 795 cm−1. The spectra recorded for different temperatures of gaseous CCl4 confirm this assumption [29].

It is reasonable to resume that the pyramidal structure relates to the cluster shape, in which the chlorine atom provides the binding between molecules. The shift of the chlorine atom in the cluster can occur in the condensed phase owing to the association of the molecules [29]. It leads to transformation of the molecular geometry to the almost planar D3h symmetry, in which A1 stretching band is forbidden in IR spectra.

In the high-frequency region of chloroform in the liquid phase, the bands at 3020 cm−1 (stretching vibration of the CH bond) and at 2401 cm−1 with a shoulder at 2435 cm−1 are detected (**Figure 2**, left side, spectrum 2) [30]. A similar spectrum is observed for bromoform: the bands at 3021 и 2256 cm−1 are shown (**Figure 2**, left side, spectrum 1). The bands in 2400–2200 cm−1 range cannot be assigned to overtones 2δ (bending of angle CHalH, where Hal-halogen), because the band have the isotopic D/H shift (**Figure 2**, right side), corresponding to the theoretical values (1.32–1.34) and their intensities relative to CH stretching band at 3020–3021 cm−1 are considerably stronger than it could be expected for anharmonic components. Besides, in bromoform molecule, the contribution of anharmonic components should decrease due to a significant increase of molecule mass. However, the opposite picture is observed in the spectra: the relative intensity of the band at 2256 cm−1 for bromoform is stronger than the band at 2401 cm−1 of chloroform [30].

The similar spectral picture is detected for water associates in liquid phase: the band of OH stretching vibration in 3400–3600 cm−1 region is combined with the band in 2200–2400 cm−1 range, assigned to stretching vibration of hydronium ion. This band manifests the hydrogen atom transfer in hydrogen bonded structures [31, 32].

The presented results indicate that the specific interaction between molecules exists in liquid haloforms due to the proton binding and its intermolecular shift, leading to the transformation of initial geometry.

The geometry of single benzene molecule is taken as a planar ring due to the conjugation of pz-orbitals and π-aromatic configuration appearance. This state corresponds to D6h symmetry point group. According to selection rules, only one stretching CH band (E-specie) should be active in IR spectra. However, in the real spectra of liquid benzene there are three bands (3092, 3071, 3036 cm−1) and in the

#### **Figure 2.**

*Fragments of FTIR spectra of chloroform and bromoform (left side) and chloroform-d in liquid phase (right side).*

#### *The Formation Mechanism and Structure of Organic Liquids in the DFT Challenges DOI: http://dx.doi.org/10.5772/intechopen.100429*

solid benzene (spectra were recorded at 20 K) there are four bands (3090, 3071, 3034, 3005 cm−1) (**Figure 3**, right side) in CH stretching region [33].

Three bands that have isotopic H/D shift (**Figure 3**, left side) close to the theoretical prediction also were detected in benzene-d6 spectrum. This spectral picture can be assigned to the existence of two molecular shapes existing in the liquid phase: planar, in which one IR band is active (D6h symmetry), and shape with two IR bands (A1- and E-species), corresponding to the C3V symmetry.

For the interpretation of the observed effect, we can assume that the benzene molecule exists in the liquid state as a cluster system. In this shape CH bonds deviate from the ring plane to the neighboring molecule (**Figure 4**). This leads to a distortion of the symmetry of the stretching vibrations of the CH bond.

In the middle IR range two bands at 1952, 1814 cm−1 were revealed (**Figure 5**, bands A), which cannot be assigned to the internal vibration modes. They also cannot be assigned to combination modes, because these bands have an isotopic shift close to the CH bond vibrations of aromatic ring, and the same bands are observed in substituted homologs of benzene as well (**Figure 5**). Besides, their amount in C6H5X spectra corresponds to the non-equivalent *o-, m-* and *p-*position in molecules [23]. These data confirm the idea that the benzene can form the π-stacks structure in the liquid phase (see section 2.3). The hydrogen atom in this system, as it was mentioned above, can interact with the carbon atom of the neighboring molecule in the stacks. Therefore, CH stretching band can shift to the middle IR interval due to the mixing of ν(C–C) and ν(C–H) stretching vibrations in the intermolecular bond (C–C–H‧‧‧C–C).

#### **Figure 3.**

*Fragment of liquid benzene-d6 FTIR spectrum at 295 K (left side) and solid benzene film at 20 K (right side) in CH/CD stretching region.*

**Figure 4.** *Hydrogen atom shift in molecular stacks of liquid benzene.*

**Figure 5.** *Fragment of liquid benzene (left side) and benzene-d6 (right side) FTIR spectrum in middle IR region.*

#### **2.5 Dynamic simulation data source**

This approach was developed in numerous works and their complete citation was beyond the scope of this study. Therefore, we would like to highlight the studies that present the molecular arrangement of organic systems as a supramolecular structure, forming not only under hydrogen bond, but also existing in the shape of identical molecular associates [34–37].

The obtained data allow to discuss the thermodynamic and kinetic aspects of the processes in liquid phase, in terms of the supramolecular organization. These investigations give the knowledge about the composition of the aggregates as well as the general approach to their transformation due to the intermolecular binding. In author's opinion, the dynamic simulation model combined with the experimental findings allows to explain the supramolecular formation mechanism under the long-ranged interactions. The model of the homogeneous molecules' association is a key to manifold conclusions of the liquid properties' nature [37].

This approach gives the complementary information about the structure of aggregates, which is not possible by other methods, especially when long-ranged molecular binding is present. The application of this method is useful to interpret some regularities of organic liquids, although it does not detail the specific interactions' appearance [34, 37].

#### **2.6 Section conclusion**

As we have shown in the short review of the modern approaches to the problem of molecular arrangement in organic liquids, all data sources agree that the structure of liquid phase can be described as a supramolecular system with non-polar or weak polar long-ranged interactions. This binding arises at the distances close to 4–5 Å, but its energy contribution is comparable with vaporization enthalpy of liquids [7, 17, 18, 20]. The thermochemical, wave scattering, X-Ray and dynamic simulation data interpret well the lability of the system, but do not explain its high stability.

Unlike the mentioned methods, IR data predict that the initial molecular geometry transformation in comparison to a single molecule occurs in the liquid system. These changes are caused by the intermolecular forces between neighbors in bulk, even if they do not correspond to the hydrogen bond [23, 29, 30].

Applied to this problem, the DFT calculation could be useful for explaining the reasons of the transformation under weak interactions in terms of general electronic distribution. Therefore, in the following section we consider the essential results of the DFT study referring to the problems discussed in this paper.

*The Formation Mechanism and Structure of Organic Liquids in the DFT Challenges DOI: http://dx.doi.org/10.5772/intechopen.100429*


*3.4. Data Sets*

The reference genome sequence and general feature format (GFF), coding sequence (CDS), protein sequence (PEP) files, and expression data are included in the GDS [20]. A summary of the genomic data currently available in the GDS is presented in Table 1. The versions of the genomes from top to bottom in the table are v1.0a2, r1.1, v1.0, r1.1, r1.1, r1.1, v4.0a2, and v1.0 [21–25].

**Table 1.** Statistics of the genome features for eight strawberry species.

#### *3.5. Completeness of the Genomes*

BUSCO provides measures for the quantitative assessment of genome assembly, gene set, and transcriptome completeness based on evolutionarily informed expectations of gene content from near-universal single-copy orthologs selected from OrthoDB [26]. On the species introduction page, we have integrated BUSCO5 results of the eight species, including complete and single-copy, complete and duplicated, fragmented, and missing orthologs. The results are presented as a bar graph (Figure 2), showing that *F. ananassa* currently has the most complete gene assembly of the eight strawberry species.

**Figure 2.** Genomic completeness of eight strawberry species evaluated using BUSCO V5 software. Fana, *Fragaria* ×*ananassa*; Fiin, *Fragaria iinumae*; Fnil, *Fragaria nilgerrensis*; Fnip, *Fragaria nipponica*; Fnub, *Fragaria nubicola*; Fori, *Fragaria orientalis*; Fves, *Fragaria vesca*; and Fvir, *Fragaria viridis*.

## *3.6. Phylogenomic Relationships among the Strawberry Genomes*

OrthoFinder [27,28] is a fast, accurate, and comprehensive platform for comparative genomics. It identifies orthogroups and orthologs, infers rooted gene trees for all orthogroups, and identifies all of the gene duplication events in the gene trees. To analyze the relationships among eight strawberry genomes, we constructed a phylogenetic tree using OrthoFinder (V2.5.2) software and included two additional species, *Rosa chinensis*, and *Arabidopsis thaliana* (Figure 3a).

**Figure 3.** Phylogenomic analysis of eight strawberry species. (**a**) A phylogenomic species tree of eight strawberry species; (**b**) *F.* ×*ananassa*; (**c**) *F. iinumae*; (**d**) *F. nilgerrensis*; (**e**) *F. nipponica*; (**f**) *F. nubicola*; (**g**) *F. orientalis*; (**h**) *F. vesca*, and (**i**) *F. viridis.* The picture (**c**) is from https: //en.wikipedia.org/wiki/File:Fragaria\_iinumae\_(fruits).jpg (accessed on 10 November 2021); picture (**i**) is from https://en.wikipedia.org/wiki/File:%D0%9A%D0%BB%D1%83%D0%B1%D0%BD% D0%B8%D0%BA%D0%B0\_(Fragaria\_viridis).jpeg (accessed on 10 November 2021); picture (**f**) is from http://www.fpcn.net/uploads/allimg/131107/2-13110G30231V8.JPG (accessed on 10 November 2021).

*Fragaria virginiana* and *Fragaria chiloensis* are the genomes of the progenitor species of *Fragaria* ×*ananassa*. However, the dispute over its diploid ancestor has lasted for more than half a century and is still unresolved. In 2019, Edger et al. speculated that it has four different diploid ancestors, *F. vesca*, *F. iinumae*, *F. viridis* and *F. nipponica* [11]. Unexpectedly, just a few months later, Liston and others completed a reanalysis of the same set of data, but they came to a completely different conclusion. They believed that the octoploid strawberry has only two existing ancestors, *F. vesca* and *F. iinumae* [29]. Edger et al. insisted on the previous conclusion [30]. The structure of our phylogenetic tree [31] clearly indicates that *Fragaria vesca* is closest to *Fragaria* ×*ananassa*. However, because of the low sequencing and assembled technology or gene introgression, there is no direct evidence indicating the origin of the cultivated strawberry.

The pictures below the tree (Figure 3b–i) shows *F.* ×*ananassa*, *F. iinumae*, *F. nilgerrensis*, *F. nipponica*, *F. nubicola*, *F. orientalis*, *F. vesca*, and *F. viridis*. Photographs b, d, e, g and h were provided by our colleague, Dr. Qiao and the others were obtained from Wikipedia or Baidu.

#### *3.7. Genomic Comparison of Gene Orthogroups*

To provide an overview of the comparison among these strawberry genomes, we compared the number of gene orthogroups identified by Orthofinder in the strawberry genomes. We uploaded the orthogroup data to an online Venn diagram tool (https://www. vandepeerlab.org/?q=tools/venn-diagrams, accessed on 10 November 2021) to generate a Venn diagram showing the shared and unique gene orthogroups in *F. iinumae*, *F. nilgerrensis*, *F. nipponica*, *F. vesca*, and *F. viridis*. The gene orthogroup numbers are shown in each segment of the diagram (Figure 4); 13,766 gene orthogroups were shared among the five species, and 13,380 gene families appeared to be unique to *F. nipponica.*

**Figure 4.** Venn diagram of gene orthogroups in five diploid and wild *Fragaria* species. Comparison of the number of shared gene families among five diploid strawberries, *F. iinumae*, *F. nilgerrensis*, *F. nipponica*, *F. vesca*, and *F. viridis*.

The number of predicted genes was quite higher in F. nipponica, compared with the other four species, and the number of specific genes for F. nipponica was also extremely high (13,380). The reason why so many genes were predicted in *F. nipponica* is because there were more than 80,000 proteins annotated, only using Illumina sequencing technology, in 2014. In the future, with the improvement of the technology of sequencing, this problem will disappear.